A face-centred finite volume method for high-contrast Stokes interface problems

A new face-centred finite volume method (FCFV) for Stokes problems involving sharp interfaces is proposed. Two formulations, based on two strong forms of the Stokes problem, and using different mixed variables, are presented. Particular attention is paid to the symmetry of the resulting system of global equations, and a simple rewriting of the interface boundary condition is proposed to ensure that one of the formulations preserves the symmetry of the linear system that is usually lost when considering material interfaces. Four numerical examples are considered to test the implementation numerically by performing mesh convergence studies, in two and three dimensions. The examples account for discontinuous viscosity as well as the effect of surface tension. The results show that one of the formulations is less sensitive to the numerical stabilisation used in FCFV methods but does not preserve the symmetry of the global system, whereas the other formulation is more sensitive to the stabilisation, but preserves the symmetry of the resulting system of equations. The FCFV method appears as a promising alternative for the simulation of viscous flow involving internal boundaries on conformal meshes. The potential application of the FCFV method for the purpose of geodynamic modelling is discussed.


INTRODUCTION
The face-centred finite volume (FCFV) method was introduced for elliptic problems, 1 including the Poisson and the Stokes problems. The method can be seen as a particular case of the hybridisable discontinuous Galerkin (HDG) method, originally developed by Cockburn and co-workers, 2-5 using the lowest order of approximation for all the variables. As such, the method inherits a number of desirable properties from HDG methods. In the context of Stokes flow problems one of the most attractive properties is the ability to use the same approximation space for the velocity and the pressure, whilst circumventing the so-called Ladyzhenskaya-Babuška-Brezzi (LBB) condition.
The application of the FCFV method to Stokes flows has been considered for the so-called isoviscous case, 1 that is, without spatial variations in the viscosity field. Results on the lowest order HDG method have been also reported, comparing three different formulations of the Stokes problem, but again only for the isoviscous case. 6 Work on HDG for Stokes flow involving material interfaces is substantially more reduced. For instance, Wang and Khoo 7 derived an HDG formulation for Stokes flows involving sharp discontinuities in the viscosity field. However, the work of Wang and Khoo 7 does not consider the lowest order approximation that leads to the FCFV scheme. The current work presents, for the first time, the derivation of the FCFV method for Stokes flows involving material interfaces. Two formulations, studied for the isoviscous case, 6 are developed and compared using numerical examples. Particular attention is paid to the effect of the interface condition on the symmetry of the global system of equations. This work shows that a simple rewriting of the interface condition ensures that the resulting system of equations for one of the FCFV formulations presented remains symmetric, even in the presence of material interfaces. The examples include test cases that involve viscosity jumps of up to 12 orders of magnitude, which are particularly relevant to geodynamic modelling, as simulations typically involve heterogeneous materials characterised by large spatial viscosity variations. The optimal convergence properties of the method are demonstrated by means of numerical experiments. In addition, examples in two and three dimensions are used to demonstrate the performance of the method when the interface jump condition accounts for the effects of surface tension. Finally, a numerical experiment is performed to study the influence of the stabilisation parameter on the accuracy of the considered FCFV methods.
The remaining of the article is structured as follows. In Section 2 two variants of the Stokes flow problem are considered and a simple rewriting of the material interface condition is proposed. Section 3 derives two weak FCFV formulations for the Stokes flow problem with material interfaces and the discrete systems are presented in Section 4. A set of numerical examples are presented in Section 5, including mesh convergence studies and experiments to illustrate the effect of the stabilisation parameter. Finally, Section 7 summarises the conclusions of the work that has been presented.

GOVERNING EQUATIONS
Let us consider an open bounded domain Ω ∈ R n sd , where n sd is the number of spatial dimensions. The boundary of the domain, Ω is partitioned into the disjoint Dirichlet and Neumann boundaries Γ D and Γ N respectively. In addition, a material interface Γ I is assumed to split the domain into two disjoint subdomains Ω 1 and Ω 2 , with different material properties. Two different formulations are considered in this work, namely a formulation based on the gradient of the velocity tensor and a formulation based on the symmetric gradient of the velocity tensor. 8 They are briefly recalled in this section.

Formulation in terms of the gradient tensor
The strong form of the Stokes equation in terms of the gradient tensor can be written as n( u − pI n sd ) = g on Γ N , where u is the velocity vector, p is the pressure, > 0 is the viscosity, s is the volumetric source term, f is the imposed velocity on the Dirichlet boundary, g is the imposed pseudo-traction on the Neumann boundary and h is the interface jump condition due to the effect of surface tension. The matrix E = [n T , t T 1 , … , t T n sd −1 ] contains the outward unit normal vector n and the tangential vectors t k , for k = 1, … , n sd − 1, such that {n, t 1 , … , t n sd −1 } form an orthonormal system of vectors. The interface jump ⟦ ⟦⋅⟧ ⟧ operator is defined as the difference between the values to the left and right of the interface, namely Ω L and Ω R , that is This formulation is sometimes referred to as the velocity-pressure formulation of the Stokes problem. 8

