Interactive comment on “ Accuracy of the zeroth and second order shallow ice approximation – numerical and theoretical results ” by J . Ahlkrona

Accuracy of the zeroth- and second-order shallow-ice approximation : numerical and theoretical results


Introduction
The cryosphere is an important part of the climate system, and includes, among other features, the Greenland Ice Sheet and the Antarctic Ice Sheet.With a volume of about 3 × 10 7 km 3 , these two ice sheets represent the largest component of the cryosphere, and store about 77 % of the global freshwater.Cryospheric research, and specifically, ice sheet modelling, is a vibrant discipline which receives much scientific, political and societal attention because of its relevance for predictions of the future sea level rise in a warming world (Solomon et al., 2007).
Increasingly complex numerical ice sheet models are used to simulate the response of the Earth's ice sheets to different climate forcings, and internationally coordinated efforts are devoted to ensure inter-comparability of such modelling exercises by means of benchmarking (Pattyn et al., 2008(Pattyn et al., , 2012;;Calov et al., 2010).Participating models range from zerothorder shallow-ice approximation (SIA) models, to higherorder models and full Stokes models; the common understanding being that higher-order models are more accurate than zeroth-order models, while full Stokes models are the most accurate (but also the most costly) ones.However, as increasingly accurate numerical ice sheet models and correspondingly powerful hardware become more readily available, overall less research is focused on e.g.analyses devoted to basic foundations of ice sheet modelling, such as the validity and limitations of approximation schemes and associated scaling procedures.These schemes and scaling procedures were especially important when virtually all ice sheet modelling relied on SIA schemes because computing power was more restricted than it is today.Yet, because of challenging applications such as palaeoglacial simulations or uncertainty quantifications, they have not lost their relevance, and a large portion of ice sheet simulations is still performed with SIA codes or higher-order approximations today.
In this paper we investigate the order of accuracy and validity of a higher-order extension to the SIA; the secondorder shallow-ice approximation (SOSIA).In the process we also study the accuracy of the classical SIA.The SIA was constructed in the 1970s and 1980s by Fowler and Larson (1978), Hutter (1983) and Morland (1984).The derivation is based on scaling and asymptotic series expansion in terms of the aspect ratio, , which expresses the shallowness of an ice sheet.In the SIA, only the zeroth-order terms are kept.In the end of the 1990s, the second-order shallow-ice approximation (SOSIA) was derived by Baral (1999) and Baral et al. (2001), pushing the series expansion to second order in with the objective of including dynamics not captured by the SIA.Computing a solution with SOSIA can be viewed as two steps in an iteration.First the SIA solution is determined, and then the SOSIA solution.A fully iterative algorithm is developed in Souček and Martinec (2008), based on an asymptotic expansion in .
The SOSIA is thus a higher-order model based on the same theory as the SIA, and is almost as computationally cheap as the SIA.The scaling assumptions underlying the SOSIA (and the SIA as derived in Baral, 1999;Baral et al., 2001;Greve, 1997) do, however, neglect a high-viscosity boundary layer near the ice surface.This boundary layer is thick (Ahlkrona et al., 2013), and other scaling assumptions by Johnson and McMeeking (1984) and Schoof and Hindmarsh (2010) are more appropriate than the classical SIA scalings (Ahlkrona et al., 2013).This calls for a proper analysis of the accuracy of the SOSIA equations, and also an investigation of the true order of accuracy of the SIA.The boundary layer is associated with the non-linear rheology of ice.There is no boundary layer in a Newtonian fluid at the upper surface and the scaling assumptions made to derive the SIA and SOSIA are valid for the whole fluid.
Our analysis is based on analytical solutions of the SOSIA (and SIA) as well as numerical solutions of SIA, SOSIA and the full Stokes equations when there is no sliding at the ice base.Ultimately, it is of interest to assess under which circumstances the SOSIA model can be regarded as a significant improvement on the SIA at low computational costs.The SOSIA, as described in Baral et al. (2001), has to our knowledge never been implemented before this study, most likely because the second-order expressions are long and tedious to code.The SOSIA was applied by Mangeney and Califano (1998) for Newtonian, anisotropic ice, and recently was applied by Egholm et al. (2011) -not, however, in its pure form, but in a depth-averaged, iterative scheme.
The outline of this paper is as follows: Sect. 2 is devoted to a summary of the general equations pertaining to ice dynamics, and to their zeroth-, first-, and second-order approximations.Recent works by Ahlkrona et al. (2013) and Kirchner et al. (2011) allow us to keep this presentation to a minimum.Sect. 3 describes the model problem which we focus on throughout the paper: a two-dimensional flow over a bumpy bed.In Sect. 4 an analytical SOSIA solution for this model problem is presented and discussed, extending previous work by Baral (1999) and Baral et al. (2001).The solution is compared to the SIA and SOSIA solutions of a Newtonian fluid.In Sect. 5 we implement the SOSIA numerically and compute the accuracy of both the SIA and the SOSIA by comparing their solutions with the solution obtained with Elmer/ICE (Gagliardini et al., 2013), for varying .The same computations are repeated for a Newtonian fluid to verify that the comparison technique is reliable.In both Sects.4 and 5 we compare our results with the theory in Schoof and Hindmarsh (2010) and investigate the effect of an extra parameter, σ res , which is necessary due to the neglect of the boundary layer at the surface.The paper concludes with a discussion in Sect.6.
2 Derivation of the zeroth-and second-order shallowice approximation

The exact equations
Ice sheet flow is commonly described using concepts from continuum mechanics, materials science, and thermodynamics, which allow for the formulation of the spatio-temporal evolution of ice masses as an initial boundary value problem, with free boundaries.Ice flow is momentum, mass, and internal energy conserving.As we will only study isothermal flow, the equations regarding energy are not described here.
The equations for balance of mass and momentum are where ρ is the density, v the velocity field, T D the deviatoric stress tensor and g gravitational acceleration.The deviatoric stress tensor, T D , and the Cauchy stress tensor, T are related by T = −pI +T D , where p is the pressure.The acceleration term, v -which is the material time derivative of the velocity in Eq. ( 2) -is very small and is therefore neglected in glaciological applications.The resulting equations are called the Stokes equations, or in glaciology rather the full Stokes equations.Velocity and stress are related by the constitutive equation, in glaciology called Glen's flow law.The strain rate tensor D is defined as , where * denotes transpose and A(T ) accounts for coupling the viscosity, η, to the pressure melting point corrected temperature T .For isothermal conditions A is merely a constant.The so-called creep response function, f , is defined by f (σ ) = σ n−1 , where we let n be equal to the standard value 3. Its argument, σ , the effective stress, is the square root of the second invariant of T D , and is defined by Here, t D ij (i, j = x, y, z) are the components of T D in a Cartesian coordinate system where the z axis is pointing in the opposite direction of gravity.We refer to t ) as the horizontal plane shear stress.As the creep response function depends on the effective stress, ice is a non-Newtonian fluid with viscosity, η, given by 1/η = 2A(T )f (σ ).This nonlinearity makes the ice flow simulations a computationally heavy task.The computations are simpler in a Newtonian fluid with n = 1.Then f = 1 and the relation between D and T D in Eq. ( 3) is linear for a constant T .
To complete the system, boundary conditions are imposed at the ice base and ice surface.At the impermeable ice base, the velocity satisfies no slip conditions on a rigid bedrock, giving the condition v = 0. (5) The ice surface is assumed to be stress-free, Here n is the outward-pointing normal vector of the ice surface.In the time-dependent case, a transport equation for the ice surface elevation is solved, where the velocity field enters as coefficients and the accumulation or ablation enters as a forcing.The equation for the height of the free ice surface h(x, y, t) is where v x , v y and v z are the velocity components and a s is the accumulation/ablation function.

