Equal-order finite element approximation for mantle-melt transport

Mantle convection and melt migration are important processes for understanding Earth’s dynamics and how they relate to observations at the surface. Recently it has been established that melt migration can be modelled by coupling variable-viscosity Stokes flow and Darcy flow, where Stokes flow generally captures the long-term behaviour of the mantle and lithosphere, and Darcy flow models the two-phase regime. It is known that approximating the solution by finite element methods requires the use of mixed inf-sup stable elements or additional stabilization terms. Here, we propose a formulation with a coercive non-symmetric linear operator which allows the use of simple equal-order elements.


Introduction
The need to solve systems of coupled Navier-Stokes and Darcy flow arises is various fields, such as modelling the interaction of surface water and groundwater aquifers [1,2], blood flow problems [3,4], and fuel cell dynamics [5,6]. Approximating the solution of coupled Stokes-Darcy systems in most methods requires solving Stokes or Darcy on adjacent sub-domains coupled with appropriate interface boundary conditions [7][8][9][10], and [11]. In melt migration modelling the partially molten rock cannot be handled using the approach above because we cannot identify the boundary between the solid and fluid phases. The model derived by Mc McKenzie [12] for melt migration assumes a dual continuum mixture of solid matrix and fluid melt. The mixing parameter is the porosity ϕ, i.e. the volume fraction of fluid melt, which is assumed to be much smaller than 1 and may be zero in parts of the domain where there is no fluid melt. The motion of the solid is governed by Stokes flow, and the melt is transported according to Darcy's law. In addition, he included a compacting relation that relates the solid and fluid pressures.
Most of the previous studies of geodynamics do not consider that melt migration can be modelled by coupling variable-viscosity Stokes flow and Darcy flow, where Stokes flow generally captures the long-term behaviour of the mantle and lithosphere, and Darcy flow models the two-phase regime. Their studies either do not consider melting or treat it in a simplified way [13][14][15][16]. Others have simplified the model by approximating the two-phase flow and the effects of compaction [17][18][19], and [20]. The correct setting was later proposed by [21]. However, the required computational effort was limited to 2D problems, and although these models take into account the compaction of the solid matrix to allow melt to be expelled or to flow in, they treat both individual phases as incompressible and often assume their densities to be constant.
The discretization of a similar (but not identical) complex melting problem is suitable described by Schimenz et al. [33] by using a mixed discontinuous Galerkin method in vertical direction and a Fourier method in vertical direction for the elliptichyperbolic system of equations.
In [23], the authors derive such a model for 2D and 3D simulations. However, the model they used require the so-called compaction pressure as an additional unknown in the system which not only increases the dimension of the system, but also makes the underlying operator non-coercive (see, [21][22][23]). Thus, approximating the solution using the finite element method requires the use of mixed inf-sup stable elements or additional stabilization terms. Here, we propose a formulation with a coercive nonsymmetric linear operator which allow the use of simple equal-order elements.

Governing equations
Let Ω ⊂ R d , d ∈ {2, 3}, be a Lipschitz-domain. The model we consider here is a two-phase flow consisting of melt and a deformable matrix. The melt fraction ϕ defines averaged quantitiesX out of solid (matrix) X s and fluid X f quantities (the subscripts f and s stand for fluid and solid, respectively): The governing equations for a two-phase flow were given by McKenzie in [12]. They include the equations for conservation of mass, momentum, and energy. The mass conservation for fluid and solid are described by the two equations Where, f and s are the fluid-and solid densities, u f and u s the fluid-and solid velocities, respectively, t the time, and Γ the melt production rate following a volume element of matrix. The conservation of momentum of the melt is given by with the constant K D = k ϕ /η f , where k ϕ is the permeability and η f is the melt viscosity, p f the pressure within the melt, g the gravity acceleration vector. For the case of vanishing solid velocity, u s = 0, Eq. (4) reduces to the well known Darcy's law. Further quantities and associated symbols used in this paper are given in Table 1.
The conservation of momentum of the matrix is given by where T is the strain tensor given by with the shear viscosity of the solid η, the bulk viscosity of the solid ξ , and the identity matrix I .
Thus, the two-phase PDE system reduces to the original two-phase PDE system given by McKenzie [12]: with α, the weighted average of the flux of compressibility of solid and fluid densities with respect to the corresponding pressures, and transport-scaling parameter β given by