Formulation in terms of the symmetric gradient tensor
The strong form of the Stokes equation in terms of symmetric gradient tensor can be written as where g denotes now the imposed traction in the Neumann boundary. This formulation is sometimes referred to as the Cauchy stress formulation of the Stokes problem. 8 It is worth noting that two main differences are observed between the two formulations of the Stokes problem considered, namely the momentum equation and the Neumann boundary condition. The numerical examples considered in this work will only feature Dirichlet boundary conditions, and therefore the differences observed between the two formulations are only due to the difference in the momentum equation.

Interface condition
The interface boundary condition in Equations (1) and (3) imposes a net zero stress on the interface. 9 As it will be shown, the discretisation of the problem with the interface boundary condition written as in Equations (1) and (3) leads to a nonsymmetric linear system of equations for both formulations of the Stokes problem. Here a simple rewriting, not usually considered in the literature, 9,10 is proposed to ensure that, after the spatial discretisation is performed, the resulting linear system of equations, corresponding to the Stokes formulation in terms of the symmetric gradient tensor, is symmetric. As the matrix E is orthonormal, the interface condition can be rewritten as whereh ∶= hE T L , and the jump operator ⟦⋅⟧ is defined as the sum between the values to the left and right of the interface, namely Ω L and Ω R , that is

FCFV WEAK FORMULATIONS
This section details the FCFV weak forms for the Stokes formulations in terms of the gradient tensor and the symmetric gradient tensor. In the remaining of the article, the two FCFV methods are referred to as the gradient-based FCFV and the symmetric gradient-based FCFV, respectively.

Mixed formulation
The domain is partitioned into n el disjoint subdomains and the strong form is written in the broken domain after splitting the momentum equation in two first-order equations, leading to the mixed strong form in Ω e , and for e = 1, … , n el , in Ω e , and for e = 1, … , n el , in Ω e , and for e = 1, … , n el , The last two equations are introduced to enforce the continuity of the velocity and the normal flux across Γ ⧵ Γ I , where Γ denotes the so-called mesh skeleton, that is, the set of internal faces (edges in two dimensions).

Local and global problems
Following the rationale of previous work on FCFV 1,[11][12][13][14] and HDG, 6,15-17 the problem is split into a set of n el local problems with Dirichlet boundary conditions and a global problem to compute the hybrid velocity defined on the mesh skeleton and the mean value of the pressure in each element.
In the local problems, for e = 1, … , n el , the solution (L e , u e , p e ) is written in terms of the unknown hybrid variablê u, namely As only Dirichlet boundary conditions are considered in the local problems, an extra condition is required to ensure that the pressure is completely determined. In this work the extra condition is used, where e denotes the mean pressure on the boundary of element Ω e and ⟨p, q⟩ Ω e ∶= ∫ Ω e pq dΓ. The global problem is introduced to compute the velocity on the mesh skeleton,û, and the mean value of the pressure in each element, e , namely It is worth noting that the first continuity condition in Equation (9) is satisfied automatically because the hybrid variable, u, is unique on each face and the Dirichlet boundary condition u =û is imposed in the local problems. In addition, the free divergence condition in Equation (7) induces the compatibility condition 3.1.3 Weak form of the local and global problems The discrete weak formulation of the local problems reads: for e = 1, … , n el , find In the above expressions,  h (Ω e ) denote the space of constant functions in each element and ) the weak form of the local problems is: In a similar fashion, the following global problem accounting for the transmission conditions, the Neumann boundary condition and the interface condition, given in Equation (9) for allŵ ∈ [ h (Γ ∪ Γ N )] n sd . In the above expressions, h (S) denote the space of constant functions in S ⊂ Γ ∪ Ω and denote the inner  2 (Γ i ) products for scalar and vector valued functions respectively. To simplify the notation in the next sections, the super-index h , denoting discrete approximations, is omitted.

Symmetric gradient-based FCFV weak formulation
The process to obtain the weak form of the symmetric gradient-based formulation is very similar to the one described above for the gradient-based FCFV formulation. First, the mixed form for the symmetric gradient-based formulation is written in the broken domain as in Ω e , and for e = 1, … , n el , in Ω e , and for e = 1, … , n el , in Ω e , and for e = 1, … , n el , Without detailing all the intermediate steps, the resulting weak form of the local problems is: for e = 1, … , n el , find In a similar fashion, the following global problem accounting for the transmission conditions, the Neumann boundary condition and the interface condition is:

