Long-period variability in ice-dammed glacier outburst floods due to evolving catchment geometry

We combine a glacier outburst flood model with a glacier flow model to investigate decadal to centennial variations in outburst floods originating from ice-dammed marginal basins. Marginal basins ::: can form due to ::: the retreat and detachment of tributary glaciers, a process that often results in remnant ice being left behind. The remnant ice, which can act like an ice shelf or break apart into a pack of icebergs, limits the basin : a :::::: basin’s ::::: water : storage capacity but also exerts pressure on the underlying water and promotes drainage. We find that during glacier retreat there is a strong, nearly linear relationship between 5 flood water volume and peak discharge for individual basins, despite large changes in glacier and remnant ice volumes that are expected to impact flood hydrographs. Consequently, peak discharge increases over time as long as there is remnant ice remaining in a basin , the :: and : peak discharge begins to decrease once a basin becomes ice free, and : . :::: Thus : similar size outburst floods can occur for very different glacier volumes : at :::: very :::::::: different ::::: stages ::: of :::::: glacier :::::: retreat. We also find that the temporal variability in outburst flood magnitude depends on how the floods initiate. Basins that connect to the subglacial hydrological 10 system only after reaching flotation ::::: depth yield greater long-term variability in outburst floods than basins that are continuously connected to the subglacial hydrological system (and therefore release floods that initiate before reaching flotation :::: depth). Our results highlight the importance of improving our understanding of both changes in basin geometry and outburst flood initiation mechanisms in order to better assess outburst flood hazards and :::: their impacts on landscape and ecosystem evolution.

Outburst flood theory dictates that flood characteristics, such as event timing and peak discharge, depend on glacier and basin geometry, both of which evolve as glaciers advance or retreat. Consequently, outburst floods can be viewed as semiperiodic disturbances to glaciated landscapes that switch on/off and evolve in response to climate change. We are motivated by a desire to understand the evolving hazard of outburst floods as well as the impacts of these extreme events on landscape and 55 ecosystem evolution. ::::: Thus, ::: our ::::: work ::::::::::: complements :::::: efforts :: to ::::::::: understand ::::::::: long-term :::::::: variations :: in :::::: glacier :::::: runoff :::::: during :::::: glacial ::::::: recession ::::::::::::::::::::::::::::::::::::::: (e.g., Milner et al., 2017;Huss and Hock, 2018) : . In situ observations of outburst floods from individual glaciers over multiple years or decades are limited to a few sites. Due to a lack of observations, no previous work has tried to develop a theoretical understanding of the impact that glacier retreat has on outburst flood hydrographs. We address this problem with a one-way coupled glacier-basin-outburst flood model and focus on quantifying the long-period variability in outburst floods 60 that arise due to changes in catchment geometry. Our primary objective is to investigate changes in outburst flood hydrographs as a glacier retreats by exploring different basin geometries and flood onset mechanisms. In addition we account for remnant ice left behind in a basin, which reduces basin storage capacity ::: the :::::: storage :::::::: capacity :: of ::::: water :: in ::: the ::::: basin : but also acts like a gravity piston that pushes water out of a basin. We do not attempt to address the significant year-to-year variability in outburst flood hydrographs that has been observed at some glaciers (e.g. Huss et al., 2007;Neal, 2007;Kienholz et al., 2020) At peak storage capacity Model schematic illustrating the glacier and basin geometry (for a wedge-shaped basin).

