Instructions to participate

We distinguish between "physical constants" (first table) and "parameters" (second table). The former should not be used for tuning whereas the latter can be used if needs be. We suggest that most poorly constraint parameters should be used for tuning, i.e. probably the ones directly related to water flow.

The given parameters are based on the ones used by the GlaDS-model (Werder et al. 2013), therefore, please report the different parameters your model uses and we can include them here. Script files containing the parameters are located in the parameters/ folder in the bitbucket repository.

Note that the ice flow constant \(A\) is for the usual channel closure relation of the form \(\frac{\partial S}{\partial t}=\frac{2A}{n^n} SN^n\) with channel cross-sectional area \(S\), effective pressure \(N\) and Glen's \(n\). This is how \(A\) is usually defined, e.g. in Cuffey & Paterson 2010. In the GlaDS paper a closure relation of the form \(\frac{\partial S}{\partial t}=A' SN^n\) is used, and the value of the constant is then \(A'=\) 2.5e-25 Pa^-3 s^-1.

To obtain comparable results for different models we suggest to use test case A3 and A5 as a "common ground", which is achieved by tuning models to the output of GlaDS (see figures below and provided as ascii files for A3 and A5). The focus should lie on tuning to the effective pressure \(N\), but sheet \(q\) and channel discharge \(Q\) can also be used (please indicate in the questionnaires).

The parameters thus obtained are then used for the subsequent experiments. A3 is chosen as a case where GlaDS (for the used parameters) produces a distributed system only, whereas A5 features channelisation up to around mid-domain.

Models which do not include both inefficient and efficient drainage may need to use a different set of parameters to tune to A3 and A5. Or may choose to only apply the model to part of the test cases.

If you feel that your model should not or cannot be tuned to A3 and/or A5, then run the model with the most appropriate set of parameters. However, in such a case it is encouraged to submit two sets of results: one tuned as best as possible and one using the most natural set of parameters.

This geometry is used for the suites A-D and mirrors a synthetic land-terminating ice sheet margin as was used in Werder et al. 2013, but with a few modifications. The domain is 100km long in the \(x\) direction and 20km wide in the \(y\) direction with the terminus along \(x=0\). The bedrock is a flat surface at 0m elevation and the surface is defined by a square root function¹:

surface(x,y) = 6*( sqrt(x+5e3) - sqrt(5e3) ) + 1
bed(x,y) = 0

This means that the maximal ice thickness is 1521m, and minimal 1m.

One dimensional models (1D) should use the following functions:

# 1D models
surface(x) = surface(x,0)
bed(x) = 0
width(x) = 20e3

The two suites E & F are performed on a 2D valley geometry to investigate the impact of the smaller glacier-size and of the bedrock shape on the behavior of the models. The geometry is based on the shape of Bench Glacier (Alaska, approximately 6km by 1km). The bed-geometry has one parameter para which determines whether the bed has an overdeepening or not, with para=0.05 mimicking Bench Glacier with no overdeepening.

The defining functions are

surface(x,y) = 100(x+200)^(1/4) + 1/60*x - 2e10^(1/4) + 1
bed(x,y, para) = f(x,para) + g(y) * h(x,para)

with the helper functions f(x,para), g(x,para), and h(x,para). f determines the flow-line geometry and it is constructed such that f(6e3, para)==surface(6e3,0), i.e. the ice thickness is always 0 at the upper end of the glacier (=6e3m). g determines the cross-sectional geometry which is modified by h. The function h is chosen such that the glacier outline is independent of para. They are given by

para_bench = 0.05
f(x,para) = (surface(6e3,0) - para*6e3)/6e3^2 * x^2 + para*x
g(y) = 0.5e-6 * abs(y)^3
h(x, para) = (-4.5*x/6e3 + 5) * (surface(x,0)-f(x, para)) /
               (surface(x,0)-f(x, para_bench)+eps())

The half-width is given by

outline(x) = ginv( (surface(x,0)-f(x,0.05))/(h(x,0.05)+eps()) )
ginv(x) = (x/0.5e-6).^(1/3) # the inverse of g

which gives a maximum full-width of approx. 1080m.

1D simulations with flowline and flowband models should be performed using the centre-line of the valley as their flowline. The geometry is given by

# 1D models
surface(x) = surface(x,0)
bed(x,para) = f(x,para)
width(x) = 2*outline(x)

Test cases A is performed using different steady and spatially uniform water inputs into the sqrt geometry. The aim is to show simple steady state configurations and to allow for model tuning, see Model tuning. The water input values are as follows:

To obtain the source in m/d multiply by 86400, and for m/a by 31536000 (i.e. one year is defined to have 365 days). For flowline models, multiply the source by 20km to obtain the total for the width.

The importance of input localisation is investigated in test B. To reach this goal, the spatially uniform input which was used in the preceding simulation is replaced by a moulin input described by an input flux and an input location. The simulations are run to steady state from the given water input (which has the same total as the one of simulation A5 plus distributed basal input equivalent to A1). The varying parameter here is the number of moulins that is used and so the amount of water that is injected into each moulin decreases with increasing number of moulins. Moulins positions are specified on a regular 1km grid to allow the participation of models based on regular or unstructured grids. The moulins need not be placed at the exact specified location but as close as possible. The location and input of each moulins is given in a .csv with four columns (moulin index, X position [m], Y position [m], input value [m³/s]). For each simulation, the position of each moulin is defined at random excluding the boundary grid points and the ones in the first five kilometres from the glacier terminus. Additionally to the moulin input, use distributed input as in model run A1 (i.e. representing a small basal melt contribution).

1D flowband models should collapse all moulins onto the flow-line and sum the input of colocated moulins as needed.

Test C is designed to investigate the effect of the diurnal melting cycle on the response of the subglacial drainage system, i.e. short time scale dynamics. The starting point for this experiment is the steady state achieved in simulation B5 (i.e. a steady input into moulins). The amplitude of the runoff is changed for the different simulations of the test. The runoff into each moulin is given by

runoff(t, ra) = max(0, moulin_in * (1 - ra*sin(2*pi*t/day)))

with time t (s) and day is seconds per day. The background moulin input (moulin_in, the runoff for experiment B5: B5_M.csv) is modulated by a sine function, set to zero if negative. The value of the relative amplitude (ra) of the signal is given on the table bellow for the different experiments of the test. Added to this is a distributed basal input as in A1:

runoff_basal(t) =  7.93e-11 # m/s

The model should be run until a periodic state is reached, then one day is submitted with output interval 1h.

Again, 1D flowband models should collapse all moulins onto the flow-line.

Test case D simulates the seasonal evolution of the drainage system, i.e. the long time scale evolution. This test uses initial conditions from test case A1 which represent the water input during winter. From this starting point, a seasonal cycle is applied to the water input. The model should be run for enough years until it settles into a periodic state. Once this state is reached, please provide output for one year at daily resolution. The forcing is computed from a simple degree day model driven by a temperature parameterization.

The temperature at 0m elevation is given by

temp(t) = -16*cos(2*pi/year*t)- 5 + DT

Where the temperature (temp) is function of the time (t) and year=31536000 is the number of second per year. The different runs of this suite are achieved by modifying the value of delta-temperature DT as presented in the table bellow.

The runoff (distributed) is then computed from the following degree day model formulation

runoff(z_s,t) = max(0, (z_s*lr+temp(t))*DDF) + basal

where z_s is the surface elevation, lr=-0.0075 K/m is the lapse rate and DDF=0.01/86400 m/K/s is the degree day factor. basal=7.93e-11 m/s is a basal melt rate equal to the source of scenario A1.

Scripts are available for the computation of these forcings on the Bitbucket repository: Julia (reference), matlab

Test case E is designed to investigate the effect of bedrock slope on the models. The common base for this experiment is the synthetic valley geometry modelled after the Bench Glacier geometry. The different simulations of this test are achieve by altering the shape of the bedrock to define a more or less pronounced overdeepening, see section "valley": Bench Glacier-like geometry. The water input on this experiment is uniformly distributed at the bed of the glacier with a value of 1.158e-6 m/s (twice the rate of scenario A6).

It is suggested that models which include a pressure-melt term, run this experiment with and without this term.

Test F runs a seasonal forcing for the synthetic Bench Glacier using topography parameter para=0.05 as in E1 (and all other parameters as in the Suite E). The water forcing mirrors Test D. First run your model to a steady state with water input as in A1 (m=7.93e-11); use this steady state to start all the model runs from. The model should be run for enough years until it settles into a periodic state. Once this state is reached, please provide output for one year at daily resolution. The forcing functions are as in Suite D but with different temperature forcings: