Sensitivity of isochrones to surface mass balance and dynamics

Abstract

The interior of an ice sheet consists of layers of accumulated snow, which contain important information on accumulation and ice dynamics that are imprinted on layer shapes over time. This work describes how changes in accumulation influence the stratigraphy of an ice sheet. The thickness of each layer at present day depends both on accumulation and on the effect of dynamic thinning after its deposition. An isochronal numerical model is used to simulate the evolution of a 2-D, idealized ice sheet while explicitly representing the layers. A series of simulations was carried out to quantify the changes that anomalous accumulation at different locations and times has on the stratigraphy. These simulations form the basis of a linear response function. A second set of simulations with more sustained changes in accumulation is then used to describe large-scale and long-term impacts on the layering of the ice sheet as well as to test the quality of the linear approximation. The aim is to examine whether long-term effects can be extrapolated from small differential changes. The result confirms a certain degree of linearity between changes in accumulation and layer thickness that may be exploited for future inverse modeling applications.

1 Introduction

The interior of major ice sheets is characterized by a stratigraphic record of distinct layers of equal age, also called isochrones, which, in the field can be identified in several ways, including the extraction of ice cores (e.g. Mojtabavi and others, 2020) and the use of ice penetrating radar (e.g. Legarsky and Gao, 2006; MacGregor and others, 2015; Winter and others, 2019). The present day thickness of these layers is a result of past accumulation and the cumulative effect of ice flow, two processes that in reality are hard to separate. In the past, attempts have been made to describe a direct relation connecting surface mass balance (SMB) to layer thickness, by finding ways to simplify the effect of ice flow. For shallow layers and far from regions of active ice flow, the dynamic effect on layer thinning can be neglected and an average precipitation rate can be found by dividing the current depth of a layer with its age (e.g. Pinglot and others, 2001; Spikes and others, 2004; Medley and others, 2013). A more sophisticated approach introduces a correction factor for the thinning of the layers (Nye, 1963), thus connecting SMB with layer thickness via a natural logarithm (Dansgaard and Johnsen, 1969). This involves a 1-D flow model, neglects horizontal advection and introduces a constant strain rate in order to calculate a mean accumulation rate (Cuffey and Clow, 1997; Leysinger Vieli and others, 2004; Siegert and Payne, 2004; Huybrechts and others, 2009).

These approaches require considerable simplifications of the ice dynamics and are only appropriate for shallow layers. At the same time however, attempts to simulate deeper isochronal surfaces have also been made with thermomechanical ice-sheet models. The inclusion of an Eulerian tracer for the age since deposition is relatively easy to do, but requires the introduction of artificial diffusion that negatively affects the results (Greve, 1997; Born and Robinson, 2021). A Lagrangian approach is also possible (Rybak and Huybrechts, 2003; Sutter and others, 2021) and produces more accurate results, but requires interpolation procedures between gridpoints, which over time accumulate errors and add complexity. To circumvent this complication, Tarasov and Peltier (2003) introduced a semi-Lagrangian transport scheme that back-tracks tracer trajectories onto the Eulerian grid at every time step. This approach may be further improved by using depositional provenance markers instead of individual tracers, because the provenance marker field is relatively smooth in space and therefore less prone to interpolation errors (Clarke and Marshall, 2002; Clarke and others, 2005; Goelles and others, 2014).

Given that complex 3-D models that simulate ice tracer flow exist, the question now is whether these can be used in order to go beyond the 1-D strain rate formula of Dansgaard and Johnsen (1969) and generalize the direct link between layer thickness and SMB for deeper layers as well, while also considering horizontal flow. Baldwin and others (2003) calculated mean accumulation patterns by tracing backward ice particles within a given field of balance velocities. Waddington and others (2007) formulated a least squares inversion problem in order to link directly layer thickness with the smoothest SMB that minimized the misfit between data and model. Comparisons between 3-D and 1-D models show that horizontal ice flow is indeed important and preferable to be taken into consideration (Leysinger Vieli and others, 2011; Nielsen and others, 2015). All these studies featured the limitation of using given velocity fields, meaning that the ice sheet was considered to be in a steady state, with unchanging flow.

This study adds to this body of work by using an isochronal layer model (Born, 2017). The model treats the discretization of the vertical axis not as a grid that is fixed in space, but rather as individual layers corresponding to different times. The explicit simulation of each layer allows for a continuous calculation of ice-sheet thickness and surface slope, meaning that ice flow is constantly updated and the ice sheet is not in a steady state. This means that a perturbation of SMB will trigger a response from ice dynamics and will leave a footprint in the thickness of the isochronal layers. This study focuses on using the isochronal model in order to establish a linear relation between SMB and layer thickness. By doing so, ice flow is parameterized inside this linearization formula and thus changes in layer thickness can then be approximated directly from changes in SMB. This is very useful for a future inversion of the problem, which will allow for the reconstruction of the SMB based on available data of the internal layers of ice sheets. This is not an easy problem, because it requires the use of regularization methods in order to find a least squares solution, however, an additional method for simplifying the relation between SMB and layer thickness is worth analyzing, since the reconstructed SMB will serve as supplementary information on variability in the hydrological cycle, complementing water isotope tracers (e.g. Noon and others, 2003; Lasher and others, 2017).

Section 2 gives a brief description of the model used for this analysis and the simulations performed in order to establish the linearization between SMB and layer thickness. Section 3 puts this relation into test in order to check its accuracy. We find that, within certain limits, the linear reconstruction is representative of the behavior that SMB perturbation has on layer thickness. Section 4 attempts to use this linear function in order to reconstruct SMB for a given layer thickness. Lastly, Section 5 describes the conclusions of this analysis and gives some remarks of how to expand on the linearization in future work.

2 Methods

2.1 Model description