Shallow-ice approximations
Here we briefly describe how the SIA and the SOSIA are derived from the exact, full Stokes equations by scaling and perturbation expansions.We exemplify the procedure by the momentum balance Eq.(2).We follow the scalings presented in Baral (1999), Baral et al. (2001), andGreve (1997), as these are the ones most commonly used today for deriving the SIA, and also those used to go further and arrive at the SOSIA equations.These scalings are where the aspect ratio, , has been introduced.The dimensionless quantities are denoted by tilde and are assumed to be of the order of magnitude O(1).[L] so that the aspect ratio, , is small.The scaling reflects that the vertical shear stresses are assumed to dominate over the normal deviatoric and normal shear stresses, and that this dominance is stronger the more shallow the ice sheet is.Also, the horizontal velocity is assumed to dominate over the vertical velocity, The scalings Eq. ( 8) are inserted into the equations and a perturbation expansion is performed, i.e. the dimensionless variables, q, are expanded in a power series as q = q(0) + q(1) + 2 q(2) + . . .
Collecting terms of equal order in gives rise to a hierarchy of models, called the SIA for zeroth order in , the FOSIA (first-order shallow-ice approximation) for first order in , and SOSIA for second order.The momentum balance Eq.(2) (in component form), for the SIA model is For the FOSIA, it is For the SOSIA, it is Note that there are no other stress components besides the vertical shear stress in the zeroth-order equations.This is the main reason why higher-order models are becoming increasingly popular: t D xx(0) , t D yy(0) , t D zz(0) and t D xy(0) are important in the coupling to ice shelves and in ice streams and ice stream shear margins.
The boundary conditions are also expanded.The stressfree condition that T • n = 0 corresponds to in zeroth order, and in second order.The no-slip condition will for the problem studied in this paper reduce to v x = v y = v z = 0 for both zeroth and second order.The stress components t D xx(0) , t D yy(0) , t D zz(0) and t D xy(0) do not require explicit boundary conditions but will be determined at the boundaries implicitly through the stress-strain relation.
Integrating Eq. (10c) in the vertical direction and using the boundary condition that the pressure is zero at the ice surface gives an explicit expression for the pressure.Inserting the pressure in Eqs.(10a) and (10b), integrating and using the stress-free condition again yields simple expressions for the vertical shear stresses.The shear stresses are sufficient to obtain the zeroth-order (SIA) velocities by integration of the stress-strain rate relation Eq. (3); see Sect. 4 for the outcome of this standard procedure.Once the zeroth-order solution is available, it can be used to solve the FOSIA, and subsequently the SOSIA equations by the same method.
Just like the SIA, the FOSIA only contains the vertical shear stresses, but no normal deviatoric stress or horizontal plane shear stress.However, the FOSIA does account for first-order boundary effects and forcings that cannot be represented in SIA (Baral, 1999).For the model problem that we study, the FOSIA solution is trivially zero, and the FOSIA will hence not be further discussed.It is in the SOSIA that the normal deviatoric stresses and the horizontal plane shear stresses are present for the first time, suggesting that more complex ice dynamical behaviour can be captured by the SOSIA than by the SIA.The SOSIA is computationally inexpensive compared to many other higher-order models, because the equations can be solved without a coupling of variables in a system of linear equations.

Boundary layer treatment
The SIA and the SOSIA are both based on the same scaling arguments in Eq. ( 8), where it is assumed, in addition to the ice body being shallow, that the dominating stress component is vertical shear stress and that the velocity components in the horizontal plane dominate over the vertical velocity.These assumptions are not valid in a number of well-known situations, including fast sliding at the ice-bedrock interface (as in for example ice streams), or at the ice divide where the ice flows mainly downwards.Where fast sliding occurs it is common to use other scaling arguments, such as e.g. the shelfy stream approximation in MacAyeal (1992).Different traction numbers representing different sliding speeds combined with other scalings are introduced in Schoof and Hindmarsh (2010).
Another region where the SIA scalings break down is in a boundary layer near the entire ice surface, which develops when there is a bumpy bed due to the non-linear rheology of ice (Ahlkrona et al., 2013).Johnson and McMeeking (1984) made the first attempt to derive a solution for the boundary layer by matched asymptotics, rescaling the pressure and stress components in the boundary layer.This allows for all stress components to influence the dynamics near the ice surface.By theoretical analysis they found the boundary layer thickness to be O( 13 ).Schoof and Hindmarsh (2010) extended the boundary layer theory by including the degree of slip at the bed in the expansion, and pushing it to second order.In the case of slow sliding, the rescaling of the pressure and stress components in the boundary layer in two dimensions is as in Johnson and McMeeking (1984); Ahlkrona et al. (2013) and in Schoof and Hindmarsh (2010).Using their traction number λ = 1/3 for the conditions at the bedrock: . By numerically solving the full Stokes equations, Ahlkrona et al. (2013) essentially confirmed the appropriateness of these rescalings for the problem described in Sect.3. In fluid dynamics, boundary layers are usually assumed to be thin, but as found in Ahlkrona et al. (2013), they may be thick and indistinct at the ice surface (depending on ), and this needs to be considered in model development.
Matched asymptotics, as used in Johnson and McMeeking (1984) and Schoof and Hindmarsh (2010). is unfortunately rather involved and has, to our knowledge, never been implemented into a numerical code for practical use.Baral (1999) and Baral et al. (2001) indeed dismiss it as too complicated and use the scalings in Eq. ( 8) for the entire ice sheet in the derivation of the SIA and SOSIA.This, in combination with the fact that the rheology of ice is singular and nonlinear, does introduce complications that need to be remedied.The singularity arises from the fact that the creep response function is zero for zero effective stress, i.e. f (0) = 0.This means that the viscosity is infinite where the effective stress is zero.By the shallow-ice scalings Eq. ( 8), the normal and horizontal plane shear stress t D xx , t D yy , t D zz , t D xy in Eq. ( 4) are neglected such that f (σ (0) ) = σ 2 (0) = t 2 xz(0) + t 2 yz(0) .Thus the (zeroth-order) effective stress and the creep response function are zero wherever the zeroth-order vertical shear stresses are zero, which is at the entire ice surface (Baral, 1999;Baral et al., 2001;Greve, 1997).In reality, normal deviatoric stresses (and depending on the situation, horizontal plane shear stresses) develop, implying non-zero effective stress at the ice surface except for a few points.
As SIA is the zeroth-order expansion, it can be derived using several different scaling arguments.Therefore, it is not very sensitive to the neglect of normal and horizontal plane shear stress in the boundary layer.Indeed, the computation of the SIA solution does not include the reciprocal of the creep response function (except for in the normal deviatoric stress), and the fact that the zeroth-order effective stress is zero at the entire surface does not imply singularities in the velocity field, pressure, or shear stress.The SOSIA is however bound to be more sensitive.When pushing the asymptotic expansions to second order, the creep response function of the zeroth-order effective stress does occur in the denominator.To remedy this, an extra parameter, σ res , is used in Baral (1999) and Baral et al. (2001), to regularise the problem, following earlier suggestions by Lliboutry (1969) and Colbeck and Evans (1973).This parameter, which we will call the finite viscosity parameter σ res , is added to f as: + σ 2 res (see e.g.Colbeck and Evans, 1973).Nonsingular creep functions of non-additive structure have been proposed by e.g.Lliboutry (1969).The question of the appropriateness of modifying the material law instead of rescaling the variables in the boundary layer immediately arises.As already mentioned, the singularities do not affect SIA.Hence, in practice σ res is not needed in the SIA, but the solution will not be accurate to the order predicted by the theory in Baral (1999) and Baral et al. (2001).It is unclear whether the SOSIA will be an improvement on SIA, because of the neglect of the special boundary layer dynamics.Moreover, there is no obvious way to choose the value of the finiteviscosity parameter.In Baral et al. (2001), σ res = √ 10 9 is used for the Greenland Ice Sheet, but with no explanation of this choice.In Sect.5, we will investigate how to choose σ res , how accurate the SIA is, and whether the SOSIA really is an improvement on the SIA.
J. Ahlkrona et al.: The second order shallow ice approximation 5 As SIA is the zeroth order expansion it can be derived using several different scaling arguments.Therefore it is not very sensitive to the neglect of normal and horizontal plane shear stress in the boundary layer.Indeed, the computation of the SIA solution does not include the reciprocal of the creep response function (except for in the normal deviatoric stress), and the fact that the zeroth order effective stress is zero at the entire surface does not imply singularities in the velocity field, pressure, or shear stress.The SOSIA is however bound to be more sensitive.When pushing the asymptotic expansions to second order, the creep response function of the zeroth order effective stress does occur in the denominator.To remedy this, an extra parameter, σ res , is used in Baral (1999) and Baral et al. (2001), to regularize the problem, following earlier suggestions by Lliboutry (1969) and Colbeck and Evans (1973).This parameter, which we will call the finite viscosity parameter σ res , is added to f as: res (see e.g.Colbeck and Evans, 1973).Nonsingular creep functions of non-additive structure have been proposed by e.g.Lliboutry (1969).The question of the appropriateness in modifying the material law instead of rescaling the variables in the boundary layer immediately arises.As already mentioned, the singularities do not affect SIA.Hence, in practice σ res is not needed in the SIA, but the solution will not be accurate to the order predicted by the theory in Baral (1999) and Baral et al. (2001).It is unclear whether the SOSIA will be an improvement of SIA, because of the neglect of the special boundary layer dynamics.Moreover, there is no obvious way to choose the value of the finite viscosity parameter.In Baral et al. (2001), σ res = √ 10 9 is used for the Greenland ice sheet, but with no explanation of this choice.In Sect.5, we will investigate how to choose σ res , how accurate the SIA is, and whether the SOSIA really is an improvement of the SIA.

