1. Introduction
Ice sheets significantly contribute to the global mean sea level (GMSL) changes which have immediate impacts on the socioeconomic activities especially in coastal regions. Nevertheless, considerable uncertainties remain regarding the future contributions of ice sheets to the GMSL (IPCC, 2021). In particular, because ice sheets are highly sensitive to global climate change, their predicted contributions to the GMSL vary greatly among different emission scenarios (Edwards and others, Reference Edwards2021). Therefore, it is of considerable importance to constrain the uncertainties in the near-future evolution of ice sheets under a warming climate. Nevertheless, while ice sheets can be viewed as an incompressible non-Newtonian fluid (Blatter and others, Reference Blatter, Greve and Abe-Ouchi2010; Kirchner and others, Reference Kirchner2016) and thus can be fully described by the Stokes equations (Greve and Blatter, Reference Greve and Blatter2009), it is difficult to obtain an analytical solution due to the high nonlinearity of ice flow. Accordingly, numerical models are generally used to predict the dynamics of ice sheets. As shown in the Intergovernmental Panel on Climate Change (IPCC) Sixth Annual Report (AR6) (IPCC, 2021), however, notable uncertainties remain in the use of ice-sheet models to predict sea level rise due to limited knowledge of both physical dynamic processes and their numerical implementation.
To predict the dynamic evolution of ice sheets, a number of ice-sheet models have been proposed, including ‘full’ Stokes models (FSMs) and other approximation models of varying numerical complexity. FSMs consider all stress tensor components and tensor gradients and thus can describe the dynamics of ice sheets most accurately (Kirchner and others, Reference Kirchner2016). Unfortunately, FSMs are generally notorious for being computationally expensive therefore it is difficult to use such models to calculate the evolution of ice sheets over a large spatiotemporal domain. Consequently, various methods for approximating FSMs have been proposed (Hindmarsh, Reference Hindmarsh2004; Kirchner and others, Reference Kirchner2016), such as the shallow ice approximation (SIA) (Hutter, Reference Hutter1983; Morland, Reference Morland1984), shallow shelf approximation (SSA) (Morland, Reference Morland1987), DIVA approximation (Goldberg, Reference Goldberg2011), hybrid approximation (SIA+SSA) (Pollard and DeConto, Reference Pollard and DeConto2012) and first-order approximation model (FOM) (Blatter, Reference Blatter1995; Pattyn, Reference Pattyn2003). Note that SIA and SSA are zeroth-order model, and the FOM is one of the higher-order models. More details about other higher-order models were described in Hindmarsh (Reference Hindmarsh2004). Over the past few decades, several numerical ice flow models based on the abovementioned approximation methods have been developed and applied to simplified benchmark studies or real ice-sheet simulations. For example, Bueler and Brown (Reference Bueler and Brown2009) combined the features of the SIA and SSA to capture the fast ice stream features observed in real ice sheets. Moreover, so-called higher-order approximation methods have been widely used to study the dynamic processes in the region where all the components of the stress tensor become equally important, such as the processes occurring near grounding lines, in the transition zones between marine ice sheets and ice streams where the ice flow is complex (Pattyn, Reference Pattyn2003; Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012; Perego and others, Reference Perego, Gunzburger and Burkardt2012; Lipscomb and others, Reference Lipscomb2019). Additionally, Perego and others (Reference Perego, Gunzburger and Burkardt2012) proposed a finite element method (FEM) model for ice sheets based on the FOM. More recently, Hoffman and others (Reference Hoffman2018) developed a new, three-dimensional (3D) variable-resolution thermomechanically coupled FOM, while Wang and others (Reference Wang2018, Reference Wang, Zhang, Xiao, Ren and Wang2020) applied a two-dimensional(2D) first-order ice flow model to study the dynamics of mountain glaciers, and Dias dos Santos and others (Reference Dias dos Santos, Morlighem and Brinkerhoff2022) presented a higher-order ice flow model using a depth-integrated formulation. Moreover, to more accurately simulate the dynamics of the ice sheets, various FSM approaches have also been utilized in different ice-sheet modeling efforts (Jarosch, Reference Jarosch2008; Sargent and Fastook, Reference Sargent and Fastook2010; Seddik and others, Reference Seddik, Greve, Zwinger, Gillet-Chaulet and Gagliardini2012; Leng and others, Reference Leng, Ju, Gunzburger, Price and Ringler2012, Reference Leng, Ju, Gunzburger and Price2014; Zhang and others, Reference Zhang2017; Räss and others, Reference Räss, Licul, Herman, Podladchikov and Suckale2020).
To better understand the differences among different ice-sheet models, several ice-sheet model intercomparisons have been conducted. For example, Le Meur and others (Reference Le Meur, Gagliardini, Zwinger and Ruokolainen2004) compared the SIA and FSM for glacier simulations, Pattyn and others (Reference Pattyn2008) designed numerical experiments to investigate the abilities of higher-order ice flow models, and Pattyn and others (Reference Pattyn2012, Reference Pattyn2013) applied ice flow models of varying numerical complexity to study the evolution of marine ice sheets. However, the above studies did not consider the coupling of ice temperature. Payne and others (Reference Payne2000) and Saito and others (Reference Saito, Abe-Ouchi and Blatter2003, Reference Saito, Abe-Ouchi and Blatter2006) included the thermomechanical coupling process, but they considered only the SIA and FOM approaches, respectively. Zhang and others (Reference Zhang, Ju, Leng, Price and Gunzburger2015) performed a comparison between the thermomechanically coupled 2D FOM and 3D FSM, but under different numerical frameworks. Recently, Rückamp and others (Reference Rückamp, Kleiner and Humbert2022) compared the 3D FOM and FSM behaviors at the Northeast Greenland Ice Stream, and found that basal undulations, basal drags and sliding laws had large and important impacts on model disagreements. Thus, considering the increasing number of 3D thermomechanically coupled FOM applications in ice-sheet modeling, there is an urgent need to further verify the 3D thermomechanical abilities of FOM.
Here, we compare the 3D thermomechanically coupled FOM and FSM under the same numerical framework based on FEM. The details of our numerical model structure can be found in Leng and others (Reference Leng, Ju, Gunzburger, Price and Ringler2012, Reference Leng, Ju, Gunzburger and Price2014). We use the output of the FSM as the ‘true’ (reference) value to analyze the numerical (velocity and temperature) biases of the FOM in transient modeling. In Section 2, we review the governing equations and numerical details of the 3D FSM and FOM. In Section 3, we introduce the experimental design employed in this study. Then, in Section 4, we analyze the simulation results and discuss the possible reasons for the biases of the FOM. Finally, we present our conclusions in Section 5.
2. Numerical model descriptions
2.1. The ‘full’ Stokes flow model (FSM)
The incompressibility of ice can be described by
where u = (u 1, u 2, u 3) represents the ice flow velocity vector and the subscripts 1, 2 and 3 represent the x, y and z directions, respectively, in Cartesian coordinates. The momentum balance equation is given as
where σ is the full stress tensor, ρ is the density of ice and g is the acceleration due to gravity. The velocity and stress tensors of ice are related by the following constitutive equation (Glen, Reference Glen1955),
where τij = −pI + σij(i, j = 1, 2, 3) represents the components of the deviatoric stress tensor; $\dot {{\boldsymbol \epsilon }}$
represents the strain-rate tensor, where $\dot {{\boldsymbol \epsilon }}_{ij} = ( 1 / 2) ( \partial u_{i} / \partial x_{j})$
; p = tr(σii)/3 represents the isotropic ice pressure; and I is the identity matrix. The viscosity of ice can be calculated (Nye, Reference Nye1957; Paterson, Reference Paterson1994) as
where $\dot {\epsilon }_{{\rm e}}$
represents the effective strain rate in the FSM and can be expanded as
the temperature-dependent rate factor A(T) in Eqn (4) is described by the Arrhenius law,
where T is the ice temperature in Kelvin (K), A 0 is the pre-exponential constant, Q is the activation energy and R is the universal gas constant. At the ice surface, we use a stress-free boundary condition,
where n is the outwards unit normal vector at the ice surface. At the lateral surfaces, we apply no-slip boundary conditions.
2.2. The first-order approximation model (FOM)
The momentum balance equations for the first-order approximation of 3D ice flow can be expressed as in Hoffman and others (Reference Hoffman2018),
where x and y are the horizontal coordinate vectors in a Cartesian reference frame, s(x, y) is the ice surface elevation, and $\dot {{\boldsymbol \epsilon }}_{1, 2}$
are given by
Different from the FSM, the effective strain rate $\dot {\epsilon }_{{\rm e}}$
in the FOM is defined as
where $\dot {\epsilon }_{ij} ( i,\; \, j = 1,\; \, 2,\; \, 3)$
are the corresponding strain-rate components. The stress-free boundary condition at the upper surface can be derived from Eqns (3) and (7) as
where n is the outward-pointing normal vector at the ice surface. As in the FSM, the lateral surfaces are not allowed to slide in our experiments. Note that the vertical velocity in the FOM is recovered from the incompressibility of ice flow (Eqn (1)).
2.3. Ice temperature model
The conservation of energy can be described by the evolution of temperature in the following form,
where c represents the heat capacity of ice and k is the heat conductivity (Table 1). At the ice surface, a Dirichlet-type boundary condition is applied, i.e. T = T s, where T s is the ice surface temperature. At the ice base, we use a Neumann-type condition,
where G represents the geothermal heat flux at the ice–bedrock interface. Additionally, for all of our model experiments, the ice temperature is not allowed to exceed the pressure-melting temperature of ice. Note that the temperature diffusion is accounted for in three dimensions.
Table 1. Model parameters used in our experiments
3. Experimental design
The ice-sheet geometry we employ in this study is similar to that in experiment A of the Ice Sheet Model Intercomparison Project for Higher-Order Ice Sheet Models (ISMIP-HOM) (Pattyn and others, Reference Pattyn2008). The ice-sheet surface and basal topographies are given by
where H s and H b are the surface and bed elevations, respectively, α is the slope of the ice sheet, $\bar {H}$
is the mean thickness across the model domain, L is the length of the ice sheet, and M is the amplitude of the sinusoidal fluctuation of the ice base. In all experiments, we set the same upper thermal boundary condition as follows,
where T t is the ice surface temperature at the terminus, H t is the elevation at the ice-sheet terminus and L r is the lapse rate.
In order to show the geometrical conditions of the ice sheet more clearly, a 3D diagram of our model geometry is shown in Figure 1. In addition, the time step in our models is set to 1 year. To obtain an accurate numerical solution, we use a few steps of Picard iteration first to get a good initial state and then switch to Newton iterations, the stopping criterion for the nonlinear iteration is that, either the maximum difference between the velocity solution of the current step and the last step is smaller than 10−4 m a−1, or the difference of pressure solution is smaller than 10−3 Pa, or the iteration reaches the maximum allowed nonlinear steps. More details can be found in Leng and others (Reference Leng, Ju, Gunzburger, Price and Ringler2012, Reference Leng, Ju, Gunzburger and Price2014). Note that in our experiments, FOM uses P1-P1 basis functions and FSM uses P2-P1 basis functions.
Fig. 1. The 3D model geometry used in this study. x and y axis represent horizontal directions, and z is the vertical direction. Basal terrain is perturbed by a sinusoidal function. In the horizontal and vertical direction, we use 21 and 6 gridpoints, respectively.
The details regarding the parameter configuration are shown in Table 1. As mentioned above, we take the FSM results (both velocity and temperature) as the ‘true’ (reference) values and quantify the biases of the FOM results by defining a mean relative error (%) as follows,
where $\bar {r}_{{\rm m}}$
represents the mean relative error, $v_{i}^{fo}$
and $v_{i}^{fs}$
denote the velocity or temperature results for the FOM and FSM, respectively, and n is the number of degrees of freedom that we use for the calculations.
Here, we test two sets of experiments, one is the thermomechanically decoupled experiments (see Appendix A1), where the rate factor A is set to be constant (Table 1) and the other is thermomechanically coupled experiments, where we test the sensitivity of the thermomechanically coupled FOM and FSM to various parameters, i.e. the ice slope, geothermal heat flux, basal topography and sliding condition. Furthermore, we assume a geometric steady-state for all of our model experiments, i.e. fix the ice geometry in time to simplify the model experiments. The ice temperature fields for the FOM and FSM are initialized in a purely diffusive manner using the same initial thermal boundary conditions. In this way, we guarantee that the only factor impacting the transient evolution is the difference between the thermomechanical physics and model complexities of the FOM and FSM, a case that is often encountered in real simulations of transient ice sheets and glaciers after initializing the model (when it is difficult to achieve a true thermodynamic steady-state).
In our thermomechanically coupled experiments, we prescribe the basal topography as exhibiting sinusoidal perturbations. For the sinusoidal ice base, the mean thickness $\bar {H}$
is set to 1000 m, and we change the fluctuation amplitude M from 300 to 500 m at an interval of 100 m with three different ice slopes (0.02, 0.03 and 0.04). To study the impacts of the thermal status on the velocity field of the ice sheet, we adopt three different geothermal heat fluxes (10, 15 and 40 mW m−2) at the ice–bedrock interface. The surface temperature can be obtained from automatic weather stations (AWS) (Fausto and others, Reference Fausto2021) or some reanalyses datasets (Delhasse and others, Reference Delhasse2020), whereas the geothermal heat flux is difficult to directly obtain and large uncertainties remain, a reason we focus on the sensitivity of the geothermal heat flux in our numerical experiments. In addition, to discuss the effect of basal sliding, we conduct additional sliding experiments with a linear friction law,
where τb is basal drag, β2 is friction coefficient (Pa a m−1), ub is basal velocities, and we set β2 as:
Note that in ISMIP-HOM (Pattyn and others, Reference Pattyn2008), the biases between the FSM and non-full-Stokes models (NFSM) decrease with the extent of ice sheet. Thus, for all of numerical experiments in this study, the ice sheet length and width were fixed as 10 km (Table 1), which showed relatively larger model biases in ISMIP-HOM, in order to further our studies for the thermomechanically coupled behavior in the FOM and FSM. Furthermore, in Appendix A2, we discuss the details of the differences of strain rate and strain heating between the FOM and FSM. The details of coupled experimental parameters are shown in Table 2 and the thermomechanically decoupled experiments are discussed in Appendix A1.
Table 2. Configurations of our numerical thermomechanically coupled experiments. C1, C2 and C3 represent the experiments performed to test the sensitivity of the model to the ice slope, geothermal heat flux and basal topography, respectively.
EXP-C1~C3 are conducted with a frozen ice base, and EXP-S1 is performed to test the impact of basal sliding condition
4. Results and discussion
4.1. Impacts of the ice slope (EXP-C1)
In this experiment, we use a geothermal heat flux G of 15 mW m−2 and a basal fluctuation amplitude M of 500 m with three different slopes α of 0.02, 0.03 and 0.04 (Table 2). As shown in Figure 2, as the slope increases, the ice-sheet flows faster due to the increased gravitational driving force. It should be noted that all velocities in this paper are referred to by the mean magnitude of the model outputs (modulus of u 1 and u 2). In addition, the differences in the mean velocity (hereafter referred to as MVD) and mean temperature in the domain (hereafter referred to as MTD) between the FSM and FOM after 100 years of evolution vary among three slopes. Nevertheless, the distribution of the relative error of depth-averaged velocity (hereafter referred to as DAV) is in good agreement with relative error of the depth-averaged temperature (hereafter referred to as DAT), and a clear pattern emerges: as the slope increases, the velocity and temperature biases of the FOM both increase at year 100.
Fig. 2. Plan views of the depth-averaged relative error distributions for the velocity and temperature fields at year 100 with three different ice slopes of 0.02 (top), 0.03 (middle) and 0.04 (bottom), respectively: (left, a, d, g) modeled horizontal velocities simulated by the FOM (m a−1); (middle, b, e, h) relative velocity errors between the FOM and FSM (%); (right, c, f, i) relative temperature errors between the FOM and FSM (%).
In addition, the relative errors of the MVD and MTD increase at each time step as the slope increases. When α = 0.02 (blue dotted lines in Figs 3a, b), the relative errors of the MVD and MTD are ~5.6 and 0.07%, respectively, at year 100, whereas when α = 0.04 (red solid lines in Figs 3a, b), the relative errors of the MVD and MTD increase to ~6.7 and 0.15%, respectively. The trends in the relative errors of MVD and MTD over time are approximately identical; that is, as the ice slope decreases, the error curves change relatively smoothly, but as the ice slope increases, the change rate of mean relative errors in the FOM exhibits fluctuations (Fig. 3).
Fig. 3. Time series of the mean relative errors in velocity (a) and temperature (b) with different ice slopes in the domain.
Moreover, Figure 2 shows that the changes of the relative error in DAV are generally consistent with the changes in the DAT. Both the DAV and DAT in the FOM seem to be overestimated for higher velocities but underestimated for smaller velocities. Thus, we argue that the thermomechanical coupling of velocity and temperature can cause the model biases generated at each time step to accumulate over time. More detailed discussion about the difference of strain rate and strain heating can be found in Appendix A2.
4.2. Impacts of the geothermal heat flux (EXP-C2)
As shown in Table 2, we use an ice slope α of 0.02 and a basal fluctuation amplitude M of 500 m to test the sensitivity of the FOM model to different geothermal heat fluxes G of 10, 15 and a more realistic value of 40 mW m−2. By changing G, we alter the thermal status of our model domain. After running the model for 100 years, the DAV of the FOM increases with increasing geothermal heat flux (Figs 4a, d, g). In addition, the relative errors of both the DAV and DAT increase when we apply a larger geothermal heat flux (Fig. 4). Additionally, with the increasing geothermal heat flux, the relative error of MVD shows an increasing trend: when the geothermal heat flux is 40 mW m−2, the relative error of MVD reaches ~6.5% at year 100 (red solid line in Fig. 5a). Nevertheless, when the geothermal heat flux is as low as 10 mW m−2, the relative error of MVD is only ~4.8% and remains almost proportional over the 100-year period (blue dashed line in Fig. 5a).
Fig. 4. Plan views of the depth-averaged relative error distributions for the velocity and temperature fields at year 100 with three different geothermal heat fluxes of 10 (top), 15 (middle) and 40 mW m−2 (bottom): (left, a, d, g) modeled horizontal velocities simulated by the FOM (m a−1); (middle, b, e, h) relative velocity errors between the FOM and FSM (%); (right, c, f, i) relative temperature errors between the FOM and FSM (%).
Fig. 5. Time series of the mean relative errors in velocity (a) and temperature (b) with different geothermal heat fluxes in the domain.
Furthermore, we detect a similar trend for the MTD. As the geothermal heat flux increases, the mean relative error of the MTD between the FSM and FOM increases accordingly. When the geothermal flux is 10 mW m−2, the relative error of MTD changes only slightly over time (only ~0.02% at year 100, blue dashed line in Fig. 5b). In contrast, when the geothermal heat flux reaches 40 mW m−2, the relative error grows from ~0.03 to 0.12% over the 100-year period (red solid line in Fig. 5b).
Similar to the results for EXP-C1, Figure 4 clearly reveals that the distributions of relative errors of DAV and DAT are closely related to the FOM velocity. That is, the mean relative errors are overestimated in the FOM regions with larger ice velocities and underestimated in the regions with smaller velocities.
4.3. Impacts of basal topography (EXP-C3)
The sensitivity of the model to basal topography is determined by changing the amplitude M of the sinusoidal fluctuation at the ice base. Here, we use M values of 300, 400 and 500 m while keeping an ice slope α of 0.02 and a geothermal heat flux G of 15 mW m−2 (Table 2). As the amplitude increases, the basal topography becomes increasingly undulatory; at the same time, the concavity of the ice sheet (the area stretching from the upper left corner to the lower right corner in Figs 6a, d, g) deepens as the ice sheet thickens, increasing the FOM velocity and consequently increasing the mean relative errors in the velocity and temperature.
Fig. 6. Plan views of the depth-averaged relative error distributions for the velocity and temperature fields at year 100 with three different basal fluctuation amplitudes of 300 (top), 400 (middle) and 500 m (bottom): (left, a, d, g) modeled horizontal velocities simulated by the FOM (m a−1); (middle, b, e, h) relative velocity errors between the FOM and FSM (%); (right, c, f, i) relative temperature errors between the FOM and FSM (%).
As shown in Figure 7, for different basal fluctuation amplitudes, the mean relative errors of both the MVD and MTD exhibit similar trends over time, with both the relative errors increasing over the 100-year period as the basal fluctuation amplitude grows. When M is 300 m (blue dotted lines in Figs 7a, b), the relative errors of MVD and MTD are ~3.3 and 0.06%, respectively, at year 100, whereas the corresponding errors increase to ~5.7 and 0.07% when M is 500 m (red solid lines in Figs 7 a, b). Hence, under the influence of complex basal terrain, the FOM may induce considerable errors compared to the FSM. Furthermore, similar to the results of EXP-C1 and EXP-C2, Figure 6 shows that the FOM overestimates the DAV and DAT in regions with larger velocities.
Fig. 7. Time series of the mean relative errors in the depth-averaged velocity (a) and temperature (b) with different basal fluctuation amplitudes in the domain.
4.4. Impacts of sliding boundary condition (EXP-S1)
It is noted that the in all above-mentioned experiments we use a no-slip basal boundary condition. In this subsection we add two additional experiments to discuss the impacts of basal boundary condition on MVD and MTD. We set ice slope α = 0.02, geothermal flux G = 15 mW m−2 and amplitude M = 500 m (Table 2). As ice slides at the base, the FOM velocity increases significantly compared to the case of frozen bed, and the relative error of MVD and MTD increases accordingly.
Similar to EXP-C1, C2, C3, the same patterns appear under the sliding condition, i.e. we find larger relative errors of DAV and DAT where the FOM velocity is higher (Fig. 8). From Figure 9, when ice slides at the base, the relative error of MVD and MTD becomes ~22.5% (red solid line in Fig. 9a) and 0.5% (blue dotted lines in Fig. 9a), respectively. However, for the frozen-bed condition, the corresponding values are approximately only 5.2% (red solid line in Fig. 9b) and 0.06% (blue dotted lines in Fig. 9b).
Fig. 8. Plan views of the depth-averaged relative error distributions for the velocity and temperature fields at year 100 with two different boundary conditions of sliding (top), and frozen one (bottom): (left, a, d) modeled horizontal velocities simulated by the FOM (m a−1); (middle, b, e) relative velocity errors between the FOM and FSM (%); (right, c, f) relative temperature errors between the FOM and FSM (%).
Fig. 9. Time series of the mean relative errors in the depth-averaged velocity (a) and temperature (b) with different basal boundary conditions in the domain.
5. Conclusions
To achieve a compromise between accuracy and efficiency, the 3D FOM has been widely used in ice-sheet modeling. To quantify the possible causes and impacts of the biases in the FOM, in this paper, we compare its 3D thermomechanically coupled behaviors with those of a 3D FSM.
To do so, we test the sensitivity of both ice-sheet models to four boundary conditions: the ice slope, geothermal heat flux, the fluctuation amplitude of the topography and sliding condition at the ice base. By changing these conditions, we analyze their impacts on the mean relative errors in the FOM velocity and temperature fields compared to the values for the FSM. All experiments are conducted under the same numerical framework based on the FEM to minimize errors caused by other sources (e.g. the discretization method and meshing approaches).
In our thermomechanical experiments, the distributions of the relative errors in MVD and MTD generally show good consistency. Compared to the FSM, the temperature bias increases as the ice slope, geothermal heat flux, basal fluctuation amplitude increase. In addition, when ice slides at the base, the relative error of MVD and MTD will increase significantly. Furthermore, we find that the FOM bias tends to be overestimated in regions with larger velocities and underestimated in regions with smaller velocities. Finally, the relative error between the FOM and FSM increases over time, indicating that the error accumulates as the FOM runs. Therefore, modelers must be wary of the results running a thermomechanically coupled FOM for a long time, especially in regions with complex ice geometries. Thus, on the basis of a series of numerical experiments, we suggest that the FOM should cautiously be used for dynamic simulations in regions with very complex ice-sheet geometries and large aspect ratios.
It should be noted that all of our numerical experiments are designed under simplified geometric and model conditions, whereas the real dynamics of ice sheets are much more complicated. Our experiments are limited to idealized geometries, and the conclusions here are not generalized but show some possible scenarios that may cause increased biases between the FOM and FSM. Further efforts should be dedicated in the near future to more realistic comparisons between the FOM and FSM. In addition, model uncertainties still remain if ice surface geometry evolves, due to different approaches of calculating vertical velocities for FOM and FSM. The FOM, considered as ‘posteriori model’ (Pattyn, Reference Pattyn2003; Pattyn and others, Reference Pattyn2008), calculates vertical velocity from the mass incompressibility condition, while in FSM the vertical velocity is part of its own solution. It is still unknown exactly how this difference impacts glacier dynamics during a prognostic simulation with an evolving ice geometry.
Acknowledgements
This work was supported by the National Key Research and Development Program of China (No. 2018YFC1406104), the State Key Laboratory of Earth Surface Processes and Resource Ecology (2021-TS-06, 2021-KF-06, 2022-ZD-05) and the Fundamental Research Funds for the Central Universities (2021NTST16). We thank Fuyuki Saito and an anonymous reviewer for their helpful comments. We also thank the editor Frank Pattyn for his constructive inputs.
Data availability
The FOM and FSM results in this study can be downloaded at https://doi.org/10.5281/zenodo.6523531.
Appendix A
A.1. Appendix A1: thernomechanically decoupled experiments
This section mainly introduces the thermomechanically decoupled experiments. In the thermomechanically decoupled experiments, we set the rate factor A to be constant (Table 1); because the whole ice sheet is viewed as an isothermal body, the rate factor does not change with the ice temperature. In addition, the mean velocity is kept in a steady-state for a fixed ice-sheet geometry and a constant rate factor, resulting in a constant relative error at each time step. The geometric parameters in the thermomechanically decoupled experiments are the same as those in the thermomechanically coupled experiments (Table 3). However, since the thermomechanically decoupled experiments do not include the evolution of temperature, we do not consider the impacts of heat flux, and we present only the relative error of MVD. At the lower and lateral surfaces, we apply no-slip boundary conditions. The corresponding results are shown in Table 4.
Table 3. Configurations of the thermomechanically decoupled experiments
Table 4. Relative errors of the mean velocity in the domain in the thermomechanically decoupled experiments
By comparing the relative errors under different geometric conditions (Table 4), we find that the ice slope does not obviously impact the relative error of MVD. That is, under different ice slopes of 0.02, 0.03 and 0.04, the relative errors remain almost the same, i.e. ~4.5%. In contrast, the basal topography seems to have a greater impact. For all three different slopes, the mean relative error increases to ~0.3% when the basal fluctuation amplitude increases from 300 to 400 m and to ~0.6% when basal fluctuation further increases from 400 to 500 m (Table 4). The maximum relative error of 4.526% occurs at a slope of 0.02 and an amplitude of 500 m.
In addition, we calculate the ice temperature in a decoupling manner with Eqn (12) and plot the FOM temperature distribution for three different slopes in Figure 10. As the slope increases, the FOM velocity becomes larger, which brings more cold ice from upstream to downstream, and accordingly produces more strain heat, inducing more temperate ice gather at the terminal by advection. Therefore, in this case we see the upstream becomes colder while the downstream becomes warmer (Fig. 10).
Fig. 10. Plan views of the depth-averaged temperature distributions of the FOM at year 100 with three different ice slopes conditions of 0.02 (left), 0.03 (middle) and 0.04 (right) (K).
A.2. Appendix A2: strain rate and strain heating
In this section we focus on the reasons for the differences in velocity and temperature fields between the FOM and FSM. We investigate the absolute difference among strain rate $\dot {\epsilon }_{11}$
, $\dot {\epsilon }_{22}$
, $\dot {\epsilon }_{13}$
and $\dot {\epsilon }_{23}$
(Fig. 11) and corresponding error of strain heating E (Fig. 12) between the FOM and FSM at different slopes (EXP-C1). The strain heating can be calculated as (Greve and Blatter, Reference Greve and Blatter2009),
Fig. 11. Plan views of the difference of depth-averaged strain rate distributions between the FOM and FSM at year 100 with three different ice slopes conditions of 0.02 (top), 0.03 (middle) and 0.04 (bottom): (1st column, a, e, i) depth-averaged difference of $\dot {\epsilon }_{11}$
; (2nd column, b, f, j) depth-averaged difference of $\dot {\epsilon }_{22}$
; (3rd column, c, g, k) depth-averaged difference of $\dot {\epsilon }_{13}$
; (4th column, d, h, l) depth-averaged difference of $\dot {\epsilon }_{23}$
(a−1).
Fig. 12. Plan views of the difference of depth-averaged strain heating distributions between the FOM and FSM at year 100 with three different ice slopes conditions of 0.02 (left), 0.03 (middle) and 0.04 (right) (%).
It can be found that the difference of strain rate between the FOM and FSM increases as the slope increases (Fig. 11) – this pattern is especially obvious for $\dot {\epsilon }_{13}$
(Figs 11 c, g, k) and $\dot {\epsilon }_{23}$
(Figs 11 d, h, l). The difference of $\dot {\epsilon }_{13}$
between the FOM and FSM is approximately between 0 and 56 a−1 at the ice slope of 0.04 (Fig. 11 k). However, when ice slope decreases to 0.02, the corresponding values decrease to 0 and 5 a−1 (Fig. 11 c). The difference of $\dot {\epsilon }_{23}$
between the FOM and FSM is approximately between 0 and 35 a−1 at the ice slope of 0.04, while when ice slope decreases to 0.02, corresponding values decrease to 0 and 6 a−1, respectively (Fig. 11 d, l). In addition, the spatial distribution of the difference of strain rate is similar to DAV of the FOM as well as the relative error of DAV and DAT, i.e. the overestimation in strain rate tends to occur where FOM velocity is larger, similar to the case in the strain heating distribution (Fig. 12). The depth-averaged relative error of strain heating between the FOM and FSM also increases with ice slope, which is ~50% across most of the ice domain when ice slope is 0.04. However, when ice slope decreases to 0.02, the relative error of strain heating between the FOM and FSM decreases to ~10%.
Therefore, we argue that there may be a positive feedback in the differences between FOM and FSM in our experiments. As the slope increases, the ice velocity and its bias between the FOM and FSM become larger, leading to more biased strain rate, strain heating and ice temperature, which in turn causes even larger biases of ice velocity field in time.
Open access