The model used for the present analysis is an isochronal layer model that simulates the evolution and movement of layers of accumulated snow in the interior of an ice sheet (Born, 2017). This means that all variables, including layer thickness, are advected horizontally, but only within each layer and not across isochronal horizons. The vertical dimension consists of the isochronal layers. The model grid consists of layers that are not equidistant in space, but instead they are in time. Every 200 years a new layer is added on top of the ice sheet, thus increasing the computational domain of the model during run time. Due to lateral advection, the layers’ thickness changes with time at each location, and is directly connected to the mass that moves horizontally within the layer. None of the variables are advected vertically through the isochronal grid. Vertical movement is solely due to thinning of older layers below. Horizontal velocities are calculated by using the shallow ice approximation and Glen's flow law. All variables, including layer thickness, are advected using a first order implicit scheme. In the ablation zone, where SMB is negative, the layers melt and their thickness is reduced accordingly. Horizontal velocities depend primarily on the slope of the ice-sheet surface, which is calculated anew at each time step of the model. This means that ice flow is also changing during the simulation and the ice sheet is not in a steady state. The time dependence of the model allows us to emphasize and examine the long-term response of layer thickness to changes in SMB. The model does not represent the firn layer on top of the ice sheet and processes at the ice-sheet bed, e.g. basal freeze-on are not included although they may have a significant effect on the stratigraphy at depth (Leysinger-Vieli and others, 2018). Note that this work prioritizes the effect of boundary conditions, specifically accumulation, in layer thickness and their uncertainty and deliberately omits the also important uncertainty in other model parameters. While both are equally important and influence each other, we choose to separate our work in studies of parameter uncertainty (Born, 2017; Born and Robinson, 2021) and the present study.

For the current study, the model domain is chosen to represent a 2-D cross section of the Greenland ice sheet (GrIS) at the ice divide spanning 133 points laterally with a spacing of 10 km (Fig. 1d). The reference point 0 km will be set at the location of the ice divide. The model has recently been updated to a 3-D version that represents the isochrones of the entire GrIS (Born and Robinson, 2021). This new version improves the numerical efficiency by decoupling the layer tracing scheme from the simulation of ice physics, where the latter are carried out on the much coarser grid of the host model. However, since the 2-D advection equation still has to be solved for every isochrone, the 3-D simulation for the entire Greenland domain is more expensive than the zonal section presented in Born (2017). This simplification does not impact the validity of our results.

Fig. 1. (a) Increased, reduced and parabolic distribution of SMB. All simulations use the parabolic SMB between 0 and 150 000 years. After that, between 150 000 and 200 000 years, CTRL continues with the parabolic, but SMB+, SMB− with the increased and reduced SMB distributions respectively. OSC oscillates between the two with a period of 5000 years. (b) Same as (a) but with the full domain of the SMB, including the melting regions. (c) Evolution of total ice volume in CTRL. (d) The isochronal layers at 160, 170, 180, 190 and 200 (same as the surface) ka for CTRL, SMB+ and SMB− in the ice sheet's domain at 200 ka.

2.2 Description of simulations

Two sets of simulations are conducted which differ in the prescribed SMB and the response of the bedrock to changing ice load (Table 1). The simulations are idealized with a flat bedrock and a constant flow factor of 1.397 × 10⁻²⁴ Pa⁻³ s⁻¹, which corresponds to a temperature of − 5^°C (Cuffey and Paterson, 2010; Born, 2017). We chose a relatively high value to represent the region where most deformation takes place, near the ice-sheet bed. Since temperatures are closer to the melting point near the bedrock (Johnsen and others, 1995), we chose this temperature in order to increase the ice flow of the simulated ice sheet within realistic boundaries. Temperature and therefore ice viscosity are constant everywhere in an effort to limit the number of free parameters. Note that this idealized setup is chosen to approximate an ice sheet with the physical and geographical characteristics of the GrIS, but does not attempt a faithful reproduction of all details. At the upper boundary, we apply an idealized SMB, representing a net effect of accumulation minus melting. The bedrock subsides under the weight of the ice sheet as described by the local lithosphere, relaxing asthenosphere model with a timescale of 10 000 years (Le Meur and Huybrechts, 1996). At the beginning of the simulations, no ice is present, so the bedrock is relaxed. All simulations span 200 000 years so that the vertical isochronal dimension contains 1000 layers at the end of the simulation (a new layer every 200 years). The first 150 ka are common for all simulations and the SMB follows a simple parabolic distribution to reach a steady state before the onset of the SMB anomaly, at which point the ice sheet is not in a steady state anymore, and this is the period that mostly interests us. The SMB is positive between − 180 and 180 km, where precipitation dominates (Fig. 1a), while outside this region the parabolic has negative SMB, and melting dominates (Fig. 1b). After 150 ka, the total ice volume is stable (Fig. 1c) indicating that a steady-state condition has been reached. A control simulation (CTRL) is carried out, in which the reference isochronal state of the ice sheet is established. CTRL is run with constant parabolic SMB forcing for 200 ka.

Table 1. SMB distribution of all simulations on each time period

The first set of simulations aims to establish a sensitivity matrix, comprising parameters that quantify the sensitivity of the thickness of the isochronal layers to very small SMB anomalies. These sensitivity simulations consist of very small perturbations of the SMB around CTRL, both local and short-lived, meaning they only affect one layer at one particular location, a single grid box. Only locations of positive SMB are perturbed, of which there are 40 equidistant ones between − 180 and 180 km. Perturbations are applied for all layers between 150 and 200 ka, layers 751–999 for a total of 9960 simulations (layer 1000 is not considered because it is created as soon as the simulation ends and has no impact on layer thinning). At each location, the perturbation is defined as an increase of SMB by 1/1000th of CTRL at the same location. Since perturbations are applied to all locations and layers, all results are accumulated in a sensitivity matrix. The 9960 sensitivity simulations are used for the calculation of the sensitivity parameter σ, which is defined as the increase in layer thickness caused by the infinitesimal increase of SMB by 1/1000th and is given by:

(1)$$\sigma_{i_0 \cdot N + j_0, i \cdot N + j} = {\partial d_{i, j}\over \partial M_{i_0, j_0}} = {{d_{i, j}-d_{{\rm CTRL}, i, j}}\over ( {1}/{1000}) M_{i_0, j_0}}$$