Model problem -ice flow over a bumpy bed
Throughout this paper we will consider the model problem described in this section.It is a slight modification of the problem studied in the ISMIP-HOM benchmark experiment B (Pattyn et al., 2008).As in Pattyn et al. (2008), we investigate a diagnostic, isothermal 2-D-problem for ice flow over an inclined, bumpy bed.The ice surface is fixed, periodic boundary conditions are applied and no-slip conditions are imposed at the base, see Fig. 1.The rate factor A is 10 −16 Pa −3 a −1 , ice density is the standard value of 910 kg m −3 , and accumulation and ablation are neglected.The mean ice thickness, [H], is 1000 m and the ice surface, h, and ice base, b, are given by h The ice base is smooth, thus avoiding additional difficulties with SIA and SOSIA for a bedrock with less regularity.The amplitude of the bumps is µ[H] and the inclination angle of the surface slope is α.The typical horizontal extent of the problem equals the wavelength of the sinusoidal bumps, i.e.L = [H]/ m.The wavelength L of the bumps is varied while [H] is kept constant, which corresponds to varying .As shown in Ahlkrona et al. (2013), and as can be seen for instance by non-dimensionalising the ice surface and ice bed, the surface slope angle α should be proportional to e.g.arctan instead of α = 0.5 o which was used in the ISMIP-HOM benchmark.The bump amplitude µ should be independent of .In our numerical experiments we will use µ = 0.5 and µ = 0.1.
4 Analytical solutions for the SIA and the SOSIA For a deeper understanding, we now compute analytical solutions for the second order field variables t xz(2) , p x(2) and v x(2) , and also give the zeroth order expressions for completeness.Note that for the 2-D case v y = t D yy = t D xy = 0.The second order z velocity is excluded for brevity of presentation, since it follows from the mass balance in the same way for all orders.

General solution
The following expressions (Eqs.19-24) hold, not only for our model problem, but for all isothermal, 2-D problems with no-slip conditions at the base and without higher order boundary terms.Generalising to 3-D, including higher order boundary terms, and sliding, is straightforward.The well-known zeroth order expressions, denoted by subscript

Model problem -ice flow over a bumpy bed
Throughout this paper we will consider the model problem described in this section.It is a slight modification of the problem studied in the ISMIP-HOM benchmark experiment B (Pattyn et al., 2008).As in Pattyn et al. (2008), we investigate a diagnostic, isothermal 2-D problem for ice flow over an inclined, bumpy bed.The ice surface is fixed, periodic boundary conditions are applied and no-slip conditions are imposed at the base (see Fig. 1).The rate factor A is 10 −16 Pa −3 a −1 , ice density is the standard value of 910 kg m −3 , and accumulation and ablation are neglected.The mean ice thickness, [H ], is 1000 m and the ice surface, h, and ice base, b, are given by h The ice base is smooth, thus avoiding additional difficulties with SIA and SOSIA in modelling a bedrock with less regularity.The amplitude of the bumps is µ[H ] and the inclination angle of the surface slope is α.The typical horizontal extent of the problem equals the wavelength of the sinusoidal bumps, i.e.L = [H ]/ m.The wavelength L of the bumps is varied while [H ] is kept constant, which corresponds to varying .
As shown in Ahlkrona et al. (2013), and as can be seen, for instance, by non-dimensionalising the ice surface and ice bed, the surface slope angle α should be proportional to e.g.arctan instead of α = 0.5 o which was used in the ISMIP-HOM benchmark.The bump amplitude µ should be independent of .In our numerical experiments we will use µ = 0.5 and µ = 0.1.
For a deeper understanding, we now compute analytical solutions for the second-order field variables t xz(2) , p x(2) and v x(2) , and also give the zeroth-order expressions for completeness.Note that for the 2-D case, v y = t D yy = t D xy = 0.The second-order z velocity is excluded for the sake of brevity, since it follows from the mass balance in the same way for all orders.

General solution
The following expressions (Eqs.19-24) hold not only for our model problem, but also for all isothermal, 2-D problems with no-slip conditions at the base and without higher-order boundary terms.Generalising to three-dimensionality, including higher-order boundary terms and sliding, is straightforward.The well-known zeroth-order expressions, denoted by subscript (0), for shear stress (cf.Eq. 10) and velocity are with v x(0) = 0 at b (0) , cf.Eq. ( 5).In order to compute from Eq. ( 12) we need the zeroth-order normal deviatoric stresses (and depending on the situation, horizontal plane shear stress), which are not computed when applying only the zeroth-order approximation, as they are not needed for the zeroth-order velocities.The normal deviatoric stress in the x direction is given directly by the stress-strain relation in Baral et al. (2001) and Greve (1997): Note that t D xx = −t D zz in the 2-D case.The shear stress, pressure and the horizontal velocity are obtained from the horizontal momentum balance, vertical momentum balance and stress-strain relation, respectively, where they occur in derivatives.
It is customary to vertically integrate the SIA equations in order to avoid having to solve for the variables numerically.This is also convenient for the SOSIA.In Baral et al. (2001), the vertical integration is carried out for the shear stresses and the pressure only, and the resulting expressions are presented in Eqs. ( 22) and ( 23) (with misprints in Baral et al., 2001, corrected).Going beyond the results of Baral et al. (2001), we derive here also the second-order velocity, v x(2) , in Eq. ( 24), computed from the stress-strain relation and boundary conditions at the base: Equations ( 22)-( 24) yield explicit expressions for secondorder variables.

