Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM)

. We present the results of the ﬁrst ice sheet model intercomparison project for higher-order and full-Stokes ice sheet models. These models are compared and veriﬁed in a series of six experiments of which one has an analytical solution obtained from a perturbation analysis. The experiments are applied to both 2-D and 3-D geometries; ﬁve experiments are steady-state diagnostic, and one has a time-dependent prognostic solution. All participating models give results that are in close agreement. A clear distinction can be made between higher-order models and those that solve the full system of equations. The full-Stokes models show a much smaller spread, hence are in better agreement with one another and with the analytical solution.


Introduction
According to the recent IPCC report (IPCC, 2007), dynamical ice-flow processes not included in current models but suggested by recent observations could increase the vulnerability of the ice sheets to warming, increasing future sea level rise. Understanding of these processes is limited, and there is no consensus on their magnitude. It was also stressed that a net ice mass loss could occur if dynamical ice discharge dominates the ice sheet mass balance (IPCC, 2007). Although the viscous flow of ice is rather well understood on a theoretical level, confidence in models is low because processes at the ice base and the seaward margin are poorly understood and have not been represented in models. Several stress components come into play in regions of high variability in basal topography and/or basal slipperiness.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Despite the lack of comprehensive predictive ice sheet modeling, the ice sheet modeling community has evolved considerably over the last decade. Increasing computational power has led to the development of more complex ice sheet models, with varying degrees of approximations to the Stokes equations. However, progress has been hampered by the lack of a universal verification framework and in particular by a lack of full-Stokes analytical solutions. While methods exist for constructing exact solutions to the Newtonian full-Stokes equations in the flowline case (e.g. Ladyzhenskaya, 1969), obtaining analytical solutions for the 3-D case is less straightforward. Nevertheless, solutions based on perturbation analysis exist (e.g. Gudmundsson, 2003), and were used for this study. Benchmark validation exercises were carried out on large-scale ice-sheet and ice-shelf models in the 1990s MacAyeal et al., 1996;Payne et al., 2000), but these tests were largely restricted to zeroorder (SIA) models and solutions. Here we present an intercomparison exercise which involves 28 higher-order flow models of varying complexity. The experiments described in this paper are designed to evaluate the conditions under which different higher-order solutions are viable and to determine whether numerical issues affect the result.
During the first and second EISMINT 1 model intercomparison exercises, a number of benchmarks were proposed specifically for ice sheet models (Huybrechts et al., , 1998Payne et al., 2000) and ice shelf models (MacAyeal et al., 1996). These ice sheet models were based on the zeroth-order shallow-ice approximation (SIA; Hutter, 1983), incorporating only vertical shear stresses in the force balance. The ISMIP-HOM exercise focuses on so-called higherorder models, i.e. models that incorporate further mechanical effects, principally longitudinal stress gradients, as well as those that solve the full system of equations of the Stokes problem.
The six experiments which comprise this benchmark exercise are designed to be universally accessible to many different types of models, i.e. flowline (2-D), vertically integrated planform (2.5-D) and full 3-D models. Furthermore, the experiments are defined as well-posed continuum problems so that their application is not limited to any specific numerical methodology. The equations have been solved using well established finite-difference and finite-element methods, in addition to more esoteric spectral techniques. The latter hold particular promise for providing high-quality results in the absence of analytical solutions, but since only one participant provided spectral model results, they were not isolated for comparison.
With the exception of Exp. F, all experiments are diagnostic; i.e. time evolution is not considered. This means that for a given ice geometry, a Glen-type flow law, and appropriate boundary conditions, the stress and velocity fields can be calculated. Exp. F considers the time-dependent response (the experiment is run until the free surface and velocity field reach a steady state) for a constant viscosity (linear flow law). Constant viscosity is assumed because in this case there exist analytic solutions derived from a first-order perturbation analysis of flow down an inclined plane (Gudmundsson, 2003). In all experiments, thermomechanical effects are neglected and ice is considered to be isothermal and isotropic.