where d is the thickness of layer j at horizontal location i of the sensitivity simulations, d _CTRL is the thickness of layer j at location i of the CTRL simulation, (1/1000) M is the amount of net SMB perturbation per layer (200 years) at location i ₀ for the layer j ₀ according to the parabolic distribution, and N = 249 are the total amount of layers. σ is essentially a metric of how a perturbation conducted at each of the 40 locations and each of the 249 layers affects the thickness of all layers at all locations. For each one of these 9960 sensitivity simulations, σ contains the normalized values of the thickness difference at all locations and all layers, which are also 9960, making σ a matrix with dimensions 9960 × 9960. If $\sigma _{i_0 \cdot N + j_0, i \cdot N + j}$ is positive then a positive perturbation at i ₀,  j ₀ creates a thicker layer at i,  j for the perturbed simulation over the CTRL, while if the value of σ is negative, then a positive perturbation at i ₀,  j ₀ creates a thinner layer at i,  j.

The second set of simulations consists of more sustained changes in SMB from CTRL. The emphasis is put again on analyzing how these affect the thickness of the layers. The changes in SMB last longer, span larger areas and entail either an increased or a decreased anomaly of the SMB. In a third case the SMB oscillates between these increased and decreased anomalies. The anomalies are limited to one section of the ice sheet, in order to examine if the thickness of the layers is affected only downstream or whether there is an impact on the isochrones across the ice divide as well.

The simulation with increased SMB (SMB+) is based on CTRL for the first 150 ka, followed by the increased SMB distribution between 150 and 200 ka. The increase is centered approximately at 80 km and follows a cosine function, with an amplitude of 0.01 m a⁻¹ (Fig. 1a). In the simulation with reduced SMB (SMB−) the identical spatial and temporal anomaly pattern as in SMB+ is applied, but this time with a negative SMB anomaly of the same size. A third simulation uses an oscillating SMB (OSC) after the spin-up. It oscillates between the increased and reduced SMB distributions with a period of 5 ka, following a sine function so that between 150 and 200 ka there are a total of ten full oscillations. In contrast to SMB+ and SMB− where the total amount of precipitation differs from CTRL, OSC has the same total SMB as CTRL.

Two additional simulations were conducted with SMB distribution identical to the CTRL and SMB+ whose bedrock, instead of being deformable, remains flat for the whole period of the simulations. The simulations were performed in order to differentiate the effects that the bedrock itself has on the internal layers from the effects of dynamics and SMB. These two simulations are only performed in order to explain a very specific signal of the age difference located in the bottom layers of the ice sheet. Since they do not contribute to anything more in the results described below, there is no reason to expand these simulations for the cases of SMB− and OSC.

Lastly, we also examined briefly the impact of sliding (SLID). We ran one more simulation, identical to the CTRL when it comes to the distribution of SMB, but between 170 991 and 171 000 years, we activated a constant sliding at 130 km of 10⁻⁶ m s⁻¹.

3 Results

3.1 Sensitivity simulations

The σ matrix containing the values of the sensitivity parameter can be visualized in two ways: first, by showing the effect that the perturbation of SMB of a particular layer at a particular location has over the whole ice sheet. Second, by showing how the perturbation at all possible locations and layers of the ice sheet affects one particular location and layer. It is important to understand that when a perturbation of SMB occurs at a very specific location and layer, it still affects the layer thickness of all layers at all locations. This is a result of the ice sheet responding to the change of dynamic thinning caused by the small perturbation. As an example, the value of the sensitivity matrix when the perturbation occurs at location 80 km and layer 850, corresponding to an SMB perturbation at 170 ka, is visualized at the end of the simulation (200 ka, Fig. 2a). All the anomalies shown here are caused by a combination of the direct reaction due to the increase in SMB and an indirect due to changes in dynamic thinning. A complex pattern of layer thickening and thinning results from this simple small SMB perturbation, with most changes occurring below the perturbed layer. These anomalies can also be seen in the time domain, i.e. the model grid, which we will use for further discussions (Fig. 2b). Unlike the spatial domain, where results get distorted by different elevation changes in different simulations, the time domain allows for a direct comparison of the isochronal layers between the different simulations.

Fig. 2. (a) Sensitivity parameter of all regions of the ice sheet, as affected by a perturbation at location 80 km and layer 850 (point where dashed lines intersect) at 200 ka. The domain is split into four areas 1, 2, 3 and 4 in order to better explain the phenomena observed. (b) Same as (a) but the y-axis shows each one of the isochronal layers instead of the elevation. (c) Same as (b) but zoomed around the perturbed layer 850.

The direct effect of the increase in SMB is very local and can only be found directly downstream of the perturbed layer (layer 850) (Fig. 2c). The signal of the increased thickness extends to the whole downstream section of the perturbed layer (locations > 80 km), a result of the increased ice mass being transported toward the margins by the flow, thus increasing the thickness of the layer at all downstream locations.

The layers below the perturbed layer (areas 3 and 4 in Fig. 2b) show a dual signal of increased and decreased thickness. The border between the two (located in area 3) is not constant and for younger layers it is moved more to the right than for older layers. The dual signal of this region is a result of the alteration in surface gradient created by the increase in SMB at layer 850. Because of the increase in SMB of the perturbed layer, the total thickness of the ice sheet at location 80 km as well as the gradient downstream, increases. This causes an increase in the downstream velocities that forces more mass transport. Consequently, the layers below the perturbation become thinner due to more ice loss. The negative sensitivity extends on the whole downstream section of the ice sheet up to the right margin. Upstream of location 80 km (area 4) layer thicknesses increase. Because of the local increase in elevation at 80 km, the gradient from the ice divide to 80 km decreases, reducing the horizontal velocities there. The layers there lose less ice mass and remain thicker. The reason why the border separating the two signals is not at the same location for all layers but instead is more to the right for the younger layers is related to the flow of ice. Since horizontal velocity is greater for younger layers, the signal is advected faster and moves faster toward the right margin when compared to the layers closer to the bed.