Solution for the model problem
For the geometry in Sect.3, we have computed the integrals and derivatives in these expressions and thus obtained analytical solutions to the SOSIA.To avoid infinite viscosity, we have followed the suggested regularisation in Baral et al. (2001), i.e. adding a constant to the creep response function, f (σ ) = σ 2 +σ 2 res .The solutions are expressed in terms of the inclination angle of the ice surface α, relative amplitude µ, finite-viscosity parameter, σ res , and wavelength L. In the Appendix A the same solutions are expressed in a more general form.By Eq. ( 19), the zeroth-order shear stress for the model problem is The zeroth-order velocity in Eq. ( 20) is given by The maximum magnitude of the velocity in Eq. ( 26) does not depend on L, but it does depend on the relative amplitude µ.Also, there is an extra term stemming from the finite-viscosity law with σ res > 0. As mentioned in Sect.2.3, a finite-viscosity law is technically not needed for the zerothorder model, since the creep response function never appears in the denominator in the calculations of the velocity field or shear stress (see Eqs. 25 and 26), which are usually the variables of interest.However, for reasons of consistency, it should also be used in the zeroth-order model when continuing to second order.To calculate second-order shear stress and velocity, the zeroth-order normal deviatoric stress in Eq. ( 21) is needed.It is given for the model problem by Having calculated t D xx(0) , the second-order shear stress is computed from Eq. ( 23) as An explicit dependence on L is introduced in Eqs. ( 27)-( 28), and hence on the aspect ratio .Also, Eqs. ( 27) and (28) depend on α, µ and σ res .In t D xx(0) , σ res appears in the denominator, preventing singularities from occurring at the ice surface when z is equal to h = −x tan(α); see Eq. ( 27).In the second-order shear stress in Eq. ( 28), σ res is both in the numerator and the denominator.The second-order shear stress is dominated by the last term, i.e. by the last five lines in Eq. ( 28).
Having calculated t D xx(0) and t D xz(2) , the second-order velocity can be derived from Eq. ( 24).The expression for the second-order solution is very long and is included in the Appendix A for the interested reader.In fact it behaves similarly to the second-order shear stress, where σ res appears both in the numerator and the denominator.Remember that to get the full second-order solution, the zeroth-and second-order contributions should be added together as v x = v x(0) + 2 v x(2) (the first-order solution is zero for our model problem).The extra term arising from the finite-viscosity law in v x(0) in Eq. ( 26) is quite large, and thus the SOSIA velocity does not allow too large a σ res .The SOSIA shear stress is not as sensitive to large σ res , since the zeroth-order shear stress in Eq. ( 25) does not contain any terms from the finite-viscosity law.Both the second-order shear stress correction and velocity correction are, however, very sensitive to too small a σ res .
Note that all terms involving σ res in Eqs. ( 27), (28), and (A6) in the Appendix A are pre-multiplied with µ; thus, the importance of σ res decays when µ decreases.This is consistent with the fact that the boundary layer near the ice surface disappears when µ = 0 (Ahlkrona et al., 2013).

Choices of σ res and impact on scalings
In order for any scaling relations (Eq.8 or 17) to hold, tan(α) should vary linearly with , e.g.tan(α) = , and µ should be independent of (Ahlkrona et al., 2013).Note, however, that the scaling relations were derived in a context where the creep response function was not modified by an additional finite-viscosity parameter.We now discuss how σ res , introduced in an a posteriori fashion, can be chosen such that compatibility with either the classical SIA scalings in Eq. ( 8) or the ones in Eq. ( 17) is achieved.

Choosing a σ res consistent with the classical SIA-scalings
If we choose the finite-viscosity parameter, σ res , to vary in the same way as the effective stress is assumed to (linearly with ), the SIA and SOSIA solutions fulfil the scaling relations Eq. ( 8) that they are derived from.Inserting tan(α) = , σ res = C σ ρg[H ] and the scaled variables in Eq. ( 8) into Eqs.( 25)-(28) yields where C σ is a constant and x and z are the nondimensionalised x and z coordinates in Eq. ( 8).The expressions in Eqs. ( 29)-( 32) only depend on geometry, material constants, and C σ .In line with Eq. ( 8), is pre-multiplied by ρg[H ] 3 .Note that the multiplying fac- Blatter (1995) and Schoof and Hindmarsh (2010).In the same manner, the analytical solution for v z(0) (which is not given here for brevity of presentation) is multiplied by

Choosing a σ res consistent with boundary layer theory
We know that the scaling relations Eq. ( 17) from Schoof and Hindmarsh (2010) (slightly rearranged) are more correct than the classical scaling relations in Eq. ( 8) (Ahlkrona et al., 2013;Schoof and Hindmarsh, 2010;Johnson and McMeeking, 1984), and that σ res is merely a parameter introduced to address this problem when using Eq. ( 8).However, we show now that instead of setting σ res = C σ ρg[H ] , with C σ constant, we can choose C σ such that the field variables fulfil the scaling relations in Eq. ( 17) rather than Eq. ( 8).
Near the ice surface, σ is dominated by t D xx (Ahlkrona et al., 2013;Schoof and Hindmarsh, 2010) 31), inserting σ res = t D xx and x + z = 0, solving for C γ and finally ignoring terms depending on , we find that C γ is approximated by Thus C γ decreases with decreasing bump amplitude only depend on the geometry and is of O(1).Due to the assumptions made when obtaining Eq. ( 33), it cannot be used directly to determine C γ , but gives an understanding of its behaviour which we will recognise in our numerical experiments.
The first correction term applied to the SIA solution t D xz(0) in Eq. ( 29) is 2 t D xz(2) in Eq. ( 32).With C σ = C γ 1/3 , at the ice surface and outside the boundary layer

A Newtonian fluid
For comparison, we derive the asymptotic expansions for a Newtonian fluid with n = 1 in Glen's flow law in Eq. ( 3) using the same techniques as for the ice model with n = 3 and the same model problem in Sect.3. The shear stress is independent of the flow law in Eq. ( 19), and the zeroth-order term is the same as in Eq. ( 25): Since f (σ (0) ) = 1 in Eq. ( 20) there is no need to introduce σ res , and the expression for the zeroth-order velocity is simplified compared to Eq. ( 26): This expression is equal to the contribution by the constant part of the creep function in Eq. ( 26).Also the lowest order term in the normal deviatoric stress in Eq. ( 21) is simplified when n = 1: This formula is similar to the formula for t D xx(0) when n = 3 in Eq. ( 27).The second-order shear stress is obtained by Eq. ( 23): The first term is the same as in Eq. ( 28) in the non-Newtonian case, while the rest of the expression is simpler and does not include any terms with σ res or arctan.Neither Eq. ( 37) nor ( 38) is singular for any x and z, as t D xx(0) and t D xz(2) are with σ res = 0 in Eq. ( 27) and (28).
Since x tan(α) + z ∼ [H ] and with tan(α) ∼ we have through Eqs. ( 35) and (36) that t D xz(0) ∼ [H ] and v x(0) ∼ [H ] 2 .Hence, for the zeroth-order velocity in the z direction, and for the second-order velocity in x direction, The scaling in Eq. ( 37) is such that t D xx(0) ∼ 2 [H ], and in Eq. ( 38), such that 2 t D xz(2) ∼ 3 [H ].These scalings are all in agreement with the assumptions made in Eq. ( 8) in the derivations of the shallow-ice approximations.