FCFV DISCRETISATION
To simplify the presentation, let us introduce the following notation for different sets of faces:  e ∶= {1, … , n e fa } is the set of indices for all the faces of element Ω e ;  e ∶= {j ∈  e | Γ e,j ∩ Γ D ≠ ∅} is the set of indices for all the faces of element Ω e on the Dirichlet boundary Γ D ;  e ∶= {j ∈  e | Γ e,j ∩ Γ N ≠ ∅} is the set of indices for all the faces of element Ω e on the Neumann boundary Γ N ;  e ∶=  e ⧵  e = {j ∈  e | Γ e,j ∩ Γ D = ∅} is the set of indices for all the faces of element Ω e not on the Dirichlet boundary Γ D ;  e ∶= {j ∈  e | Γ e,j ∩ Γ I ≠ ∅} is the set of indices for all the faces of element Ω e on the interface Γ I .

Gradient-based formulation
Assuming a constant degree of approximation for all the variables in Equation (14), that is for L e , u e and p e in each element and forû in each face (edge in two dimensions), the discrete local problems are for e = 1, … , n el . It is worth noting that Equation (20c) becomes redundant when using a constant approximation. First, it does not contain any of the unknowns of the local problem and, second, it coincides with the discrete version of Equation (15b).
The local problem leads to three uncoupled equations that express L e , u e and p e as functions of the global unknownŝ u and e , namely where the following auxiliary quantities can be precomputed as they only depend on the data, the stabilisation parameter and the mesh Similarly the global problem of Equation (15) becomes where  e is the indicator function of the set  e , that is, By inserting the expressions of Equation (21) into Equation (23), the global problem, only in terms of the global unknownsû and , is obtained, namely The vector û contains the value of the hybrid variable on Γ ∪ Γ N and the vector contains the values of the mean pressure on each element Ω e . The blocks composing the matrix and vector of the global linear system are obtained by assembling the elemental contributions given by With the gradient-based formulation, the occurrence of an interface with discontinuous viscosity requires specific modifications of both the left and right-hand sides of the discretised momentum equations (Equation (26a) and (26c)). It is worth emphasising that for a pure Dirichlet problem, an extra condition must be imposed in order to remove the indeterminacy of the pressure. Following the work of Cockburn et al., 6,16 the usual strategy of imposing mean zero pressure on the whole domain, that is, is considered here.

Symmetric gradient-based formulation
Following the same procedure described above for the gradient-based formulation, the local problem of the symmetric gradient-based formulation leads to three uncoupled equations that express L e , u e and p e as functions of the global unknownsû and e , namely where the following auxiliary quantities can be precomputed as they only depend on the data, the stabilisation parameter and the mesh Similarly, the global problem leads to a system of equations where the vector û contains the value of the hybrid variable on Γ ∪ Γ N and the vector contains the values of the mean pressure on each element Ω e . The blocks composing the matrix and vector of the of the global linear system are obtained by assembling the elemental contributions given by A notable difference with the gradient-based formulation is the occurrence of interface contributions only on the right-hand side of the discretised momentum equations, as shown in Equation (31c).

