An Arctic delta reduced-complexity model and its reproduction of key geomorphological structures

Arctic river deltas define the interface between the terrestrial Arctic and the Arctic Ocean. They are the site of sediment, nutrient, and soil organic carbon discharge to the Arctic Ocean. Arctic deltas are unique globally because they are underlain by permafrost and acted on by river and sea ice, and many are surrounded by a broad shallow ramp. Such ramps may buffer the delta from waves, but as the climate warms and permafrost thaws, the evolution of Arctic deltas will likely take a different course, with implications for both the local scale and the wider Arctic Ocean. One important way to understand and predict the evolution of Arctic deltas is through numerical models. Here we present ArcDelRCM.jl, an improved reduced-complexity model (RCM) of arctic delta evolution based on the DeltaRCM-Arctic model (Lauzon et al., 2019), which we have reconstructed in Julia language using published information. Unlike previous models, ArcDelRCM.jl is able to replicate the ramp around the delta. We have found that the delayed breakup of the so-called “bottom-fast ice” (i.e. ice that is in direct contact with the bed of the channel or the sea, also known as “bed-fast ice”) on and around the deltas is ultimately responsible for the appearance of the ramp feature in our models. However, changes made to the modelling of permafrost erosion and the protective effects of bottom-fast ice are also important contributors. Graph analyses of the delta network performed on ensemble runs show that deltas produced by ArcDelRCM.jl have more interconnected channels and contain less abandoned subnetworks. This may suggest a more even feeding of sediments to all sections of the delta shoreline, supporting ramp growth. Moreover, we showed that the morphodynamic processes during the summer months remain active enough to contribute significant sediment input to the growth and evolution of Arctic deltas and thus should not be neglected in simulations gauging the multi-year evolution of delta features. Finally, we tested a strong climate-warming scenario on the simulated deltas of ArcDelRCM.jl, with temperature, discharge, and ice conditions consistent with RCP7–8.5. We found that the ramp features degrade on the timescale of centuries and effectively disappear in under 1 millennium. Ocean processes, which are not included in these models, may further shorten the timescale. With the degradation of the ramps, any dissipative effects on wave energy they offered would also decrease. This could expose the sub-aerial parts of the deltas to increased coastal erosion, thus impacting permafrost degradation, nutrients, and carbon releases.

1 Introduction

Arctic deltas are key components of the permafrost landscape, connecting the permafrost areas upstream and the Arctic Ocean. They act as records and filters of the particulate and dissolved matter, such as sediments and nutrients, that originated from the Arctic and sub-Arctic regions – regions that contain a substantial portion of the Earth's soil organic carbon (Hugelius et al., 2014; Schuur et al., 2015) – which could potentially exacerbate climate warming through positive feedback. As permafrost thaws and the Arctic Ocean trends towards being free of ice, especially under amplified warming in polar regions (Stocker, 2014; Rantanen et al., 2022), Arctic deltas face an uncertain future. On a local level, the ecosystems surrounding Arctic deltas will also face significant impacts as a result of climate change (Pisaric et al., 2011). The largest Arctic delta is the Lena Delta, which has fan-shaped morphology with multiple channels leading discharge from its epicentre to the delta edge. Although complicated by neotectonics (Are and Reimnitz, 2000), its large-scale structure is probably determined by regional relative sea level effects on fluvial sediment budgets within the delta and at its edge (Whitehouse et al., 2007).

An important feature ubiquitous to Arctic deltas is the 2 m ramps, which we will interchangeably refer to as “ramp features” or simply “ramps” (Are and Reimnitz, 2000). Figure 1 shows such a ramp feature surrounding the Lena Delta, which dips gradually from below the sea surface to roughly 2 m depth but with localized variations on the order of 1 m (Fuchs et al., 2022). These ramps extend from the above-water shoreline of Arctic deltas over tens of kilometres towards the open ocean (Reimnitz, 2002). The shallow inclination of the ramp serves to diffuse wave energy offshore and may protect the delta from direct wave impact. The shallow water depth of the ramp corresponds to the range of maximum winter ice thickness (Wegner et al., 2017), meaning that ice freezes to the seabed at some point during winter, which is called “bottom-fast” or “bed-fast” ice (Dammann et al., 2018, 2019; Eicken et al., 2005), and that a seasonally frozen layer can develop in the sediment. When seasonal freezing exceeds thawing in subsequent years, permafrost can develop below the seabed (Osterkamp et al., 1989). Land- or bottom-fast ice and seasonal or permanent freezing may all act to stabilize the ramp (Overeem et al., 2022). Therefore, the ramp features, aside from being an integral part of Arctic deltas, may also play an important role in protecting Arctic shorelines from coastal erosion (Dean and Dalrymple, 2002) and could enhance carbon sequestration (Overeem et al., 2022). Moreover, the shallow-water platform provided by the ramp could play an important role in the surrounding ecosystems (Lopez et al., 2006). The origin, evolution, and stability of the ramp features are therefore an important feature of Arctic deltas. A better understanding of all three is required to predict how changing sea ice cover and air temperatures will impact Arctic delta morphology.

In order to explore the possible origin of these ramps, numerical models form an important tool in the toolbox. Due to the complexity of the system (which involves permafrost, flow on low-slope environments, ice cover, spring floods, and more), delta models are typically divided into two classes. The first is built on physically based equations, simplified to be computationally tractable (e.g. models involving Delft 3D; Lesser et al., 2004). They are able to simulate more directly delta dynamics but at a cost of computational resources and the ability to cover timescales of years or longer. To address these issues, the second class of models – reduced-complexity models (RCMs) – simulate phenomenological processes of arctic delta evolution using rule-based trajectories of cellular automata (works originating this method in this field include Murray and Paola, 1994, 2003; Murray, 2007). The rules governing the automaton units are typically informed by physical equations under specific sets of conditions and by empirical observations. Due to the simplified construct, there are greater flexibilities in choosing the spatial and temporal step sizes, resulting in much greater spatial and temporal coverage whilst keeping computational requirements feasible. As Bokulich (2013) pointed out, there are ongoing debates about the predictive power and explanatory insights offered by these classes of models, and their use by researchers is typically with “division of cognitive labour” – to use a term coined by Bokulich (2013) – in mind. Within this context, we take the RCM approach in order to explore what physical-process rules could be favourable to the formation of a ramp feature.

Figure 1 Map of the Lena Delta and the bathymetry of the southern Laptev Sea region (Fuchs et al., 2022). The orange-to-green relief shows the sub-aerial portion of the Lena Delta and its surrounding land. The blue-to-white colour scale shows the bathymetry. Dark-blue channels within the delta show deep channels that do not freeze in winter (Juhls et al., 2021). Light blue within the delta shows the maximum channel area during the spring flood. The hatched area displays the shallow-water 2 m ramp feature. Some deeply incised channels are visible within the ramp feature. The land area relief is visualized in a false colour Landsat-8 mosaic, courtesy of the United States Geological Survey, processed in Google Earth Engine. The red box in the inset shows the location on the world map; the inset's world-map background is from OpenStreetMap (© OpenStreetMap contributors 2022. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.).

One representative RCM for understanding deltas is the DeltaRCM (Liang et al., 2015b), the Arctic extension of which (called DeltaRCM-Arctic; Lauzon et al., 2019) will serve as the starting point of our work herein. DeltaRCM and its Arctic extension have been demonstrated to efficiently reproduce numerous observed features of natural and experimental deltas (Liang et al., 2015b, a; Lauzon et al., 2019; Piliouras et al., 2021). Our goal is to build upon previous work to reproduce the ramp features in order to explore their physical-process origin.

We start by rewriting the DeltaRCM base and reconstruct its Arctic extension in Julia, which has comparable performance to C and FORTRAN but retains the syntactical convenience of MATLAB and Python. We then explored the effects of modifying the rules to include additional processes and physics, with the hope of identifying and understanding the circumstances that favour the formation of the ramp feature. Moreover, motivated by the high flow rate during summer in large Arctic deltas, we also address the importance of summer months in Arctic delta evolution when examining growth and evolution over timescales of years or longer. We call the modified model ArcDelRCM.jl for succinctness and in keeping with the conventions of Julia code packages and to signify that it is a reconstruction of the Arctic extension of DeltaRCM based on published articles (Lauzon et al., 2019; Piliouras et al., 2021) and not a direct translation of the original DeltaRCM-Arctic source code (which is not publicly available).

The purpose of this article is two-fold. First, we present and motivate the physical-process rule changes we made to our reconstructed version of DeltaRCM-Arctic to arrive at the ramp-producing ArcDelRCM.jl version.¹ Second, we present the model outputs, including ones that are intended to simulate the evolution of large-scale Arctic deltas. Through these outputs, we demonstrate the model's capability for reproducing the 2 m ramp, identify the processes that led to their formation in the model, and make exemplary comparisons of these outputs with a real-world case, the Lena Delta. Aside from being the largest Arctic delta, the Lena Delta is chosen for the real-world case example because of its fit to the modelled geometry. Many of the other Arctic river mouths (e.g. Ob, Yenisei, Mackenzie) are estuarine or geologically confined and thus match the modelled geometry poorly, which may confound the analyses herein. Other Arctic deltas such as Olenyok and Colville are much smaller and therefore more likely to be influenced in their growth by processes other than those modelled here (we discuss some of these processes at the end of Sect. 4). Through the exemplary case of the Lena Delta, we also argue for the importance of the summer months in modelling delta evolution on timescales of years or longer. Finally, we gauge the possible fate of these ramps under a warming climate.

2 Methods

2.1 Description of DeltaRCM(-Arctic)

In this section, we provide an overview of DeltaRCM and its Arctic extension (referred to as “DeltaRCM-Arctic” by Piliouras et al., 2021) based entirely on the model-describing publications (Liang et al., 2015b, a; Lauzon et al., 2019; Piliouras et al., 2021), the source codes of the (non-Arctic) DeltaRCM (Liang, 2015; Perignon, 2018), and our observations during the process of reproducing these models prior to extending them. These descriptions will facilitate the introduction of all of our modifications, left to Sect. 2.2, and the specifications of our simulations in Sect. 3.

2.1.1 Simulation domain