Numerical computation of the accuracy of SIA and SOSIA
In this section, we compute the SIA and SOSIA solutions for the problem described in Sect.3, and compute the solutions' accuracy by comparing them with a full Stokes solution.All our results presented in this section regard the accuracy for the x velocity, v x , and shear stress, t D xz .The normal deviatoric stresses t D xx and t D zz are not calculated to second order, since this is not necessary in order to obtain the velocity field.The vertical velocity is also excluded, as it follows directly from the mass balance, which is the same for all orders.

Method
We are interested in the order of accuracy, i.e. how the error in SIA and SOSIA varies with .For this purpose we perform repeated simulations using L = 10, 20,40,80,160,320,640,1280,2560,5120, and 10 240 km while keeping H constant; this is equivalent to varying the aspect ratio between 9.77 × 10 −5 and 0.1.We do this in order to investigate the accuracy of shallow-ice approximations in the limit → 0.
The SOSIA Eqs. ( 19)-( 24) are implemented in MATLAB and our implementation follows the standard in SImulation COde for POLythermal Ice Sheets (SICOPOLIS) (Greve, 1995).Finite differences are used on a staggered grid in order to avoid having oscillatory solutions with the same wavelength as small multiples of the grid size.The velocities, horizontal volume fluxes, vertical shear stresses and the horizontal derivatives of the bedrock topography and the icesurface topography are defined in between the grid points.All other quantities are defined at grid points.When a quantity is needed at a point where it is not defined, linear interpolation is used.To ensure that the grid points coincide with physical boundaries, a σ transformation is used (Greve and Blatter, 2009;Greve, 1995).Central differences are applied when possible; otherwise one-sided differences are used.Integrals are computed by the trapezoidal method if the integrand and the integral are defined at the same points.Otherwise, the midpoint rule is used.
The full Stokes solution that we use for comparison is obtained using the finite-element code Elmer/Ice (Gagliardini et al., 2013).The same mesh with the nodes on vertical lines was used for SIA, SOSIA and in Elmer/Ice.We use a mesh fine enough to keep the relative numerical error below 10 −4 for both the velocity and shear stress.This error can be seen in some of the figures below, but as the mesh is refined it decreases even more.Even though the singular behaviour of the viscosity does not introduce singularities in the field variables in the full Stokes setting, an extra parameter, the critical shear rate γo , is introduced in Elmer/Ice in order to treat numerical instabilities at low stress (Råback et al., 2013), occurring e.g.near the ice surface in the shear stress.The critical shear rate is a lower bound for the shear rate γ , which is related to the effective stress by γ = 2tr(D 2 ) = 2d = 2Aσ 3 .We have used γo = 10 −10 throughout the simulations.For large aspect ratios there is no γo which suppresses the numerical instabilities in the shear stress without altering the velocity.
The accuracy is measured in terms of the relative error defined by where q is v x or t D xz and || • || 2 denotes the L 2 norm defined by where V is the area of .The integral in Eq. ( 42) is computed on a discrete grid using the trapezoidal rule.
Since the assumptions behind the SOSIA are valid only below the boundary layer, we measure the accuracy in horizontal layers, viz. at 0.01[H ], 0.5[H ], 0.75[H ], 0.9[H ] and 0.95[H ] mean height above the ice base, to see if the errors due to the boundary layer assumptions in SOSIA spread down into the basal ice.

Results
We know from our analytical solutions in Sect.4.3 that the SOSIA solution is very sensitive to the parameter σ res .Therefore, we experiment with different values of and ways of setting σ res .
Our analytical solutions suggest that if the scaling relations used to derive the SOSIA were correct, σ res should vary with as C σ ρg[H ] , C σ being a constant.This relationship is used in Fig. 2a and b, showing the accuracy of the SIA and SOSIA solution for a bumpy bed with relative bump amplitude µ = 0.5.The SIA solution for the horizontal velocity component is v x(0) in Eq. ( 26) with C σ = 0, since a finiteviscosity law is unnecessary in this case.The SOSIA solution is computed for C σ equal to 0.35.

Accuracy of the SIA
We start by analysing the error in SIA.According to classical SIA theory in Baral (1999), Baral et al. (2001), andGreve (1997), the SIA relative error should be of O( 2).However, we know from Ahlkrona et al. (2013) that the scalings used in Schoof and Hindmarsh (2010) are more correct than the classical SIA scalings.The asymptotic expansions in Schoof and Hindmarsh (2010) yield the same zeroth-order solution as the SIA, but the correction terms are different.The relative error in the velocity is estimated by the first neglected term in the expansions of v x in Eqs.(3.73) and (3.108) in Schoof and Hindmarsh (2010) and is of order 4/3 in the boundary layer and 5/3 outside the layer.The estimated slope, log(error)/ log( ), of the relative SIA error using all values in Fig. 2a (thick red line with nodes) is 1.43, which is between the theoretical rates in Schoof and Hindmarsh (2010).Similarly, the relative error in the SIA solution of t D xz(0) should, according to Eqs. (3.108) and (3.73) in Schoof and Hindmarsh (2010) and Eqs. ( 29) and ( 34), be 4/3 in the boundary layer and 5/3 outside it.This is in good agreement with the slope of the relative error in Fig. 2b, which is 1.38.

Accuracy of the SOSIA: σ res = C σ ρg[H ]
We now move to analysing the SOSIA error in Fig. 2a and b (thick blue line with nodes).If the classical SIA theory were correct, SOSIA would in principle be much more accurate than SIA for sufficiently small .Clearly this is not the case in Fig. 2a and b.In addition to computing the SOSIA solution with C σ = 0.35, we also tried C σ = 0.11 and C σ = 1.12, with no improvement in accuracy.
The SOSIA is thus not a correction to the SIA with σ res =C σ ρg[H ] and C σ constant.This is because the scaling relations in Eq. ( 8) are not correct for the thick boundary layer near the ice surface (which will dominate in the global error).Since the scaling relations do hold below this layer (Ahlkrona et al., 2013;Schoof and Hindmarsh, 2010;Johnson and McMeeking, 1984), one would like to know if the SOSIA solution is more accurate further down in the ice.To investigate this, we compute the accuracy of SIA and SOSIA at horizontal layers at 0.01[H ], 0.5[H ], 0.75[H ], 0.9[H ] and 0.95[H ] mean height above the ice base; the result is included in Fig. 2a and b (the thin red and blue lines).The accuracy of the SOSIA velocity is slightly higher-deeper in the ice for small , but even there the SOSIA is less accurate than SIA.The shear stress is more accurate deeper in the ice for both SIA and SOSIA.