Layers younger than the perturbation (areas 1 and 2) are also affected and generally show a reduction in ice thickness. This is a result of the increase in surface slope caused by the anomalous accumulation in older layers. The resulting anomalies in dynamic thinning outlast the direct effect of higher accumulation rates and therefore also impact layers that are deposited after the perturbation period. Since the downstream part below the perturbed layer has thinner layers while the upstream part has thicker layers, the slope of the surface is expected to increase. This increases the dynamic thinning and produces thinner younger layers. The ice mass is then advected toward the right margin, slightly increasing the thickness of the layers at this location.

It is interesting to note that the perturbation at 80 km, layer 850, also affects the locations < 0 km with a similar pattern as for the region between the ice divide at 0 km and up to 80 km. This can be explained if one considers the relation between the two regions. Between 0 and 80 km, older layers remain thicker (area 4). Thicker layers on one side of the ice divide mean that more mass is transported to the opposite side to locations < 0 km. As a result, the layers of these locations remain thicker as well. For the younger layers of locations 0–80 km (area 1), which are thinner, the layers accumulate less mass near the ice divide, and thus less mass also flows at locations < 0 km. It becomes clear that a perturbation on one side of the ice sheet, despite how small it might be and despite the fact that it does not create a horizontal movement of the location of the ice divide, still affects the layering of the ice sheet as a whole. Because layer thickness, horizontal velocities and surface slope are all factors that interact at all times, a small change of SMB affects the physical mechanism of ice flow not only around the point of the perturbation, but instead leaves its impact in the whole ice sheet, even including the section at the opposite side.

The sensitivity matrix can also provide information about another question: How do the perturbations at all possible horizontal spatial locations and all layers of the ice sheet affect one particular location and layer? As an example, the effect that perturbations at all locations and layers have on location 80 km and layer 850 are shown in Figure 3. Since only layers younger than 750 were perturbed, the sensitivity to SMB changes before that (1–750) is unknown and therefore not shown here. The vertical axis in Figure 3 is labeled as the perturbed layer, which directly corresponds to the accumulation time when the perturbation occurred. The best way to read this slice of the sensitivity matrix is by associating each entry with the SMB of a certain point in space and time. So for example, at location 100 km and layer 950, we see how a perturbation at 100 km and 190 ka into the simulation affects the thickness at location 80 km, layer 850, at the end of the simulation. Layer 850 increases in thickness by perturbations that occur upstream from location 80 km in the same layer, because these changes are advected toward location 80 km. Perturbations that occur on older layers (751–849) have a primarily negative effect on the layer thickness at 850, due to the general increase in slope. Regarding accumulation that occurs after the deposition of the layer in question (851–999) perturbations that occur upstream of location 80 km cause a negative sensitivity because they increase the effect of dynamic thinning. Perturbations that occur downstream (> 80 km) or on the opposite section of the ice sheet (< 0 km) area cause a positive sensitivity, primarily because of the reduction of the surface gradient that we discussed previously.

Fig. 3. Sensitivity parameter at 80 km and layer 850 (point where dashed lines intersect) as affected by perturbation in all regions of the ice sheet. Layers 1–750 are not shown because no perturbation of the SMB is applied for those layers.

3.2 Increased and reduced SMB

We can now use our understanding of single-point perturbations to analyze variations in SMB with a broader spatial and temporal reach. Simulation SMB+ has the increased SMB curve between 150 and 200  ka. At 155 ka, after 5 ka of anomalous SMB, layers above the 150 ka isochrone (> 750) are thicker between 50 and 110 km due to the higher accumulation of the anomalous forcing (area 5, Fig. 4a). Older layers are not affected directly by an increase in SMB but only indirectly by changes in ice flow (area 6). These layers are thinner because the higher SMB increases the elevation of the ice sheet between 50 and 110 km and steepens the surface gradient downstream. Velocities increase and more ice mass is lost, resulting in thinner layers. The mass that is lost from area 6, gets transported to area 3 and the layers there increase in thickness. This effect is stronger in younger layers that are closer to the surface and therefore subject to higher velocities (areas 2 and top of area 3).

Fig. 4. Simulations SMB+-CTRL. The increased SMB applies inside the area marked by the two vertical dashed lines (50–110 km), and between the dashed horizontal and the ice surface (layer 750–last layer). Thickness difference at (a) 155 ka and (b) 200 ka.

The region from the left margin up to 50 km (areas 1 and 4) also shows an increase in layer thickness. This is explained by the increase in elevation at location 50–110 km, that decreases the slope upstream, thus reducing the velocities there. Smaller velocities mean less ice mass being lost and the layers remain thicker. It is notable that the entire section of the ice sheet of areas 1 and 4 is also affected with a positive thickness difference due to this change in surface slope at 50–110 km, even locations close to the left margin.

At the end of the simulation, at 200 ka, the negative thickness anomaly between SMB+ and CTRL can be seen from 110 km all the way to the right ice-sheet margin (area 3, Fig. 4b). At this point, the negative thickness anomaly from area 6 advects downstream to area 3 and overwhelms the initially positive signal. The same pattern is identified in the sensitivity simulation at 200 ka (area 3, Fig. 2b). In both cases, the older layers, which became thinner as a result of dynamic thinning, advected their signal downstream thus covering the whole section with a negative thickness difference. On the contrary, the younger layers (areas 2 and 5, Fig. 4b) show both an increase and a decrease of their thickness. The two signals coexist because these layers are subject to both increased thickness due to the direct SMB surplus and also thin more quickly due to the dynamic thinning created by younger layers. In addition, the thick layer anomaly is moved downstream by advection. The accumulation anomaly dominates where the cumulative effect of increased flow did not yet have enough time to act, which also explains the boundary between the two regions. It is important to note that the direct effect of the increase in SMB is much more prominent in the case of SMB+ than it was for the sensitivity simulation. In the latter one, the direct effect was only noticeable on the perturbed layer and did not really affect the general thickness of the layers. In the case of SMB+, the fact that SMB is increased for all layers > 750 combined with the increase in the magnitude of SMB anomaly, make the direct effect of SMB change a significant and influential factor for the final shape of the layers. As for areas 1 and 4 of the ice sheet, these preserve their positive thickness difference.

