A comparison between three-dimensional, transient, thermomechanically coupled first-order and Stokes ice flow models

Abstract

In this study, we investigate the differences between two transient, three-dimensional, thermomechanically coupled ice-sheet models, namely, a first-order approximation model (FOM) and a ‘full’ Stokes ice-sheet model (FSM) under the same numerical framework. For all numerical experiments, we take the FSM outputs as the reference values and calculate the mean relative errors in the velocity and temperature fields for the FOM over 100 years. Four different boundary conditions (ice slope, geothermal heat flux, basal topography and basal sliding) are tested, and by changing these parameters, we verify the thermomechanical behavior of the FOM and discover that the velocity and temperature biases of the FOM generally increase with increases in the ice slope, geothermal heat flux, undulation amplitude of the ice base, and with the existence of basal sliding. In addition, the model difference between the FOM and FSM may accumulate over time, and the spatial distribution patterns of the relative velocity and temperature errors are in good agreement.

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, 2021). 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, 2010; Kirchner and others, 2016) and thus can be fully described by the Stokes equations (Greve and Blatter, 2009), 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, 2016). 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, 2004; Kirchner and others, 2016), such as the shallow ice approximation (SIA) (Hutter, 1983; Morland, 1984), shallow shelf approximation (SSA) (Morland, 1987), DIVA approximation (Goldberg, 2011), hybrid approximation (SIA+SSA) (Pollard and DeConto, 2012) and first-order approximation model (FOM) (Blatter, 1995; Pattyn, 2003). 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 (2004). 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 (2009) 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, 2003; Larour and others, 2012; Perego and others, 2012; Lipscomb and others, 2019). Additionally, Perego and others (2012) proposed a finite element method (FEM) model for ice sheets based on the FOM. More recently, Hoffman and others (2018) developed a new, three-dimensional (3D) variable-resolution thermomechanically coupled FOM, while Wang and others (2018, 2020) applied a two-dimensional(2D) first-order ice flow model to study the dynamics of mountain glaciers, and Dias dos Santos and others (2022) 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, 2008; Sargent and Fastook, 2010; Seddik and others, 2012; Leng and others, 2012, 2014; Zhang and others, 2017; Räss and others, 2020).

To better understand the differences among different ice-sheet models, several ice-sheet model intercomparisons have been conducted. For example, Le Meur and others (2004) compared the SIA and FSM for glacier simulations, Pattyn and others (2008) designed numerical experiments to investigate the abilities of higher-order ice flow models, and Pattyn and others (2012, 2013) 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 (2000) and Saito and others (2003, 2006) included the thermomechanical coupling process, but they considered only the SIA and FOM approaches, respectively. Zhang and others (2015) performed a comparison between the thermomechanically coupled 2D FOM and 3D FSM, but under different numerical frameworks. Recently, Rückamp and others (2022) 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 (2012, 2014). 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

(1)$$\nabla \cdot {\boldsymbol u} = 0,\; $$

where u = (u ₁,  u ₂,  u ₃) 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

(2)$$\nabla \cdot {\boldsymbol \sigma} + \rho {\boldsymbol g} = 0,\; $$

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, 1955),

(3)$${\boldsymbol \tau} = 2 \eta \dot{{\boldsymbol \epsilon}},\; $$

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, 1957; Paterson, 1994) as

(4)$$\eta = {1\over 2} A^{-{1\over n}} \dot{\epsilon}_{{\rm e}}^{\left({1\over n}-1\right)},\; $$

where $\dot {\epsilon }_{{\rm e}}$ represents the effective strain rate in the FSM and can be expanded as

(5)$$\eqalign{\dot{\epsilon}_{{\rm e}} & = \sqrt{{1\over 2} \displaystyle\sum_{i j} \dot{\epsilon}_{i j}^{2}} \cr & = \sqrt{{1\over 2}\left(\left(\dot{{\epsilon}}_{11}\right)^{2} + \left(\dot{\epsilon}_{22}\right)^{2} + \left(\dot{\epsilon}_{33}\right)^{2} + 2\left(\dot{\epsilon}_{12}\right)^{2} + 2\left(\dot{\epsilon}_{23}\right)^{2} + 2\left(\epsilon_{13}\right)^{2}\right)},\; }$$

the temperature-dependent rate factor A(T) in Eqn (4) is described by the Arrhenius law,

(6)$$A( T) = A_{0} {\rm e}^{-Q/RT},\; $$