In the most basic setup, the simulation domain (depicted in Figs. 1 and 2 of Liang et al., 2015b) consists of a rectangular grid of N_x by N_y cells (typically, N_y=2N_x). Each cell is a square with a width of δc. Along the y dimension, the first N_wall (typically three) cells in the x direction are defined as the inlet wall, which are impermeable and static. Centred around the $\left(\frac{\mathrm{1}}{\mathrm{2}}{N}_y\right)$th cell, there is an opening in the inlet wall N₀ cells wide from where the water (volume) discharge, Q_w, and the sediment (volume) discharge, Q_s, enter the simulation domain. This opening is the inlet channel.

The domain is initialized with a water-surface elevation of H (typically 0 m, taken to be the sea-surface height), a water depth of h_B (“B” for ocean basin), and a corresponding bed elevation of $\mathit{\eta }=H-h_{\mathrm{B}}$. Within the inlet channel, an initial surface slope, S₀, is added to mimic the backwater slope. An inlet flow depth, h₀, is given and reimposed at the start of each time step, such that the inlet flow speed, u₀, can be determined in conjunction with Q_w and N₀×δc.

2.1.2 Flow field

To build the flow field within each time step, the input discharge is divided into n_w (which is 2000 in all examples) packets and sent through the simulation domain using a weighted random-walk scheme. The weights of each of the eight neighbouring grid cells are determined by a linear combination of two factors: (i) water-surface gradient (as a proxy for gravity),

$$\begin{array}{}\text{(1)}& w_{i,\mathrm{surface}}=\max \left(\mathrm{0},-\mathrm{\nabla }{H}_i\right),\end{array}$$

where ∇H_i is the gradient of water-surface elevation from the current cell towards its ith neighbouring cell; (ii) flow depth, expressed as a resistance measure (as a proxy for inertia), scaled by both the projection of the local flow field to the eight neighbouring cells and by the distances to those neighbouring cells,

$$\begin{array}{}\text{(2)}& w_{i,\mathrm{inertia}}={\displaystyle \frac{\mathrm{1}}{R_i}}{\displaystyle \frac{\max \left(\mathrm{0},{\mathit{q}}_{\mathrm{w}}\cdot {\mathit{d}}_i\right)}{{\mathrm{\Delta }}_i}},\end{array}$$

where $R_i^{-\mathrm{1}}=(h_i-h_{\mathrm{ice},i})\left(\mathrm{1}-\frac{h_{\mathrm{ice},i}}{h_i}\right)$ is the flow resistance measure of the ith neighbouring cell, taking into account the full water depth, h_i, and the portion of it that is ice, h_ice,i (Lauzon et al., 2019). q_w is the unit-discharge vector at the current cell (serving as the flow-direction vector), d_i is the unit vector pointing from the current cell towards the ith neighbouring cell, and Δ_i is the distance between the centres of the current cell and that of its ith neighbour. Cells with water depth shallower than 0.1h₀ (up to a maximum of 0.1 m) are classified as “dry” and thus assigned a weight of 0. Discharge is conserved throughout, with each packet's contribution (q_w) to each visited grid cell recorded in an additive manner. The total random-walk weight is combined through a “partitioning coefficient”, γ:

$$\begin{array}{}\text{(3)}& w_i=\mathit{\gamma }{w}_{i,\mathrm{surface}}+(\mathrm{1}-\mathit{\gamma })w_{i,\mathrm{inertia}}.\end{array}$$

Liang et al. (2015b, a) described γ as a free parameter that controls the lateral spread of water. By increasing the importance of the water-surface gradient, its cross-channel component is also emphasized. The value of γ is typically small (e.g. 0.05 according to Liang et al., 2015b) and is also given by Liang et al. (2015a) as

$$\begin{array}{}\text{(4)}& \mathit{\gamma }={\displaystyle \frac{g\mathit{\delta }c{S}_{\mathrm{0}}}{u_{\mathrm{0}}^{\mathrm{2}}}}\end{array}$$

as a guideline to choose an appropriate value, where g is the gravitational acceleration. The latter expression may have arisen from taking the ratio between the pressure gradient and inertia terms (without local acceleration) of the shallow-water equations. The latter expression is the version for γ implemented in the source codes of DeltaRCM (Liang, 2015; Perignon, 2018), which would result in, for example, a value of 0.098 in the demonstration cases with a 50 % sand fraction² instead of 0.05 in Liang et al. (2015b). We note that this remains a free parameter that has been given a range of values (e.g. 0.02 to 0.15) in various tests of its influence on delta planforms (Liang et al., 2015b, Supplement). In our own experiments mimicking the size and conditions of the Lena Delta, with much larger δc, we have found that γ=0.135 (±0.02) works best in producing the planform of deltas similar to the Lena Delta.

The slope S₀ is imposed along each water packet's path, so that an averaged, approximate backwater slope is formed in H in each time step. The details of this procedure can be found in Liang et al. (2015b, Sect. 3.2.2).³

Given the bed-elevation field at the time, η (either from initialization or from the previous time step), the water depth is determined from $h=H-\mathit{\eta }$. Using the accumulated q_w from all passing water packets and the flow depth h (minus any portion that is h_ice), the flow speed u is calculated. Finally, the direction of flow at each grid cell is the average entry and exit directions of all water packets that passed through that cell (Fig. 4 of Liang et al., 2015b). We note that, in practice, the unit-discharge vector field serves as the flow-direction field and is thus computed during the routing of the water packets. Its under-relaxation by the free parameter γ is done directly when computing the routing weights (Eq. 3).

The flow-field determination process is iterated multiple (n_iter) times per time step to suppress any instability due to the randomized nature of the scheme (Liang et al., 2015b). Note that there is one more step that is not documented in the original article but exists in the non-Arctic source codes (Liang, 2015; Perignon, 2018): the flow field q_w also undergoes an “under-relaxation” identical to that undergone by H (Liang et al., 2015b, Eq. 12), except that it is with a different coefficient and is applied between each iteration instead of merely across time steps. The under-relaxation coefficient for the q_w field is 0.9 for the first iteration (in each time step) and $\mathrm{2}/n_{\mathrm{iter}}$ for all subsequent iterations.

2.1.3 Ice cover

Lauzon et al. (2019) extended DeltaRCM to simulate Arctic deltas. They simulate only the spring-flood period (assumed to be 10 d). At the beginning of each flood, the maximum ice cover is defined with a pair of parameters: maximum ice thickness, h_ice,max, and maximum ice extent, f_ice, which is a fraction between 0 and 1. At its maximum, ice thickness is h_ice,max everywhere except on dry grid cells and except in cells within (1−f_ice) of the mean distance between the inlet channel and the (average) coast of the delta. To improve simulation stability, a taper from 0 to h_ice,max over the equivalent of 10δc is imposed from the edge of the ice-free zone towards the ocean. For similar purposes, h_ice may not be equal to the local water depth (e.g. h_ice<0.9999 h).

Two processes contribute additively to the melting of the ice: discharge-based heat flux and atmospheric heat. The former is proportional to the flow-speed field, u, of the current time step and is given by an expression from Searcy et al. (1996, Eq. 2 and the subsequent unnumbered equation) (also in Lauzon et al., 2019, Supplement, and Piliouras et al., 2021). The latter (i.e. atmospheric contribution) is given by

$$\begin{array}{}\text{(5)}& {\left({\displaystyle \frac{\mathrm{d}{h}_{\mathrm{ice}}}{\mathrm{d}t}}\right)}_{\mathrm{atmospheric}}=-a{\displaystyle \frac{h_{\mathrm{ice},\max }}{t_{\mathrm{melt}}}},\end{array}$$

where t_melt is the period during which the entire ice thickness would melt due to atmospheric heat alone and a is a scaling factor, between 0 and 1, to tune the contribution of atmospherically induced melting towards the total melt. In the model result of Searcy et al. (1996), a≈0.58, and in Piliouras et al. (2021), a=0.5. In DeltaRCM-Arctic, t_melt is taken to be 10 d (i.e. the entire simulation period for each model year) (Lauzon et al., 2019; Piliouras et al., 2021).

2.1.4 Sediments and permafrost

The sediment entering the delta in DeltaRCM(-Arctic) is split into two categories: sand (representing bed load) and mud (representing suspended load). Their relative fraction of the total is given by the “sand-fraction” parameter, f_sand. The routing of sediment packets is similar to the routing of water packets and is described in Liang et al. (2015b, Sect. 3.2.4).

Deposition and erosion of sediments are determined by a set of threshold relations (detailed in Liang et al., 2015b, Sect. 3.2.5). Deposited sediments are tracked and stored in a grid with the same lateral spatial coverage as the simulation domain (N_x by N_y) and with a depth-wise dimension of δz per cell. After all sediment packets have passed through during one time step, any bed-elevation gains are added to the sediment grid. Each of these added grid elements records the sand fraction (i.e. the relative fraction of sand in the total volume of sand and mud) of all sediments deposited by all passing packets. Grid cells corresponding to eroded sediments (i.e. the volume picked up by passing packets) are simply removed.

In DeltaRCM-Arctic, Lauzon et al. (2019) assume a constant active layer (or thaw depth) of 0.5 m. A sediment cell is considered a permafrost cell if it has remained below the thaw depth for at least 2 years. If a column of sediment contains 75% permafrost cells or if the permafrost cells amount to ≥75 % of the inlet-channel depth h₀, the corresponding planar grid cell (i.e. in the x, y dimensions) is flagged as a “permafrost cell”.

To simulate the erosion of a permafrost bed, DeltaRCM-Arctic (Lauzon et al., 2019) used a multiplicative erodibility factor, E≤1, to scale the erosional thresholds (given in Liang et al., 2015b, Sect. 3.2.5) in (planar) permafrost grid cells, such that erosion is harder to achieve.

2.1.5 Bed diffusion and flood correction

Immediately after each round of sediment packet routing, a bed-diffusion process is applied “to take into account the influence of topographical slope on sediment flux”, and to introduce “lateral erosion by allowing sediment on the bank to be removed and added to the channels” (Liang et al., 2015b). This is achieved by calculating the “diffusive sand flux”:

$$\begin{array}{}\text{(6)}& q_{\mathrm{sand},\mathrm{diff}}=\mathit{\alpha }\left|\mathrm{\nabla }\mathit{\eta }\right|q_{\mathrm{sand}},\end{array}$$