Given the changes that occur in layer thickness due to dynamic thinning, the elevation of the layers above the bedrock is also going to differ between SMB+ and CTRL. The result is a notable difference in the vertical age profiles (Fig. 5a), where the year 0 is defined as the surface layer. At 155 ka there is a column of negative age difference between 50 and 110 km. The negative difference shows the presence of younger layers in SMB+ than in CTRL, at the same depth. This is directly related to the thickness anomaly (Fig. 4a). Since the layers of area 6 are thinner, they occupy deeper depths than in the CTRL experiment and thus have a smaller age difference. This can also be seen on the layer contours (Fig. 1d). The layers of SMB+ at the location where the increased anomaly occurs have sunk when compared to CTRL. Similarly, the section downstream of the 50–110 km column shows a positive age difference, because the thickness of the older layers is larger and thus they occupy shallower depths. The section upstream of the 50–110 km column shows an increase in the age difference, again due to the increase in layer thickness. The distribution of the age of the layers in the ice sheet and the corresponding depth at which each layer is found depends heavily on their thickness.

Fig. 5. Simulations SMB+-CTRL. The increased SMB applies inside the area marked by the two vertical dashed lines (50–110 km), and between the dashed horizontal and the ice surface (layer 750–last layer). Age difference at (a) 155 ka and (b) 200 ka.

One anomaly in the age difference that does not have a direct correspondence in the perturbed layer thicknesses appears in the lower parts of the ice sheet. Figure 5a shows negative age difference located on the deeper layers of the ice sheet. This can be explained by the bedrock deformation in combination with the steep increase in isochrone age near the glacier bed. Since SMB+ has an ice sheet of a larger mass due to the increased SMB, the bedrock is subjected to greater weight. The ice sheet sinks due to the bedrock deformation and the layers move downward, indicating a negative age difference. The slight difference due to the sinking of the bedrock is mostly noticeable on the thinner bottom layers. The same difference on two simulations equivalent to SMB+ and CTRL but with a flat and non-deformable bedrock, shows that no age difference appears at the bottom of the ice sheet (Fig. 6a). With the exception of these very old layers near the bedrock, all other layers have similar age difference in both the solid and the deformable bedrock, indicating that the presence of bedrock deformation does not affect the final thickness of the layers significantly.

Fig. 6. Simulations SMB+-CTRL but with no bedrock deformation. The increased SMB applies inside the area marked by the two vertical dashed lines (50–110 km), and between the dashed horizontal and the ice surface (layer 750–last layer). Age difference at (a) 155 ka and (b) 200 ka.

At 200 ka (Fig. 5b), the section downstream of the 50–110 km column, for elevations lower than 850 m, has a negative age difference, following similar patterns as with area 3 in Fig. 4b. Thinner layers move downward and thicker remain in higher elevation points. Near the bedrock, there is a larger presence of thin layers close to the bottom at 200 ka when compared to 155 ka, because dynamic thinning has acted for a longer period. This explains the reason the negative age difference effect appears much stronger. In the simulation without any bedrock deformation, there is no negative difference on the layers close to the bottom (Fig. 6b).

The differences in layer thickness between SMB− and CTRL at 200 ka follow a pattern that is opposite to SMB+ (Fig. 7). Because of the reduction of SMB, the slope of the ice sheet downstream of 50–110 km decreases. This reduces the velocities and weakens the dynamic thinning. Old layers have increased thickness because of the reduced dynamic thinning. Younger layers experience both an increase and a decrease on their thickness according to whether the decrease in SMB or the reduced dynamic thinning prevails. When the reduction of SMB is stronger, the difference in thickness is negative, thus creating thinner layers.

Fig. 7. Simulations SMB−-CTRL The reduced SMB applies inside the area marked by the two vertical dashed lines (50–110 km), and between the dashed horizontal and the ice surface (layer 750–last layer). Thickness difference at 200 ka.

The inverted pattern of layer thickness difference anomalies in SMB+ and SMB− indicate that even a relatively large SMB anomaly approximately causes a linear response, something that cannot necessarily be expected given the nonlinearity of Glen's flow law which governs dynamic thinning. We will test the linearity further in Section 4.

3.3 Oscillatory SMB

The oscillating SMB anomalies in OSC start with the increased SMB distribution. After 155 ka the SMB distribution completes its first oscillation. Layers 750–775 show one positive and one negative horizontal signal, indicating that the effect of the SMB oscillation is immediately recorded on the thickness of these layers (Fig. 8a). The information is primarily recorded downstream (area 2), but the region upstream (area 1) is also affected and the anomalies are consistent with the early phase of SMB+. The older layers (area 6) are only affected by dynamic thinning and to a lesser degree than SMB+ and SMB− (notice the difference in scale between the figures of OSC and SMB+ and SMB−) because, in the case of OSC, the average SMB anomaly is zero. The thickness difference of area 6 has both negative and positive signals. The negative signal is a result of the initial increased SMB that dynamically thins older layers as seen in SMB+. It gets transported downstream and affects also area 3. The positive signal of area 6 is a result of the following decreased SMB. Since the dynamic thinning occurs fast, it quickly affects all older layers between 50 and 110 km. These rapid changes in the thicknesses of the older layers are observed until the end of the simulation. Area 6 and the part of area 4 around the ice divide change sign according to the different phases of the oscillation. Area 6 shifts from a negative thickness difference when the SMB is increased at 197 ka (Fig. 8b) to a positive one when the SMB is reduced at 200 ka (Fig. 8c). These older layers are affected from the current signal of the SMB distribution and the response is fast. The signal alternates with the change in SMB and is always replaced with the current signal. The history of the oscillations is not present and no memory of previous signals appears. This fast reaction is in contrast with the downstream area 3. There, the anomalous ice thickness of the older layers preserves the same constant signal of both positive and negative thickness difference in both Figures 8b and 8c. The fact that the signal does not change together with the oscillations shows a long-term and slow effect of dynamic thinning. Lastly, at 200 ka, all younger layers (area 2) show the transition between the oscillations and all ten oscillations (50 000/5000) are clearly visible as long horizontal signals. The thickness of the younger layer is affected by the direct change in SMB. These horizontal stripes are visible on the opposite side of the ice sheet (area 1) albeit with much less intensity.