Model physics, parameters and constants
Higher-order models are ice-sheet or glacier models that incorporate effects not present in the shallow-ice approximation. In most cases this implies the inclusion of longitudinal stress gradients apart from the two horizontal-plane shear components (Hindmarsh, 2004). Longitudinal stresses have recently been termed "membrane stresses" when considered in three dimensions (Hindmarsh, 2006). The suite of models is based on conservation laws of mass and momentum, i.e.
where ρ is the ice density, g is gravitational acceleration, v is the velocity vector, and T is the stress tensor. Values for parameters and constants are given in Table 1. Generally, acceleration terms in Eq.
(2) are neglected. Ice incompressibility is more easily described if the stress tensor is split into a deviatoric part and an isotropic pressure P , where P =− 1 3 tr(T). The constitutive equation for ice then links deviatoric stresses to strain rates: where T andė are the deviatoric stress and strain-rate tensor, respectively, and η is the effective viscosity. Both linear and nonlinear ice rheologies are considered. In the latter case (Glen's flow law), η is strain-rate dependent and defined by whereε e is the effective strain rate. For the case of linear rheology, Eq. (5) reduces to η=(2A) −1 , where A is spatially uniform (i.e., the ice is isothermal). Neglecting acceleration terms, the momentum balance is written as: The Cryosphere, 2, 95-108, 2008 www.the-cryosphere.net/2/95/2008/ Solving Eqs. (7)-(9) leads to the full-Stokes solution. Simplifications of these equations lead to the various higherorder approximations discussed below.

Boundary conditions
In Exps. A, B, E1 and F1 the ice is frozen to the bed (v b =0).
For the other experiments, basal sliding is introduced through a friction law, characterized by a friction coefficient β 2 . This friction law has the form of where n b is the unit normal vector pointing into the bedrock, t is the unit tangent vector, and β 2 (Pa a m −1 ) is a positive scalar (MacAyeal, 1993). Basal shear stress τ b is not equal to the driving stress but is part of the solution. The stress is negligible at the upper ice surface, implying that n s ·(Tn s )=P atm ≈0. Kinematic boundary conditions apply at the upper and lower surfaces of the ice mass, i.e.
for i=(s, b). Since the vertical velocity field must obey the incompressibility condition Eq. (1), and the surface accumulation/ablation is zero, the vertical velocity at the surface contains the local imbalance as well and becomes a model output. not predefined, and any type of discretization scheme can be used. The number of grid points in the horizontal and vertical directions can be chosen freely. The basic parameter for the experiments is the length scale L of the domain. Exps. A-D are carried out for L=160, 80, 40, 20, 10 and 5 km, respectively, which results in aspect ratios =H /L varying from 0.006 to 0.2. Periodic boundary conditions are applied in the horizontal, so that the domain is surrounded by an infinite number of copies of itself.

Experiment description
3.1 Exp. A: Ice flow over a bumpy bed Exp. A considers a parallel-sided slab of ice with a mean ice thickness H =1000 m lying on a sloping bed with a mean slope α=0.5 • . This slope is maximum in x and zero in y. The basal topography is then defined as a series of sinusoidal oscillations with an amplitude of 500 m. The surface elevation is defined as The basal topography is then given by z b (x, y) = z s (x, y) − 1000 + 500 sin (ω x) · sin (ω y) , where x ∈ [0, L] and L=160, 80, 40, 20, 10 and 5 km, respectively. The basal bumps have a frequency of ω=2π/L. The bed topography is shown in Fig. 1  The only difference with Exp. A is that the basal topography does not vary with y, so that the experiment is suitable for 2-D flowline models as well. The basal topography is thus formed by a series of ripples with an amplitude of 500 m: z b (x, y) = z s (x, y) − 1000 + 500 sin (ω x) .
3.3 Exp. C: Ice stream flow I The experiment setup is similar to Exp. A, except that the bedrock topography is flat, so that the ice thickness is spatially uniform (H =1000 m): where x ∈ [0, L] and L=160, 80, 40, 20, 10 and 5 km, respectively, and where α=0.1 • . The basal friction coefficient is prescribed as The β 2 -field is shown in Fig. 2. The basal friction oscillations have a frequency ω=2π/L.

Exp. D: Ice stream flow II
The only difference with Exp. C is that the basal friction coefficient does not vary with y, so that the experiment is suited for 2-D flowline models as well. The basal friction field is thus formed by a series of ripples defined as Note that in Exps. C and D the basal friction coefficient β 2 goes to zero within the domain.

Exp. E: Haut Glacier d'Arolla
Exp. E is a diagnostic experiment along the central flowline of a temperate glacier in the European Alps (Haut Glacier d'Arolla), based on earlier experiments by Blatter et al. (1998) and Pattyn (2002). Model input consists of the longitudinal surface and bedrock profiles of Haut Glacier d'Arolla, Switzerland, according to the Little Ice Age geometry (Fig. 3). The longitudinal profile of this glacier has a very simple geometry, hence the resulting stress field is not influenced by geometrical perturbations such as the presence of a steep ice fall. In a first experiment (E1), a zero basal velocity is considered (β 2 =∞), and the width of the drainage basin is kept equal to 1 along the entire flowline. The flow-law rate factor A=10 −16 Pa −n a −1 is assumed constant. Upstream and downstream boundary conditions imply zero ice thickness and velocity. The input file has a resolution x=100 m, but the authors were free to choose any grid resolution.
A second experiment (E2) considers a narrow zone of zero traction, similar to the experiment described in Blatter et al. (1998): The Cryosphere, 2, 95-108, 2008 www.the-cryosphere.net/2/95/2008/ 3.6 Exp. F: Prognostic experiment Exp. F is a prognostic experiment for which the free surface is allowed to relax until a steady state is reached for zero surface mass balance: where v h is the horizontal velocity vector (m a −1 ). Basic model setup differs from the setup in Exps. A and C. A slab of ice with mean ice thickness H (0) =1000 m is considered, resting on a sloping bed with a mean slope α=3.0 • (Fig. 4). This slope is maximum in x and zero in y. The bedrock plane is parallel to the surface plane and is perturbed by a Gaussian bump. Initial bedrock B (0) and unperturbed surface S (0) elevation are thus governed by where σ =10000m=10H (0) and where x, y (m) are the horizontal coordinates with respect to the center of the Gaussian bump. The basal perturbation has a maximum height of one-tenth of the mean ice thickness, i.e. a 0 =100=0.1H (0) (Fig. 4). The domain size L is taken to be 100 H (0) in x and y. The horizontal coordinates for output are scaled bŷ Periodic boundary conditions are applied in the horizontal. The major difference with the previous experiments is that n=1 in Eq. (5), so that the effective viscosity is constant and becomes η=(2A) −1 . Therefore, the unperturbed velocity field at the surface is defined by where τ b =ρgH (0) sin α is the unperturbed basal shear stress, and A=2.140373×10 −7 Pa −1 a −1 , so that Experiments are carried out for different values of the slip ratio c, which determines the relation between the basal velocity and basal drag. The basal velocity is written in terms of a basal friction coefficient β 2 , or Following the scalings given by Gudmundsson (2003), the basal friction coefficient is related to the slip ratio c by b =cU (0) . Table 2 lists the main constants used for Exp. F. Using these settings, the model should run until a steady state of the free surface is reached.

Model classification
A total of 27 numerical models and one analytical model from 20 contributors participated in the intercomparison exercise. Table 3 summarizes model characteristics and contributions. The different Stokes approximants extend the shallow-ice approximation in various ways (Hutter, 1983;Hindmarsh, 2004). We will follow here the classification scheme for higher-order models by Hindmarsh (2004), who gives a detailed description of the different "longitudinal stress schemes" widely used in ice sheet modeling. The most common longitudinal stress approximations introduce the two horizontal velocity components as field variables. This leads to an elliptic system with two rather than four variables of the full system at points in three-dimensional space (Pattyn, 2003;Hindmarsh, 2004), and the resulting linear systems are generally better conditioned than those resulting from the numerical analysis of the full system (Hindmarsh, 2004). These models are termed "multilayer models". A number of these models solve an elliptic system at one www.the-cryosphere.net/2/95/2008/ The Cryosphere, 2, 95-108, 2008  (2008)  oso1 SIA Hindmarsh (2004)  The L1L1 approximation is a one-layer longitudinal stress scheme using τ xx at the surface computed by solving elliptic equations and is identical to the approximation used by MacAyeal (1989). An alternative approximation is the L1L2 approach, or one-layer longitudinal stress scheme, usingε xx at the surface computed by solving elliptical equations with a vertical correction of τ xx . Here, the surface velocities used in computing the non-horizontal plane stresses are computed using the vertical shear stresses in the shear strain relationship and in the sliding relationship.
The most common approximation is the LMLa or multilayer longitudinal stress scheme. This is the classic longitu-dinal stress scheme used by Blatter (1995) and Pattyn (2003). Compared with L1L2, the longitudinal stresses use the velocity at the corresponding elevations rather than at the surface, and the stress-invariant calculations are self-consistent rather than using the SIA stress (Hindmarsh, 2004). Finally, there is the LTSML or multilayer longitudinal stresses scheme with horizontal shear stress gradients approximated by SIA. Here, horizontal gradients of the vertical velocity are neglected. Horizontal plane shear stresses, when needed to calculate the horizontal gradient of such shear stresses, are approximated by SIA values. This approach is similar to LMLa, but with inclusion of the vertical resistive stress, as in Van der Veen and Whillans (1989

Results
A graphical representation of all the results for each of the contributing authors as well as the submitted data files are found in the supplemental files http://www.the-cryosphere. net/2/95/2008/tc-2-95-2008-supplement.zip. The detailed description of the experimental setup is given in http://www. the-cryosphere.net/2/95/2008/tc-2-95-2008-supplement.zip. An analysis on the CPU performance of each of the experiments is presented in an accompanying paper (Gagliardini and Zwinger, 2008  The numerical results are displayed in both numerical and graphical format. A distinction is made between full-Stokes (FS) and non-full-Stokes (NFS) models, such as LMLa, LTSML, L1L2 and L1L1. All parameters refer to horizontal surface velocities, which are determined as the norm of the horizontal components of the velocity field, defined by ||v s ||= v 2 x +v 2 y . Velocities are displayed along a section y=L/4 for Exps. A and C and along a central line y=L/2 for Exp. F. For the 2-D experiments we show results along the flowline. The graphs show the mean velocity and its standard deviation along each section or flowline for both FS and NFS models. The tables present two parameters for each flowline or section: the maximum value of the norm of the surface velocity and the mean value of the velocity. For both parameters the mean and standard deviations are given for both FS and NFS models.

Experiments A and B
The results for Exps. A and B are shown in Figs. 5 and 6, respectively, and the numerical results are given in Tables 4  and 5. Each graph displays the norm of the surface velocity across the bumps at y=L/4 (for Exp. A) and along the central flowline (for Exp. B), for the different length scales and model groups (FS and NFS). The experiments were set up such that for this longitudinal profile the SIA gives a solution independent of L, which is not the case for higher-order models. The surface velocity according to the SIA is given by where v x (z b )=0 is the basal velocity (Fig. 7). The maximum surface velocity according to the SIA remains constant for all length scales (119.69 m a −1 ). However, whenever topographic differences occur, longitudinal stress gradients must develop which tend to smooth out the velocity field. For high aspect ratios =H /L (hence low values of L) this leads to a more or less constant surface velocity field as the ice sheet does not "feel" the individual bedrock undulations. Rather, it feels the fast sequence of large bed undulations as a viscous drag. The aspect ratio determines the amplitude of the horizontal surface velocity field, and the surface velocity decreases from ∼100 to ∼10 m a −1 .
Full-Stokes models closely agree with each other when calculating the velocity field for different length scales, compared to the larger spread of solutions for the higher-order approximations (Tables 4 and 5). Several factors could be responsible for the larger spread among NFS models: (i) More models are participating; (ii) These models are solving different continuum equations (LMLa, LTSML, L1L1, and L1L2); (iii) At the highest aspect ratios, the different approximations are not valid, so that the full-stress model needs to be solved; and (iv) There are numerical errors relative to the unknown exact solutions of the continuum equations. The FS results are only subject to a spread from the latter cause. For the smallest length scales the standard deviation for the full-Stokes models reduces to <0.2 m a −1 (Tables 4 and 5). The coarsest grid used by the participating models had dimensions 40×40, leading to a numerical mismatch of the order of 2.5%, which is far smaller than the standard deviation shown in Fig. 5 for the higher-order models. It is therefore very likely that the large spread associated with the higherorder models is due to the invalidity of the approximations compared to the full-Stokes solution.
The flowline experiments (Exp. B) show similar results compared to the 3-D experiments (Exp. A). There are differences, however, associated with the absence of transverse stress gradients. A peculiarity is observed in the surface horizontal velocity for the smallest length scale L=5 km. The surface velocity according to all full-Stokes models is larger over the bump than over the trough, hence anti-correlated with ice thickness (Fig. 6). All the other approximations (LMLa, LTSML, L1L2 and L1L1) predict that ice velocity and thickness are positively correlated for all length scales. The marked difference can be attributed to mass conservation, as at such high aspect ratios the horizontal ice flux cannot be balanced by the vertical flux at the free surface since the vertical velocity would be too large for the given depth. The horizontal ice flux is therefore more or less constant, inducing larger velocities for smaller depth and vice versa. This behavior is noticeable only for the flowline experiments, because in 3-D the ice can flow around the bumps. The flow inversion is an artifact stemming from the diagnostic nature of the experiments and would disappear if the free surface were allowed to respond to the applied stress field. Higherorder models fail to produce the velocity inversion, since the stress field is determined solely from horizontal strain rate components and vertical velocities are an a posteriori model result (e.g. Pattyn, 2003

Experiments C and D
In this series of experiments, variations in basal conditions (slipperiness) determine where longitudinal stress gradients must develop. Due to the importance of basal sliding, the ice behaves as in an ice stream, in which vertical shearing is present, though minimal. Ice flow in this experiment can be considered as ice shelf flow with minimal basal traction. The invalidity of the SIA solution is shown in Fig. 7, where the analytical SIA solution is plotted for a simplified basal sliding relationship in which the basal shear stress is supposed to balance the driving stress without longitudinal stress gradients, so that in Eq. (27). However, a singularity in the velocity occurs in Eq. (28) for β 2 =0. Not all velocities are therefore plotted in Fig. 7. As for Exps. A and B, the SIA solution is independent of L. Again, full-Stokes models show a smaller spread than the higher-order approximations (Figs. 8-9 and Table 4). The SIA is definitely not suited for simulating this type of ice flow where longitudinal stresses dominate over vertical shear stresses (compare with Fig. 7). As in Exps. A-B, the amplitude of the surface velocity field decreases with increasing aspect ratio . At the smallest length scales the surface velocity field is almost constant as the transition between high and low friction regions is smoothed out by longitudinal stress transmission. Similar to Exp. B, a surface velocity field anticorrelated to basal friction is observed for full-Stokes models at L=5 km (Exp. D). For the other higher-order approximawww.the-cryosphere.net/2/95/2008/ The Cryosphere, 2, 95-108, 2008 tions this is not observed (but due to the larger disparity in solutions, this effect is unnoticeable in Fig. 9).
In general, the spread in results of the modeled velocity field is higher than for the experiments over the bedrock bumps. The smallest spread is obtained with full-Stokes models, and this spread reduces with increasing , contrary to the results from Exps. A-B. L1L1 and L1L2 models have larger spreads than the LMLa models, despite the fact that they were designed for coping with such type of ice flow in the first place.

Experiment E: Haut Glacier d'Arolla
Although the input file lists the bedrock and surface data along the flowline of Haut Glacier d'Arolla with a fixed grid spacing of x=100 m, most participants interpolated this dataset at a higher resolution (Fig. 10). The effect of subsampling is captured in Fig. 11, where the oscillations in the basal shear stress along the flowline are either jagged when undersampled or smoother when a sufficiently small grid size is chosen. Again, the spread in surface velocity for full-Stokes models is smaller than for the other approximants, albeit that for the no-slip case, the standard deviation is small for all models.
Inclusion of a sliding zone (an area of zero basal friction) leads to larger differences among the different participating models. Also the full-Stokes models show a much larger spread of solutions. Here, increasing resolution results in other complexities compared to the no-sliding case, such as the occurrence of oscillations in the basal shear stress. The slip/no-slip boundaries are very sensitive to model resolution, as they can be regarded as singularities where the friction parameter β 2 suddenly jumps from zero to infinity and vice versa. In particular, the linear interpolation leads to break points in basal and surface topography that influence the result. The results of the sliding experiment underline the difficulty of simulating end-member behavior in basal sliding (slip/no-slip). Benchmarking of numerical ice sheet models is possible when analytical solutions exist for a particular problem. The analytical solutions used here are derived from first-order perturbation analysis of flow down a uniformly inclined plane (Gudmundsson, 2003). It is inherent in this type of analysis that the resulting flow perturbations are linear functions of basal amplitudes. Numerical solutions are usually not limited by this assumption and can therefore, for any finite amplitude perturbation in basal properties, be better approximations to the Stokes equations than the analytical solutions given by Gudmundsson (2003). For an accurate comparison with the analytical solutions, numerical solutions must be calculated for a number of different amplitudes and then scaled by forming the ratio between each solution and the respective basal amplitude. If this ratio is found to be independent of amplitude for small amplitudes, the scaled numerical solutions can be compared to the analytical ones. This kind of test was done by Raymond and Gudmundsson (2005). The exact error estimate depends on wavelength, amplitude, and slip ratio. (The reader is referred to Raymond and Gudmundsson (2005) and Gudmundsson (2008) for a detailed discussion.) For example, the analytical solutions were found to be generally accurate to within ∼1% for sinusoidal pertubations with wavelength larger than ∼10H and amplitudes less than ∼0.1H . For Exp. F, we expect a similar degree of agreement between the analytical solutions and exact Stokes solvers. Only a few models participated in this experiment, since not all of them treat an evolving free surface. The objective of the test was to run the models forward in time until a steady state was reached. The definition of steady state was left to the interpretation of each participant. The resulting steadystate surface elevation and velocity are shown in Figs. 12-13.
We present the analytical solution here for the case without basal sliding. The two full-Stokes models show a very good agreement with the analytical solution. Higher-order model solutions also fit also well but show a slightly larger vari-The Cryosphere, 2, 95-108, 2008 www.the-cryosphere.net/2/95/2008/  of the full-Stokes equations. However, the analytical solution for the velocity is well away from the tightly-clustered FS numerical results. Here, the FS solution might be more accurate, as the analytical solution remains a solution based on a perturbation expansion and therefore is an approximate method for solving the full-Stokes equations.

Conclusions
In this paper we present the results of the first intercomparison exercise of higher-order and full-Stokes ice sheet models. A total of 27 different numerical models participated in this benchmarking effort. Six experiments were designed to evaluate complex ice flow with high spatial variability in basal topography and slipperiness. All experiments were designed in such a way that the Shallow-Ice Approximation (SIA) is not valid, especially at high aspect ratios. Although the SIA is valid for large parts of ice sheets, higher-order models are necessary for describing ice flow in areas of high basal topographic variability and slipperiness, which are generally the most dynamic regions of ice sheets. Compared to previous benchmark experiments MacAyeal et al., 1996;Payne et al., 2000), a significantly higher number of ice sheet models participated in this benchmark. Despite the greater complexity of the problem, all models produce results that are in approximate agreement, even for high aspect ratios. This shows that numerical ice sheet models have improved considerably during the past decade, and are capable of simulating ice flow in regions where longitudinal stress gradients are important.
Full-Stokes models produce reliable results in the sense that (i) their spread of results is very low (<0.2 m a −1 ) and (ii) they give a result in concordance with analytical solutions based on perturbation theory. Results of the higherorder approximants show a significant larger spread that cannot be attributed solely to numerical issues, such as the discretization error. The greatest spread is found for high aspect ratios where all stress components (not only membrane/longitudinal stresses) are equally important, and higher-order approximations are insufficient. Also, coding errors could be present, since all higher-order models are coded by the authors (which is less true for the full-Stokes models).
Only a limited number of L1L1 and L1L2 models participated. They show the same spread in solutions for basal topographic perturbations (Exps. A and B) as the LMLa models. However, they produce a larger spread for the ice stream simulations at relatively low aspect ratios (Exps. C and D) compared to LMLa models, although they were designed for such flows. The LMLa results differ from the full-Stokes results The Cryosphere, 2, 95-108, 2008 www.the-cryosphere.net/2/95/2008/ for high aspect ratios and in regions of variable basal sliding. All models (including full-Stokes) agree poorly when sudden variations in basal friction are considered, such as the slip/no-slip jumps in Exp.E. Finally, the full-Stokes models give the most consistent results. Results from different models show a very small spread and are in good agreement with the analytical solution based on perturbation theory. For most experiments a clear distinction can be made between results from full-Stokes models and from higher-order approximants.
Future model intercomparisons should definitely focus on transient problems, as preliminarily explored in Exp. F, and slip/no-slip boundary problems, such as Exp. E. Analytical solutions to specific (simplified) problems are something to look into as well, which would lead to more reliable benchmarks and verification.