Time dependent melt fraction
Although we previously assumed that the flow field is in temporal equilibrium, this is not necessarily the case for the melt fraction ϕ. The corresponding equation is given by Assuming s > 0 a.e. in Ω, this equation can also be written in the form

Boundary conditions
The system of Eqs. (10), (11) and (15) has to be supplemented with suitable boundary conditions. We impose Dirichlet conditions for u s and the lithostatic pressure gradient for p f : The boundary condition for ϕ are of Dirichlet type on the inflow part of the boundary ∂Ω − := {x ∈ ∂Ω : u s · n < 0}.

Compaction pressure
In order to solve the PDE system (10), (11) numerically several authors treat it as a system of three-field saddle point problems [21,23], and [22] by introducing a secondary pressure p c , called compaction pressure, defined as Then, Eq. (10) reduces to the form Where, the reduced strain tensorT u s is given bỹ The resulting system (11), (18) and (19) now has three unknowns u s , p f , and p c . The corresponding system is of saddle point structure, so that its finite element discretization is not possible without considering inf-sup stable elements or special stabilization techniques, for instance the pressure stabilized Petrov-Galerkin method (PSPG) or the local projection pressure stabilization (LPS). We refer to [24,26,29], and [25] for details. In this work we will follow a simpler approach by dealing with the original equations. So, the compaction pressure p c will not be a primary variable. In case this quantity is needed for special purposes, it can be recovered from u s by solving the discrete system corresponding to (18), as we shall explain below.

Weak formulation and discretization for constant melt fraction
Let us first consider the case that the melt fraction ϕ is given and constant in time. Then the density becomes a given constant and we only have two unknown variables p f and u s , determined by a linear stationary system of equations, i.e. without any time dependency. Moreover, let us assume for ease of presentation that we have homogeneous Dirichlet data for the velocities, u 0 = 0. The extension to non-homogeneous Dirichlet data is straight forward.

Weak formulation
The associated function spaces are

As bilinear form
The right hand side is given by In order to have this bilinear form A and the right hand side F[¯ ] to be well-defined, it is sufficient to make certain regularity assumptions on the involved coefficients. This will be discussed below. We seek u s ∈ V and p f ∈ Q s.t.
Assumption 1 We assume the following regularities for the coefficients and the partial densities are bounded from below as The regularity assumptions for s , f , κ s and κ f are sufficient to ensure that α determined by (12) has the regularity α ∈ L 3/2 (Ω). This can be used to show the well-posedness of the term (αu · g, q) in the bilinear form A: Here we used the Sobolev embedding H 1 (Ω) → L 6 (Ω), for d = 2, 3.

Lemma 1 Under the regularity assumption (Assumption 1) the linear form F
Proof Let us now check the individual term of F[¯ ] in view of appropriate definition. The boundary integral is well-defined due to the assumed continuity of K D , g, ρ s and the trace inequality: The regularity¯ ∈ H −1 (Ω) is obtained by Eq. (1) and ϕ ∈ L 3/2 ⊂ H −1 (Ω) and . It remains to show an upper bound for the dual pairing of forcing By the trace inequality we can bound ||q|| L 2 (∂Ω) by the H 1 -norm of q. The Hölder inequality and Sobolev embedding H 1 (Ω) → L 6 (Ω) yields Remark 1 Even a stronger assumption ϕ ∈ L ∞ (Ω) is not a severe restriction, because, from the physical point of view, 0 ≤ ϕ ≤ 1 is reasonable. The lower bound for f and s are necessary for ensuring that the right hand side f is properly defined.