where α is a coefficient set to 0.1 in all demonstrated cases, |∇η| is the absolute slope of the bed, and q_sand is the sand flux into and out of the grid cell concerned. Both |∇η| and q_sand are calculated across the boundary between the grid cell concerned and each of its neighbouring cells. Contributions of each of the neighbours are summed to give q_sand,diff, which is then used to determine the diffusion-induced bed-elevation change Δη of the grid cell concerned. In DeltaRCM-Arctic, α is further scaled by the erodibility factor, E.

After each full update of the water surface and each full update of the bed elevation, some grid cell may contain a water surface that is higher than the dry bed in a neighbouring cell, causing an unphysical shoreline or riverbank location. In the source codes of DeltaRCM, this is referred to as “flood correction” (Liang, 2015; Perignon, 2018) but not explicitly described in the article (Liang et al., 2015b). This means that any dry grid cells are refilled with water when one or more neighbouring cells have higher water surfaces. We assume that DeltaRCM-Arctic also inherited this mechanism and describe our interpretation of it in Sect. 2.2.4.

2.2 Modifications and interpretations in ArcDelRCM.jl

We completely rewrote the model in the Julia language, first according to the model descriptions of Liang et al. (2015b) and Lauzon et al. (2019) summarized above. Due to the non-availability of the source code of DeltaRCM-Arctic, whenever we present simulations results of DeltaRCM-Arctic, we mean the DeltaRCM-Arctic configuration of ArcDelRCM.jl. Aside from reconstructing DeltaRCM-Arctic, we also made significant changes to the model itself to improve its ability to account for processes that are climate-sensitive, with the goal of exploring the physical processes that give rise to ramp features. We describe them in turn in this section.⁴ Since the ramps appear to be related to winter ice, their depth being the same as the thickness of winter sea ice, we begin by examining ice-related processes.

2.2.1 Bottom-fast ice protection and shielding

Due to the weighted random-walk scheme and the limit of h_ice to 99.99 % of the water depth (for numerical stability), water packets can still go through grid cells where the entire water depth is effectively in the form of ice (albeit with a decreased probability, since the un-frozen water depth plays an important role in determining random-walk weights). This can generate unrealistic flow pathways with anomalously high speeds (due to ice constriction) and consequent ice melting. To eliminate this unrealistic behaviour, we prohibit flow-speed-induced melting when ice is effectively in contact with the bed (i.e. h_ice≈h). We call this “bottom-fast ice protection”. In the same locations, the cell is considered entirely blocked by bottom-fast ice and no erosion or deposition can occur. We call this “bottom-fast ice shielding” of the bed.

2.2.2 Time-dependent thaw depth

During the winter months, bottom-fast ice can bond with the sediments below through conductive heat transfer. Since we introduced ice protection and ice shielding, we deemed it logically necessary to also introduce a time-dependent thaw depth to avoid an abrupt jump (in time) from ice-bonded top sediment layers to a finite thaw depth (or active-layer depth, as defined in DeltaRCM-Arctic) when bottom-fast ice ceases to be bottom-fast. We do so according to the Stefan Model (Riseborough et al., 2008; Lunardini, 1981): $X=\sqrt{\frac{\mathrm{2}\mathit{\lambda }I}{L}},$ where X is the thaw depth, λ≡2.22 W m⁻¹ K⁻¹ is the thermal conductivity of ice near 0 ^∘C, $L\equiv \mathrm{3.3355}\times {\mathrm{10}}^{\mathrm{8}}$ J m⁻³ is the volumetric latent heat of fusion of water, and I is the “positive-degree-day index”, which is the integrated number of days times the positive temperature since winter. For I in our simulations, we use a mean temperature of 4 ^∘C to get I (see Appendix B for the reason for this choice and the sensitivity of the model output to this value). Where there is bottom-fast ice at the start of the simulation, the thaw depth starts at 0 and I only begins to increment when the ice in the pixel is no longer bottom-fast. Otherwise, I starts at 10 d in our standard simulation (see Appendix B for the reason behind this choice and the sensitivity of the model output to this value). As a result of the time-variable thaw depth, we forgo the “permafrost” label for horizontal (or planar) grid cells. In uses where the classification of vertical sediment cells as permafrost is relevant, it can be defined as the vertical cells that stayed below the maximum thaw depth (instead of a static depth of 0.5 m in DeltaRCM-Arctic) for at least 2 years, although we do not use such labels in this study. To avoid confusion, we will use the terms “frozen” or “thawed” (as opposed to “permafrost”) in the context of the erosional rules in ArcDelRCM.jl.

2.2.3 Erosion of frozen ground

To simulate the erosion of the permafrost bed, the original model used a scaled erodibility factor, E≤1, for horizontal grid cells with over 75 % permafrost content in the sediment column. With the introduction of time-dependent thaw depth, we find it more self-consistent to check the calculated erosional depth of the grid cell against its corresponding sediment column: if the calculated erosion reaches deeper than the available thawed layers, the erosional depth is limited to the thawed (vertical-cell) layers only. Whilst the sediment column is immediately updated as a sediment packet passes through, which prevents duplicate erosion/deposition by successive packets, the value of the bed elevation is kept unchanged in the memory during each time step, thus tracking the exact layer that is at the bottom of the thawed section. As the model proceeds to the next time step, more layers are thawed as the thaw depth deepens as described in the previous subsection. This approach is taken in both erosion by packets routed through individual pixels and in the bed-diffusion process explained in Sect. 2.1.5. We thus forgo the use of the erodibility factor, E (and, as mentioned in Sect. 2.2.2, the labelling of any pixel as permafrost), and consider only the depth of the boundary between frozen and thawed grounds.

Conceptually, both DeltaRCM-Arctic and ArcDelRCM.jl allow erosion of sediment columns with frozen (vertical) cells. Consider a sediment column with a significant fraction of its cells frozen during some number of time steps (i.e. with a shallow thaw depth in ArcDelRCM.jl or with the corresponding horizontal grid cell labelled as permafrost in DeltaRCM-Arctic). DeltaRCM-Arctic would allow for erosion with a scaled-up flow-speed threshold through E (and a similarly scaled-down bed diffusion); the cumulative erosional effect on the corresponding column can thus be thought of as a fraction of the equivalent no-permafrost case (i.e. E=1). By allowing only thawed (vertical) sediment cells to be eroded, the same sediment column in ArcDelRCM.jl also undergoes erosion that is a fraction of an equivalent case without any frozen (vertical) cells. However, the main difference between the two is the timing of the erosional events. ArcDelRCM.jl allows the thawed vertical cells to be eroded unhindered and thus earlier (considering the chance of having a lower versus higher flow speed across the domain at any given moment) but then delays the erosion of the vertical cells frozen at the time, which would become thawed in the next time step and freely erodible again.

2.2.4 Flood correction

A flood correction mechanism is built into the original, non-Arctic source codes of Delta RCM (Liang et al., 2015a; Perignon, 2018) but not explicitly described in the article of Liang et al. (2015b). This mechanism ensures that dry grid cells are refilled with water when one or more neighbouring cells have higher water surfaces. Without access to the source code of DeltaRCM-Arctic (Lauzon et al., 2019), how they handled the presence of ice in this context is unclear. In ArcDelRCM.jl, we interpret the water-surface elevation in this context as the below-ice surface. Therefore, only dry cells with at least one neighbouring cell that has a higher liquid water-surface elevation are rewetted. The water-surface of the rewetted cell is calculated from the average of all surrounding wet cells.

2.2.5 Time-step size

To keep the simulation numerically stable, the original model determines the time-step size based on the volume of sediments entering the simulation domain. Specifically, $\mathrm{\Delta }t=\frac{N_{\mathrm{0}}^{\mathrm{2}}{h}_{\mathrm{0}}{\mathit{\delta }}_c^{\mathrm{2}}}{\mathrm{10}{Q}_{{\mathrm{s}}_{\mathrm{0}}}}$ is given as a guideline by Liang et al. (2015b) to prevent too much sediment from entering in each time step relative to the accommodation space for deposition in the grid cells. In our simulations intended to mimic the Lena Delta, we discovered that, where the grid-cell dimensions (which are terms in the numerator) are several times larger than in Liang et al. (2015b) and Lauzon et al. (2019), the 10 in the denominator of the expression for Δt needed to be increased by a factor of a few (e.g. to 40–80, depending on the volume of sediments entering the domain in each time step in the specific simulation). Otherwise, the model becomes numerically unstable. In order to facilitate the use of time-dependent input discharges (described in the next subsection), we let Δt be user-determined but implemented an internal checking procedure to warn users of potential numerical instability. For this check, we introduce a quantity that we call “scale-height measure”: $\mathit{\zeta }\equiv \frac{Q_{{\mathrm{s}}_{\mathrm{0}}}\mathrm{\Delta }t}{{\mathit{\delta }}_{\mathrm{c}}^{\mathrm{2}}{h}_{\mathrm{B}}}$. ζ is a rough scaling between the volume of sediment entering through the inlet channel in each time step and the “available” volume of an average single ocean-basin grid cell. This is in the same spirit as the “reference volume” introduced by Liang et al. (2015b, Eq. 1) that led to their guideline expression for Δt. Based on experiments with both the dimensions of Liang et al. (2015b), Lauzon et al. (2019), and our Lena-like dimensions, we found that ζ should be ⪅3 for the model to be numerically stable. Due to the random-walk nature of the RCMs, this threshold is not sharp.

2.2.6 Input discharges as time series

In large drainage basins such as the Lena watershed, discharge beyond the spring flooding season remains significant through the summer. In order to capture these deltas' summer evolution, we modify the model to take in time series of input discharges (both water and sediments) and extend the simulation model year to include the summer months. Under this setup, the inlet flow speed, u₀, the inlet flow depth, h₀, and the reference water-surface slope (used as an approximation of the backwater slope), S₀, all become time series themselves and are dependent on the water discharge $Q_{{\mathrm{w}}_{\mathrm{0}}}$ time series.⁵ To simplify the process for the user and to reduce the chance of mistakes, ArcDelRCM.jl users input the minima of u₀, h₀, and S₀, corresponding to the minimum $Q_{{\mathrm{w}}_{\mathrm{0}}}$. The simulation then uses the time series of $Q_{{\mathrm{w}}_{\mathrm{0}}}$ to calculate the corresponding time series of u₀ and h₀ based on a scaling derived from the Gauckler–Manning formula (Gauckler, 1867; Manning et al., 1890) and of S₀ based on simple geometric arguments. We describe these in order.