Accuracy of the SOSIA: σ res as a lower threshold
The SOSIA solution is very sensitive to σ res , and the different choices of C σ yield very different results.To limit the sensitivity of σ res , we can choose to use it only where the effective stress is too small.We therefore experiment with σ res as a lower bound on the effective stress, viz.: Figure 2c and d show the error using this approach when C σ = 0.11, 0.35 and 1.12.The SOSIA velocity error decreases considerably, and there are now aspect ratios for which SOSIA is more accurate than SIA when C σ = 0.35.The stress t D xz is not largely affected.We have also computed the accuracy of SOSIA in layers throughout the ice, in the same way as in Fig. 2a and b.For the velocity there is a significant change for small aspect ratios.At 0.01[H ] mean height over the ice surface the SOSIA velocity solution is more accurate than the SIA solution for all aspect ratios smaller than 10 −2 .For the deviatoric shear stress, the layer-wise error is very similar to what is shown in Fig. 2b.
The error of the velocity for small aspect ratios is mainly due to the term on the third line of Eq. ( 26).It is an extra term in the zeroth-order velocity in SOSIA introduced by the use of a finite-viscosity law.This error is reduced by using σ res as a lower threshold and it is the most influential near the ice surface where the scaling relations do not hold and where the zeroth-order shear stress and therefore the zeroth-order effective stress is zero.This explains why the error decreases further down in the ice.Since there is no such zeroth-order term in the vertical shear stress, it is not affected by this type of error.In fact, for small aspect ratios, even an extremely large σ res does not cause a large error in the second-order shear stress.The stress is thus less sensitive to the handling of σ res , which can be seen in Figs.2a-d.The zeroth-order term involving σ res is dominant if C σ is too large and results in the velocity being too high overall (see Fig. 3a) where C σ = 1.12 and = 1/160.
For larger aspect ratios, the error is of a different character.It arises from an excessive second-order correction, resulting in a dip in the velocity and shear stress at x = 3L/4; see Fig. 3b for the velocity.This type of error is dominant when is large and σ res is too small (see Fig. 3a for C σ = 0.11 and Fig. 3b), and it is difficult to avoid for the largest .The character of the dip can be understood by studying the dominant     The relative error in v x decays as 2.81 in Fig. 2c for ≥ 1/320.This is the order of the next term in the expansion.The reduction of the error ends when the term proportional to the constant σ res with C σ = 0.35 in Eq. ( 26) becomes important.The same observation is valid also for the SOSIA error in the vertical shear stress in Fig. 2d.The decay of the relative error for ≥ 1/320 is here 2.51 before the reduction is damped by the constant σ res .With the classical theory in Baral (1999) and Baral et al. (2001), the relative SOSIA error should be O( 3 ).In practice, the larger aspect ratios are of interest.Numerical errors are commonly of the order 10 −2 , which is higher than the model error in the SIA solution for aspect ratios smaller than 3×10 −3 .For ice sheet flow, aspect ratios smaller than 10 −3 are seldom applicable.

Accuracy of the SOSIA:
and further adjustments Figures 2c and 2d indicate that the optimal choice of C σ might not be independent of the aspect ratio.Indeed, we found in Sect. 4 that if we choose σ res as C γ ρg[H] 4/3 (where C σ = C γ 1/3 ) the SOSIA correction terms are remedied so that they are consistent with the scalings in John-son and McMeeking (1984), Schoof and Hindmarsh (2010), and Ahlkrona et al. (2013), see Eq. ( 34).We can even go further in limiting the influence of σ res by only using it as a lower threshold in the computations where the creep response function is in the denominator (that is in the computation of the normal deviatoric, and normal shear stresses, see Eqs. ( 31) and ( 32)).The combined results of these two measures are shown in Figs.4a and 4b.The SOSIA velocity and shear stress are now more accurate than SIA for all aspect ratios smaller than 10 −2 with C γ chosen to be 3.The slope of the errors in Fig. 4 is almost equal for both SIA and SOSIA indicating that the order in in the remaining error is the same for both approximations.Using a different parameter e.g.C γ = 4 results in SOSIA being more accurate than SIA for even larger aspect ratios, but on the other hand the improvement is no longer as significant.The parameter C γ depends on the geometry (see Eq. ( 33)) and is thus problem dependent and difficult to determine beforehand.Interesting to note is that there is no difference in accuracy of the velocity in the different layers through the ice.The accuracy of the deviatoric shear stress through the layers is distributed similarly to Fig. 2b.term in the second-order shear stress, i.e. the last five lines in Eq. ( 28), or the last four lines in Eq. (32).As this term depends on (1 − µ sin (2πx/L)), it is largest at x = 3L/4, and as it is pre-multiplied by µ/C σ , it will in general decrease with increasing C σ and decrease with decreasing bump amplitude.The error is attributed to the improper handling of the singularity in Glen's flow law.At the point x = 3L/4 the effective stress is zero in both the SOSIA as well as in the full Stokes setting.
The relative error in v x decays at 2.81 in Fig. 2c for ≥ 1/320.This is the order of the next term in the expansion.The reduction of the error ends when the term proportional to the constant σ res with C σ = 0.35 in Eq. ( 26) becomes important.The same observation is also valid for the SOSIA error in the vertical shear stress in Fig. 2d.The decay of the relative error for ≥ 1/320 is here 2.51 before the reduction is damped by the constant σ res .With the classical theory in Baral (1999) and Baral et al. (2001), the relative SOSIA error should be O( 3 ).In practice, the larger aspect ratios are of interest.Numerical errors are commonly of the order 10 −2 , which is higher than the model error in the SIA solution for aspect ratios smaller than 3 × 10 −3 .For ice sheet flow, aspect ratios smaller than 10 −3 are seldom applicable.

and further adjustments
Figure 2c and d indicate that the optimal choice of C σ might not be independent of the aspect ratio.Indeed, we found in Sect. 4 that if we choose σ res as C γ ρg[H ] 4/3 (where C σ = C γ 1/3 ) the SOSIA correction terms are remedied so that they are consistent with the scalings in Johnson and McMeeking (1984), Schoof and Hindmarsh (2010), and Ahlkrona et al. (2013); see Eq. ( 34).We can even go further in limiting the influence of σ res by only using it as a lower threshold in the computations where the creep response function is in the denominator (that is in the computation of the normal deviatoric, and normal shear stresses, see Eqs. 31 and 32).The combined results of these two measures are shown in Fig. 4a and b.The SOSIA velocity and shear stress are now more accurate than in SIA for all aspect ratios smaller than 10 −2 with C γ chosen to be 3.The slope of the errors in Fig. 4 is almost equal in both SIA and SOSIA, indicating that the order in in the remaining error is the same for both approximations.Using a different parameter, e.g.C γ = 4, results in SOSIA being more accurate than SIA for even larger aspect ratios, but on the other hand the improvement is no longer as significant.The parameter C γ depends on the geometry (see Eq. 33) and is thus problem-dependent and difficult to determine beforehand.Interesting to note is that there is no

Low bump amplitude
Here we investigate how the accuracy changes as the bump amplitude is decreased.There is a common perception that shallow ice approximations are more accurate for lower bump amplitudes.Also, we found in Sect. 4 that the terms involving σ res in the stresses and (Eqs.31 and 32) and velocity Eq. (A8) were pre-multiplied with the relative bump amplitude µ, suggesting that the influence of σ res decreases with decreasing bump amplitude.We have applied the SIA and SOSIA for µ = 0.1 (a bump amplitude of 10 % of mean ice thickness), and the result is illustrated in Fig. 4c and d.We set σ res in the same way as in Fig. 4a and b but with C γ = 1.5 instead of C γ = 3.Indeed, Eq. ( 33) shows that C γ decreases with µ.There is a small improvement of the accuracy of SIA and SOSIA compared to the case when µ = 0.5.For large aspect ratios of almost 0.1 there is a notable improvement in both SIA and SOSIA, but at these large aspect ratios the errors are still very large and neither SIA nor SOSIA is a good model.An explanation of the lack of significant improvement for lower amplitudes is that even if the classical scalings in Eq. ( 8) do hold for a flat bed, the thick boundary layer where variables rescale develops very rapidly as a small bump is introduced (Ahlkrona et al., 2013).Note that for very small aspect ratios numerical errors are present, so that the error for SIA and SOSIA is not smaller than 10 −4 .This causes a bend in the curves in the lower left parts of Figs.4c and 4d, but also in Figs.2a-2b and 4a-4b.