Solution algorithm and symmetry of the global system
For both FCFV approaches, the solution algorithm can be summarised in three stages. First, the global system is built by assembling the elemental contributions given by Equations (26) and (31), respectively. Second, the global system is solved to obtain the hybrid velocity on each face, û, and the mean pressure on each element, . Finally, the local problems of Equations (21) and (28) are used to obtain the velocity, the velocity gradient or symmetric velocity gradient and the pressure, L e , u e and p e on each element Ω e . Both FCFV approaches were implemented into an open source Julia package (https://github.com/tduretz/FCFV_NME23), which provides useful details on the practical implementation of the considered FCFV schemes. The resulting global system of equations for the gradient-based FCFV is nonsymmetric, due to the term containing the exterior product n j ⊗ n i , in Equation (26a), only for faces that belong to the material interface. The reason for the loss of symmetry is due to the original formulation of the Stokes problem in Equation (1) that features the gradient tensor on the momentum equation but the Cauchy stress tensor on the interface condition.
In contrast, the symmetric gradient-based formulation leads to a symmetric system of linear equations due to two main reasons. First, the momentum equation and the interface condition in Equation (3) both feature the Cauchy stress tensor. Second, the simple rewriting of the interface condition, proposed in Section 2.3, ensures that the contribution to the left-hand side of the system of equations for faces on the interface is the same as for all the other faces of the mesh skeleton.
Without this simple rewriting, the resulting global system of equations for the symmetric gradient-based FCFV formulation is where the blocks composing the matrix and vector are obtained by assembling the elemental contributions given by It is worth noting that without the rewriting of the interface condition, not only is the diagonal block of the global systemK ûû nonsymmetric but also, the off diagonal blocksK û andK û do not satisfy thatK û =K T û . To further analyse the difference between the gradient-based and symmetric gradient-based FCFV formulations, let us consider the sparsity pattern of the global matrices in Equations (25) and (30), respectively. For the gradient-based formulation, the number of nonzero entries in the blockK ûû is where n Γ fa is the number of interior faces (i.e., not on the boundary), n I fa is the number of faces on a material interface and n fa is the number of boundary faces. Similarly, for the symmetric gradient-based formulation, the number of nonzero entries in the blockK ûû is n sg nz (K ûû ) = (2n e fa − 1)n 2 sd n Γ fa + n e fa n 2 sd n fa .
To better compare both formulations, let us assume that the number of boundary faces and the number of faces on a material interface are negligible, when compared to the number of interior faces. In addition, let us consider the case of simplicial meshes, where the number of element faces is n e fa = n sd + 1. Under these assumptions, the number of nonzero entries in the blockK ûû is approximately and for the gradient-based and symmetric-gradient based formulations, respectively. It is worth noting that for the symmetric gradient formulation, only half of the entries are required due to the symmetry of the blockK ûû , hence the two in the denominator in Equation (37). The ratio between the nonzero entries of the blockK ûû corresponding to the symmetric gradient and gradient formulations is approximately n sd ∕2. Therefore, in two dimensions, both formulations lead to a blockK ûû with approximately the same number of nonzero entries, whereas in three dimensions, the blockK ûû of the symmetric formulation contains 1.5 more nonzero entries. Despite the extra nonzero entries in the blockK ûû , the symmetric gradient formulation only requires to store one of the off-diagonal blocks of the global matrix and more efficient solution algorithms can be applied.
The numerical examples in the next section will also compare the accuracy provided by both schemes for the same spatial discretisation.

NUMERICAL EXAMPLES
Two examples usually employed by the geodynamic community are presented here, namely the SolCx and SolVi benchmarks. [18][19][20][21][22] Both tests provide an exact full flow-field solution of 2D Stokes flow in the presence of a discontinuous viscosity field. In the case of SolCx, the flow is driven by buoyancy and the viscosity jump occurs along a linear segment. In the SolVi case, the flow is driven by far field kinematics and the viscosity jump occurs along a circular or elliptical interface. The examples are used to verify the implementation of the proposed FCFV scheme and to compare the gradient-based and symmetric gradient-based formulations of the FCFV methods described above. Two further numerical examples in two and three dimensions are used to show the ability of the approach to handle problems that involve a nonvanishing surface tension. Unless otherwise stated, the stabilisation parameter is selected as wherex e denotes the barycentre of the element Ω e and is a numerical parameter selected as = 10 in two dimensions and = 1 in three dimensions. Numerical experiments are included to demonstrate the suitability of this choice for a variety of problems and for different viscosity contrasts. The accuracy of the numerical solution is assessed by measuring the relative error of the velocity, pressure, gradient tensor in the  2 (Ω) norm and the hybrid velocity in the  2 (Γ) norm, namely where the superscript h denotes the numerical solution, the superscript ⋆ denotes the analytical solution and the  2 (Ω) norm for scalar, vector and tensor valued functions is the norm induced by the scalar products in Equation (12).
To further assess the implementation of the methods proposed, the rate of convergence for the error of all the variables is also included in the numerical examples. The rate of convergence of the error for the velocity u is defined as where h i denotes the characteristic mesh size of i-th mesh employed in the convergence study. In a similar fashion the rate of convergence for the errors of the pressure, the velocity gradient and the hybrid velocity are denoted by r p , r L , and r û .
To compute the domain integrals that are involved in the evaluation of u , p , and L , a quadrature with six integration points in each triangle is employed. Similarly, to compute the integrals on the mesh skeleton that are involved in the evaluation of û , a quadrature with three integration points is used. This is done to ensure that the variation of the analytical solutions is properly captured.
The presented results were obtained using an in-house code based on the open-source HDGlab, 23 which is a MATLAB serial implementation of the HDG method. An open-source Julia implementation of the FCFV method applied to the viscous inclusion benchmark is provided online (https://github.com/tduretz/FCFV_NME23).

5.1
The SolCx test: Discontinuous viscous profile in a square The first example considers the domain Ω = [0, 1] 2 and the viscosity given by The flow is driven by a source term s(x) = (x)g, with (x) = cos( x 1 ) sin( x 2 ) and g = (0, −1) T . Dirichlet boundary conditions, corresponding to the analytical solution, derived by Zhong,18 are imposed in whole boundary of Ω.
Four triangular meshes, shown in Figure 1, are considered to perform a convergence study. The meshes are generated by splitting into four triangles a Cartesian grid of quadrilateral elements. The number of elements in the four meshes is 1024, 4,096, 16,384, and 65,536 respectively, leading to a global system of equations of dimension 4032, 16,256, 65,280, and 261,632, respectively. It is worth noting that further numerical examples, not reported here for brevity, demonstrated that the mesh is insensitive to mesh distortion, so the errors vary only marginally when nonregular meshes are used. This confirms the conclusions reported in previous works on FCFV. 1,[11][12][13][14] The numerical solution computed with the symmetric gradient-based FCFV formulation is shown in Figures 2 and  3 for the first and last mesh, respectively. The results, corresponding to a discontinuous viscosity field with 1 = 1 and 2 = 10 6 , clearly show the increased resolution of the complex features, induced by the large viscosity jump, as the mesh is refined. It can be clearly observed how the strong jump in the viscosity field induces a disparate behaviour of the velocity and pressure fields in the two regions. The discontinuous character of the FCFV approximation and the alignment of the mesh with the interface results in a sharp resolution of the pressure discontinuity, with no artefacts or oscillations displayed in the numerical solution.
The results of the mesh convergence study are displayed in Table 1. The error and rate of convergence for the velocity, pressure, gradient tensor and hybrid velocity are reported for both the gradient-based FCFV and the symmetric gradient-based FCFV formulations. The results show the optimal, (h), rate of convergence for both FCFV formulations. The most accurate results are obtained for the pressure field, followed by the hybrid velocity, the gradient tensor and the velocity.  TA B L E 2 SolCx: Wall clock (in seconds) for the solution of the linear system. CPU g and CPU sg denote the wall clock for the gradient-based and symmetric gradient-based formulations, respectively. The ratio is also reported.  When comparing the two FCFV formulations, it is observed that they produce almost identical results for the velocity and hybrid velocity. However, the symmetric gradient-based formulation provides a lower error for the pressure and gradient tensor. To provide a comparison of the computational cost of the two formulations, Table 2 shows the wall clock, in seconds, required to solve the global system of equations for the two FCFV methods. Only the time required to solve the linear system is reported as this is the dominant cost, especially for fine meshes. For the gradient-based formulation the multifrontal LU factorisation algorithm in UMFPACK 24 is applied for solving the global system, whereas for the symmetric-gradient formulation the multi-frontal MA57 25 algorithm is employed. The computations were performed in an Intel ® Core CPU @ 3.70GHz and 32GB main memory available.

Gradient-based FCFV
The results clearly show the superior performance of the symmetric gradient-based formulation. For a given mesh, not only the time required to solve the linear system is lower, but also the accuracy obtained is higher. The superior performance is significantly more evident as the mesh is refined.
To further analyse the behaviour of the FCFV formulations considered, Table 3 shows the error on the last mesh and the rate of convergence for an increasing value of the viscosity contrast. In all cases, the optimal rate of convergence is exhibited, but the value of the viscosity contrast influences the accuracy of some of the computed quantities. For both FCFV methods, the errors on the velocity and gradient tensor show the greatest dependence with respect to the viscosity contrast, with an error approximately three times higher when the viscosity contrast exceeds 10 2 , compared to the isoviscous case. It is worth noting that the error in all the variables with 2 = 10 3 is almost identical to the error in velocity with = 10 6 , showing the robustness of the proposed schemes for high contrast problems. The error of the hybrid velocity and pressure are less influenced by the values of the viscosity contrast. For the gradient-based FCFV, the error on the pressure for high values of the viscosity contrast is less than two times higher than the error of the isoviscous case, whereas for the symmetric gradient-based approach a marginally lower error is actually observed for high contrast cases. The error on the hybrid velocity for high values of the viscosity contrast is approximately two times higher than the error of the isoviscous case for both the gradient-based and the symmetric gradient-based FCFVs. Overall, the error of the symmetric gradient-based formulation seems to be always lower or the same provided by the gradient-based formulation. When combined with the fact that the symmetric gradient-based approach requires the solution of a symmetric linear system of equations, indicates that the symmetric gradient-based FCFV is the most efficient alternative.
To conclude with this example, the effect of the stabilisation parameter is studied. Figure 4 shows the error of velocity, pressure, gradient tensor and hybrid velocity as a function of the stabilisation parameter, which is selected independently in the two subdomains Ω 1 and Ω 2 , in the range [10 −2 , 10 2 ]. This study considers the case with equal viscosity 1 = 2 = 1 and the gradient-based FCFV. The results show that low values of the stabilisation parameter destroy the accuracy of the velocity field, whereas a high value destroys the accuracy of the other variables. The results indicate that choosing = 10 provides the most accurate results for all the variables. Further simulations, not reported here, for the symmetric gradient-based formulation are almost identical.
The case with a high contrast is considered in Figure 5, corresponding to 1 = 1 and 2 = 10 6 and by using the gradient-based FCFV. The behaviour in this example is completely different, when compared to the isoviscous case. The most accurate results correspond to a value of = 10 in Ω 1 and are almost independent of the stabilisation parameter used in the second subdomain, except for the velocity field. The velocity error shows a higher dependence on the selected value of the stabilisation and very low values of the stabilisation also induce a higher error. The optimal value of the stabilisation in the first subdomain is again = 10 whereas in the second subdomain a value of between 10 and 10 8 does not produce a visible change in the accuracy obtained. Once more, further experiments, not reported here, for the symmetric gradient-based formulation are almost identical.
This study also confirms that the choice adopted here for the stabilisation parameter, given in Equation (38), provides the most accurate results for both FCFV formulations. This factor 10, appearing in Equation (38), was also identified in previous FCFV studies 1 as an optimal value of the stabilisation for isoviscous cases.

The SolVi test: Circular inclusion
The second example considers the domain Ω = [−3, 3] 2 with a circular inclusion of radius 1. The viscosity is given by Dirichlet boundary conditions, corresponding to the analytical solution, derived by Schmid and Podladchikov, 19 are imposed in whole boundary of Ω. Four triangular meshes generated using the advancing front method, shown in Figure 6, are considered to perform a convergence study. The number of elements in the four meshes is 2004, 8130, 32,840, and 132,022 respectively, leading to a global system of equations of dimension 7896, 32,280, 130,880, and 527,128, respectively.
Two cases with different viscosity contrast are considered, namely 1 = 1 and 2 = 10 6 and 1 = 10 6 and 2 = 10 −6 . The numerical solution with 1 = 1 and 2 = 10 6 computed with the symmetric gradient-based formulation on the finest mesh is shown in Figure 7. This case exhibits a moderate jump in the pressure, with a maximum absolute value jump of the pressure of approximately 4 outside the inclusion and a constant value of the pressure equals to zero inside the circular inclusion.
The numerical solution with 1 = 10 6 and 2 = 10 −6 computed with the symmetric gradient-based formulation on the finest mesh is shown in Figure 8. This case exhibits a much larger jump in the pressure, six order of magnitude higher than in the previous case, due to the much larger jump in the viscosity across the interface, of 12 orders of magnitude in this example.
The results of the mesh convergence study are displayed in Table 4. As in the previous test case, optimal convergence is observed for all the variables, for both high contrast values considered and for both FCFV formulations. The accuracy of both FCFV formulations is comparable, with minor differences observed on the pressure and gradient tensor for the case with lower viscosity contrast and on all the variables for the case with higher viscosity contrast. The results also show that for higher viscosity contrasts the accuracy of the velocity decreases, whereas the FCFV methods considered offer higher accuracy for the pressure field. This was also observed in the previous example.
To corroborate the conclusions of the previous example, a comparison of the computational cost of the two formulations is also performed. Table 5 shows the wall clock, in seconds, required to solve the global system of equations for the two FCFV methods. Similar to the previous example, the symmetric gradient-based formulation clearly shows the superiority when compared to the gradient-based formulation. For the same number of degrees of freedom the accuracy is slightly better and the wall clock is more than two orders of magnitude lower.
As in the previous example, a further experiment is performed to study the effect of the viscosity contrast on the accuracy and convergence rate of the FCFV formulations considered here. Table 6 shows the error on the last mesh and the rate of convergence for an increasing value of the viscosity contrast. The conclusions are slightly different when compared to the previous example. Here, an increase in the viscosity contrast results in a less accurate computation of the pressure and the gradient tensor, whereas the velocity error seems to be independent on the viscosity contrast. It is worth noting that optimal rates of convergence are exhibited by both FCFV formulations for all the viscosity contrasts considered.
Further computations, not reported here for brevity, show that the errors reported in Table 4 are almost identical if an exact description of the interface is adopted, for instance using the NEFEM rationale. 26 Given that the error is expected to converge with (h), the geometric error is lower than the spatial discretisation error. This shows that the different behaviour observed in both examples is not due to the curvature of the interface.
As in the previous example, the effect of the stabilisation parameter is studied. Figure 9 shows the error of velocity, pressure, gradient tensor and hybrid velocity as a function of the stabilisation parameter, which is selected independently in the two subdomains Ω 1 and Ω 2 . This study considers the case with 1 = 1 and 2 = 10 6 and the gradient-based FCFV. The results are qualitatively similar to the ones reported in the SolCx test with the same jump in viscosity. However, there are some important differences. The errors on all variables show that a low value of the stabilisation parameter in the first subdomain offers the most accurate results, whereas almost no influence is observed by the selected value of the stabilisation in the second subdomain. A low value of the stabilisation parameter does not destroy the accuracy of the velocity field as it was observed in the previous test case. It is worth noting that, as in the previous example, the pressure is the variable that shows the larger variation of the error in terms of the stabilisation parameter.
The results for the same study applied to the symmetric gradient-based FCFV formulation are displayed in Figure 10. Similar to the previous example, the error seems to be not influenced by the choice of the stabilisation in the second subdomain. However, in this case a low value of the stabilisation parameter in the first subdomain seems to be not the optimal choice. The results indicate that a value of = 10 in the first subdomain is required to obtain the maximum TA B L E 4 SolVi: Error of the velocity, pressure, gradient tensor and hybrid velocity and rate of convergence for the two FCFV formulations and for different viscosity contrasts. TA B L E 5 SolVi: Wall clock (in seconds) for the solution of the linear system. CPU g and CPU sg denote the wall clock for the gradient-based and symmetric gradient-based formulations, respectively. The ratio is also reported. accuracy. Despite that the results suggest that a value of = 10 in the first subdomain might lead to better results, the error on the pressure is again highly influenced by an increase in the stabilisation, whereas the other variables only show a minor improvement. Finally the case with higher viscosity contrast is considered. Figure 11 shows the error of velocity, pressure, gradient tensor and hybrid velocity as a function of the stabilisation parameter, for the gradient-based FCFV. The results show that the gradient tensor is the variable most affected by the choice of the stabilisation parameter, specially the stabilisation used in the second subdomain. The other variables show that a value of the stabilisation below 10 7 provides little differences in the results. TA B L E 6 SolVi: Error of the velocity, pressure, gradient tensor and hybrid velocity and rate of convergence, computed on the finest mesh, for the two FCFV formulations with 1 = 1 and an increasing value of 2 . The results for the same study applied to the symmetric gradient-based FCFV formulation are displayed in Figure 12. The results confirm the higher dependency of the accuracy on the stabilisation parameter, previously observed for this FCFV formulation. A value of the stabilisation parameter equal to 10 7 in the first subdomain is crucial to obtain accurate results, whereas the results are less influenced by the choice of the stabilisation in the second subdomain.

Gradient-based FCFV
It is worth noting that the choice recommended in this work, given in Equation (38) provides the maximum accuracy and the optimal rate of convergence in all the examples studied and for all the different viscosity contrasts studied.

Manufactured 2D solution with nontrivial interface jump
The last example, taken from Wang and Khoo, 7 considers the domain Ω = [0, 2] × [−0.5, 1, 5] and the viscosity given by The analytical solution is given by where = 1 2 − √ 1 4 2 + 4 2 . As in previous examples, the pressure and the gradient tensor exhibit a discontinuity across the interface, but contrary to the previous examples, the right hand side of the interface condition is nonzero, that is, h ≠ 0. Therefore, this test is selected to show that ability of the proposed FCFV scheme to deal with problems in which surface tension is taken into account.
Dirichlet boundary conditions are imposed in the whole boundary and the meshes considered for this example are the same used in the SolCx test, but scaled and translated to the domain of interest, here [0, 2] × [−0.5, 1, 5], instead of [0, 1] 2 .
The numerical solution with 1 = 1 and 2 = 10 −4 , computed with the symmetric gradient-based FCFV formulation on the finest mesh is shown in Figure 13. The different variation of the pressure field in the horizontal direction, due to the different value of the viscosity, is clearly observed. For = 1, the pressure presents a sharp gradient near the left part of the boundary, whereas for = 10 −4 the variation is almost not noticeable. For x 1 = 0 the pressure is continuous across the interface, whereas the jump in pressure across the interface grows as x 1 grows. For the velocity, the higher variations are observed in the region with lower viscosity, whereas in the region with = 1, the velocity is almost constant. The results also show that the maximum error in the second component of the velocity is observed in the region where the horizontal velocity presents a sizeable variation. It is worth remarking that the numerical approximation of the second component of the velocity is purely numerical error as the analytical value is zero in this example. The results of the mesh convergence study are displayed in Table 7. Both formulations, the gradient-based and the symmetric gradient-based FCFV, show an optimal rate of convergence in all the variables. In addition, it can be observed that both methods provide almost identical accuracy for the velocity and pressure whereas the stress based formulation seems to offer some marginal extra accuracy on the gradient tensor. As mentioned earlier, when combined with the fact that the symmetric gradient-based approach leads to a symmetric linear system of equations, indicates that the symmetric gradient-based FCFV is the most efficient alternative.

Manufactured 3D solution with nontrivial interface jump
The last example, taken from Kirchhart et al., 27 considers a three dimensional problem in [−1, 1] 3 and the viscosity given by where r = ||x|| 2 and the interface is located at r = r I with r I = 2∕3. The analytical solution is given by where Dirichlet boundary conditions, corresponding to the analytical solution, are imposed in whole boundary of Ω. As in the previous example, not only do the pressure and gradient tensor exhibit a discontinuity across the interface, but also h ≠ 0.
Four tetrahedral meshes, shown in Figure 14, are considered to perform a convergence study. The number of elements in the four meshes is 5025, 40,660, 326,589, and 2,626,862 respectively, leading to a global system of equations of dimension 34,005, 279,598, 2,265,387, and 18,304,694, respectively.
The numerical solution with 1 = 1 and 2 = 10 2 , computed with the gradient-based FCFV formulation on the finest mesh is shown in Figure 15. The continuity of the velocity field is clearly observed, with the first and second component TA B L E 8 Manufactured 3D solution: Error of the velocity, pressure, gradient tensor and hybrid velocity and rate of convergence for the symmetric gradient-based FCFV formulation with 1 = 1 and 2 = 10 2 . exhibiting the main variation in the vertical and horizontal directions respectively. The third component mainly shows the discretisation error, as the exact solution is exactly zero. The norm of the velocity field, shows the dependency on the radial position. Finally, the pressure shows the dependency on the horizontal direction only, which is different in the two subdomains. To better appreciate the discontinuity of the pressure on the interface and the variation of the pressure in the two subdomains, a nonlinear colour scale is used for the pressure plot. Inside the sphere, the pressure varies between 9.7 and 10.3, whereas outside, the pressure varies between −1 and 1. It is worth noting that the values used to display the solution on the interface are taken from the adjacent elements in the first subdomain, that is, inside the sphere. The results of the mesh convergence study are displayed in Table 8. The optimal rate of convergence is observed in all variables, confirming the correct implementation of the FCFV for three dimensional problems with nontrivial interface jumps that account for surface tension.

PERSPECTIVES
The presented results demonstrate that the FCFV method may readily be used to simulate viscous flow with internal viscosity discontinuities and surface tension. The occurrence of discontinuous viscosity interfaces is a common feature in geodynamic models which involves lithological heterogeneities (e.g., layers, inclusions) characterised by large and sharp effective viscosity jumps. Additionally, the effect of interfacial tension can be important for the simulation of flows involving multiple phases (solid, liquid, gas) or microstructural processes (e.g., grain boundary migration). Geodynamic modelling tools are often based on Cartesian representations which fail at precisely capturing interfaces 21 and, more particularly, viscosity jumps. The resolution of viscosity discontinuities has been addressed by using either immersed interface or ghost-fluid methods 10 or boundary-fitted finite-element (FE) method. 28 Nevertheless, due to LBB requirements, robust FE solutions generally require high-order elements (e.g., Crouzeix-Raviart elements 22,28 ). By providing robust Stokes flow solution on conformal meshes while circumventing LBB requirements, the FCFV method could represent an attractive low-order alternative to typical FE. By preserving the symmetry of the velocity block, the symmetric gradient-based FCFV formulation render possible the use of efficient solving strategies based on Schur-complement reduction 29 and Cholesky factorisations. 28 After carefully addressing the issue of discontinuous viscosity, the next development step of the FCFV will be to address other prototypical geodynamic problems which should account for smooth and large spatial variations of viscosity, 29 free surface flow, 30 material nonlinearities and compressibility. 31

CONCLUDING REMARKS
A new FCFV method for the solution of Stokes problems with sharp interfaces have been presented. Two FCFV formulations, based on two different strong forms of the Stokes problems and using different mixed variables, were presented. The first formulation is based on the Stokes problem written in terms of the gradient tensor, whereas the second formulation is based on the Stokes problem written in terms of the symmetric gradient tensor. In addition, a simple rewriting of the interface condition was used to ensure that the symmetric gradient-based FCFV formulation led to a global symmetric system of equations. The proposed approach, inherits a number of attractive properties from the hybridisable discontinuous Galerkin method, including the ability to use the same approximation space for the velocity and the pressure, whilst circumventing the LBB condition.
Four numerical examples were presented, including two and three dimensional examples, problems of interest to the geodynamics community, problems with nontrivial interface jumps accounting for surface tension and viscosity contrasts of up to 12 orders of magnitude. In all the examples considered, the two FCFV formulations showed the optimal rate of convergence under mesh refinement. A comparison of both formulation in terms of accuracy was also presented, showing that the accuracy of both methods was similar, with the error in pressure and gradient tensor being smaller for the gradient-based formulation. Together with the fact that the symmetric gradient-based formulation leads to a symmetric system of equations, results in this formulation being the preferred choice in this work. Numerical experiments were also presented to show the suitability of the recommended stabilisation parameter for the FCFV methods.

ACKNOWLEDGMENTS
The first author acknowledges the support of the Engineering and Physical Sciences Research Council (Grant Number: EP/T009071/1).

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in FCFV_NME23 at https://github.com/tduretz/ FCFV_NME23.