Model description
We build on the outburst flood modeling work of Nye (1976), Fowler (1999, and Kingslake (2013) by accounting for changes in glacier and basin geometry (Fig. 2), both of which are expected to affect the magnitude and duration of outburst floods.

70
We first use an idealized glacier flow model to quantify changes in glacier geometry, ice dam thickness, and the amount of remnant, floating ice in a basin as a glacier retreats. For each year of the glacier flow model we extract the glacier geometry and remnant ice volume, which we then feed into the glacier outburst flood model. In the following subsections we describe the outburst flood model, the hypsometry and evolution of the marginal basin, and the glacier flow model. A list of model variables is included in Table 1. The outburst flood model consists of four coupled equations that conserve mass, momentum, and energy as water flows from a marginal basin and through a semi-circular conduit to the glacier terminus, assumed to be open to the atmosphere (Nye, 1976;Fowler, 1999). The ice dam seal is assumed to be located immediately adjacent to the basin. The cross-sectional area of the 80 conduit, S, evolves by melting and creep closure, and consequently discharge Q, effective pressure N (ice-overburden minus water pressure), and melt rateṁ (expressed as mass per unit length per unit time) vary temporally and spatially. We define the densities of ice and water as ρ i = 917 kg/m 3 :::: m −3 and ρ w = 1000 kg/m 3 ::::: m −3 , gravitational acceleration as g, and the latent heat of fusion as L f = 3.34×10 5 J/kg :::: kg −1 : (Cuffey and Paterson, 2010). Following Fowler (1999), we use the basic hydraulic gradient ψ = ρ w g sin θ − ∂P i /∂s, where θ is the conduit slope (assumed to equal the bed slope), P i = ρ i gH is the ice pressure,

85
H is the glacier thickness, and s is the along-flow coordinate parallel to the bed. The conduit length L, glacier thickness, and glacier thickness gradient evolve as the glacier thins and retreats (Section 2.2).
The assumption that the channel walls enlarge by melt and shrink due to creep closure results in an expression for the rate of change of conduit area given by

95
where M represents additional water flux supplied to the conduit per unit length. We prescribe a small value of M = 10 −5 m 2 s −1 to ensure that the conduit always remains open :::::::::::: (Fowler, 1999). We use Manning's equation to describe conservation of momentum, yielding an expression relating the discharge and conduit area to the basic hydraulic gradient and effective pressure,
Finally, conservation of energy requires that Two boundary conditions are required to solve this system of equations. We set the effective pressure at the terminus equal to 0. At the basin outlet, the effective pressure is where H b and h w are the glacier thickness and water depth at the basin outlet and h i is the thickness of floating ice in the basin.

115
where A b is the mapview area of the basin at different elevations, z b is the elevation relative to the lowest point in the basin, and a and p are constants that describe the basin shape. For reference, p = 1, p = 2, and p = 3 describe box-, wedge-, and semicircular-cone-shaped basins, respectively. We define W b , L b , and θ b as the basin width, length, and bed slope (Fig. 2).
For a box-shaped basin a = W b L b , for a wedge-shaped basin a = W b cot θ b , and for a semicircular-cone-shaped basin a = (π/2) cot 2 θ b .

120
The basin is assumed to be completely filled with ice at year 0, at which point the tributary glacier detaches from the trunk glacier and leaves behind remnant ice. Initially the remnant ice may be attached to the trunk glacier and act like a floating ice tongue, but ultimately it breaks into a pack of icebergs. We assume that the remnant ice thins at a rate given by the specific surface mass balance rate. Thus we neglect replenishment of ice into the basin via glacier flow or iceberg calving. We further assume that the remnant ice is sufficiently mobile and fractured to form a horizontal layer of thickness h i as the basin fills. We 125 therefore assume that drainage proceeds quickly enough that the floating ice thickness does not change during the course of the outburst flood and consequently ice is stranded on the basin walls (see Fig. 2). The floating ice volume at time t is given by where subscript , 0 refers to initial conditions,Ḃ b is the specific surface mass balance rate (see Section 2.2) at the basin's elevation, and we apply the mass balance rate to the surface of the remnant ice.
The simulations are run until the glacier terminus retreats past the basin, which is initially located 75% of the way from the head of the glacier to its terminus. We ran additional simulations with different parameter values for bed slope, climate, and basin location. Although these parameters affect the details of how outburst floods change from year to year, they do not affect the overall pattern of how outburst floods evolve.

180
The glacier flow model is based on conservation of momentum, which requires that the glaciological driving stress is balanced by gradients in longitudinal stress, lateral drag, and basal drag (van der Veen, 2013), such that where ν is the depth-and width-averaged viscosity, U is the depth-and width-averaged velocity, W is glacier width, τ is the basal shear stress, and h s is the glacier surface elevation. The viscosity depends on the strain rate according to Glen's Flow 185 Law: We assume a simplified ad hoc parameterization of the basal shear stress, in which τ = τ max (U/ max(U )) with τ max = 2.5 × 10 5 Pa. This parameterization results in shear stresses on the order of 10 5 Pa, which are typical values for valley glaciers (e.g., Braedstrup et al., 2016), and produces realistic glacier geometries and velocities across a wide range of bed slopes and 190 climates. Importantly, the parameterization ensures that the resistive stresses never exceed the glaciological driving stress.
where ELA is the equilibrium line altitude. We use an initial ELA of 1500 m, a balance gradient of dḂ/dz = 0.005 a −1 , and 200 a maximum balance rate ofḂ max = 2 m a −1 . The ELA increases at a rate of 5 m a −1 to approximate expected changes under climate warming scenarios (Huss and Hock, 2015).

Simulations
2.3 :::::::::: Simulations In the glacier flow model it takes about 300 years for the glacier terminus to retreat from its initial position to the location of the 210 marginal basin, a distance of ∼12 km . : 5 ::: km :::: (Fig. ::: 3). For each year of the glacier model output, we extract (i) the distance from the basin to the terminus, which we take to equal the conduit length, (ii) the glacier thickness profile and ice dam thickness, and (iii) the specific mass balance rate of the ice dam. Then (i) and (ii) are fed directly into the outburst flood model and (iii) is used to calculate the volume of floating ice remaining in the basin (Fig. 4).
To demonstrate how remnant ice affects outburst floods, we first run simulations in which we use the glacier geometry from 215 one time step in the glacier flow model, assume a box-shaped basin, and run the outburst flood model with different starting water volumes. We run the simulations both without ice and with enough ice to force the ice dam to be at :: ice :::: dam : flotation. between ::: the ::::: glacier :::: flow, :::: basin :::::::: evolution, ::: and :::::: outburst :::: flood :::::: models.
Thus these initial simulations are similar to those that we run in the flotation scenario (next paragraph) except that here the basin is not necessarily at flotation ::::: depth unless it contains remnant ice.
We then use the evolving glacier and basin geometries to model long-period variations in outburst floods using the flotation 220 and filling scenarios described in Section 2.1.1. In the flotation scenario, we assume that the initial water pressure at the basin outlet equals the overburden pressure of the ice dam. In this scenario, the initial conduit area is 1 m 2 . Thus we assume that the basin is not connected to the subglacial drainage system until the onset of the outburst flood. To test the effect of basin geometry and floating ice on outburst flood evolution, we run simulations with (i) box-shaped, wedge-shaped, and semicircular-coneshaped basins and (ii) both with and without floating ice . :::: (Fig. ::: 5). For the box-shaped basin we used a value of a = 8.5×10 5 m 2 , 225 for the wedge-shaped basin we used a basin width of 1910 m and a bed slope of 15 • , and for the semicircular-cone-shaped basin we used a bed slope of 10.6 • . These values were chosen so that the basins would initially have the same :::: basin : storage capacity V s , which we define as the water volume when the ice dam is at flotation, if they were ice free.

Results
For glaciers with a fixed geometry, floating ice in a basin causes outburst floods to have higher peak discharge and shorter duration than might otherwise be expected based solely on the consideration of flood water volume (Fig. 6)   simulations started with the basins filled to flotation ::::: depth (Fig. 7) or if the basins were connected to the subglacial hydrological system as they filled (Fig. 8). The floods that occur in the years immediately following the formation of a marginal basin have low peak discharge on account of the basin's small storage capacity. As the climate warms, the remnant ice thins more quickly than the ice dam, which is partially replenished by the delivery of ice from upstream. The largest outburst floods occur when the basin becomes ice free, after which the peak discharge decreases until the basin is no longer dammed by the trunk glacier.

245
For the simulations in which the basins were filled to flotation ::::: depth before flood onset, we considered three different basin hypsometries (semicircular-cone-, wedge-, and box-shaped) that had identical storage capacities at the time of basin formation (year 0). The cone-shaped basin produced the largest outburst floods in terms of peak discharge, although the number of years with large outburst floods was less than for the ::::::: however :::: flood ::::::::: magnitude ::::::::: decreased :::: more :::::: rapidly :::::: across ::: the :::::::::: simulations :: for ::: the ::::::::::: cone-shaped ::::: basin :::::::: compared ::: to ::: the wedge-or ::: and : box-shaped basins (Fig. 7a-c). The differences in peak discharge 250 and duration of large magnitude floods arise because, owing to their hypsometry, cone-shaped basins lose their floating ice more rapidly than wedge-or box-shaped basins and because as they drain the floating ice in the basin exerts pressure on the underlying water that helps to drive the water out of the basin. However, ice-dam thinning reduces basin capacity faster in cone-shaped basins than in wedge-or box-shaped basins and therefore the period of large outburst floods tends to be shorter :::::::: magnitude ::: of ::: the ::::::: outburst :::::: floods : in cone-shaped basins ::::::: decrease ::::: more :::::: rapidly. In the early years of the simulations, flood 255 magnitude increases in basins that are initially filled with ice (solid line) until the basin is ice free, whereas in basins that are initially ice-free (dotted line) the flood magnitude always decreases (Fig. 7d-f). For all three hypsometries we observe a nearly linear relationship between peak discharge and storage capacity , where storage capacity is defined as the water volume when the ice is at flotation (Fig. 7g-i). This relationship holds regardless of whether the basin contains ice or is ice-free.
For the simulations in which the basin is initially drained of water but remains connected to the subglacial hydrological 260 system as the basin fills, the relationship between peak discharge and peak water volume (which is often less than the storage capacity, as defined above) takes a slightly different form. First, the basin often does not reach flotation :::: depth in our simulations because the conduit enlarges at the same time as the basin is filling and consequently the outburst floods tend to be smaller in magnitude (Fig. 8). This behavior is sensitive to the model parameters though, as the basin could be made to reach or even exceed flotation ::::: depth by selecting a larger influx Q in . Second, there is a more prominent spike in the peak discharge curve 265 that occurs as the remnant ice is about to melt away completely (Fig. 8b). Similar to the flotation scenario, the relationship between peak discharge and peak water volume is approximately linear; however, in the filling scenario the peak water volume is less than the storage capacity because the basin does not completely fill (Fig. 8c). As a result the relationship between peak discharge and (total) storage capacity is not linear. During decadal to centennial scale glacier retreat, the peak discharge of ice-dammed outburst floods will tend to increase with time as long as there is remnant ice in a basin that is melting away. The peak discharge will begin to decrease only once the remnant ice is gone. This result is independent of basin geometry and the mechanism of drainage onset and is ultimately a consequence of the proportionality between peak discharge and basin storage capacity that occurs for individual basins despite large changes in glacier geometry and remnant ice. In other words, the model exhibits very little hysteresis between 285 peak discharge and :::: basin : storage capacity (Figs. 7g-i and 8c). :: As :: a ::::: result, ::: the :::: time :::: rate :::::: change :: of ::: the ::::: basin :::::: storage :::::::: capacity ::::::::: illuminates :::: how :::: peak ::::::::: discharge :::::: evolves ::::: with :::: time. : The storage capacity is found by inserting Equation 11 into Equation 8, which gives  Taking the derivative of Equation 17 with respect to time, we find that storage capacity evolves according to The term in square brackets in Equation 18 is always positive, and thus the storage capacity will always increase as long as dH b /dt > dh i /dt (i.e., the ice dam is thinning less quickly than the remnant ice). Thinning of the ice dam due to surface melting is partially offset by ice flow from upglacier and therefore the storage capacity, and by extension the peak discharge of outburst floods, will continue to increase until a basin is ice free.

295
However, in our simulations we did not account for ice flow or calving of icebergs into a basin, which would require a significantly more sophisticated ice flow model. Ice flow into a basin shortens the basin and reduces the storage capacity.
Calving changes the basin geometry but tends to have little net impact on storage capacity because it has two competing effects: it results in retreat of an ice dam away from a basin, which increases storage capacity, but it also adds to the volume of remnant ice, which reduces storage capacity. This is easiest to see for a box-shaped basin, for which the storage capacity is Since we are now allowing for :: For :: a ::::: more ::::::: detailed ::::::::: discussion ::: on ::: the ::::::: impacts ::: of ice flow and calving , the basin length is no longer treated as a constant and the remnant ice thickness varies in response to the addition of new icebergs and compaction/extension due to changes in the location of the ice dam. The rate of change of the storage capacity is The thickness of the remnant ice changes at a rate that is given by where U c is the calving rate. The three terms on the right-hand side of Equation A3 describe the changes in ice thickness due to the surface mass balance, the influx of freshly calved ice, and changes in the ice dam location. The rate of change of the basin length is simply dL b /dt = U c − U b , where U b is the rate at which ice is flowing toward the basin. By inserting these 310 expressions for dh i /dt and dL b /dt into Equation A2 and rearranging, we find that which indicates that the storage capacitywill increase as long as For a box-shaped basin the effects of calving cancel out completely and changes in storage capacity are only due to thinning of the ice dam, the surface mass balance rate, and the ice flux toward the basin. The effect of ice flow is to reduce the maximum storage capacity that occurs in 315 a basin and to increase the time that it takes for the maximum storage capacity to be reached since contraction of the remnant ice reduces the surface area that is susceptible to melting. :: on ::::::: storage ::::::: capacity, ::: see ::::::::: Appendix :: A. : Equation A4 illustrates that ice flow toward a basin may have important consequences for basin storage capacity. For example, the remnant ice in Suicide Basin, the source of recent outburst floods at Mendenhall Glacier, has a surface mass balance flux (Ḃ b L b W b ) of about −2.5 × 10 6 m 3 a −1 and the ice flux toward the basin (U b H b W b ) is roughly 3.5-7.0×10 5 m 3 a −1 320 (both expressed as ice equivalent) (Kienholz et al., 2020); thus ice flow is currently offsetting the growth in storage capacity due to melting by about 25%. Our analysis here has focused solely on basin storage capacity. The relationship between peak discharge and water volume is likely to become more complicated than presented in Figures 7 and 8 when changes in basin geometry due to ice flow and calving are accounted for. Moreover, we do not account for lateral variations in glacier thickness that may cause the seal of the ice dam to be located some distance from the basin. Flow re-direction toward a marginal basin 325 due to lateral surface gradients will affect the ice dam thinning rates and location of the seal in ways that we are unable to capture in our one-dimensional flowline model. These additional complexities should be considered in more detail in future studies.

Comparison to the Clague-Mathews relationship
Observations across a range of systems are suggestive of a power-law relationship between the peak discharge and total water 330 volume drained, ∆V w , during outburst floods (Clague and Mathews, 1973;Walder and Costa, 1996): This relationship is commonly referred to as the Clague-Mathews relationship. Ng and Björnsson (2003) examined the Clague-Mathews relationship by analyzing the equations describing flood evolution. Using a simplified version of the outburst flood model used in this study, they demonstrated that for :::::::: individual : basins that do not drain completely, (i) each flood trajectory 335 has a unique set of initial and final water levels and peak discharge, (ii) as a result peak discharge monotonically increases with water volume, and (iii) for individual basins there is a power-law relationship between discharge and water volume for floods. They focused on analyzing basins that experience incomplete drainage because some information on flood mechanics is lost if a basin drains completely. Their analysis predicts an exponent in the power-law relationship of about 1-2 for individual basins, depending on basin geometry and ice coverage. When observed flood data from multiple glaciers was scaled and placed 340 into their theoretical framework, they arrived at an exponent close to 1. They hypothesized that the difference between their theoretical exponent and the exponent in the Clague-Mathews relationship is due to confounding factors such as differences in flood initiation, basin geometry, and complete drainage.
Our simulations extend the work of Ng and Björnsson (2003). We modeled variations in outburst floods over decadal to centennial timescales, from different shaped basins, and with different drainage scenarios (flotation vs. filling). In addition, in 345 our simulations the basins always drained completely. We observe that the relationship between peak discharge and volume drained ::: peak ::::: water ::::::: volume ::::::: reached ::::: (equal ::: to :::::: volume :::::::: drained) is nearly linear in the flotation scenario (power law exponent of ∼1; Fig. 7g-i) and superlinear in the filling scenario (power law exponent >1; Fig. 8c). These trends occur regardless of whether the basins contain remnant ice, in which case peak discharge and storage capacity increase with time, or are ice free.
We also find that the slopes of the discharge-volume curves depend on basin geometry, where basins that contain less volume 350 near their outlets produce a steeper slope. This is likely a result of cone-and wedge-shaped basins being able to maintain high water pressures as they drain, thus favoring more rapid conduit growth.
The explanation for the lower, 2/3 exponent in the Clague-Mathews relationship remains elusive. Ng and Björnsson (2003) suggest that the lower exponent is due to differences in flood initiation across different basins (implying that flood initiation may depend on basin hypsometry). Flood initiation could also depend on some time-varying property such as ice dam thickness 355 or the size of a previous year's flood, both of which could influence the state of the subglacial hydrological system at the onset of a flood (see also Kingslake, 2015). One possibility is that floods :: For :::::::: example, :: in ::::: some :::: flood :::::: events involving large volumes of water :::: (i.e., ::::: when :: the ::::: basin :::::: storage :::::::: capacity :: is ::::: large) have a persistent impact on the subglacial system, so that when a basin refills it does so while slowly draining, whereas floods involving small volumes of water may have a less persistent impact and as a result subsequent floods will only initiate after the basin reaches flotation ::::: depth.
A second consequence of floating ice is that it affects the duration of outburst floods (Fig. 11). The role of floating ice 380 is again most clear in the flotation scenario. Early in the simulations, when the storage capacity is small, outburst floods can occur that have higher peak discharge than might be expected because the pressure from the floating ice helps to drive water out of the basin. However, the small amounts of water (relative to the size of the basin) in these events are not able to melt the conduit walls as rapidly as later floods and the overburden pressure from the ice dam, which favors creep closure, is high. Consequently the floods tend to be slower building and the basins may take about a week to drain. Later, when the  Figure 11. Comparison of peak discharge and time ::: from :::: lake ::::::: highstand : to peak discharge for the box-shaped basin under the (a) flotation and (b) filling scenarios.
Our simulations predict differences in how outburst floods will evolve with time, depending on whether a basin begins to 390 drain once it has filled or if the basin remains connected to the subglacial hydrological system and begins to drain as it is filling. Furthermore, our model does not address the large year-to-year variability in peak discharge and total water volume of outburst floods, which may vary by a factor of two or more in subsequent years (e.g., Huss et al., 2007;Neal, 2007;Kienholz et al., 2020) and is also likely related to the onset mechanism. A deeper understanding of the onset of outburst floods is therefore critical to improving our ability to assess both the short-and long-term risk associated with outburst floods.

Conclusions
We modeled the effect of changes in glacier and basin geometries on the magnitude and duration of ice-dammed glacier outburst floods. In our simulations we accounted for remnant, floating ice that is left behind in marginal basins during the retreat of tributary glaciers. The remnant ice acts to exerts pressure on the underlying water and ::: thus : helps to increase discharge by enlarging the subglacial conduit. Because the remnant ice is not replenished by ice flow from upglacier, it thins more quickly 400 than the adjacent ice dam. As a result the basin storage capacity increases with time as long as the basin contains remnant ice, regardless of basin hypsometry. Despite complex relationships between glacier and basin hypsometry, remnant ice thickness, and discharge, we find nearly linear relationships between peak outburst flood discharge and total water volume for individual basins. This is regardless of whether the basin does not start to drain until it is full ::::: drains :::: once :: it :::::: reaches :::::::: flotation ::::: depth or if the two different drainage scenarios that we considered highlight the importance of improving our understanding of drainage onset. Basins that are continuously connected to the subglacial drainage system tend to produce similar outburst floods from one year to the next, except during the years immediately before and after the loss of remnant ice, whereas basins that do not begin to drain until full produce much larger variability in the magnitude and timing of outburst floods.
In our simulations we made a number of simplifying assumptions in order to garner a fundamental understanding of the 410 long-period variability of outburst floods in an evolving catchment. In particular, we :: (1) assumed that the seal of the ice dam was immediately adjacent to the basin and did not account for changes in the hydraulic potential gradient that could drive water from the glacier into the basin as it is filling, :: (2) : treated remnant ice as fluid that spreads out as a basin fills, instead of accounting for the granular nature of the icebergs, ::: (3) did not consider the state of the glacier's hydrological system at the time of drainage, which may impact flood evolution, or changes in ice flow due to the evolving subglacial hydrology, :: (4) : did not 415 allow for ice flow into the basin from the trunk glacier, and :: (5) : did not account for interannual variability in climate and its affects on glacier geometry and basin filling rates. Year-to-year variability in the timing, duration, and magnitude of outburst floods (e.g., Huss et al., 2007;Neal, 2007;Kienholz et al., 2020) may mask the longer period changes in outburst floods due to changes in glacier and basin geometry that we modeled here. Additional and more sophisticated modeling studies will be needed to elucidate the impact of these processes on the decadal and centennial evolution of outburst floods and to connect 420 outburst floods to landscape and ecosystem evolution.
Author contributions. JA and EH conceived the study, AJ ran the simulations with assistance from JA and JK, AJ and JA prepared the manuscript with contributions from JK and EH.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This project was supported by funding from the Alaska Climate Adaptation Science Center and the US National Science