Fig. 8. Simulations OSC-CTRL. The oscillatory SMB applies inside the area marked by the two vertical dashed lines (50–110 km), and between the dashed horizontal and the ice surface (layer 750–last layer). Thickness difference at (a) 155 ka (after one full oscillation), (b) 197 ka and (c) 200 ka.

In summary, the older layers for the two regions downstream and upstream of 750 km behave differently: on the upstream section (areas 6 and 4 around the ice divide), the thickness and age difference of the layers changes with the oscillations in SMB, while on the downstream (area 3) it remains constant. The periodic change of the total thickness difference for layers 1–500 at the upstream section (e.g. at 60 km, Fig. 9a) is of course explained because, after the initial negative offset disappears, the mean of the oscillations is very close to 0, and thus the signal switches from positive to negative values. At the downstream section (e.g. at 180 km, Fig. 9b) the offset appears with a delay since it is advecting from upstream, and the amplitude of the oscillations is not enough for the thickness difference to cross the 0 axis, thus remaining steadily on the negative for the whole examined period of 200 ka.

Fig. 9. Simulations OSC-CTRL. Evolution of the total thickness difference of all layers 1–500 at (a) 60 km and (b) 180 km. The thick line is a 5-point running average.

3.4 Sliding

For the simulation SLID, sliding is activated between 170 991 and 171 000 years at 130 km and is an increased horizontal velocity of all the layers of the ice sheet by 10⁻⁶ m s⁻¹. The relatively short duration of 10 years was chosen analog to the small perturbation in SMB in the sensitivity simulation above. It is long enough to excite a measurable response. The location was chosen in one side of the ice sheet where horizontal ice flow is significant. Immediately after sliding is activated, two columns are formed due to this increase in velocity, one upstream with negative thickness difference between SLID and CTRL and one downstream with positive difference (Fig. 10a). The increased velocity of all layers enhances ice movement at 130 km and these two columns are formed, one that loses mass faster, making the layers thinner in SLID than CTRL, and one that receives this mass, making the layers thicker. At 200 ka, this initial sliding effect does not affect newer layers, while in the old layers, the columns have advected downstream and are not vertical anymore (Fig. 10b). We have the appearance of three signals: an increased layer thickness difference at the lower right near the margin, a decreased layer thickness difference exactly on top and to the left, and again an increased thickness difference on the top and to the left. The first two are the columns of 171 000 years that have advected completely downstream. The third signal is created after 171 000 years and is a result of the two columns. Since the layer thickness of the column upstream of 130 km is decreased, the surface's elevation decreases. Similarly, downstream of 130 km the surface's elevation increases. This creates a flatter surface at approximately 130 km, which reduces the long-term dynamic thinning, meaning that all layers upstream of this perturbation will have more mass. This explains the appearance of the positive third signal.

Fig. 10. Simulations SLID-CTRL. The sliding applies at 170 991–171 000 years, at 130 km (point where dashed lines intersect). Thickness difference at (a) 171 ka and (b) 200 ka.

Since sliding increases ice velocity, its dynamical impact is identical to the increase in velocity due to changes in SMB. The emphasis of this analysis is to find a way to circumvent the effect of dynamic thinning by establishing a linear relation between SMB and layer thickness. As a result, we will not examine sliding further but instead we will isolate changes in SMB as the only factor that affects the stratigraphy of the ice sheet. Yet, at least we could confirm that the impact of sliding gives results that are well aligned and follow similar patterns to the effects of dynamic thinning that we have already covered.

4 Testing the quality of linear approximation

The analysis so far has shown how sustained SMB anomalies may be similar to the short perturbation, but also that non-linearities and asymmetries arise. This section will test the linearity of the three sustained SMB anomalies (SMB+, SMB− and OSC) further. Here we will use the linearized version of our model as represented by the sensitivity matrix. The main question is if a linear superposition of the individual elements of the sensitivity matrix, scaled by the known amount of SMB, captures the thickness anomalies of the SMB+, SMB− and OSC cases, in spite of the known differences in timescales. The reconstructed RecSMB+, RecSMB− and RecOSC which will be shown in this section are not independent simulations, but linear summations of the 9960 sensitivity simulations, scaled by the known values of SMB. The agreement between SMB+-CTRL (ΔSMB+), SMB−-CTRL (ΔSMB−), OSC-CTRL (ΔOSC) and RecSMB+, RecSMB−, RecOSC respectively will be the objective of the following analysis.

We define as $\Delta M_{i_0, j_0}$ the deviation of SMB from CTRL, where i ₀,  j ₀ are the horizontal location and layer where this deviation takes place. As shown in Figure 1a, this occurs for i ₀ between 72 and 77 (the six nodes corresponding to the 50–110 km area, with 10 km spacing) and, since it applies only for the years 150 000–200 000 (with 200 years per layer), j ₀ is between 751 and 999. This change in SMB creates an ice sheet with a new set of layer thicknesses (SMB+, SMB−, OSC), whose thickness difference from CTRL is Δd _i,j, where i,  j are the location and layer of the whole ice sheet. i varies between 0 and 133 (spatial model domain) while j between 1 and 999 (temporal domain). The relation between $\Delta M_{i_0, j_0}$ and Δd _i,j is found from Eqn (1) and the Taylor series for multiple variables. By eliminating all second order derivatives, the series is written as:

(2)$$\Delta d_{i, j} \approx \sum_{i_0, j_0}^{} \bigg({\partial d_{i, j}\over \partial M_{i_0, j_0}} \Delta M_{i_0, j_0}\bigg) = \sum_{i_0, j_0}^{} ( \sigma_{i_0, j_0, i, j} \Delta M_{i_0, j_0}) $$

where the derivative is centered on CTRL, and it can thus be substituted with the sensitivity matrix σ. Note, that Eqn (2) is an approximation. The left-hand side of the approximation represents the thickness difference ΔSMB+, ΔSMB− and ΔOSC. The right-hand side represents a linear reconstruction RecSMB+, RecSMB− and RecOSC.