where T is the ice temperature in Kelvin (K), A ₀ 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,

(7)$${\boldsymbol \sigma} \cdot {\boldsymbol n} = 0,\; $$

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 (2018),

(8)$$\matrix{& -\nabla \cdot ( 2 \eta \dot{{\boldsymbol \epsilon}}_1) + \rho g{\partial H_{{\rm s}}\over \partial x} = 0,\; \cr & - \nabla \cdot ( 2 \eta \dot{{\boldsymbol \epsilon}}_2) + \rho g{\partial H_{{\rm s}}\over \partial y} = 0,\; \cr }$$

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

(9)$$\matrix{ \dot{{\boldsymbol \epsilon}}_1 = \left(\matrix{ 2\dot{\epsilon}_{11} + \dot{\epsilon}_{22},\; & \dot{\epsilon}_{12},\; & \dot{\epsilon}_{13}}\right)^{T},\; \cr \dot{{\boldsymbol \epsilon}}_2 = \left(\matrix{\dot{\epsilon}_{12},\; & \dot{\epsilon}_{11} + 2\dot{\epsilon}_{22},\; & \dot{\epsilon}_{23} }\right)^{T}. \cr }$$

Different from the FSM, the effective strain rate $\dot {\epsilon }_{{\rm e}}$ in the FOM is defined as

(10)$$\dot{\epsilon}_{\rm e} \equiv \left(\dot{\epsilon}_{11}^{2} + \dot{\epsilon}_{22}^{2} + \dot{\epsilon}_{11} \dot{\epsilon}_{22} + \dot{\epsilon}_{12}^{2} + \dot{\epsilon}_{13}^{2} + \dot{\epsilon}_{23}^{2} \right)^{{1 \over 2}},\; $$

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

(11)$$\dot{{\boldsymbol \epsilon}}_1 \cdot {\bf n} = \dot{{\boldsymbol \epsilon}}_2 \cdot {\bf n} = 0,\; $$

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,

(12)$$\rho c\left({\partial T\over \partial t} + {\boldsymbol u} \cdot \nabla T\right) = k \nabla^{2} T + 2 \eta \dot{{\boldsymbol \epsilon}}\, \colon \, \dot{{\boldsymbol \epsilon}},\; $$

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,

(13)$$G = -k{\partial T\over \partial z},\; $$

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, 2008). The ice-sheet surface and basal topographies are given by

(14)$$\eqalign{& H_{{\rm s}} = 3000-x \times \alpha \times L,\; \cr & H_{{\rm b}} = H_{{\rm s}}-\bar{H} + M \times \sin ( 2 \pi \times x) \times \sin ( 2 \pi \times y) ,\; }$$

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,

(15)$${{T}_{{\rm s}}} = T_{\rm t} + ( H_{\rm s}-H_{\rm t}) \times L_{\rm r},\; $$

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⁻⁴ m a⁻¹, or the difference of pressure solution is smaller than 10⁻³ Pa, or the iteration reaches the maximum allowed nonlinear steps. More details can be found in Leng and others (2012, 2014). 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,

(16)$$\bar{r}_{{\rm m}} = {1\over n} \sum_{i = 1}^{n} {\left\vert v_{i}^{\,f o}-v_{i}^{\,f s}\right\vert \over \left\vert v_{i}^{\,f s}\right\vert } \times 100,\; $$

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⁻²) at the ice–bedrock interface. The surface temperature can be obtained from automatic weather stations (AWS) (Fausto and others, 2021) or some reanalyses datasets (Delhasse and others, 2020), 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,

(17)$${\boldsymbol \tau_{{\rm b}}} = \beta^{2} {\bf u_{{\rm b}}},\; $$

where τ_b is basal drag, β² is friction coefficient (Pa a m⁻¹), u_b is basal velocities, and we set β² as:

(18)$$\beta^{2} = 1000 + 1000 \times \sin( 2 \pi \times x) \times sin( 2 \pi \times y) .$$

Note that in ISMIP-HOM (Pattyn and others, 2008), 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⁻² 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 ₁ and u ₂). 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⁻¹); (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⁻². 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⁻², 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⁻², 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⁻² (bottom): (left, a, d, g) modeled horizontal velocities simulated by the FOM (m a⁻¹); (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⁻², 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⁻², 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⁻² (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⁻¹); (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⁻² 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⁻¹); (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, 2003; Pattyn and others, 2008), 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.