Figure 2 Satellite imagery of the Lena Delta from 2018 showing the delayed breakup of bottom-fast ice on the ramp feature. The upper panels (a–c) show the imagery of the Lena Delta. The hash pattern in red marks the location of the ramp, and (a), (b), and (c) correspond respectively to the start, middle, and completion of ice breakup in the delta itself. The corresponding dates and discharge (measured at Kyusyur station and corrected for the distance to the Lena Delta; Juhls et al., 2020) are shown in (d), marked by the red arrows. Satellite imagery is acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) and obtained from NASA Worldview.

We first assume open-channel flow and that the (overall) channel bed slope is approximately constant. The average flow velocity is given by the Gauckler–Manning formula (Gauckler, 1867; Manning et al., 1890), which, for our purpose and under our assumptions, can be expressed as a proportionality:

$$\begin{array}{}\text{(7)}& \overline{u}\propto R_{\mathrm{h}}^{\frac{\mathrm{2}}{\mathrm{3}}},\end{array}$$

where R_h is the hydraulic radius:

$$R_{\mathrm{h}}={\displaystyle \frac{A}{P}},$$

where A is the cross-sectional area of the channel and P is the “wetted perimeter” of the channel. Assuming a simplistic, rectangular channel shape, with width, w, and flow depth, h, we have A=hw and $P=\mathrm{2}h+w$. Then, using the expression for discharge, $Q=A\stackrel{\mathrm{‾}}{u}$, we calculate the flow depth under some new discharge, Q_n, in relation to an “original” discharge, Q_o, as follows:

$${\displaystyle \frac{Q_{\mathrm{n}}}{Q_{\mathrm{o}}}}={\displaystyle \frac{A_{\mathrm{n}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{n}}}{A_{\mathrm{o}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{o}}}}={\displaystyle \frac{h_{\mathrm{n}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{n}}}{h_{\mathrm{o}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{o}}}}={\displaystyle \frac{h_{\mathrm{n}}{\left(\frac{h_{\mathrm{n}}w}{\mathrm{2}{h}_{\mathrm{n}}+w}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}}{h_{\mathrm{o}}{\left(\frac{h_{\mathrm{o}}w}{\mathrm{2}{h}_{\mathrm{o}}+w}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}}}\approx {\left({\displaystyle \frac{h_{\mathrm{n}}}{h_{\mathrm{o}}}}\right)}^{\frac{\mathrm{5}}{\mathrm{3}}},$$

where we have assumed a constant w and that w is much greater than both h_o and h_n. Similarly, the transformation of $\stackrel{\mathrm{‾}}{u}$ from Q_o to Q_n can be calculated, via Eq. (7), by

$${\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{n}}}{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{o}}}}={\left({\displaystyle \frac{h_{\mathrm{n}}}{h_{\mathrm{o}}}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}.$$

Given an original slope, S_o, corresponding to h_o under the discharge Q_o (with ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{o}}$), one can define a “baseline reach”, $L_{\mathrm{B}}\equiv h_{\mathrm{o}}/S_{\mathrm{o}}$, such that the new slope under the discharge Q_n (with h_n and ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{n}}$) can be approximated by

$$S_{\mathrm{n}}\approx {\displaystyle \frac{h_{\mathrm{n}}}{L_{\mathrm{B}}}}.$$

In the manner described above, the inlet flow depth, the inlet flow speed, and the reference slope corresponding to the input discharge time series, Q_w, can be calculated as long as h₀, u₀, and S₀ corresponding to the minimum Q_w are supplied.

The sediment discharge, $Q_{{\mathrm{s}}_{\mathrm{0}}}$, can also be input as a time series. However, many large Arctic deltas have sediment fluxes that scale approximately with water discharge (Holmes et al., 2002). Therefore, following DeltaRCM, a simple multiplicative factor is used to translate $Q_{{\mathrm{w}}_{\mathrm{0}}}$ to $Q_{{\mathrm{s}}_{\mathrm{0}}}$.

2.2.7 Controls on the melting of ice

We have added the capability for the users to shift the timing and duration of the ice cover's melt. This modification is motivated by the observations that bottom-fast ice just offshore from the delta remains intact longer than the ice cover in the delta and does not disintegrate until the delta itself is almost ice-free or only contains the so-called “serpentine (floating) ice” in the channels (Fig. 2). This ice remains bottom-fast for longer perhaps due to its additional thickness from congelation from below during surges or from water overflowing its top and depressing it during floods (Reimnitz, 2002). Moreover, based on examinations of satellite imagery, the duration of ice-cover melting in the Lena Delta is closer to 20 d instead of 10 d. During this time, both the water temperature and air temperature remain close to the freezing point (Reimnitz, 2002; Yang et al., 2005; Juhls et al., 2020). An example observation of this breakup sequence and duration is shown in Fig. 2. In ArcDelRCM.jl, users can specify the length of the ice-cover melting period, during which atmospheric heat contributes to the melting. Note that flow-speed-induced melting is always active, subject to the protection described in Sect. 2.2.1.

Figure 3 Exemplary overview of workflow to extract the hydrological graph from the imaged delta simulation. (a) Bed elevation, (b) binarization of bed elevation to differentiate between channels and non-channels, (c) skeleton (in red) of the channels overlaid on the binary image (b), and (d) graph of the channel network derived from the skeleton in (c) with edges in red and nodes in black overlaid on the delta elevation.

2.3 Graph analyses on ensemble runs

To assess quantitatively and statistically the effects of the various modifications described in Sect. 2.2, we apply graph theory to derive metrics on a collection of ensemble runs, each set with a different process (corresponding to Sect. 2.2.1 to 2.2.4) or a combination thereof enabled. Previous work by, for example, Smart and Moruzzi (1971), Edmonds et al. (2011), Tejedor et al. (2015), and Nesvold et al. (2019) has shown that the topologies of deltas can be described with quantitative graph metrics, such as the “loopiness” and the structural overlapping of the subnetworks, the “recombination factor” describing the ratio between the number of junctions and the number of forks in the delta systems, and the fractal dimension characterizing a delta's self-similarity. While these metrics give interesting insights into the environmental properties of the real-world deltas, we sought a holistic approach that would quantify the differences between simulation results with the simplest yet meaningful descriptors. We thus made use of the approach introduced by Rettelbach et al. (2021), which provides an end-to-end approach starting from the extraction of the graph from the delta images to providing the quantitative metrics of interest for the comparison of the ensemble runs with different parameters and forcings. While the approach by Rettelbach et al. (2021) was initially developed for characterizing hydrological networks in polygonal permafrost landscapes, the methodology of automated graph extraction from underlying imagery remains exactly the same. In the original publication, the authors used binarized digital elevation models to compute the skeleton of the channel network, while we here extracted the graph from the deltas binarized based on the location of the channels (see Fig. 3). In this context, channel pixels are defined as those having a water depth of 0.1 m or over, which is the threshold value for “dry” versus “wet” pixel labels used in DeltaRCM (Liang et al., 2015b) and inherited by both DeltaRCM-Arctic and ArcDelRCM.jl. “Open-ocean” pixels without any depositions (i.e. with depth equalling the ocean-basin depth, h_B) are excluded. From the derived graph (as seen in Fig. 3d in red and black), we then calculated the following metrics:

• number of nodes and edges, giving us an idea of the size and complexity of the underlying hydrological network.

• the number of connected components, where each connected component represents a hydrological subnetwork. One – and only one – component will always be connected with the apex of the delta, and thus any edges in this component (channels in this subnetwork), are considered “active”. Any other component will represent an abandoned channel (or cluster of channels) that is no longer fed by the upstream river. Evaluating this number in combination with the graph density (see below), we can gather similar information as Tejedor et al. (2015) do with their metric of resistance distance.

• the total length of all channels combined, which makes quantification of the number of all potential waterways possible.

• the graph’s density, which lets us estimate the network flow effectiveness. It is defined via the ratio of the number of existing edges over the number of edges that could exist based on the number of nodes, n. For planar graphs, this equals 3(n−1). It can be seen as a simpler alternative to the self-similarity measure of fractal dimension introduced by Edmonds et al. (2011).

• the graph's diameter, representing the longest of all shortest path lengths between vertices. In combination with the total length of channels, this gives us an idea of the asymmetry of the delta.

These metrics from basic graph theory are efficient to calculate and can already provide valuable insights into the properties of deltaic networks and accurate parameters to quantify and compare the results of the different ensemble runs, thus providing a simple yet powerful alternative to the more complex ones introduced and described by Smart and Moruzzi (1971), Edmonds et al. (2011), Tejedor et al. (2015), and Nesvold et al. (2019).

3 Results

We present here the outputs of ArcDelRCM.jl and its DeltaRCM-Arctic configuration (for comparison), with particular attention paid to the ramp feature. All simulations share the following parameters: (N_x, N_y) = (150, 300); N_wall=3; the number of water and sediment packets (sand and mud separately) are $n_{\mathrm{w}}=n_{\mathrm{s}}=\mathrm{2000}$; the coefficient for bed diffusion α=0.1; and the scaling factor for atmospherically induced melting a=0.5. Further parameters applicable to individual cases are specified in the respective subsections below. Any remaining parameters and conditions not explicitly listed take on values and specifications given in Liang et al. (2015b). Note that the colour-blindness-friendly colour scheme, uniform across all the filled-contour figures in this section, is chosen to highlight the per-metre gradation of elevations below the water surface, with a focus on the shallower depths where the existence or absence of the ramp can be seen. Interpretations and speculations are left to Sect. 4.

3.1 Analogous setup to DeltaRCM-Arctic demonstrations

In this subsection, we present comparisons of simulations, run with identical parameters and identical random seeds, in ArcDelRCM.jl and its DeltaRCM-Arctic configuration. Specifically, we adopt δc=50 m; N₀=5; h₀=5 m; h_B=5 m; u₀=1 m s⁻¹; $S_{\mathrm{0}}=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{4}}$; Q_w=1250 m³ s⁻¹; $Q_{\mathrm{s}}/Q_{\mathrm{w}}=\mathrm{0.001}$; a sand fraction (of the total sediment volume) of 25 %, a maximum ice extent of 40 %, $h_{\mathrm{ice},\max }=\mathrm{2}$ m; and γ=0.0735. With the exception of h_ice,max, these values are chosen after the demonstrated cases of DeltaRCM-Arctic in Piliouras et al. (2021) and Lauzon et al. (2019). In the case of DeltaRCM-Arctic, we use an erodibility factor E=0.8, which is the “middle” value used in Lauzon et al. (2019) and Piliouras et al. (2021). For ArcDelRCM.jl, we specify Δt=25 000 s, matching the Δt given by the DeltaRCM(-Arctic) guideline. All cases are run for 5000 time steps, with an additional 300 “lead-up” steps under non-Arctic conditions to build up a seed flow field (following Lauzon et al., 2019; Piliouras et al., 2021).⁶