Newtonian Fluid
As seen in Sect.4.4, the classical SIA-scalings hold for a Newtonian rheology (n = 1).For a Newtonian fluid we therefore expect the error in the SIA and SOSIA to behave difference in the accuracy of the velocity for the different layers through the ice.The accuracy of the deviatoric shear stress through the layers is distributed similarly to how it is shown in Fig. 2b.

Low bump amplitude
Here we investigate how the accuracy changes as the bump amplitude is decreased.There is a common perception that shallow-ice approximations are more accurate for lower bump amplitudes.Also, we found in Sect. 4 that the terms involving σ res in the stresses and Eqs. ( 31 and 32) and the velocity in Eq. (A8) were pre-multiplied with the relative bump amplitude µ, suggesting that the influence of σ res decreases with decreasing bump amplitude.We have applied the SIA and SOSIA for µ = 0.1 (a bump amplitude of 10 % of mean ice thickness), and the result is illustrated in Fig. 4c and d.We set σ res in the same way as in Fig. 4a and b but with C γ = 1.5 instead of C γ = 3.Indeed, Eq. ( 33) shows that C γ decreases with µ.There is a small improvement of the accuracy of SIA and SOSIA compared to the case when µ = 0.5.For large aspect ratios of almost 0.1, there is a notable improvement in both SIA and SOSIA, but at these large aspect ratios the errors are still very large and neither SIA nor SOSIA is a good model.An explanation of the lack of significant improvement for lower amplitudes is that even if the classical scalings in Eq. ( 8) do hold for a flat bed, the thick boundary layer where variables rescale develops very rapidly as a small bump is introduced (Ahlkrona et al., 2013).Note that for very small aspect ratios, numerical errors are present, so that the error for SIA and SOSIA is not smaller than 10 −4 .This causes a bend in the curves in the lower left parts of Fig. 4c and d, but also in Figs.2a-b and 4a-b.

Newtonian fluid
As seen in Sect.4.4, the classical SIA scalings hold for a Newtonian rheology (n = 1).For a Newtonian fluid we therefore expect the error in the SIA and SOSIA to behave as predicted in Baral et al., i.e. we expect the SIA error to be O( 2) and the SOSIA error to be O( 3 ).Verifying this confirms that the methodology of all our numerical experiments     as predicted in Baral et al., i.e. we expect the SIA error to be O( 2) and the SOSIA error to be O( 3 ).Verifying this confirms that the methodology of all our numerical experiments in this paper is accurate.We measure the error for linear rheology in the same way as for the non-Newtonian case.The results are shown in Fig. 5.
The relative error decreases much faster with in Fig. 5 than in the non-Newtonian case.The error does not decrease below about 10 −4 , but as in the non-Newtonian case the reason is the numerical errors in the Elmer/Ice and the SIA and SOSIA solutions which decrease with mesh size.These error are of the same order as in Sect.5.2.A polynomial fit reveals that the error in the SIA velocity is O( 1.91 ) and in the      as predicted in Baral et al., i.e. we expect the SIA error to be O( 2) and the SOSIA error to be O( 3 ).Verifying this confirms that the methodology of all our numerical experiments in this paper is accurate.We measure the error for linear rheology in the same way as for the non-Newtonian case.The results are shown in Fig. 5.
The relative error decreases much faster with in Fig. 5 than in the non-Newtonian case.The error does not decrease below about 10 −4 , but as in the non-Newtonian case the reason is the numerical errors in the Elmer/Ice and the SIA and SOSIA solutions which decrease with mesh size.These error are of the same order as in Sect.5.2.A polynomial fit reveals that the error in the SIA velocity is O( 1.91 ) and in the in this paper is accurate.We measure the error for linear rheology in the same way as for the non-Newtonian case.The results are shown in Fig. 5.The relative error decreases much faster with in Fig. 5 than in the non-Newtonian case.The error does not decrease below about 10 −4 , but, as in the non-Newtonian case, the reason is the numerical errors in the Elmer/Ice and the SIA and SOSIA solutions which decrease with mesh size.These errors are of the same order as in Sect.5.2.A polynomial fit reveals that the error in the SIA velocity is O( 1.91 ) and in the SIA shear stress it is O( 1.93 ).The error in the SOSIA velocity is O( 3.14 ) and in the SOSIA shear stress it is O( 3.23 ).
The polynomial fit was computed with 1/640 ≤ ≤ 1/10 for SIA and 1/160 ≤ ≤ 1/10 for SOSIA in order to avoid the influence of the numerical errors at small .Our numerical results are in good agreement with theory, supporting our suspicion that the reduction in accuracy for non-Newtonian ice is due to the boundary layer near the surface.We also note that for a Newtonian fluid, σ res is no longer needed (it is set to zero), and the SOSIA error is significantly smaller than the SIA error, even for quite large .
We have solved the SIA and the SOSIA equations both analytically and numerically, and determined and analysed the accuracy of the velocity, v x , and the shear stress, t D xz , compared to the full Stokes equations.This was done for a model problem representing ice flow on a bumpy, sloping bed with no-slip conditions at the base, slightly modified from the ISMIP-HOM B benchmark set-up (Pattyn et al., 2008).The purpose was to determine whether the SOSIA is an improvement on the SIA which could be used in practice.We also wanted to show that a thick, high-viscosity boundary layer at the ice surface cannot be overlooked in higher-order asymptotic expansions, and that it results in the SIA not being second-order accurate.
We know from Ahlkrona et al. (2013), Johnson and McMeeking (1984), and Schoof and Hindmarsh (2010) that the scaling arguments behind the SIA and SOSIA as stated in Baral et al. (2001) are not valid in the boundary layer near the surface.Consequently, the numerically computed accuracy of the SIA is O( 1.43 ) for the velocity and O( 1.38 ) for the shear stress instead of O( 2), as expected from the classical SIA theory.Our results rather agree with the analysis in Schoof and Hindmarsh (2010), which predicts the accuracy of the SIA velocity and shear stress to be O( 4/3 ) in the boundary layer and O( 5/3 ) outside it.Our accuracy was measured over the whole domain including the boundary layer.
The neglect of the special boundary-layer dynamics in the classical SIA scalings results in singularities in the SOSIA field variables.To remedy this, the finite-viscosity parameter, σ res , is added to the effective stress in the creep response function when computing the SOSIA.Both our analytical solutions and numerical experiments show that the SOSIA solution is very sensitive to this parameter.It can neither be chosen too large nor to small, since it appears in both the numerator and the denominator of the expressions.Too large a σ res yields too high a velocity, especially for small aspect ratios.This is due to an extra zeroth-order term in the zerothorder velocity.This error is more important near the surface than further down in the ice.If σ res is simply added to the creep response function as f (σ ) = σ 2 + σ 2 res , the error is so large that the SOSIA velocity is less accurate than that of the SIA for all aspect ratios.If σ res is instead used as a lower threshold for the effective stress, as in Eq. ( 43), the error is decreased so that there is a range of (quite small) aspect ratios for which SOSIA is an improvement on SIA.Further reducing the effect of σ res -so that it is only used in the calculations where singularities would arise otherwise -improves the accuracy even more.The shear stress is not as sensitive to too large a σ res as the velocity is.Too small a σ res yields errors in both the velocity and shear stress for large aspect ratios, due to excessive second-order corrections caused by the singularity in Glen's flow law.
Depending on how we choose the value of σ res , the analytical expressions for the field variables conform with either the classical shallow-ice scalings that were used in the derivation, or with the scaling relations in Schoof and Hindmarsh (2010) which we know are accurate (Ahlkrona et al., 2013).For agreement with the classical SIA scalings we conclude that σ res should be C σ ρg [H ] where C σ is a constant.In order for the solution to scale with the aspect ratio as in Schoof and Hindmarsh (2010) (both inside and outside the boundary layer), σ res should be C γ ρg[H ] 4/3 , where C γ is a constant.In our numerical experiments we find that σ res = C γ ρg[H ] 4/3 yields more accurate results than σ res = C σ ρg[H ] , in the sense that the former choice, combined with using σ res as a lower threshold where singularities occur, makes SOSIA more accurate than SIA even for aspect ratios almost as large as 0.1.
Our analytical solutions suggest that the sensitivity to σ res decreases when the amplitude of the bumps at the ice base decreases.In our numerical experiments, lowering the bump amplitude from half of the mean ice thickness to 10 % of the mean thickness results in a slight improvement of the accuracy of both SIA and SOSIA for small aspect ratios.For large aspect ratios of around 0.1 the improvement is significant, but for these aspect ratios the relative error is too large, about 0.1, for the models to be applicable.
Contrary to the case in which n = 3 in Glen's flow law, the scaling assumptions in the SIA and SOSIA are valid everywhere in a Newtonian fluid with n = 1.The analytical solutions all scale in as expected.The differences between the SIA and SOSIA solutions and the computed solutions with Elmer/Ice also show the expected order of decay with diminishing for a Newtonian fluid.
Since the SIA is accurate for very small aspect ratios, and as the aspect ratio of a real ice sheet would not be less than 10 −3 , the potential applicability of SOSIA is for aspect ratios larger than that.SOSIA is an improvement on SIA for aspect ratios as large as 5 × 10 −2 depending on how the finiteviscosity parameter σ res is set and on the amplitude of the bedrock bumps.A relative error at 1 % is quite good and SOSIA achieves that for larger than SIA does.A crucial issue to overcome is how to determine σ res or C γ and C σ .These parameters are problem-dependent and the accuracy of the results is sensitive to their values.It is unclear how to choose them in practice for more complicated geometry and forcings than those treated in this paper.
As SOSIA is almost as computationally cheap as the SIA, it could be a convenient tool and replace SIA where it is applicable in e.g.palaeoglacial simulations or to provide an initial guess in an iterative solution for a more advanced ice model.However, considering the increased complexity of the model, the sensitivity to the σ res parameter, and the often marginal improvement in accuracy it may not be worth the effort to implement the SOSIA model.