Existence and uniqueness of solutions
Throughout this section we presuppose that Assumption 1 is valid. In order to prove the continuity and the coercivity of the bilinear form A(u, p; v, q), and the continuity of the right hand side linear form, we equip the product space X with the norm: In order to show the coercivity of this bilinear form, we first show the non-negativity of the strain tensor when it is tested diagonally: Proof Because of I : ∇u = ∇ · u we obtain a. Let us first assume, that u ∈ (H 1 0 (Ω) ∩ C 2 (Ω)) d . Due to this regularity of u the terms in the sum above can be integrated by parts: Using this in the equation derived above, yields the assertion. b. For less regular u, the same equation is obtained by a density argument: We take the limit of a smooth limiting sequence

Lemma 4 (Coercivity)
We assume that with a certain constant c = c(η, K D ). Then the bilinear form A is continuous and X-coercive; in particular Proof The assertion follows by diagonal testing, use of the previous Lemma, and the fact that the convective term vanishes when tested diagonally i.e. (g ·∇ p, p) = 0. The only critical (not necessarily positive) term is the zero order reaction term (αu · g, p). However, due to the assumption on α, this term can be bounded by (23) and Young's inequality as follows This leads to the desired lower bound for A(u, p; u, p).

Remark 2
The constant c in (25) Proof Taking into account that the quantities¯ and f are functions of ϕ and Γ , respectively, we see that the H −1 -norms of¯ and f are well-defined. Now, the assertion is an immediate consequence of the previous Lemma and the theorem of Lax-Milgram.

Equal-order finite elements
The coercivity of the problem allows us to use several standard equal order elements for the discrete velocity u h and the discrete pressure p h . Let T h be a shape regular partition of Ω into d-dimensional simplices, quadrilaterals or hexahedra [30]. The diameter of a cell K ∈ T h will be denoted by h K and the mesh parameter h = max{h K |K ∈ T h } represents the maximum diameter of all cells. Let S h ⊂ H 1 (Ω) be a finite element space of continuous, piecewise polynomial functions defined over T h , with degree of interpolation order r ≥ 1. We consider triangulations T h of shape regular elements, so that an interpolation operator i h : H 1 (Ω) → S h exists (Scott-Zhang interpolation [31]) with the following properties for all 1 ≤ l ≤ r + 1, all K ∈ T h , and all u ∈ H l (ω K ): Here, ω K denotes a union of cells in the neighbourhood of K , and the expression a b means a ≤ cb with the constant c independent of the mesh parameter h.
Possible choices on shape regular meshes are P r or Q r elements, or finite element spaces containing these spaces, e.g. locally enriched spaces. Here, we consider equalorder finite element spaces for velocity and pressure: Now, the linear system to be solved reads as follows: The main advantage of the equal-order discretization of the problem is that the degrees of freedom of all variables can be assigned to the same geometrical identities, e.g. to vertices. Therefore, the discrete systems can be designed in a block-wise manner, so that an efficient block preconditioner or block smoother for multigrid solvers can be utilized.

Theorem 2 We make the same assumptions as in Theorem 1. The discrete solution of Eq. (28)
Proof Due to coercivity (Lemma 4) we can apply Cea's lemma, see e.g. [30]: with arbitrary interpolation i h : H l (Ω) d → V h and i h : H l (Ω) → Q h . For l = 1 we choose the Scott-Zhang interpolant, and for l ≥ 2 we take the nodal interpolant. This leads e.g. to with a constant c depending only on K D and Ω. The other parts in the norm ||| · ||| are obtained analogously.

Post-processing of the compaction pressure
Working with equal-order elements on simplices, i.e. with P r elements, the divergence of the discrete velocity, ∇·u h , is piece-wise polynomial of order r −1 but discontinuous across element faces/edges. Hence, for constant (or cell-wise constant) bulk viscosity, the discrete compaction pressure p h c := −ξ ∇ · u h is of the same type. Therefore, the evaluation of this quantity inside the cells is straight forward. The situation is different for its evaluation on vertices or edges/faces due to the discontinuity. If the discrete compaction pressure is required on the vertices, denoted by p h c , a common strategy is to define it by the L 2 -projection onto Q h ,

Lemma 5 Under the same assumptions as in Theorem
Proof Letp h c ∈ Q h be the solution of the problem

