Face‐centred finite volume methods for Stokes flows with variable viscosity

Six face‐centred finite volume formulations are derived and compared for the simulation of Stokes flows with spatially varying viscosity. The main difference between the methods derived is the mixed variable used in the mixed formulation and the use of a weak or strong form in each element using integration by parts. A brief discussion about the properties of the different methods is provided, including comments on the computational cost and the symmetry of the resulting global system of equations. Finally, numerical examples in two and three dimensions are used to compare the accuracy of all the formulations presented. The examples include a problem where the methods are employed to simulate a steep variation of the viscosity, showing the ability to perform these simulations without using a mesh conforming to a material interface. The performance of different element types and different choices of the stabilisation is also discussed.

][30] The face-centred finite volume (FCFV) method, initially proposed in Reference 31, can be seen as a particular case of the hybridisable DG method with constant degree of approximation for all the variables.3][34][35][36] In this context, the method allows to build a first-order accurate approximation of the velocity and the pressure fields as well as the gradient of the velocity, without the need of using a reconstruction of the gradient.As a result, the method has shown that the accuracy is not sensitive to mesh distortion, contrary to other finite volume techniques.Recently, the authors extended the original FCFV formulation for Stokes problems to handle sharp material interfaces with viscosity contrasts of up to 12 orders of magnitude. 37his work presents a family of FCFV methods for the simulation of Stokes flows with spatially varying viscosity.Six formulations are considered, differing on the choice of the mixed variable and the use of a weak or strong form in each element by performing integration by parts.The six formulations are discussed and it is found that some of the formulations lead to an identical global system of equations and the only difference between them is the local problem to compute the mixed variable.In addition, the different choices of the mixed variable are found to lead to important implications in terms of the solution of the global problem.Some formulations require the computation of integrals of the spatially varying viscosity field over elements and over element faces (edges in two dimensions), whereas other formulations only require computing integrals of the viscosity field over the elements.Furthermore, some formulations are found to lead to a global system of equations that is not symmetric, whereas other formulations are able to maintain the symmetry of the global system.The formulations are then tested numerically using problems in two and three dimensions with rapidly varying viscosity fields.The examples also demonstrate the possibility to compute accurate solutions of problems involving viscosity fields that present steep gradients.The numerical studies presented include mesh convergence analysis and a discussion about the influence of the numerical integration error, induced by the computation of integrals involving the viscosity field, on the accuracy of the simulations.
The remainder of the paper is organised as follows.Section 2 briefly recalls the governing equations for a Stokes flow.In Section 3 the weak formulation of four FCFV methods is presented and the discretisation is detailed in Section 4. Section 5 discusses the differences between the formulations presented and two variations of the so-called non-scaled FCFV methods are described.Numerical examples are presented in Section 6.The two-dimensional examples include mesh convergence studies, a comparison of the accuracy of the formulations presented and a numerical study showing the influence of the numerical integration error.The conclusions of these studies are used to select the one of the formulations that is then applied to a three dimensional test case to show the optimal approximation properties of the method.Section 8 summarises the conclusions of the work that has been presented.Appendix A presents further numerical results using quadrilateral elements and briefly compares the performance of triangular and quadrilateral elements.Finally, Appendix B presents some further numerical results that illustrate that the selected stabilisation parameter for the HDG formulation does not have a noticeable influence on the accuracy of the results.

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.The strong form of the formulation of the Stokes equation can be written as where u(x) is the velocity vector, p(x) is the pressure, (x) is the viscosity function, s(x) is the volumetric source, f(x) the imposed velocity on the Dirichlet boundary, g(x) the imposed traction in the Neumann boundary and n is the outward unit normal vector to the boundary.This formulation is sometimes referred to as the Cauchy stress formulation of the Stokes problem. 38t is worth noting that the so-called velocity-pressure formulation 38 is only equivalent to the Cauchy stress formulation when the viscosity is piecewise-constant.Therefore, contrary to previous work by the authors in Reference 37, here only the Cauchy stress formulation is considered.