Appendix A Analytical solutions to the SIA and SOSIA
This Appendix contains analytical solutions for v x(0) , v z(0) , t D xx(0) , t D xz(2) , p (2) and v x(2) for an isothermal, diagnostic, 2-D problem with no-slip conditions at the impermeable ice base, an inclined plane as ice surface, and a zeroth-order ice base.A finite-viscosity law of the form f (σ ) = σ 2 + σ 2 res is used, as discussed in Sect.2.3.For the model problem in Sect.3, we also express the solutions in terms of L, α, µ, and σ res .

A1 Zeroth-order x velocity, v x(0)
Computing the integral in the second equation in Eq. ( 19) gives Inserting the expressions for h, b (see Eq. 18) and H (= h−b) yields Eq. ( 26).
A5 Second-order pressure, p (2) The calculation of the second-order pressure from Eq. ( 22) is straightforward, knowing that t D zz = −t

Fig. 1 .
Fig. 1.Model set-up showing the basal topography and the ice surface.The ice flows down-slope in the positive x-direction.

Fig. 1 .
Fig. 1.Model set-up showing the basal topography and the ice surface.The ice flows downslope in the positive x direction.
and the creep response function satisfies f (σ ) ≈ (t D xx ) 2 in 2-D.In the SOSIA model, t D xx is neglected in the creep response, and instead f (σ ) = σ 2 res at the ice surface.This is a consequence of the asymptotic expansion of the solution and of equating terms of equal order.Hence, C σ in σ res = C σ ρg[H ] should be such that σ res = t D xx in the upper boundary layer.Let C σ = C γ γ and insert it into Eq.(31).Since x + z = 0 at the surface, we have t D xx ∼ 2 / 2γ ∼ σ res ∼ 1+γ and consequently that γ = 1/3.Then t D xx ∼ 4/3 in the boundary layer as derived in Eq. (3.108) in Schoof and Hindmarsh (2010).Outside the boundary layer in Eq. (31), when x + z = O(1), t D xx ∼ 2 as expected from the SOSIA equations.By replacing C σ by C γ γ in the expression for t D xx in Eq. ( in agreement with the second terms in the expansions in Eqs.(3.108) and (3.73) inSchoof and Hindmarsh (2010).The velocity component v x(0) in Eq. (4.11) is of O( 3 ) for everywhere in the ice.The first correction term v x(2) in the Appendix A depends in the same way as t D xz(2) on σ −1 res outside the boundary layer and σ −2 res close to the ice surface.In our numerical experiments we will apply both σ res = C σ ρg[H ] and σ res = C γ ρg[H ] 4/3 , where C σ and C γ are constants.
. t D xz .σres as a lower threshold.

Fig. 2 .
Fig. 2. Relative error of horizontal velocity vx and vertical shear stress t D xz with σres = Cσρg[H] for both SIA (red) and SOSIA (blue).Thick lines show the error measured over the whole domain and the thin lines in the upper two panels show the errors measured over horizontal layers.

Fig. 2 .
Fig. 2. Relative error of horizontal velocity v x and vertical shear stress t D with σ res = C σ ρg[H ] for both SIA (red) and SOSIA (blue).Thick lines show the error measured over the whole domain and the thin lines in the upper two panels show the errors measured over horizontal layers.

Fig. 3 .
Fig. 3. Surface x velocity for the full Stokes solution, the SIA and the SOSIA.For SOSIA, σres = Cσρg[H] is used as a lower threshold for the effective stress.

Fig. 3 .
Fig. 3. Surface x velocity for the full Stokes solution, the SIA and the SOSIA.For SOSIA, σ res = C σ ρg[H ] is used as a lower threshold for the effective stress.

Fig. 4 .
Fig. 4. Relative error of horizontal velocity vx and vertical shear stress t D xz for both SIA (red) and SOSIA (blue) for µ = 0.5 and µ = 0.1.σres = Cγρg[H] 4/3 is only used when needed in the denominator of the expressions.

Fig. 5 .
Fig. 5. Relative error of zeroth order and second order horizontal velocity vx and vertical shear stress t D xz for Newtonian fluid (n = 1).

Fig. 4 .
Fig. 4. Relative error of horizontal velocity v x and vertical shear stress t Dfor both SIA (red) and SOSIA (blue) when µ = 0.5 and µ = 0.1.σ res = C γ ρg[H ] 4/3 is only used when needed in the denominator of the expressions.

Fig. 4 .
Fig. 4. Relative error of horizontal velocity vx and vertical shear stress t D xz for both SIA (red) and SOSIA (blue) for µ = 0.5 and µ = 0.1.σres = Cγρg[H] 4/3 is only used when needed in the denominator of the expressions.

Fig. 5 .
Fig. 5. Relative error of zeroth order and second order horizontal velocity vx and vertical shear stress t D xz for Newtonian fluid (n = 1).

Fig. 5 .
Fig. 5. Relative error of zeroth-and second-order horizontal v x and vertical shear stress t D xz for Newtonian fluid (n = 1).
D xx , t D yy and t D zz as normal deviatoric stresses, to t D xz (= t xz ) and t D yz (= t yz ) as vertical shear stresses and to t D xy (= t xy