The anomalous layer thicknesses of RecSMB+ (Fig. 11a) are very similar to ΔSMB+ (Fig. 4b) and differences between the linearized and the full model are relatively small (Fig. 11b), indicating that the linearized reconstruction is a reasonably faithful representation of SMB+. Older layers downstream of the 50–110 km perturbation zone (area 3) show almost no mismatch, thus ΔSMB+ and RecSMB+ are almost identical.

Fig. 11. RecSMB+. (a) Thickness difference at 200 ka. (b) The difference between Figures 4b – 11a. Note that the scale is 5 times smaller than Figure 11a.

The remaining differences can help to understand the shortcomings of the linearization. Newer layers downstream of the SMB perturbation at 50–110 km (area 2) show an inconsistency between ΔSMB+ and RecSMB+. The thickness difference is positive (Fig. 11a), and ΔSMB+ seems to have a stronger positive difference on the upper right corner of the ice sheet and a weaker positive difference on the layers immediately below and to the left (areas 5 and 2, Fig. 11b). This is a result of the increase in surface slope. Because of the continuously increased SMB distribution, the SMB+ simulation creates a glacier of higher elevation than CTRL, increasing the surface slope and thus the effect of dynamic thinning. This is not captured by the linear superposition of sensitivity simulations where the individual SMB perturbations are fully independent, do not add up, and therefore only minimally alter the elevation. Since ΔSMB+ transports ice mass faster toward the right margin, the positive thickness difference is attenuated in comparison with RecSMB+. At the same time, since the same process pushes more mass downstream, a stronger positive signal appears on the upper right corner. In order to quantify the accuracy of the linearization, we define the deviation from linearization metric as:

(3)$${\Sigma \vert RecSMB + \vert - \Sigma\vert \Delta SMB + \vert \over \Sigma\vert \Delta SMB + \vert } \times 100 \% $$

The smaller the number, the more accurate the linearization. We exclude from the computation the gridpoints of the domain located near the margins, because these are prone to computational errors since they are very sensitive to small horizontal movements of the ice sheet. For the rest of the internal layers of the ice sheet, we get a total deviation of 5.22%, meaning that ~95% of the SMB+ has been reconstructed via the linearization.

We can reach similar conclusions by looking at the results of RecSMB−. Comparison of Figures 7 and 12a shows that the reconstructed thickness differences are similar in pattern. The difference between the two gives a quite accurate reconstruction for the older layers because the inconsistency between ΔSMB− and RecSMB− is very low (area 3, Fig. 12b). Since the SMB− simulation has a lower elevation than CTRL, it is expected that the dynamic thinning effect will be weaker in ΔSMB− than in RecSMB−, making the younger layers, which already have a negative thickness difference (Fig. 12a) have slower horizontal velocities for the case of ΔSMB− over CTRL. As a result of the weaker dynamics and horizontal movement, the negative thickness difference is enhanced, and this explains the negative sign in the comparison between ΔSMB− and CTRL (areas 5 and 2, Fig. 12b). In addition, since the negative signal moves less toward the right margin, it makes the layers of the upper right column appear thicker, thus explaining the positive sign at the upper part of area 2. The deviation from linearization of SMB− is found equal to 7.16%.

Fig. 12. RecSMB-. (a) Thickness difference at 200 ka. (b) The difference between Figures 7 – 12a. Note that the scale is 5 times smaller than Figure 12a.

Examining the case of RecOSC (Fig. 13a) we also find a similar pattern to ΔOSC (Fig. 8c). The difference between the two (Fig. 13b) indicates that there is no inconsistency in old layers (area 3) while the new layers (areas 5 and 2) differ. Given the fact that during the oscillations the elevation of the ice sheet continuously alternates from higher to lower than in CTRL, we expect changes in dynamic thinning from stronger to weaker that are not represented in the sensitivity simulations and thus produce alternating mismatches. The deviation from linearization of the reconstruction is only 0.23%, an order of magnitude smaller than the cases of RecSMB+ and RecSMB−. This smaller number shows that the linearization is an even better approximation for the case of OSC. Since OSC has no net SMB difference with CTRL, the change in elevation is very small, effects of differences in the strength of dynamic thinning are less impactful and thus a linearization can be considered more valid. This is the case even though the effect of dynamic thinning is fast and changes sign with every oscillation, as discussed above. However, since these anomalies are also short-lived, their long-term effect is negligible.

Fig. 13. RecOSC. (a) Thickness difference at 200 ka. (b) The difference between Figures 8c – 13a. Note that the scale is 30 times smaller than Figure 11b, indicating a much smaller inconsistency than RecSMB+ and RecSMB−.

In conclusion, the reconstruction based on the linear sensitivities yields good results in these idealized cases. The main differences in layer thickness between SMB+, SMB−, OSC-CTRL were well captured on the reconstructed RecSMB+, RecSMB−, RecOSC respectively, indicating that extrapolating the linearized equation (2) gives an accurate approximation of the dynamical behavior of the ice sheet. This indicates that a linearized parameterization can largely account for the effect of ice flow. The primary reason for disagreements is sustained SMB that gives rise to surface elevation and slope changes and eventually dynamic thinning that is not accounted for by the linearization.

Given that, despite the differences, a linearization of the relationship between $\Delta M_{i_0, j_0}$ and Δd _i,j yields defensible results, we will now briefly explore the possibility of solving an inverse problem and calculating the original SMB. The focus will be on the OSC simulation for two reasons: first because, as shown, the linearity is more accurate for the case of OSC, and second, the fact that SMB changes through time makes OSC a more realistic simulation. The main question is whether the full knowledge of the sensitivity matrix $\sigma _{i_0, j_0, i, j}$ and of the thickness difference Δd _i,j are enough in order to calculate the anomalous SMB, and compare whether the calculation is the same as its actual value. Equation (2) can be rewritten in matrix notation:

$${\boldsymbol D} = {\boldsymbol \sigma} \cdot {\boldsymbol M} ,\;$$

where D is a vector with dimensions (i · j) and contains the values Δd _i,j, σ is a matrix with dimensions (i · j) × (i ₀ · j ₀) and contains the values $\sigma _{i_0, j_0, i, j}$ and M is a vector with dimensions (i ₀ · j ₀) and contains the values $\Delta M_{i_0, j_0}$. Since σ is known, and D is the thickness difference OSC-CTRL, we can transpose for M, the only unknown:

