Fig. 1. Our model ice-dammed reservoir system. (a) A subglacial channel is hydraulically coupled to an ice-dammed lake that receives meltwater at a rate Q; The channel receives a constant and uniform supply of water along its length M. (b) The prescribed basic hydraulic gradient ~ (blue curve) is negative near the lake - a topographic seal. Also displayed is an example pair of ice surface and bed profiles (z, and zp; black and brown curves) that would produce this hydraulic gradient. In this example the bed slope Yb = 0.01. Table 2. Model parameters and physical constants Fig. 2. Synthetic air temperature and lake input time series defined by Eqns (6) and (7). Fig. 3. Simulated flood cycles. Results from three simulations using seasonally varying air temperatures to drive lake input, with (a, b) Tm = 10°C, (c, d) Ty = 12.7°C and (e, f) Tm = 15°C. (a, c, e) Time series of discharge at the lake outlet Q(0,t) (blue curves) and lake input Qi() (green curves) for 20 model years after transients have ended. (b, d, f) Solution orbits in Q(0,t)-h.(t) phase space from the complete 120-year-long simulations. Fig. 4. T,, parameter space. (a) Peak discharge of floods recorded during 120 year simulations plotted against the value of the midsummer air temperature 7 used in each simulation. The first ten flood peaks of each simulation are omitted to remove transients. The black box indicates the region shown in more detail in (b). Horizontal quantization visible in (a) (7 > 16°C) is the result of a coarse search through this region. The vertical lines in (b) correspond to orbits plotted in Figure 9. Lines that are one point in vertical extent (e.g. 6 < Tm, < 11°C and T,>19°C) correspond to limit cycles, where all the recorded floods have the same peak discharge; slightly thicker lines (e.g. 14 < T,, < 17°C) correspond to simulations during which transients lasted longer than ten flood cycles; and large blocks of points indicate dense, chaotic orbits. The colour of each point indicates the day of the year on which the flood peak occurs. In (a) the green arrows indicate two regions where floods occur progressively later in the year as 7, decreases. Fig. 5. Sensitivity to initial conditions. (a) Time series of lake depth from two simulations which used T,, = 12.7°C. Simulations are identical except for a slight difference in their initial lake depth h, (0). The blue curve corresponds to h,(0) = 40.00 m, and the red curve corresponds to h,(0) = 40.01 m. (b, c) Corresponding orbits in Q(0,t)—h_(0)—S(0,t) phase space, with the initial and final conditions shown in green and rec respectively. (d) Time series of the normalized separation between the orbits in (b, c). Normalized separation is calculated by differencing each component of the orbits’ positions normalized with the maximum value of the corresponding variable. The Pythagorean sum of the normalized separation is plotted. Note the different time axes in (a) and (d). Fig. 7. Orbits in Q(0,t)-h.(t) phase space of solutions with (a) Tm = 13.888000°C and (b) T,,=13.888002°C, which lie on each side of an abrupt bifurcation in 7,, parameter space (Fig. 4) for 20<t<90 years. Points are colored to indicate simulation time t and are separated by ~4.5 model days. When it is forced with constant meltwater inputs, Nye’s (1976) jékulhlaup model simulates floods whose magni- tudes increase with the input to the lake; the faster a lake is supplied with water, the larger the floods it produces (Clarke, 1982; Fowler, 1999; Ng and others, 2007; Kings- lake, 2013). We have demonstrated a more complex relationship between meltwater input and flood size by prescribing a more realistic, time-varying lake input. Inter- actions between the timing of lake drainage and a seasonal Fig. 6. The Nye attractor and topological mixing of phase space. (a) A solution orbit in Q(0,t)-h.(0—-t phase space of a simulation that used Tm = 12.7°C and h,(t=0)=40m. The blue curve’s distance from the vertical green line denotes the lake depth /,(0), its position along the vertical axis denotes the discharge at the lake Q(0,t) and its rotation around the green line denotes the fractional part of time t multiplied by 2m (one rotation corresponds to 1 year). (b) Poincaré sections taken perpendicular to the t-dimension at nine different stages of the year, denoted by the day of the year d. The points locate the intersection of the orbit with each section. Each section’s location in Q(0,t)—-h,()-t space is indicated by the boxes in (a). The colour of the points indicates simulation time t. Fig. 8. Time series of lake input Q; and discharge through the lake outlet Q(0,t) from two simulations that lie close to the abrupt bifurcation at T,+13.88°C. (a) Initial Q;,(t) and QV(O,t) time series from two simulations that used T,,=13.888000°C and Tm = 13.888002°C. (b, c) Long-time time series of the same variables from simulations that used 7,,=13.888000°C and Tm = 13.888002°C respectively. Fig. 9. Period-doubling bifurcations. Long-time solution orbits (transients have been removed) in Q(0,t)—h,(t) phase space of four simulations that used (a) Tm=13.701°C, (b) Tm=13.401°C, (cC) Tm= 13.201°C and (d) Ty, = 13.151°C. (a-d) correspond respect- ively to the positions in T,,-parameter space indicated in Figure 4b by the vertical lines A-D.