Then we have by Cea's Lemma and standard interpolation results
By stability of the discrete equations and Theorem 2, we have with a constant C depending on ξ and η. With the triangle inequality we arrive at the desired estimate.
For more regular compaction pressure, p c ∈ H l (Ω), we see that the obtained accuracy for the compaction pressure is of one order less (with respect to the mesh size h) than the optimal interpolation error. However, there are well known methods to increase the accuracy by special gradient recovery techniques as an alternative to solving Eq. (29). We refer to the classical work of Zienciewicz and Zhu [32].

Variational formulation and discretization for variable melt fraction rate
In this section we propose a second order explicit Runge-Kutta scheme to solve the hyperbolic PDE for the melt fraction ϕ. However, we start with the first order forward Euler method, because it will be an intermediate step in the higher-order time stepping scheme.

First-order forward Euler method
We first discretize the hyperbolic Eq. (15) for ϕ in time by using a forward Euler with time step k := t n − t n−1 > 0. The semi-discrete equation for ϕ n ≈ ϕ(t n ) reads Taking this equation in variational form and applying Q r elements results in the discrete system for ϕ h n ∈ Q h : with Solving Eq. (30) requires basically inversion of the mass matrix. Afterwards, the mean density can be updated:¯ Finally, the equation for (u h n , p h n ) ∈ X h is solved: with the right hand side F n [¯ n ] as given in Eq. (21). The entire algorithm for variable melt fraction rate looks now as follows: Algorithm (forward Euler) 1. Initialize u h 0 , p h 0 and ϕ h 0 , set n := 0 and t n := 0. 2. Increase n → n + 1 and set t n := t n−1 + k. 3. Make one forward Euler step (30) to determine ϕ h n . 4. Update mean density¯ n by (31). 5. Solve the linear problem (32) to determine p h n and u h n . 6. If t n ≤ T goto 2.
Note that we have a time step restriction (CFL condition) of the form k ≤ h for stability issues. An alternative without such a time step restriction would be the usage of an implicit time stepping scheme for ϕ h n . However, in this case, the corresponding equation includes u h n . On the other hand, the quantities α(t n ) in (12) and¯ (t n ) include ϕ h n , so that the equations for ϕ h n and (28) become mutually coupled. Solving this system is much more numerically expensive than the semi-explicit algorithm presented above.

Second-order forward Heun method
In order to obtain a second order method, me apply a second-order explicit Runge-Kutta method for ϕ, for instance the Heun method. This time stepping scheme consists of the forward Euler as a predictor step: leading to a predictor mean density¯ * ,n =¯ (ϕ h * ,n ), and a predictor velocity and pressure, u h * ,n and p h * ,n respectively, given as solution of Afterwards, the new melt fraction is obtained by solving followed by the new mean density¯ n =¯ (ϕ h n ), and new velocity and pressure by solving (32).
The numerical cost of the Heun method is just a factor of two compared to the forward Euler, but leads to a substantial increase in accuracy. This will be demonstrated in the numerical examples below. The corresponding algorithm looks as follows: Algorithm (Heun method) 1. Initialize u h 0 , p h 0 and ϕ h 0 , set n := 0 and t n := 0. 2. Increase n → n + 1 and set t n := t n−1 + k. 3. Make one forward Euler step (33) to determine ϕ h * ,n . 4. Update mean density¯ * ,n (similar to (31)). 5. Solve the linear problem (34) to determine p h * ,n and u h * ,n . 6. Update melt fraction ϕ h n by inversion of mass matrix, (35). 7, Get new mean density¯ n according to (31). 8. Update solid velocity and fluid pressure by solving (32). 9. If t n ≤ T goto 2.

Numerical example
In this section we test the methodology for two 2D examples, a stationary problem and a time-dependent problem. Both cases are designed in such a way that the exact solution is known.