FCFV WEAK FORMULATIONS
This section details six different FCFV weak forms for the Stokes problem.The formulations differ in the choice of the mixed variable and the use of a weak or strong form in each element.To simplify the presentation and due to the similarity between the different formulations, all the details are provided for the derivation of the first weak formulation whereas the other formulations omit unnecessary details.

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 The last two equations are introduced to enforce the continuity of the velocity and the normal flux across Γ, 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 31,32,34,33,35 and HDG, 19,21,20,39 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 variable û, 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 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 ( 5) is satisfied automatically because the hybrid variable, û, 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 (3) induces the compatibility condition

Weak form of the local and global problems
The discrete weak formulation of the local problems reads: for e = 1, … n el , find 1 for all (G, w, q) ∈ [ h (Ω e )] n sd ×n sd × [ h (Ω e )] n sd ×  h (Ω e ).In the above expressions,  h (Ω e ) denote the space of constant functions in each element and (p, q) Ω e ∶= ∫ Ω e pq dΩ, (p, q) Ω e ∶= ∫ Ω e p ⋅ q dΩ, (P, denote the inner  2 (Ω e ) products for scalar, vector and tensor valued functions, respectively.After integrating by parts Equation (7b) again the strong form in each element is recovered, and by introducing the trace of the numerical normal flux the weak form of the local problems is: ) Similarly, the following global problem accounting for the transmission conditions and the Neumann boundary condition, given in Equation ( 5) is: find 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.

Mixed formulation
The process to obtain the weak form of the symmetric gradient formulation is analogous to the one described above for the gradient FCFV formulation.First, the mixed form is written in the broken domain as

Mixed formulation
An alternative formulation involves selecting the mixed variable so that the viscosity function does not feature in the momentum equation.The strong form in the partitioned domain is It is worth noting that this formulation leads to an identical discrete system when the viscosity is constant, but this is not the case when the viscosity is a spatially varying function, which is the case of interest in this work.

Weak form of the local and global problems
Without detailing all the intermediate steps, the resulting weak form of the local problems is: for e = 1, … n el , find ) ) .Therefore it has an obvious physical meaning.

Weak form of the local and global problems
The resulting weak form of the local problems is: Similarly, the global problem is:

FCFV DISCRETISATION
To