(4)$${\boldsymbol M} = {\boldsymbol \sigma}^{-1} {\boldsymbol D}$$

As mentioned previously, i is between 0 and 133 and j between 1 and 999, while i ₀ is between 72 and 77 (corresponding to 50–110 km) and j ₀ between 751 and 999 (corresponding to 150–200 ka). However, in order for σ to be invertible, it needs to be a square matrix and i · j needs to be the same as i ₀ · j ₀. In order not to have an overdefined problem, with more equations than unknowns, we will solve the system of equations only for i between 72 and 77 and j between 751 and 999.

The comparison of SMB for the case of OSC at location 80 km between the reference (blue line) and the reconstruction (red line) shows that layers 900–999 are reconstructed well but the result for layers 750–900 does not have any physical meaning (Fig. 14). This means that we managed to reconstruct only the last 20 000 years. The presence of noise in the case of the older layers can be explained by examining the nature of the problem that we attempt to solve. The system of Eqns (4) is ill-conditioned, meaning that the solution is too sensitive to error, in this case computational machine error. By doing a simple inversion, this error propagates the solution to unnatural magnitudes, thus producing the noise-like result of Figure 14. In order to filter out the noise, regularization methods are usually applied. Future work will examine an implementation of these methods and how much information about SMB reconstruction can be inferred from isochronal layers.

Fig. 14. Comparison between the actual SMB difference at location 80 km (blue line) and the calculated one from the inverse solution (red line) (Eqn (4)), for the cases of OSC-CTRL.

5 Discussion and conclusion

The current study examined the influence that changes of SMB have on the internal layering in an ice sheet. Formulating this relation is challenging because the effect of dynamic thinning depends on many factors including the thickness of the layers themselves, creating a complex cause–effect relation that is only imperfectly represented when assuming steady state. Previous works have described this relation with some notable assumptions, such as the calculation of a strain rate based on a 1-D model (e.g. Leysinger Vieli and others, 2004; Siegert and Payne, 2004), or when horizontal flow is included, the ice sheet is in a steady state (e.g. Baldwin and others, 2003; Waddington and others, 2007).

In this study, we used the isochronal numerical model by Born (2017) which does not require the ice sheet to be in a steady state. The surface of the ice sheet changes continuously with time, thus ice-sheet dynamics and horizontal flow are always updated according to the SMB of the new layer. SMB affects layer thickness directly via the immediate change in precipitation and indirectly via changes in dynamic thinning. Direct changes of precipitation are more local and affect only the layer created at the time the precipitation occurred, while changes in dynamic thinning have a strong influence on the whole ice sheet and also affect previous and following layers, making this the dominant factor. We examined a range of representative cases, with increased, reduced and seasonal SMB changes, and quantified the differences in layer thickness in each case. An increase in SMB increases the elevation and surface steepness, enhancing the dynamic thinning and creating thinner older layers.

In spite of the strongly non-linear response of dynamic thinning to changes in SMB, a linearized version of our model yielded a satisfactory representation of changes in layer thickness for a given alteration of SMB. The linearization was made by using a sensitivity matrix, a set of parameters calculated by forcing infinitesimal perturbations in the SMB of the model at every location and layer, and then quantifying the sensitivity of the layer thickness at each perturbation. Of the three tested cases, the linearization performed best for the simulation with oscillatory SMB anomaly. Because of the alternating nature of sinusoidal SMB, the net mass balance of the ice sheet at OSC remained closer to CTRL than for the cases of SMB+ and SMB− where the net mass balance was larger and smaller, respectively. Since the average shape of the ice sheet and the surface slope remained more similar, the dynamic behavior was also closer to CTRL, meaning that the linear approximation was more accurate since there was a smaller deviation from the original state of the ice sheet.

To isolate the effect of SMB on layer thickness, some aspects of the model were idealized. The horizontal velocities of the model are calculated using the shallow ice approximation and Glen's flow law. This means that the effect of dynamic thinning essentially depends on the slope of the surface. Vertical advection is not taken into consideration. In addition, since the aim was to focus on the influence of SMB alone, all other factors that could have affected the stratigraphy had to be neglected. Thus, the temperature of the ice sheet was taken equal to −5^°C and consequently the flow factor is constant. We argue that albeit the impact of temperature on ice deformation is substantial, uncertainty in this parameter is smaller than in SMB, in particular considering that the majority of deformation occurs near the base where temperatures are relatively stable. In addition, sliding was not taken directly into consideration for the linearization, but its effect was examined and found to have a very similar dynamical behavior with changes in horizontal velocities due to SMB perturbations. As a result, incorporating sliding in the method is a possibility for future applications. Changes of state for the ice (freezing, melting, etc.) are outside the scope of this research and do have a significant impact on layer stratigraphy (e.g. Leysinger-Vieli and others, 2018). Lastly, a transition to 3-D flow and the presence of a non-flat bedrock are expected to add computational difficulties, but because of horizontal isotropy we expect that our main findings are adaptable to 3-D flow. Overall, within the limitations examined here, we find that the linear approximation is relatively accurate, in particular if carried out for cases where the changes in SMB are not too large.

The linearized approximation offers a new possibility for the reconstruction of past SMB from isochronal layers. Layer thickness data from various sources, notably the radiostratigraphy archives of the ice sheets of Greenland and Antarctica (e.g. MacGregor and others, 2021) that not only cover layers thousands of years old, but are also quite dense with spatial information. By simplifying the relation between SMB and layer thickness and replacing the effect of ice flow with the sensitivity matrix, it becomes theoretically possible to solve for the SMB and get important reconstructions of the polar climate of the past and at different locations of the ice sheets. However, solving for SMB is not as simple as inverting the sensitivity matrix as was shown here in the idealized case of OSC, where by using the layer thickness difference and the 9960 sensitivity simulations we obtained results that featured a lot of noise, an indication that the linear system of equations is ill-posed. This issue can be addressed by introducing regularization processes, which will be examined in future work.