Stationary example
To validate the analysis and the error estimates of the proposed scheme, we choose as an example the two-dimensional problem given in [23] with a known exact solution: The problem is solved on the unit square Ω := (0, 1) 2 in the xz− plane. The Darcy coefficient K D and the melt fraction ϕ depend only on the variable z : The shear and bulk viscosities η and ξ , respectively, are given by The resulting compaction pressure according to (18) is We Moreover, we are interested in the error behavior of the post-processed compaction pressure p c . In Fig. 2 we plot the corresponding discretization errors for Q 1 and Q 2 elements. We obtain for Q 1 -elements an error || p c − p h c || = O(h 3/2 ), and for Q 2elements an error || p c − p h c || = O(h 2 ). In comparison with the theoretical expectation, these results corresponds exactly with the theory for Q 2 , and show a super-convergence behaviour for Q 1 .
In order to get an idea about the error distribution, in Fig. 3 we depict the discretization errors on a sequence of meshes for Q 1 -elements. The errors are reduced uniformly under mesh refinement.

Time-dependent example
We extend the example in the previous section to the time-dependent case. The exact solution now reads: The coefficients η, ξ and K D remain independent of time and, hence, identical to their values in the stationary case (previous subsection). The quantities Γ , f and f 2 are adapted in such a way that the solution given above solves the set of equations.
Obviously, the solution is designed in such a way that we recover for t = 0 the same solution as in the stationary example: It is straightforward to verify that (10) and (11)  The initial forcing at t = 0, i.e. f (z, 0), is identical to f (z) of the stationary example.

Forward Euler
In order to validate the temporal error we first have a look at the error in melt fraction in different norms, see Fig. 4. In the L 2 -and the L ∞ -norm we observe first order convergence with respect to the time step, The H 1 -seminorm also begins to reduce with first order, but stagnates for smaller time steps. This is not unexpected, since Eq. (15) does not enforce H 1 -regularity nor H 1 -stability.
For the pressure and velocity variables the spatial error is much larger than the temporal error. The reason for this is probably that the time step does not enter directly into Eq. (22), but the temporal discretization enters only implicitly by the mean densitȳ which itself depends on the time dependent melt fraction ϕ. However, in order to visualize the temporal impact we plot the temporal error only by considering the quantities p h − p h n and u h − u h n at t = 0.1, where p h and u h are pressure and velocity, respectively, discretized in space but with the correct mean density¯ (t n ), i.e.
In Fig. 4 (right figure) we observe first oder convergence of the error in the L 2 -norms, || p h − p h n || ∼ k and ||u h − u h n || ∼ k for fixed spatial mesh size h. The compaction pressure, obtained by post-processing, also converges with first order: ||( p c ) h −( p c ) h n || ∼ k.

Heun method
For the Heun method on a fixed spatial mesh we obtain second order convergence for the melt fraction ||ϕ h k − ϕ|| L 2 (Ω) ∼ k 2 in L 2 -, H 1 -and L ∞ -norm (Fig. 5 left) until the (fixed) spatial error dominates and leads to stagnation of the total error. The reason that this stagnation is not observed with the forward Euler method is that with Heun, the error is a factor 100 smaller than with Euler (with k = 4 · 10 −3 ). This stagnation appears earlier on in the H 1 -norm. For solid velocity, fluid pressure and compaction pressure we obtain second order convergence as well (Fig. 5 right). Here, no stagnation appears which is due to the fact that we (once more) only depict the temporal error, e.g. ||u h − u h n || L 2 (Ω) , so that spatial effects are excluded.

Summary
We propose a variational formulation for modeling mantle-melt transport with a coercive bilinear form for solid velocity u s and fluid pressure p f . The compaction pressure p c is determined by a post-processing step if needed. We obtain existence and uniqueness of solutions, and we derived an a priori error estimate for equal-order finite elements. For the time-dependent case, we propose a splitting method which consists of an hyperbolic equation for the melt fraction ϕ and an elliptic problem for solid velocity and fluid pressure. The time-discretization for the melt fraction is carried out by explicit schemes (forward Euler or Heun method), so that the equation for velocity and pressure decouples from the equation for the melt fraction, i.e. only information of previous time steps of u s and p f enters into the equation for ϕ.
In numerical examples with known exact solutions we demonstrate that the expected convergence rates with respect to the spatial mesh size h and with respect to the time step k are obtained.