The first row of Fig. 4 shows the output from DeltaRCM-Arctic (Fig. 4a) and ArcDelRCM.jl (Fig. 4b). The second row of Fig. 4 shows the output from the same simulations as the first row but with $h_{\mathrm{ice},\max }=\mathrm{3}$ m (Fig. 4c and d). The last row of Fig. 4 shows the output of ArcDelRCM.jl under identical configurations as in Fig. 4b, except h_B is increased to 7 m in Fig. 4e, and the ice extent is increased to 100 % in Fig. 4f.

Figure 4 Bed-elevation output of (a, c) DeltaRCM-Arctic and (b, d–f) ArcDelRCM.jl after 5000 time steps (following a 300-step lead-up under non-Arctic conditions). All runs have identical parameters (see text for full configuration), except the following differences: (a, b) h_ice,max is 2 m; (c, d) h_ice,max is 3 m; (e) h_B is 7 m; and (f) the ice extent is 100 %. Note the depths of the ramp features in panels (b) and (d), which correspond to h_ice,max. Tentacle-like “offshore depositions” (as described in Piliouras et al., 2021) are visible in the $h_{\mathrm{ice},\max }=\mathrm{3}$ m cases of DeltaRCM-Arctic (in panel c) as well as outside the ramp of ArcDelRCM.jl (in panel d). Also, the ramp feature has vanished in panel (e) but has become more prominent in panel (f).

The ramp features form continuous bands around all the deltas in ArcDelRCM.jl except in the case where h_B=7 m, in which the ramp is not well-formed (Fig. 4e). The ramp appears to be more prominent in the case with 100 % ice extent (Fig. 4f) and the case with $h_{\mathrm{ice},\max }=\mathrm{3}$ m (Fig. 4d). The deltas of DeltaRCM-Arctic do not show such ramps (Fig. 4a) but rather display lobes or tentacles of “offshore depositions” (as they are called in Piliouras et al., 2021) around channel outlets (Fig. 4c). The $h_{\mathrm{ice},\max }=\mathrm{3}$ m case (Fig. 4c) has prominent offshore deposits around channel outlets that are at the right depth for the ramp feature (also observed by Piliouras et al., 2021) but has less deposition in the lateral directions to join with neighbouring lobes to form a band. The lobes are also more disrupted by additional deposits to higher elevations. As we will explore in the next subsection, we find that the lateral joining appears after switching on the protection of bottom-fast ice (Sect. 2.2.1). A visual representation of the thaw-depth patterns corresponding to Fig. 4b and f is shown in Appendix C.

3.2 Individual modifications in ArcDelRCM.jl

Figure 5 shows the effects on the simulated deltas arising from the individual modifications described in Sect. 2.2.1 to 2.2.4. The model settings are identical to the $h_{\mathrm{ice},\max }=\mathrm{2}$ m cases of Fig. 4 but with a different random seed. The random seed across all cases in Fig. 5 is identical, however. The cases where the erosion is limited to thaw depth (Sect. 2.2.3) but retains DeltaRCM-Arctic rules for bed diffusion (Sect. 2.1.5) use E=0.8 to match the other cases (Fig. 5c and e).

Figure 5 Bed-elevation output of DeltaRCM-Arctic (a) and ArcDelRCM.jl (b) and with modifications described in Sect. 2.2.1 to 2.2.4 applied individually or in tandem, as indicated by the title of each panel (c–h). Note the 2 m ramp in the full model and the case in which ice shielding and ice protections are applied together.

All the cases with erosion limited to thawed layers instead of using the erodibility factor (Sect. 2.2.3) appear similar to each other (Fig. 5c–f), but the cases in which the thaw depth is time-dependent (Sect. 2.2.2) have less tentacle-like deposits reaching seaward from channel outlets (Fig. 5e and f).

The case with only ice shielding of the bed (Sect. 2.2.1) displays tentacle-like deposits (Fig. 5g), similar to those in Fig. 5c and d. The case in which bottom-fast ice is protected from flow-induced melt (i.e. ice protection; Sect. 2.2.1) in addition to ice shielding of the bed exhibits joining of deposition lobes around channel outlets to form a visible 2 m ramp around the delta, in addition to further tentacle-like deposits seaward (Fig. 5h). The cases for DeltaRCM-Arctic and the full ArcDelRCM.jl are given for visual comparison (Fig. 5a and b).

3.3 Quantitative comparison of delta ensemble runs

To assess quantitatively and statistically the individual modifications represented in Fig. 5 and their potential influence on the ramp feature, we applied graph theory analyses as described in Sect. 2.3 to an ensemble of 105 realizations for each case. Figure 6 shows the resulting delta metrics. In addition to the cases shown in Fig. 5, we also run variants of them that are not shown but the metrics of which are included here. These variants are as follows: all the runs where the erodibility factor is applicable with E=0.65 (i.e. the cases that use DeltaRCM-Arctic rules for all erosion), the non-Arctic DeltaRCM runs, and the ArcDelRCM.jl run with unhindered bed diffusion (i.e. no thaw-depth limits, equivalent to E=1, during the bed-diffusion step).

Figure 6 Plots showing delta metrics determined by the graph analysis we performed on each of the configurations described in Sect. 3.2. The blue violin plots show the kernel density estimate, the orange box and whisker plots show the central 50 %, maximum, median, and minimum values, and the black dots show the actual values of the individual simulations. The labels along the vertical axes correspond to the titles of each panel in Fig. 5. The metrics of the variant cases not explicitly shown in this text (i.e. all cases with E=0.65, ArcDelRCM.jl with unhindered bed diffusion, and the non-Arctic DeltaRCM) are included for comparison.

In terms of the number of nodes and edges, the cases divide roughly into three clusters: the cases with individual modifications turned on (Fig. 5c–f), the DeltaRCM family, and the ArcDelRCM.jl family. Most of the variabilities overlap each other, however.

The number of connected components (reflecting the number of abandoned channels or subnetworks) and the graph density (actual versus all possible connections between nodes) display an inverse relation to each other. Here the cases can also be viewed as having the same three clusters as the number of nodes and edges. However, the separation is weaker (i.e. with more overlap across variability). Across all the cases that use the erodibility factor, E, those with E=0.8 show a tendency towards less abandoned channels or subnetworks and thus are better connected (i.e. have higher graph density). The non-Arctic case (effectively meaning E=1 but also ice-free) occupies roughly the same range as the E=0.65 cases of DeltaRCM-Arctic.

The maximum diameters, measuring the longest of all the shortest paths between vertices, are not significantly different between all the cases considered. There is a tendency towards longer single channels for the cases with ice shielding and/or bottom-fast ice protection (see also Fig. 5g and h), affecting the maximum diameter. The total lengths of all channels are also longer in these cases (which all use the erodibility factor as in DeltaRCM-Arctic). These are, in turn, similar to the lengths of channels in the DeltaRCM cases (including non-Arctic). On the other hand, the cases with erosion limited to thaw depth (but without ice-related modifications), regardless of specific variants, display tendencies towards having shorter total lengths of channels. All ArcDelRCM.jl cases occupy the range of total channel lengths in between.

To gauge the size of the ramp feature, which is not straightforward to quantify due to the irregular shape and distance from the apex, we use the number of pixels with elevations ≥ −2.5 and $<-\mathrm{1.5}$ m as a proxy (Fig. 7). The non-Arctic cases have the lowest counts of such pixels, whilst the cases with both ice shielding of the bed and bottom-fast ice protection have the highest. Protection of bottom-fast ice appears to increase this count compared to ice shielding of the bed alone. The cases with erosion limited to thaw depth, including all ArcDelRCM.jl cases, show a similar range in this pixel count, with the cases using a static thaw depth giving higher counts. Notably, both DeltaRCM-Arctic cases occupy the same range as the ArcDelRCM.jl cases, even though the spatial distributions of the pixels are different (see Fig. 5a and b). Finally, between cases that use E=0.8 and E=0.65, the cases with E=0.8 show higher counts of pixels within this elevation range.

Figure 7 An analogous plot to Fig. 6 but showing the number of pixels with elevation $\in [-\mathrm{2.5}$, 1.5) m. Since the modelled ramp feature that forms under $h_{\mathrm{ice},\max }=\mathrm{2}$ m lies mostly in this elevation range, we use this value as a proxy for the size of the ramp, even though it inevitably includes delta pixels that are not part of a ramp feature but happen to be in the same elevation range. This is evident when comparing the DeltaRCM-Arctic (E=0.8) and ArcDelRCM.jl cases with Fig. 5a and b. The pixels with the relevant elevations are more concentrated around channel outlets in the former, while they are more distributed along the delta shore in the latter (forming the ramp).

3.4 Lena Delta approximants

