Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity

Introduction Conclusions References

These solutions are constrained for a number of well defined model setups, which are of potential significance for various situations which numerical codes may face during real geodynamic simulations.Availability and broad range of 2-D and 3-D benchmark solutions are, therefore, critical for the development and testing of the next generation of numerical geodynamic modeling software which aims to combine rheological complexity of constitutive laws with adaptive grid resolution to on both global and regional scales (e.g., Moresi et al., 2003;Dabrowski et al., 2008;Tackley, 2008;Stadler et al., 2010;Gerya et al., 2013).
In the present paper we aim to significantly expand availability of benchmark solutions for both 2-D and 3-D variableviscosity Stokes flows.In contrast with previous studies, we prefer not to start from any prescribed model setups but rather derive general analytical solutions, which are potentially suitable for generating a broad range of test problems.We derive generalized solutions for incompressible Stokes problems with linearly and exponentially variable viscosity.In the following we demonstrate how these generalized solutions can be converted into 2-D and 3-D test problems suitable for benchmarking numerical codes.Finally, based on the obtained benchmark problems, we show examples of numerical convergence tests for staggered-grid discretizations schemes (e.g., Gerya andYuen, 2003, 2007).

Formulation of 2-D equations with variable viscosity
Consider the plane (3) Here (v x , v y ) is the flow velocity, η = η(x, y) is the viscosity, P is the pressure, ρ is the density, (G x , G y ) is the gravitational force.Note that Eq. ( 3) is the continuity equation.
Let us change the variables v x , v y , P in such a way that The correctness conditions for such replacement are as follows: ∂ ∂y These conditions lead to the following correlations: If we consider these conditions as a partial differential equations then we obtain the same characteristic equation for all these conditions: Evidently, η(x, y) = C is an integral of the equation.It is well known that an integral of the characteristic equation gives one good new variable.Namely, η = η(x, y) is a good new variable (the second coordinate should be orthogonal to the first one).Note that the assumption that the viscosity varies is now crucial because we use the viscosity as a new coordinate (instead of the standard Cartesian spatial coordinate).Hence, the solutions of our equations, which predetermine the correctness of the replacement suggested above, are u x = (η), u y = (η), P = P (η).
After the replacement, the Stokes Eqs.(1, 2) and the continuity condition (3) transform to the following form: Inserting the expressions for u x , u y into Eqs.(7-9), one obtains the following equations:

Linearly varying viscosity
To obtain the solutions of Eqs.(10-12) we make some assumptions.In particular, in this section we assume, first, that the viscosity η is a linear function of the Cartesian coordinates, where a, b, c are non-zero constants.Second, the restriction concerning the gravitational terms takes place: ρ(aG y − bG x ) and ρ(bG y + aG x ) are assumed to be functions of one variable -viscosity.It means that ρG y and ρG x are functions of η only.Introduce functions f, f 1 : Under these assumptions one gets the following system of equations: One can simply solve this linear algebraic system with respect to , , P : That gives us Correspondingly, one obtains v x , v y , P : Finally, Analogous transformation takes place for v y : The expression for the pressure is as follows:

I. Yu. Popov et al.: Benchmarking of geodynamic Stokes problems with variable viscosity
In particular, in the case of constant gravitational terms, i.e., for f (η Consider a more complicated case when the density is a linear function of the viscosity: ρ = β 1 η + β 2 .Then, where constants a 1 , a 2 , b 1 , b 2 are as follows: Note that in this case the continuity Eq. ( 3) should be rewritten in a more general form: Equations ( 16)-( 18) gives us The continuity Eq. ( 19) gives us the relation between c 2 , c 3 :

Exponentially varying viscosity
Let us construct the second benchmark solution.Now we assume that the viscosity is the exponential function of the Cartesian coordinates: General consideration up to Eqs. (10-12) is the same as earlier.By inserting Eq. ( 20) into Eqs.(10-12) and taking into account that one obtains the following system of equations: Using the last relation, we exclude from the first two equations: One can see that we obtain a linear algebraic system with respect to ( η 2 + η) and P .The solution is as follows: Remark.It is interesting that these formulas contain the same functions f (η), f 1 (η) as in the previous section.Equation ( 23) is a well-known Euler ordinary differential equation.One can get its solution for arbitrary function f : Taking into account the relation in Eq. ( 21), one obtains u y : Taking into account Eqs.(4, 5) one obtains v x , v y : Hence, we get the expression for v x and analogously, for v y : As for the pressure, we obtain it from Eq. ( 22) by taking into account Eq. ( 6): Hence, One can compare Eqs.(26-28) with Eqs.(16-18).For a simple particular case (constant gravitational term) when For more a complicated case when the density is a linear function of the viscosity ρ where constants a 1 , a 2 , b 1 , b 2 are the same as in the previous section.The continuity equation should be written in a more general form Eq. ( 19).It is simple to evaluate integrals in Eqs.(26)(27)(28).In such a way one obtains The continuity Eq. ( 19) gives us the same relation as earlier: 3 Three-dimensional case

Formulation of 3-D equations with variable viscosity
The situation in the 3-D case is similar to that for 2-D.In particular, we can realize the same procedure as in the 2-D case with some additional restrictions.The initial system of equations is as follows: The replacement of the variables v x , v y , v z , P by u x , u y , u z , P is analogous to that in the two-dimensional case: The characteristic equations for each triplet of equations, i.e., for u x , u y , u z , P are the same (analogously to 2-D case).The solutions for these conditional equations are u x = (η), u y = (η), u z = (η), P = P (η).(32) These relations give us correctness of the above introduced replacement for arbitrarily varying viscosity.However, it allows one to obtain the solution of the Stokes continuity equations only under some assumptions concerning the viscosity.

Linearly varying viscosity -3-D case
Below we will assume (analogously to 2-D case, Eq. 13) that where a, b, c, e are non-zero constants.Moreover, in the 3-D case we need the following restriction for the gravitational terms (compare with Eq. 14): ρG x , ρG y , ρG z are functions of one variable -viscosity (for example, these terms may be constant).
Making the substitution (Eq.31) in the system (Eq.29), one transforms it to the form, Substitution of Eq. ( 32) into Eq.(??), Eq. ( 35) leads to the following system for , , , P : Hence, The first three equations gives one a linear algebraic system with respect to , , P which can be solved without difficulties.Then, one gets from the last equation of the system.To be correct, the obtained expressions should be functions of one variable -viscosity.It is ensured by our assumption concerning the gravitational terms.The result is as follows: Here we defined four functions: f x (η), f y (η), f z (η), f P (η).
The obtained expressions are analogous to that in Eq. ( 15).Correspondingly, the integration is analogous, and one obtains In particular, in the case of constant gravitational terms, i.e., for one obtains a result analogous to the 2-D case: The continuity condition gives one the following correlation between the coefficients: The condition, is identically valid (see the expressions for f x , f y , f z ).Consider a more complicated case when the density is a linear function of the viscosity: ρ = β 1 η + β 2 .Then, where constants a 1 , a 2 , b 1 , b 2 , d 1 , d 2 , p 1 , p 2 are as follows: The continuity Eq. ( 30) in this situation has more general form: In this case the formulas for the velocity and the pressure preserve their form from Eqs. (44-47) and give us The continuity Eq. ( 49) leads to the relation ac 2x + bc 2y + cc 2z = 0.

Exponentially varying viscosity -3-D case
Consider the second benchmark solution (for exponential dependence of the viscosity on the Cartesian coordinates): The assumption concerning the gravitational terms is the same as earlier.General 3-D consideration is the same.By inserting Eq. ( 32) into the equations for u x , u y , u z , P and taking into account that one obtains the following system of equations: We can note that the first three equations give us the same linear algebraic system as in the case of linear viscosity (Eqs.36-38), if one takes ( η 2 + η), ( η 2 + η), ( η 2 + η), η P as variables instead of , , , P in the linear viscosity case.Hence, the solution of the system is as follows: The definition of the functions f x , f y , f z , f p was given above, see Eqs. (40)(41)(42)(43).One can solve these equations by the procedure used in the corresponding 2-D case (exponential viscosity).By this method, we come to the result For a particularly simple case (constant gravitational term) when f x (η) = A x = const, f y (η) = A y = const, f z (η) = A z = const, f P (η) = A P = const one has: Consider a more complicated case when the density is a linear function of the viscosity: ρ = β 1 η + β 2 .Then, where constants a 1 , a 2 , b 1 , b 2 , d 1 , d 2 , p 1 , p 2 have been determined in the previous section.Here we use the continuity equation in the form of Eq. ( 49).In this case Eqs.(51-54) give us The continuity Eq. ( 49) leads to the relation

Example problems and numerical convergence tests
The scheme of algorithm testing is as follows.Initially, we have obtained particular solutions of the Stokes and continuity equations for two types of viscosity variations.Let us choose a domain, e.g., a rectangle in the 2-D case.We calculate the values for velocity and pressure given by our analytical solution and take these values as the boundary conditions.

Linearly varying viscosity
Consider a simple example of such flow in a rectangle : 0 ≤ x ≤ x size , 0 ≤ y ≤ y size .We assume that η = ax +by +c.
We will mark the exact solution obtained in Sect. 2 as v x,a , v y,a , P a .It is the solution of the boundary problem in the rectangle with the following conditions at the boundary ∂ = {x = 0, x = x size , y = 0, y = y size } : Let us compute the velocity and the pressure by the finitedifference scheme.The corresponding solution is marked as v x,n , v y,n , P n .The deviation of these values from the exact  (2010).Calculations show that the convergence is better for the case of low-viscosity contrast (cf., Supplement).The results for high-viscosity contrast are presented in Figs.1-3: Fig. 1 shows prescribed viscosity and density distribution, Fig. 2 presents the pressure and the velocity components distributions and Fig. 3 contains the plot of the relative errors via the grid resolutions in logarithmic scale.The viscosity contrast, i.e., the values of the coefficients in the expression for the viscosity, is determined by the given values of the viscosity at three corners of the model rectangle. . .The value of the viscosity at the initial rectangle corner is 1, whereas η 2 and η 3 are the prescribed values of the viscosity at the upper-left and lower-right corner, respectively.In all figures "n" means "numerical solution", "a" means "analytical solution" (benchmark).
For the case of linearly varying viscosity we made calculations for the following system parameters: One can see that the numerical approach has rather high accuracy.We observe the conventional situation -L ∞ -error is the largest among the considered errors norms, and L 1 -error and L 2 -error are similar.To describe in more details the dependence of the error on the viscosity contrast, we fill tables with errors for different values of η 2 , η 3 (see Appendix, Tables A1-A3).

Exponentially varying viscosity
The case of exponentially varying viscosity is treated analogously.In order to make a comparison with the case of linearly varying viscosity, we take the same system parameters (geometrical size, gravitational terms and the dependence of the density on the viscosity) with the same viscosity contrast (i.e., the values of the viscosity at the rectangle corners): We examine two cases (low-and high-viscosity contrasts, cf.Supplement).Figures 4-6 present the results for highviscosity contrast.There are some similarities with the previous case.In particular, L ∞ error norm gives us the maximum relative error value among the three considered error norms.However, peculiarities are more interesting.Namely, the numerical scheme works essentially better for the case of exponentially varying viscosity than for linearly varying viscosity.We observe good convergence for both low-and high-viscosity contrasts (compare Figs. 3 and 6,also cf. Supplement).
Dependence of the errors on the viscosity contrast for exponentially varying viscosity is presented in the Appendix Tables A3-A6.

3-D example
One can note that benchmark solutions for the 3-D case are essentially more rare than for the corresponding 2-D situation.It is, therefore, noteworthy that the suggested approach allows us to obtain such solutions for 3-D Stokes and continuity equations.As earlier, we consider the cases of linearly and exponentially varying viscosity.

Linearly varying viscosity
The statement of the problem is analogous to the previous case.In the 3-D space we consider a parallelepiped : 0 ≤ x ≤ x size , 0 ≤ y ≤ y size , 0 ≤ z ≤ z size .We assume that η = ax + by + cz + e.We will mark the exact solution obtained in Sect. 3 as v x,a , v y,a , v z,a , P a .Due to the uniqueness theorem, it is the solution of the boundary problem in the parallelepiped with the following conditions at the boundary ∂ = {x = 0, x = x size , y = 0, y = y size , z = 0, z = z size } : Let us compute the velocity and the pressure numerically using chosen finite-difference algorithm.The corresponding numerical solution is marked as v x,n , v y,n , v z,n , P n .The deviation of these values from the exact solution (v x,n − v x,a , v y,n − v y,a , v z,n − v z,a , P n − P a ) is related to the er- ror of the numerical scheme.As in the 2-D case, we consider three error norms: L ∞ , L 1 , L 2 .We examine cases of low-and high-viscosity contrasts.Coefficients a, b, c, e in the viscosity formula are determined by given viscosity values (η 1 = 1, η 2 , η 3 , η 4 ) at four adjacent parallelepiped vertices.
We choose the following system parameters: Results are presented in the Supplement.It should be mentioned that the numerical procedure for the 3-D case takes essentially greater time than for the 2-D case.Qualitatively, the results are similar to that for the corresponding 2-D case.Better numerical convergence is again found for low-viscosity contrast.

Exponentially varying viscosity
As in the 2-D case, we consider also exponentially varying viscosity.The consideration is analogous to that in the respective 2-D case.In order to make a comparison we take here the same system parameters.Naturally, the viscosity and, correspondingly, the density distributions are now exponential.
The system parameters are chosen as follows: The results are presented in Figs.7-10 for the case of highviscosity contrast.Similarly to 2-D results, the numerical procedure converges better for the exponentially varying viscosity than for the linearly varying viscosity.The algorithm convergence, i.e., the dependence of the errors on the grid resolution is shown in Fig. 10.

Conclusions
In this paper, we developed new, specific analytical solutions for the 2-D and 3-D Stokes flows with both linearly and exponentially variable viscosity.We also demonstrated how these solutions can be converted into 2-D and 3-D test problems suitable for benchmarking numerical geodynamic codes.The main advantage of this new generalized approach is that large variety of benchmark problems can be easily generated, including relatively complex cases with open model boundaries, non-vertical gravity and variable gradients of viscosity and density fields, which are not parallel to Cartesian axes.These solutions can be very useful for testing numerical algorithms aimed at modeling variable-viscosity mantle convection and lithospheric dynamics.Examples of respective 2-D and 3-D MatLab codes are provided with the paper.
The Supplement related to this article is available online at doi:10.5194/se-5-461-2014-supplement.

Appendix A: Error decreasing with decreasing grid resolution
The dependence of the convergence on the viscosity contrast is shown in the following tables.Namely, the viscosities η 2 , η 3 run through a set of values from 2 to 10 000.For each pair of η 2 , η 3 we perform the benchmarking procedure described above and obtain the curve describing the dependence of the logarithm of the error norm on the logarithm of the grid resolution.The tangent of the slope angle of this curve gives us an input of the table.Full set of tables are given in the Supplement.Below, we show several examples of such tables.

A2 Exponentially varying viscosity
Tables A4-6 contain log rates of L 2 error norm decreasing with decreasing grid step (the curve slope in Fig. 6 in the logarithmic scale) for different viscosity contrasts in the case of exponentially varying viscosity.The system parameters are the same as in the case of linearly varying viscosity with respective changes: by Copernicus Publications on behalf of the European Geosciences Union.I. Yu.Popov et al.: Benchmarking of geodynamic Stokes problems with variable viscosity -2-D thermomechanical corner flows in subduction zones (van Keken et al., 2008); . Popov et al.: Benchmarking of geodynamic Stokes problems with variable viscosity Here c 1x , c 1y , c 1z , c 2x , c 2y , c 2z , c p are constants.The continuity equation gives one a relation between the coefficients: ac 1x + bc 1y + cc 1z = 0.One can compare Eqs.(51-54) with the results for the corresponding 2-D case (Eqs.26-28).

Table A1 .
Tangent of the slope of the curve showing the dependence of the logarithm of L 2 -error for pressure P on the logarithm of the grid step for different viscosity contrasts.Linearly varying viscosity.

Table A2 .
Tangent of the slope of the curve showing the dependence of the logarithm of L 2 -error for velocity component v x on the logarithm of the grid step for different viscosity contrasts.Linearly varying viscosity.

Table A3 .
Tangent of the slope of the curve showing the dependence of the logarithm of L 2 -error for velocity component v y on the logarithm of the grid step for different viscosity contrasts.Linearly varying viscosity.

Table A4 .
Tangent of the slope of the curve showing the dependence of the logarithm of L 2 -error for pressure P on the logarithm of the grid step for different viscosity contrasts.Exponentially varying viscosity.

Table A5 .
Tangent of the slope of the curve showing the dependence of the logarithm of L 2 -error for velocity component v x on the logarithm of the grid step for different viscosity contrasts.Exponentially varying viscosity.

Table A6 .
Tangent of the slope of the curve showing the dependence of the logarithm of L 2 -error for velocity component v y on the logarithm of the grid step for different viscosity contrasts.Exponentially varying viscosity.