Gradient formulation
Assuming a constant degree of approximation for all the variables in Equation (10a-d), that is for L e , u e and p e in each element and for û on each face (edge in two dimensions), the discrete local problems are ) where It is worth noting that Equation (22c) 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 (11b).
The local problem 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 of Equation (11a,b) becomes ∑ where and   e is the indicator function of the set  e , that is, By inserting the expressions of Equations (24a-c) into (26a,b), 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 ( f 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 References 21,20, the usual strategy of imposing mean zero pressure on the whole domain, that is, is considered here.

Symmetric gradient formulation
Following the same procedure described above, the local problem of the symmetric gradient formulation leads to three uncoupled equations that express L e , u e and p e as functions of the global unknowns û and  e , namely The blocks composing the matrix and vector of the global linear system are obtained by assembling the elemental contributions where and The blocks composing the matrix and vector of the global linear system are obtained by assembling the elemental contributions given by for i, j ∈  e .

Scaled symmetric gradient formulation
Following the same procedure described above, the local problem of the symmetric gradient 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 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

CRITICAL DISCUSSION
Two variations of the non-scaled gradient and symmetric gradient FCFV formulations are also considered.They are derived by performing only one integration by parts of the term containing the mixed variable L in Equations ( 7b) and (14b), respectively.This variation leads to the same discrete FCFV formulation but with a different definition of the coefficient  e , namely This variation is not applicable to the scaled formulations because the viscosity does not appear in the momentum equation, but in the equation associated to the mixed variable.The resulting global system of the gradient FCFV formulation, given by Equation (30a-d), is identical to the one derived for the symmetric gradient FCFV formulation, given by Equation (36a-d).Therefore, the only difference between the two formulations is the local problem associated to the mixed variable.Analogously, the global system of the scaled gradient FCFV formulation, given by Equation (42a-d), is identical to the one derived for the scaled symmetric gradient FCFV formulation, given by Equation (47a-d).
The global system of equations associated to the non-scaled FCFV formulations is not symmetric.The term  i , involving the integration of the viscosity field on the element faces, breaks the symmetry of the block Kûû defined in Equations (30a) and (36a).In contrast the global system of the scaled FCFV formulations maintains the symmetry even with a variable viscosity coefficient.
It is worth noting that the original FCFV formulation, presented in Reference 31, for Stokes flow problems with constant viscosity is recovered from the non-scaled gradient formulation of Section 3.1.When considering a constant approximation for the mixed variable, the first term of Equation (10b) vanishes, leading to the original FCFV formulation.Results, not reported here for brevity, show that the original formulation does not lead to satisfactory results when applied to problems with variable viscosity.The coupling of the velocity and the gradient of the velocity provided by the first term of Equation (10b) is crucial to obtain a formulation suitable for problems with variable viscosity.
In the context of problems with constant viscosity, there is no difference between using a scaled or non-scaled FCFV formulation.However, for problems with variable viscosity, the resulting formulations are different.Not only the scaled formulation leads to a symmetric system, but it also seems, a priori, more efficient.The scaled formulations only require the approximation of element integrals involving the viscosity field, whereas the non-scaled formulations require both the approximation of element and face integrals involving the viscosity.As it will be shown using numerical examples, a good approximation of these integrals is crucial to obtain the optimal approximation properties of the FCFV methods.

NUMERICAL EXAMPLES
Three numerical examples in two and three dimensions are considered to test the optimal approximation properties of the different FCFV formulations presented in this work.The accuracy and rate of convergence under mesh refinement are compared for the primal, mixed and hybrid variables.
The two-dimensional examples are also used to illustrate the need of using high order numerical quadratures when the viscosity present steep and localised gradients, whereas for other cases it shows that some of the FCFV formulations can lead to optimal results without the need of high order quadratures.
In all the numerical examples, the stabilisation parameter,  e , is selected as where x e denotes the centroid of the element Ω e , the parameter  is taken as 10 unless otherwise stated and  is a characteristic length, taken as 1 in the problems considered.Numerical experiments are included to demonstrate the suitability of this choice for a variety of problems.
It is worth noting that the stabilisation is limited using the maximum between (x e ) and one to ensure that for problems where the viscosity is nearly zero, enough stabilisation is introduced.
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 ⋆ 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 ( 8).
To compute the domain integrals that are involved in the evaluation of  u ,  p and  L , a high order quadrature in each triangle is employed.Similarly, to compute the integrals on the mesh skeleton that are involved in the evaluation of  û , a high order quadrature is used.This is done to ensure that the variation of the analytical solutions is properly captured.
A formal a priori error analysis of the methods proposed in this work can be performed by using the projection proposed in Reference 40.For the stabilisation selected here, independent on the element size and satisfying the properties stated in Reference 40, the main result of this work applies, meaning that the error for all the variables is expected to be (h).It is worth noting that following the analysis of Reference 40 for the methods in this work requires that the error introduced by the numerical integration of the variable viscosity is bounded by the spatial discretisation error.

The SolKz test: Exponential variation of the viscosity field
The first example considers the domain Ω = [0, 1] 2 and the viscosity given by where the parameter B is used to control the viscosity variation.In this example, B is taken equal to 6.9, as in Reference 41, to produce a variation of the viscosity in the domain of six orders of magnitude.The flow is driven by a source term s(x) = (x)g, with (x) = cos(3x 1 ) sin(2x 2 ) and g = (0, −1) T .Dirichlet boundary conditions, corresponding to the analytical solution, derived in Reference 42, are imposed in whole boundary of Ω.
Five triangular meshes are considered to perform the convergence studies.The meshes are generated by splitting into four triangles structured grids of quadrilateral elements.The number of elements in the five meshes is 256, 1024, 4096, 16,384 and 65,536 respectively.The first three meshes are shown in Figure 1.
The numerical solution computed with the scaled symmetric gradient FCFV formulation is shown in Figures 2 and 3 for the first and last mesh, respectively.The results illustrate the increase in resolution of the velocity field features near the bottom boundary as the mesh is refined.
The error and rate of convergence of the error of the hybrid variable û, the primal variable u, the pressure p and the mixed variable L is shown in Table 1 for the different FCFV formulations proposed in Section 4.
When employing the non-scaled FCFV formulations, the elemental integrals involving the viscosity are evaluated using a quadrature with a single integration point, namely  e ≈ |Ω e |(x e ), where x e is the centroid of the element.Similarly, the boundary integrals involving the viscosity are evaluated using a single integration point, namely  i ≈ |Γ e,i |(x e,i ), where x e,i is the centroid of the element edge.Analogously, when the scaled FCFV formulations are employed, the integral involving the viscosity is evaluated using a single integration point, namely  e ≈ |Ω e | −1 (x e ).
The results clearly show that the four non-scaled formulations do not provide convergence of any variable as the mesh is refined.In contrast, the scaled formulations provide optimal, (h), convergence on the hybrid velocity, the velocity and the pressure.If the symmetric gradient is considered as mixed variable, optimal convergence is also observed in the mixed variable, whereas suboptimal convergence is observed when the mixed variable is the velocity gradient.
TA B L E 1 SolKz: Error of the hybrid velocity, velocity, pressure and mixed variable and rates of convergence for the different FCFV formulations.Next, the same convergence study is repeated, but high order numerical quadratures are used to compute the integrals involving the viscosity.The results, shown in Table 2, demonstrate the importance of using an accurate numerical integration of the viscosity in the non-scaled formulations.The four non-scaled formulations show optimal convergence of the hybrid velocity, the velocity and the pressure.Similar to the scaled formulation, only when the symmetric gradient is considered as the mixed variable, optimal convergence in all variables is achieved.It is also worth noting that the scaled formulation does not benefit from using a higher order numerical quadrature as the results are almost identical to the ones shown in Table 1.

Mesh
This example shows important benefits of using the scaled formulation, in particular the scaled symmetric gradient formulation.This is the only approach that provides optimal convergence with the lowest order numerical quadrature.In addition, this formulation only requires elemental integrals of the inverse of the viscosity field, rather than integrals over faces or elemental integrals of the gradient of the viscosity.Finally, with a lower computational cost it, the scaled formulations provide a slightly more accurate results on the pressure field and the same accuracy on the rest of variables, when compared to the non-scaled formulations.
To further investigate the origin of the differences between the formulations considered, the numerical integration error is studied.The non-scaled formulation using the strong form in each element requires the computation of two integrals involving the viscosity, namely  e , given by Equation ( 23), and  i , given by Equation (27).The non-scaled formulation with one integration by parts also requires two integrals, but both on the element boundaries, namely γe , given by Equation (48), and  i .Finally, the scaled formulations only require the computation of one integral on the elements, namely  e given by Equation (38).
The numerical integration error as a function of the characteristic element size is represented in Figure 4 for quadratures of order one and two.The error is defined as the maximum relative difference, for all the elements, between the exact value of the integral and the approximation using a numerical quadrature.The results, show as expected an identical error for the computation of the element integrals  e and  e .This is due to the analytical expression of the viscosity of  e or γe .For the non-scaled gradient formulation this might be surprising as with an accurate computation of  i the dominant numerical integration error, when approximating  e , is identical to the numerical integration error of the scaled formulations, when approximating  e .However, there is an important difference between the two formulations that can be observed in the entries of the matrix Kûû of the global system.The matrix of the non-scaled formulations, given in Equation (30a), features the term  e multiplied by the stabilisation parameter  e .In contrast, the matrix of the scaled formulation, given in Equation (42a), does not contain the term  e multiplied by the stabilisation parameter.Given the high values of the stabilisation parameter, which is a multiple of the viscosity, a small variation in  e translates into a large variation in some entries of the global matrix in the non-scaled formulations, whereas this is not the case for the scaled formulation.As an example, when considering the second mesh, the maximum relative difference between the approximation of  e and  e with quadratures of order one and two is 3.4 × 10 −3 .This translates into a maximum relative difference in the entries of the global matrix of 81.6% in the case of the non-scaled formulation, but only 0.4% in the case of the scaled formulation.
Appendix A considers the same numerical example using quadrilateral meshes.For this example the results show that quadrilateral elements are able to provide optimal convergence in all the variables, even for the non-symmetric formulations.In terms of the accuracy, the performance of triangular meshes is shown to be superior for this example.
It is worth noting that the stabilisation parameter defined in Equation ( 49) induces a different stabilisation in an element edge as seen from the two elements sharing this edge.A different definition of the stabilisation parameter that ensures a unique stabilisation in an element edge has also been considered.The results, presented in Appendix B show that the accuracy and convergence are not influenced by this choice.
For purpose of reproducibility, an application of the scaled symmetric gradient FCFV formulation to the SolKz test is available via the open source Julia package (https://github.com/tduretz/FCFV_NME23).

Manufactured solution with a steep viscosity layer
The next example considers a viscosity field that presents a steep layer.The objective is to show the potential of the proposed approach to resolve an interface-type variation of the viscosity without the use of meshes that are fitted to the interface and without the need of an interface condition, as done in Reference 37.
The problem, taken from Reference 43, considers a viscosity field given by where  1 = 1 and  2 = 10 −4 .The analytical velocity field is given by and the exact pressure is The same meshes employed in the previous example are considered here.Figure 5A shows the viscosity field in logarithmic scale.A detailed view of the viscosity field and two meshes near the region that presents the steep viscosity layer is shown in Figure 5B,C, clearly illustrating that the mesh is not aligned with the very steep gradient of the viscosity field.In fact, for the first mesh considered in this study, the viscosity field varies four orders of magnitude within a single triangular element.
Figure 6 shows the numerical solution using the scaled symmetric gradient FCFV formulation on the coarsest mesh.Some artifacts in the central part of the domain can be clearly observed, denoting that the mesh is too coarse to capture the effects of the source induced by such a steep gradient of the viscosity field.
The numerical solution on the finest mesh, displayed in Figure 7, shows smooth velocity and pressure fields with no artifacts in the central part of the domain.In this mesh, the variation of four orders of magnitude on the viscosity field occurs over 14 triangular elements.The results of the mesh convergence study for all the different FCFV formulations and by using a quadrature with one integration point to compute the integrals involving the viscosity, are shown in Table 3.
In this example, not only the non-scaled FCFV formulations do not display the expected rate of convergence in all the variables, but also the scaled FCFV formulations.More precisely, for the scaled symmetric gradient formulation, sub-optimal convergence is observed in the pressure.
The same convergence study is next performed by increasing the order of the quadratures employed to compute the element and edge integrals involving the viscosity.The results, in Table 4, show that all formulations recover the expected rate of convergence in the velocity and pressure fields.The formulations with the mixed variable being the gradient of the Note: The elemental and boundary integrals involving the viscosity, namely  e , γe ,  i and  e are evaluated using a quadrature with one integration point.velocity, convergence is not observed in the mixed variable, whereas suboptimal convergence in the symmetric gradient is observed for all the non-scaled formulations.It is worth noting that, again, the scaled symmetric gradient formulation is the only one that provides optimal convergence in all the variables.A study of the numerical integration error is performed to further analyse the differences between the different formulations and the differences with respect to the previous example.The numerical integration error as a function of the characteristic element size is represented in Figure 8 for quadratures of order one and two.Contrary to the previous example, the more abrupt change of the viscosity profile poses a challenge when performing the numerical integration on meshes that do not account for the high gradients of the viscosity profile.It can be clearly observed that the numerical integration errors with quadratures of order one are between one and two orders of magnitude higher, when compared to the previous example.When a quadrature of order two is used, the numerical integration error is several orders of magnitude higher than in the previous example.
The results show that in this example, the main reason for not obtaining the optimal rate of convergence with a quadrature of order one is that the numerical integration error becomes higher than the spatial discretisation error.In this example, with a viscosity field showing a steep gradient, a quadrature of order two is enough in order to ensure that the numerical integration error is below the spatial discretisation error.
This example demonstrates the possibility of solving interface-type Stokes problems with mesh not fitted to the interface.Using the scaled symmetric gradient formulation an optimal rate of convergence is observed in all the variables when employing a quadrature of order two to compute  e .It is worth noting that the cost of computing  e with a high order quadrature is negligible when compared to the cost of assembling and solving the global system of equations.In addition, the use of a high order quadrature does not induce a higher use of memory as the value  e can be precomputed for all the elements and stored in a vector whose dimension does not depend on the order of the quadrature used.Note: The elemental and boundary integrals involving the viscosity, namely  e , γe ,  i and  e are evaluated using a quadrature with three integration point in each element and two integration points on each edge, respectively.
This approach might be particularly attractive in transient problems with moving interfaces, as it will enable the solution on a fixed mesh, whereas a formulation that requires a fitted mesh, needs to perform constant local re-meshing to ensure that the interface is captured by the mesh.With the proposed technique, only the vector containing the values of  e would need to be recomputed each time step.

Stokes flow in a spherical shell
The last example considers a three-dimensional Stokes flow in a spherical shell.The problem has an analytical solution, presented in Reference 44, which is recalled here using spherical coordinates (r, , ).The inner and outer radii of the spherical shell are taken as R 1 = 0.5 and R 2 = 1 respectively, and the viscosity field is given by The exact velocity vector field in spherical coordinates is where In the above expressions, the constants  and  are defined as Four unstructured tetrahedral meshes are considered to perform a mesh convergence study.The number of elements in the four meshes is 2223, 18,474, 148,734 and 1,203,063 respectively, leading to a global system of equations of dimension 14,781, 126,021, 1,027,779 and 8,366,694 respectively.The first three meshes are shown in Figure 9.
The numerical solution computed with the scaled symmetric gradient FCFV formulation on the finest mesh is shown in Figure 10.
Table 5 summarises the results of the mesh convergence study.The optimal rate of convergence is observed in all variables, confirming the correct implementation and the optimal approximation properties of the scaled symmetric gradient FCFV formulation for three dimensional problems with variable viscosity.

DISCUSSION
Among the different FCFV formulations presented here, the scaled symmetric gradient FCFV formulation appears as the clear choice.On the one hand, the use of the symmetrised velocity gradient tensor ensures the global symmetry of the discrete problem.On the second hand, the viscosity scaling explicitly introduces the deviatoric stress tensor in the definition of the mixed variable.The fact that this scheme preserves the symmetry of the global system of equation is a notable advantage.This enables the use of efficient sparse-direct solvers for linear systems 45 and will facilitate the design of efficient solutions strategies for non-linear problems. 46The scaled symmetric gradient FCFV formulation was shown to perform reliably on both triangular and quadrilateral elements, which are common choices in geodynamic modelling.While unstructured triangular meshes allow to honour material interfaces and dynamically adapt meshes to flow patterns, 45,37,47 structured quadrilateral meshes are popular for their relative algorithmic simplicity and efficiency on parallel supercomputers. 13,48,49The latter can also offer some degree of geometrical flexibility when combined with Arbitrary-Lagrangian-Eulerian formulations 50 or adaptative mesh refinement. 51,8Since the proposed FCFV formulations can handle both types of elements while remaining efficient and stable, it could become a method of choice for future geodynamic simulation codes. 36Moreover, the FCFV formulation also offers an adequate incorporation of internal boundaries, such as velocity discontinuities or surface tension, which often arise in the study of seismotectonics 52 or magmatic processes. 27esides achieving first order convergence of primal variables (velocity and pressure), the mixed variable (deviatoric stress tensor) also exhibits first-order convergence, with both triangular and quadrilateral elements.This property is crucial as it will allow for reliable implementations of non-linear stress-dependent rheologies such as plasticity.The latter is essential to model the emergence of tectonic plates in convection models 53 or the fault-related structures in the deforming lithosphere. 54Future developments will focus on the extension of FCFV towards the incorporation of complex rheological models 55 and the implementation in the context of GPU and multi-GPU computing. 56Besides steady geophysical flow, the findings of this work will be relevant for the simulation of turbulent incompressible flows where the viscosity is a spatially varying quantity. 57

CONCLUDING REMARKS
This work presents six different FCFV formulations to solve problems governed by the Stokes equations with spatially varying viscosity.The formulations differ in the definition of the mixed variable and the use of a weak or strong form in each element.A brief comparison of the different methods shows that only in some cases the resulting global system of equations is symmetric.Contrary to existing FCFV formulations for the Stokes problem with constant viscosity, the use of an accurate quadrature to compute the integrals involving the varying viscosity is found to be crucial to obtain the optimal approximation properties of the methods.In addition, a comparison between triangular and quadrilateral elements shows that quadrilateral elements are less sensitive to the numerical integration error.
Two-dimensional examples are considered to compare the accuracy and rate of convergence of the six FCFV formulations.The results show that when the numerical error associated to the numerical integration is below the error of the spatial discretisation, first order convergence is observed for all methods when the mixed variable is the symmetric gradient of the velocity.When considering formulations with the mixed variable being the gradient of the velocity, the optimal convergence of the mixed variable is not displayed in triangular meshes, whereas it is exhibited when using quadrilateral elements.The results also show that the most efficient and accurate approach is the so-called scaled symmetric gradient FCFV formulation, where the mixed variable is related to the deviatoric stress tensor.Not only the result are more accurate and less sensitive to numerical integration errors, but also this formulation only requires computing the integral of the viscosity within an element and not over element faces (edges in two dimensions).A three-dimensional example is also used to confirm the optimal convergence properties of the scaled symmetric gradient FCFV.
TA B L E A1 SolKz: Error of the hybrid velocity, velocity, pressure and mixed variable and rates of convergence for the different FCFV formulations with quadrilateral meshes.Note: The elemental and boundary integrals involving the viscosity, namely  e , γe ,  i and  e are evaluated using a quadrature with one integration point.

Mesh
In terms of the convergence rates, the four non-scaled formulations do not provide convergence of any variable as the mesh is refined, which is the same behaviour observed for triangular meshes.The scaled formulations provide optimal convergence, (h), convergence on all the variables.It is worth noting that when using triangular meshes, optimal convergence of the mixed variable was only observed when using the symmetric gradient formulation.
To further compare the results between triangular and quadrilateral meshes, the same convergence is repeated by employing a higher order numerical quadrature to compute the integrals involving the variable viscosity.The error and rate of convergence of the error of the hybrid variable û, the primal variable u, the pressure p and the mixed variable L is shown in Table A2.
Similar to the behaviour observed with triangular meshes, the non-scaled formulations now show optimal convergence of the hybrid velocity, the velocity and the pressure.For the mixed variable, almost optimal convergence rates are observed for both choices of the mixed variable, contrary to the case of triangular meshes, when only the symmetric formulations showed the optimal convergence rate of the mixed variable.Similarly, for the scaled FCFV formulation, all variables exhibit an optimal convergence, even for the non-symmetric formulation.
In summary, the use of quadrilateral meshes enable to provide optimal convergence of both the gradient of the velocity field and the symmetric gradient of the velocity field.In contrast only the symmetric formulations are able to provide optimal convergence of the mixed variable when using triangular meshes.
In terms of the accuracy, it can be clearly observed that triangular meshes offer a superior performance.To compare the results, it should be noted that the i-th mesh of triangular elements contains the same number of elements as the i + 1-th mesh of quadrilateral elements.As an example, comparing the results of Tables 2 and A2, it can be seen that the accuracy of the hybrid velocity computed in the first triangular mesh is comparable to the accuracy obtained in the fourth quadrilateral mesh.For the velocity, pressure and symmetric gradient of the velocity, the accuracy obtained in the first triangular mesh is similar to the accuracy obtained in the third quadrilateral mesh.This implies that the size of 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 .

5 6 7
Steep layer viscosity: Viscosity field in logarithmic scale, showing two detailed views of the variation of the viscosity in the elements of the first two meshes.(A) (x).(B) Mesh 1. (C) Mesh 2. Steep layer viscosity: Numerical solution on the first mesh using the scaled symmetric gradient FCFV formulation.(A) u 1 .(B) u 2 .(C) ||u||.(D) p. Steep layer viscosity: Numerical solution on the fifth mesh using the scaled symmetric gradient FCFV formulation.(A) u 1 .(B) u 2 .(C) ||u||.(D) p.

8 9
Steep layer viscosity: Numerical integration error, for quadratures of order 1 and 2, as a function of the characteristic element size for the four integrals involving the viscosity.(A) Order 1. (B) Order 2. First three tetrahedral meshes used for the Stokes flow in a spherical shell.(A) Mesh 1. (B) Mesh 2. (C) Mesh 3. and the exact pressure field is given by p(r, , ) = 4 r (r)F(r) cos().

10
Stokes flow in a spherical shell: Numerical solution on the finest mesh using the scaled symmetric gradient FCFV formulation.(A) u 1 .(B) u 2 .(C) u 3 .(D) ||u||.(E) p. TA B L E 5 Stokes flow in a spherical shell: Error of the hybrid velocity, velocity, pressure and mixed variable, and rates of convergence, for the scaled symmetric gradient FCFV formulation.
e I n sd + |Ω e | −1 ( e ⋅ n j )I n sd + |Ω e | −1 n j ⊗  e e |Γ e,j | e , (33a)  e ∶= |Ω e |s e + ∑ j∈ e |Γ e,j | e f j + |Ω e | −1  e Z e , (33b) Z e ∶= ∑ j∈ e |Γ e,j |(n j ⊗ f j + f j ⊗ n j ).e=1 |Γ e,i |g i   e (i) for i ∈  e , i | e  ij I n sd (36a) ( Kû ) e i ∶= |Γ e,i |n i (36b) ( fû ) e i ∶= −|Γ e,i | −1 e  e  e +  i |Ω e | −1 n i Z e − |Γ e,i |g i   e (i), (36c) ( f ) e ∶= − ∑ j∈ e |Γ e,j |n j ⋅ f j .Following the same procedure described in the previous section, the local problem leads to three uncoupled equations that express L e , u e and p e as functions of the global unknowns, û and  e , namely L e = − e Z e −  e ∑ j∈ e |Γ e,j |n j ⊗ ûj , (37a) u e =  −1 e  e +  −1 e ∑ j∈ e |Γ e,j | e ûj , (37b) j ∶= |Γ e,i | {  −1 e  e  e |Γ e,j |I n sd −  e |Γ e,j | [ (n i ⋅ n j )I n sd + n j ⊗ n i ] −  e  ij I n sd  e  e +  e |Ω e | −1 n i Z e − g i   e (i) Steep layer viscosity: Error of the hybrid velocity, velocity, pressure and mixed variable and rates of convergence for the different FCFV formulations.
TA B L E 3 Steep layer viscosity: Error of the hybrid velocity, velocity, pressure and mixed variable and rates of convergence for the different FCFV formulations.