As a test case to mimic a large delta such as the Lena Delta, we ran simulations adopting spatial scales that had never been applied to ice-dominated delta using this family of models before. (However, the non-Arctic DeltaRCM has previously been applied on a similar spatial scale on the Mississippi and Selenga deltas by Moodie and Passalacqua, 2021.) Specifically, we use δc=400 m; N₀=6; $h_{\mathrm{0},\min }=\mathrm{10}$ m; h_B=15 m; $u_{\mathrm{0},\min }=\mathrm{1}$ m s⁻¹; $S_{\mathrm{0},\min }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$; $Q_{\mathrm{s}}/Q_{\mathrm{w}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ (roughly 10 times the average volume fraction measured in the Lena Delta; Boike et al., 2019); a sand fraction of 20 % (Alekseevsky, 2007, as cited in Fedorova et al., 2015); a maximum ice extent of 100 %; $h_{\mathrm{ice},\max }=\mathrm{2}$ m; γ=0.135; and a time step of 1 d.

We would first like to find out the differences arising from assuming a 10 d model year (suitable for smaller Arctic deltas, as adopted by Lauzon et al., 2019; Piliouras et al., 2021) and from assuming a 4-month model year including the summer months (suitable for large Arctic deltas such as Lena Delta and Mackenzie Delta) due to the resulting different total sediment input per year. We ran one batch of simulations for 150 model years. Within this batch, the discharge Q_w is treated differently in order to probe the difference it makes in terms of the number of days per simulation year and thus sediment input per year, given realistic Q_w(t). In each model year, the “full” simulation cases cover the 4 months from 1 June to 30 September (122 d). In these cases, Q_w(t) is a time series constructed from daily discharge data from the Stolb station near the main channel into the Lena Delta (GRDC Station Data 2903430, 2018), with daily values averaged from 1951 to 1980 inclusive (Fig. 8).

Figure 8 Daily discharge measured at Stolb station (thick blue line; GRDC Station Data 2903430, 2018), averaged from 1951 to 1980 inclusive. Discharge remains over 20 000 m³ s⁻¹ during the months from June to September (between the thin grey dotted lines), which is the simulation period for the 4-month Lena-approximant cases. The grey band spans the peak 10 d of discharge. The light-blue, thinner line shows the same discharge pattern, except the overall discharge has been scaled up by 35 % (representing the RCP7–8.5 scenario) whilst keeping the peak value and the shape of the curve the same. The period during which discharge is over 20 000 m³ s⁻¹ is longer, at 136 d (between the thick grey dashed lines).

The “10 d” cases have two variants: the “constant averaged peak discharge” case (Fig. 9a) uses the averaged value from the peak 10 d of Q_w(t) as the constant discharge, and the “time-variable discharge during peak” case (Fig. 9b) uses the peak 10 d of the Q_w(t) time series as the input discharge. The 10 d peak period is highlighted by a grey band in Fig. 8.

Figure 9 Bed elevations of deltas produced by running ArcDelRCM.jl for 150 model years on Lena Delta-like spatial scales (see text in Sect. 3.4), with input discharge derived from daily measured values from GRDC Station Data 2903430 (2018). The top row (a, b) features deltas produced by running model years of 10 d each, which is also the ice-melt period. Discharge in these 10 d cases is taken from the peak 10 d of the time series and either (a) averaged and used as a constant value or (b) used directly as a 10 d discharge time series. The bottom row (c, d) features deltas that are produced by 4-month model years (June to September), using the full input discharge time series for the corresponding period. The case in panel (c) kept the ice-melt period of 10 d, whilst the case in panel (d) has an ice-melt period of 20 d and a delayed onset of atmospherically induced melting by 20 d from the start of each model year. Note the difference in size between the top and bottom rows and the ramp feature around the delta in panel (d).

The full 4-month cases are also divided into two variants: one in which the ice-melt period is 10 d (Fig. 9c), similar to the aforementioned 10 d cases, and one in which the ice-melt period is 20 d long and is delayed by 20 d from the start of each model year (Fig. 9d; motivation described in Sect. 2.2.7). In all cases, flow-induced melting is active (where allowed) throughout the whole of the simulation. We use “ice-melt period” to refer only to the time during which the atmospheric contribution is active (i.e. t_melt in Eq. 5).

The resulting deltas in the 10 d cases are entirely under water and similar to each other in extent (Fig. 9a and b). The 4-month cases also produced deltas that are similar in extent to each other (Fig. 9c and d) but reach twice as far from the inlet wall (or 4 times the area) as the 10 d cases. A ramp is also visible around the delta in the case with 20 d ice-melt period delayed by 20 d (Fig. 9d), dipping from just below water level towards 2 m, albeit with a rougher surface with shallower depths (∼1 m; fourth contour level from the top) than in the small-scale cases in Sect. 3.1.

As an additional demonstration of the model in approximating the scale of the Lena Delta, we run the simulation with the same parameters but with $Q_{\mathrm{s}}/Q_{\mathrm{w}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$, reflecting the measured average volume fraction of sediment to water (Boike et al., 2019). This low sediment-to-water volume fraction requires a much longer run time to produce a delta. Therefore, we ran the simulation for 1200 model years (with the same Δt=1 d). To further mimic the underlying ocean bathymetry on the Laptev Sea coast, where the Lena Delta is situated, we introduced a gradual tilt of the ocean-basin elevation h_B: from 10 m at the inlet wall to 20 m on the opposite side of the simulation domain. The extent of the tilt and the 20 m maximum h_B is motivated by inspecting the bathymetry of the Laptev Sea coast (Fuchs et al., 2022). The resulting delta is shown in Fig. 10. The extent of the delta is similar to the 4-month cases (at 150 model years) in Fig. 9, whilst the ramp feature has a clearer dip towards ∼2 m depth than in the other Lena-approximant cases.

Figure 10 A delta produced by ArcDelRCM.jl after 1200 model years with configurations identical to those in Fig. 9d but with the low sediment-to-water volume ratio observed (Boike et al., 2019) and with a tilted ocean-basin bed motivated by the bathymetry of the Laptev Sea coast near the Lena Delta (see text in Sect. 3.4 for details). The panels show (a) the filled-contour view and (b) the gradient-coloured view of the same bed elevations. The latter is to facilitate visual comparisons with Fig. 1.

3.5 Ramp feature under a warming climate

To explore how a warming climate might affect the ramp feature, we continued the simulation of the delta shown in Fig. 10 for another 1200 model years. However, for this portion, we adopt an end-member scenario of roughly Representative Concentration Pathway (RCP) 7 to 8.5 (Stocker, 2014), which corresponds to ∼4 ^∘C of global warming by the year 2100. Under this scenario, maximum ice thickness during winter is not expected to reduce drastically (Nummelin et al., 2016; Sun et al., 2018), although some thinning has been suggested (Landrum and Holland, 2020). We therefore adopt $h_{\mathrm{ice},\max }=\mathrm{1}$ m thickness instead of 2 m. The discharge at the Lena Delta has also been observed to be increasing in recent years (Fedorova et al., 2015). For this warming scenario, we also adopt an overall increase of 35 % in total discharge (Mann et al., 2022), whilst the peaks remain the same (motivated by Juhls et al., 2020, in which they found an increased winter discharge; see also Mann et al., 2022). The overall time period during which discharge is over 20 000 m³ s⁻¹ is 14 d longer. This modified discharge pattern is shown in Fig. 8. Furthermore, based on the surface-temperature increase (and the rate of temperature increase during the spring and summer months; Sun et al., 2018), the atmospheric heat-induced melting of ice cover is brought forward by 10 d. The duration is also shortened by 10 d due to the thinner ice. All other parameters remain the same as described in Sect. 3.4.

Figure 11 shows various snapshots of the continued evolution of the Lena-approximant delta. Under the warm conditions, the ramp feature has diminished by the 200-year mark (Fig. 11a), and becomes mostly disrupted by the 400-year mark (Fig. 11b). From 800 years onwards, no continuous ramp features remain (Fig. 11c and d).

Figure 11 The continued evolution of the delta shown in Fig. 10, except now under conditions possible in an end-member climate-warming scenario (based roughly on RCP7 to 8.5). As before, coloured contours reflect the bed elevations. The four panels show (a) 200, (b) 400, (c) 800, and (d) 1200 model years into this continued portion of the simulation. Note the degradation and disappearance of the ramp feature.

4 Discussion

The results (Fig. 4) demonstrates that ArcDelRCM.jl is able to reproduce the 2 m ramp around Arctic deltas (contrast with Fig. 4a), and Fig. 4d shows that the ramp is related to the maximum ice thickness (h_ice,max). As observed in Lauzon et al. (2019) and Piliouras et al. (2021), offshore deposits do occur as tentacle-like features in DeltaRCM-Arctic, especially when thick ice is imposed on a shallow domain ($h_{\mathrm{ice},\max }/h_{\mathrm{B}}=\mathrm{3}/\mathrm{5}$; Fig. 4c). However, our modifications in ArcDelRCM.jl led to similar features in addition to a continuous ramp not limited to around channel outlets (Fig. 4d). Increasing maximum ice extent to 100 % led to a more prominent ramp (Fig. 4f), further supporting the idea that ice is the driving factor behind the ramp feature (Reimnitz, 2002). However, having a deeper ocean basin (i.e. more accommodation space) appears to impede the development of the ramp, as seen in Fig. 4e (in stark contrast to Fig. 4d). This suggests that the available space under maximum ice thickness (h_ice,max) and thus the dynamics of transport under ice cover around the shore of deltas play a determining role in the formation of ramp features. The less space there is, the more flow is constricted and the faster sediments build up against the bottom of the ice cover, forcing lateral deposition, which forms a continuous band that becomes the ramp feature. Previous work by Lauzon et al. (2019) and Piliouras et al. (2021) using DeltaRCM-Arctic observed that offshore deposits increase with h_ice,max and thus with decreasing accommodation space in the ocean basin. Our results are therefore in agreement.

Figures 5 and 6 together show the various effects of individual modifications detailed in Sect. 2.2. Using the original DeltaRCM-Arctic model (Fig. 5a) as a baseline, limiting erosion to thawed layers whilst keeping a static thaw depth (Fig. 5c and d) and foregoing the erodibility factor has the effect of allowing existing channels to erode more easily down to the thaw depth. The more readily available sediments combined with some retained constriction to the flow (from the thaw-depth limitation) lead to more sediments being carried seaward along the same channel paths, resulting in tentacle-like offshore deposits at ∼2 m depth (reflected in the pixel count in Fig. 7). Further making the thaw depth time-dependent (as in Fig. 5e and f), erosion can still occur (without the erodibility factor) down to the thaw depth, but the thaw depth now stays shallower throughout (approximately between 0.2 and 0.3 m). This leads to similar characteristics to the cases in Fig. 5c and d in the central (ice-free) part of the delta, but the shallower erodible depth leads to less offshore deposits (Fig. 7) and more short channels in the ice-cover zone branching in search of unblocked pathways. The latter helps spread sediments more evenly along the delta shoreline.

The relationship between the extent of offshore deposition and erodible sediments (either through erodible thaw depth or the erodibility factor) is echoed in all the cases that use DeltaRCM rules for erosion. Cases with E=0.8 result in more pronounced offshore deposits at ∼2 m depth – potential building material for the ramp – than those with E=0.65 (Fig. 7).

The cases where the thaw-depth erosional limit is included in the bed-diffusion process (Fig. 5d and f) do not appear to be substantially different from the cases where it is excluded (Fig. 5c and e). The same is reflected in all the graph and pixel-count metrics (Figs. 6 and 7). This shows that bed diffusion has a relatively minor effect on the delta's form compared to the flow-driven erosion and deposition.

The ice-shielding case (Fig. 5g) has identical erosion mechanisms to DeltaRCM-Arctic (Fig. 5a), except erosion and deposition are blocked wherever ice is bottom-fast. This enhances flow constriction by ice, which focuses erosion on the few unblocked pathways (e.g. sub-ice channels) and leads to sediment being carried farther seaward. The result is a tentacle-like offshore deposition pattern similar to Fig. 5c and d, although the underlying mechanisms differ. The far-reaching tentacles yield a higher total length of channels and a larger maximum-diameter metric (Fig. 6).

Given the locations and elevations of the tentacle-like offshore deposition, either in DeltaRCM-Arctic or the partial ArcDelRCM.jl configurations of Fig. 4c–g, they are likely the foundation of the ramp features, as Lauzon et al. (2019) and Piliouras et al. (2021) suggested. As they extend outwards from various channel outlets, however, they are still separated from each other instead of forming a band-like ramp feature around the delta shoreline.

Figure 5h shows that the protection of bottom-fast ice from flow-induced melting (Sect. 2.2.1) is ultimately the process that gave the model the ability to produce the ramp feature. In this case, the erosion and deposition rules are the same as in the ice-shielding case (Fig. 5g), giving a similar tentacle pattern (and similarly elevated metrics of maximum diameter and total length of channels; Fig. 6). However, the bottom-fast ice survives for longer when protected from melting. This gives rise to two effects: (i) sediment deposition is more evenly spread out along the shore due to the constriction of flow by ice (similar to the cases of Fig. 5e and 5f); (ii) previously deposited material is protected, allowing depositional lobes to expand and merge subaqueously without being deposited over the top. The ramp feature is thus stabilized and sustained as a result.

Considering what we observed in the model outputs so far, the following balance may be at play in the formation of the ramp feature: the existence of ice constricts flow and promotes offshore deposition (Lauzon et al., 2019; Piliouras et al., 2021). The erosional rules of ArcDelRCM.jl further allow the erodible layers (i.e. thawed layers) to be eroded to form more substantial offshore deposits. The introduction of time-dependent thaw depth leads to shallower thaw depths within the same 10 d model year, offsetting the previous effect somewhat, but the ice constriction is enhanced. Finally, the shallower thaw depth and ice shielding of the bed work in tandem with the longer survival of bottom-fast ice to both constrict flow and cause it to spread out laterally. This leads to more even deposition along the shore, which is then protected by the longer-surviving bottom-fast ice, resulting in the ramp. As a pair of examples, Fig. C1 shows the spatial views of the thaw-depth pattern corresponding to simulations shown in Fig. 4b and f. These thaw-depth patterns provide an indication of where the aforementioned processes and balances are active.

Many of the individual modifications are made to ensure logical consistency. For instance, the protection of the bed by bottom-fast ice (Sect. 2.2.1) and limiting erosion to only thawed layers (Sect. 2.2.3) directly follow from the protection of bottom-fast ice from flow-induced melting. Time-dependent thaw depth (Sect. 2.2.2) also becomes necessary due to the fact that bottom-fast ice transfers heat conductively during the winter months and delays the progression of the thaw depth during a model year, reducing erosion even during summer months. The combined effects of the individual modifications described above are what give rise to the form of the simulated deltas with ramp features in ArcDelRCM.jl (Figs. 4b and 5b).

To assess if further insights could be gained through a quantitative measure of network structures, we examine the four graph metrics corresponding to the first and second rows of Fig. 6. The cases with ice- and erosion-related modifications of ArcDelRCM.jl switched on in isolation occupy the range between full ArcDelRCM.jl cases and the DeltaRCM family of cases. The ArcDelRCM.jl cases appear to have more overlap with either the non-Arctic DeltaRCM or the DeltaRCM-Arctic case with E=0.8, depending on the metric. This may be a manifestation of the hybrid nature of the erosion rules of ArcDelRCM.jl (i.e. erosion is unhindered for thawed vertical sediment cells but restricted for frozen vertical cells until the thaw depth deepens in the next time step) versus the more constant fractionally scaled (in flow-speed thresholds or the magnitude of bed diffusion) erosion of DeltaRCM-Arctic, as explained in Sect. 2.2.3. In the context of the ramp feature, the better-connected channels (graph density) and fewer abandoned subnetworks (number of connected components) could translate into a more consistent feeding of sediments to all segments of the delta shore, which would help build an evenly distributed ramp. Future work could apply more comprehensive analyses to further interpret the differences or similarities between these different Arctic delta models, but these are beyond the scope of this ramp-focused work.

In order to gauge the long-term (multi-year) evolution of ramp features in major Arctic deltas, we first address the duration covered by each model year. To this end, Fig. 9 demonstrates the importance of delta activities (and the associated sediment input) outside the peak flooding season. The discharge of large drainage basins such as the Lena watershed remains significant during the summer months (about 53 % of the annual discharge; Holmes et al., 2012), even though it is at a much lower level compared to the peak (GRDC Station Data 2903430, 2018, and Fig. 8). The deltas produced by taking into account the summer months are 4 times the area (and more if one considers only above-water areas) of the equivalent ones that take into account only the peak flooding period. Therefore, for our purpose, we adopt a 4-month model year rather than a 10 d one adopted by Lauzon et al. (2019) and Piliouras et al. (2021). Whether a constant discharge or a time-variable one is used during the peak period does not appear to have an impact on overall areal extent (Fig. 9a and b).

Regarding the ramp feature, it is not only affected by the under-ice depth of the ocean basin but also by the timing of the ice melt. This is demonstrated in Fig. 9c and d. With 15 m of accommodation space under $h_{\mathrm{ice},\max }=\mathrm{2}$ m, deposition that leads to ramp formation is not favoured. However, when an observation-informed timing and duration of bottom-fast ice melt is introduced (Fig. 2), the ramp feature begins to emerge (Fig. 9d). This corresponds to how bottom-fast ice resting on the ramp remains in place whilst the ice in other areas of the delta is flushed or melted away and only starts to break up after the peak flood is over (Fig. 2). However, the deeper ocean and the different discharge pattern led to slower build-up of deposits during ice cover and more deposition during ice-free summers, resulting in the ramp being more hummocky and unevenly graded than in the small-scale, “benchmark” cases in Sect. 3.1.

In reality, the Lena Delta has a much lower sediment volume discharge (roughly a 10th) than we used in our demonstration cases in Fig. 9. The ocean bed on which it formed may not have been flat but rather tilted from the coast towards the Laptev Sea. Taking these factors into account, the simulation in Fig. 10 took 8 times as long to produce a delta of similar size to the one in Fig. 9d but with a ramp feature that is slightly wider and dips closer to 2 m in depth. This is consistent with the aforementioned observation that available depth below ice plays an important role in the formation of ramp features.

Continuing the simulation of Fig. 10 under a strong warming scenario (for discharge and ice cover) akin to RCP7–8.5, we find that an existing ramp feature could degrade on a timescale of centuries (Fig. 11a and b) and effectively disappear within 1 millennium (Fig. 11c and d). Marine processes, such as warming seawater, wave attack, alongshore currents, and sea ice entrainment, may accelerate these timescales. The degradation of the ramp could affect the transport distance of sediments, impacting the release or sequestration of soil organic carbon (Overeem et al., 2022). The reduction in the shallow-water platform provided by the ramp can also impact the delta ecosystems (Lopez et al., 2006). Deltas will also lose a potentially important buffer against coastal erosion (Dean and Dalrymple, 2002).

We note, however, that important ocean-driven processes are missing in the model, resulting in differences in smoothness and outer-edge shapes between the modelled ramps and those observed in reality (Fig. 1). Surges (e.g. during winter storms) can thicken the ice cover (from underneath) that becomes bottom-fast over the ramps (Reimnitz, 2002), which could enhance the protection by bottom-fast ice during spring. Moreover, compared to the deltas produced by DeltaRCM-Arctic and ArcDelRCM.jl, the observed slopes of the sediment bed beyond the outer-edge of the 2 m ramps are much gentler, typically dipping from 2 to >20 m over 𝒪(10) km rather than 𝒪(0.1) km (Reimnitz, 2002; Are et al., 2002). The lack of the gentle dipping may have resulted from the limitation of the models having abrupt thresholds for deposition (Sect. 2.1.4) and in the classification of “on-delta” and “ocean” grid cells during the flow routing (Sect. 2.1.2), in which most of the sediments carried in a packet tend to get deposited as soon as they leave the on-delta cells. The gentle slope observed in reality could be a result of marine processes such as sediment re-suspension by waves (especially during fall storms), which can be effective to water depths exceeding 10 m (Heim et al., 2014), and transport by currents. Wave action re-suspends sediment, especially during fall storms. Sediment in the water column acts as nuclei for frazil ice formation and can be integrated into the forming ice pack (Reimnitz et al., 1992). The formation and export of ice during the fall can remove material from the ramp, as can the anchoring of ice to the seabed in winter and spring (Krumpen et al., 2020). Once anchored, continued cold temperatures create a seasonally frozen layer in the sediment beneath the bottom-fast ice (Osterkamp et al., 1989). This material can be exported once the ice is lifted in spring (Reimnitz et al., 1987). Both processes, entrainment of frazil ice and adherence of frozen sediment, lead to ice rafting (Are and Reimnitz, 2008). Sediment can also be moved by ice gouging and ice bulldozing, when thicker ice masses are ploughed into the seabed (Maznev et al., 2019; Ogorodov et al., 2018). In the context of the simulations shown above, these marine and marine-ice processes could affect not just the exact morphology but also the timescale of ramp formation (balancing between the sediment supply to build and the reworking to sculpt the overall ramp) although not prevent it. The same could, of course, also affect the timescale of degradation of the ramp under a changing climate. Future work on the model could focus on improving the capability of offshore dynamics, such that a full picture of a delta's formation and destruction can be built.

5 Conclusions

The 2 m ramp feature is an integral and ubiquitous feature in Arctic deltas. The morphology and location of the ramp feature means that it could play an important role in the surrounding ecosystems (Lopez et al., 2006), in diffusing wave energy to protect the delta from coastal erosion (Dean and Dalrymple, 2002), and in enhancing carbon sequestration (Overeem et al., 2022). Although it has been suggested that the formation and existence of these ramps may be ice-related (Reimnitz, 2002), we set out to explore from a modelling perspective the conditions and processes that could give rise to the ramps.

To this end, we took the approach of using the explanatory insights (Bokulich, 2013) of RCMs and have written the ArcDelRCM.jl model based on the published descriptions of DeltaRCM (Liang et al., 2015b; Liang, 2015; Perignon, 2018) and DeltaRCM-Arctic (Lauzon et al., 2019; Piliouras et al., 2021). Even though previous work using DeltaRCM-Arctic (Lauzon et al., 2019; Piliouras et al., 2021) has shown tentacle-like offshore deposits at the right elevation range ($\sim -\mathrm{2}$ m) under conditions where the accommodation space between winter ice cover and the ocean bottom is small enough, no continuous (band-like) ramp feature could be produced.

Therefore, we tested a series of physical-rule modifications to explore the origins of ramp formation. The resulting ArcDelRCM.jl contains the following modifications over the base DeltaRCM-Arctic: (i) the protection of bottom-fast ice from flow-induced melting; (ii) the shielding of the bed by bottom-fast ice; (iii) time-dependent thaw depth; (iv) the limiting of all forms of erosion to thawed layers only (instead of using an erodibility factor); (v) the ability for users to specify the time-step size (with internal checks for numerical stability); (vi) the ability to use a time series for input discharge and its related parameters.

We showed that ArcDelRCM.jl can produce the ramp feature. We have found that the modelled ramp feature is indeed related to the winter ice cover, with its depth determined by the maximum thickness of winter ice. We have also found that the prominence of this ramp feature is affected by three factors: (i) the thickness and extent of winter ice, (ii) the available depth in the ocean basin under the winter ice cover (i.e. accommodation space; this confirms previous work by Lauzon et al., 2019 and Piliouras et al., 2021), and (iii) the timing of the melting of bottom-fast ice from atmospheric heat. Simulations of Lena-scale deltas also suggest that a delayed breakup of bottom-fast ice (simulated by a delayed onset of atmospheric melting), which is widespread on the ramp feature, plays an important role in the formation and growth of the ramp feature. The bottom-fast ice protects the ramp from degradation during peak flow.

To further elucidate the processes responsible for the ramp feature and to ensure that the results are not the product of some isolated random events, we ran ensembles of 105 realizations for multiple configurations. These configurations range from the DeltaRCM-Arctic reconstruction with which we began, through individual modifications of ArcDelRCM.jl being enabled in isolation or in groups, to the full ArcDelRCM.jl that includes all the new physical rules.

We found that the protection of bottom-fast ice from flow-induced melting is ultimately responsible for the ability of the model to reproduce the 2 m ramp features. However, the other modifications are not only necessary for internal logical consistency with the first; they also contribute individually to the under-ice deposition pattern and its subsequent preservation required for the ramp feature to form. Specifically, we found that an enhanced constriction leads to more offshore deposits (consistent with Lauzon et al., 2019 and Piliouras et al., 2021) and to more even spreading of deposits along the delta shore as the flow searches for unblocked pathways. Such constriction results from the combination of shielding of the bed from erosion/deposition by bottom-fast ice and the different timing and magnitude of erosion resulting from the updated erosional rules. Acting together, deltas produced by ArcDelRCM.jl appear to have less abandoned parts of the channel network and the channels are better interconnected, as shown in our graph analyses of the ensemble runs.

Taking the model to explore the ramp features' evolution in large Arctic deltas such as the Lena Delta, we first demonstrated that the months outside the peak spring-flood season (and thus the sediment input floods bring) are significant contributors to an Arctic delta's evolution and cannot be neglected by using a 10 d model year (as was adopted for general, smaller deltas in Lauzon et al., 2019; Piliouras et al., 2021). In a set of Lena-like simulations, we found that including summer months (from June to September) instead of limiting the model to the peak-flood period of 10 d led to a quadrupling of the delta area under similar conditions.

When compared to bathymetry data, the simulation produced ramp features that have different elevation smoothness and outer-edge shapes (i.e. the underwater “shorelines”). This may be due to the lack of marine and marine-ice processes in the model and to the clear distinctions between “delta” and “ocean” grid cells, affecting offshore sediment transport. These limitations not only impact the exact morphology produced but may also impact the timescale of formation, growth, or (under climate change) deterioration of the ramps. Future work could focus on addressing these limitations in order to improve the model's capability for predicting the future of Arctic deltas under an increasingly warm climate.

In a sequential pair of simulations (lasting millennia in model time) meant to closely mimic the Lena Delta under present-day and future ice and discharge conditions, we found that a formed ramp feature can degrade and effectively disappear on a timescale of centuries under an extreme climate-warming scenario akin to RCP7–8.5. This timescale could be accelerated further by ocean processes not included in the current model. Such degradation and disappearance of the ramp feature can impact the transport of carbon-carrying sediments, affect the delta ecosystems, and reduce future buffering of Arctic delta shorelines against coastal erosion.

Appendix A: Other feature additions in ArcDelRCM.jl

A1 Inheritance of simulation states

A convenience-motivated feature addition is the ability to start a simulation from an output state given by another simulation. This would allow users to investigate multiple change scenarios that occur after the formation of a delta, such as a sudden increase in discharge or in the duration of spring floods. This has been utilized in the case of Sect. 3.5. It could also be used to break very long simulations into stages to mitigate the risk of a computing-system crash.

A2 Pre-existing island blocks

Users can specify islands in the ocean basin of the simulation domain where no physical processes can occur (similar to the inlet wall). This provides the ability to mimic, albeit simplistically, islands such as Arga in the Lena Delta. Geometries, such as rectangular and elliptical, are available initially; more can be added in the future.

A3 Bed geometries

The initial bed of the ocean basin can have non-uniform depths (i.e. variable h_B). Simple tilt geometries such as linear (from inlet side wall towards the ocean) and radial (from the centre of the inlet channel outwards) are available, in which users can specify the distance over which the ocean bed varies from a specified depth to the h_B specified in the simulation-domain parameters. This has been utilized in Sect. 3.4 and 3.5.

Appendix B: The choice of the initial positive-degree-day index

The simplistic initial starting count of 10 d to calculate the positive-degree-day index, I, is chosen by balancing a few factors that can vary over time and specific sites: (i) comparison between the daily temperatures at the Lena Delta extracted from the ERA-Interim reanalysis data (Dee et al., 2011) (Fig. B1) and the approximate timing of a typical onset of the spring-flood season there (around 1 June); (ii) the measured water temperature in the Lena River showing that it rises above 0 ^∘C near the onset of the spring flood (Juhls et al., 2020); and (iii) the fact that all the simulated cases start when discharge is above a certain threshold (either covering the peak 10 d or when the discharge rises above 20 000 m³ s⁻¹; details in Sect. 3), which occurs some days after the temperature becomes positive. The temperature we used to calculate I is the average daily temperature (Fig. B1) from 1 June to 30 September, which is 4 ^∘C. In addition to the graph metrics of the results presented in Sect. 3.2, Fig. 6 shows the same metrics from the cases identical to the standard ArcDelRCM.jl runs except with an I starting count from 5 and 15 d and (independently) the average positive temperature being 8 and 12 ^∘C. The metrics show that, within the ranges tested, these parameters do not affect the resulting deltas beyond the internal variability in each case. To help illustrate this visually, an example of each case is shown in Fig. B2.

Figure B1 Daily temperature extracted from ERA-Interim reanalysis data (Dee et al., 2011) and averaged over all years from 1951 to 1980 (i.e. the same period from which the average daily discharge in Fig. 8 is obtained). The lighter-shaded ribbon around the average line shows the range of daily temperatures during the same period.

Figure B2 Examples of model outputs that use identical setup and parameters as Fig. 4b, except for one of the following parameters: the positive-degree-day index, I, starts the count at (a) 5 d and (b) 15 d instead of 10; the average positive temperature used to calculate the positive-degree-day index is (c) 8 ^∘C and (d) 12 ^∘C instead of 4 ^∘C.

Appendix C: Ground freeze–thaw pattern

Figure C1 shows the time-variable thaw depth in the time steps around the final winter in the simulations of Fig. 4b and f. The time snapshots shown are from the time step before the start of the final model year through to two time steps after (i.e. from step 4970 to 4973 out of 5000, not counting the 300-step ramp-up phase). The second time snapshot (i.e. step 4971) corresponds to the imposition of maximum ice thickness (h_ice,max) and extent. Since ArcDelRCM.jl forgoes the use of the permafrost label (Sect. 2.2.2 and 2.2.3), direct comparison with distribution maps of permafrost or frozen cells from DeltaRCM-Arctic output is not meaningful. To provide a visual sense of how the ground is frozen on the delta bed, we plot the depth of thaw instead. The thaw-depth patterns provide an indication of where the various processes and balances discussed in Sect. 4 are active. The figure also shows that the thaw depth is directly influenced by winter ice cover, since bottom-fast ice leads to the top layer also being frozen in our model. The resulting step division between deeper and shallower thaw depths does not delineate the ramp feature well, with the shallower thaw depths continuing into the ramp's boundary. The ice-related protective effects are thus active on parts of the ramp itself, as expected.

Figure C1 Evolution of thaw depths in the time steps around the final ice-cover maximum in the simulations in (left column panels) Fig. 4b and (right column panels) Fig. 4f. The time steps shown here are, from top to bottom, step number 4970 to 4973 out of 5000 (not counting the 300-step ramp-up phase). These correspond to the time step before the last winter ice maximum of the simulation through to two steps afterwards. Thaw depths are defined only where there are at least one vertical grid cell of sediments accumulated. Note that these plots show only the thaw depths from the top of any accumulated sediment column, and do not reflect elevation information. The contrast between the left and right columns, corresponding respectively to ice extents of 40 % and 100 %, shows that the thaw depths are directly influenced by the winter ice cover in our model. These thaw-depth patterns provide a visual sense of how the interplay and balances of processes in ArcDelRCM.jl, as discussed in Sect. 4, are spatially distributed.