the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An efficient partialdifferentialequationbased method to compute pressure boundary conditions in regional geodynamic models
Anthony Jourdon
Dave A. May
Modelling the pressure in the Earth's interior is a common problem in Earth sciences. In this study we propose a method based on the conservation of the momentum of a fluid by using a hydrostatic scenario or a uniformly moving fluid to approximate the pressure. This results in a partial differential equation (PDE) that can be solved using classical numerical methods. In hydrostatic cases, the computed pressure is the lithostatic pressure. In nonhydrostatic cases, we show that this PDEbased approach better approximates the total pressure than the classical 1D depthintegrated approach. To illustrate the performance of this PDEbased formulation we present several hydrostatic and nonhydrostatic 2D models in which we compute the lithostatic pressure or an approximation of the total pressure, respectively. Moreover, we also present a 3D rift model that uses that approximated pressure as a timedependent boundary condition to simulate farfield normal stresses. This model shows a high degree of noncylindrical deformation, resulting from the stress boundary condition, that is accommodated by strikeslip shear zones. We compare the result of this numerical model with a traditional rift model employing freeslip boundary conditions to demonstrate the firstorder implications of considering “open” boundary conditions in 3D thermomechanical rift models.
 Article
(10162 KB) 
Supplement
(4280 KB)  BibTeX
 EndNote
In Earth sciences and geodynamic modelling, computing the pressure can be essential. Specifically, numerous regional thermomechanical studies use the lithostatic pressure or a reference pressure based on some density structure as a normal stress boundary condition (e.g. Baes et al., 2018; Brune, 2014; Brune et al., 2012, 2014, 2017; Chertova et al., 2012, 2014; Glerum et al., 2018; IsmailZadeh et al., 2013; Popov and Sobolev, 2008; Quinteros et al., 2010; Yamato et al., 2008). By imposing only the normal stress, material is permitted to flow in and out of the domain in response to the other boundary conditions and or deformation in the domain interior. This is inherently closer to the reality of the dynamics within a regional segment of the Earth,compared to a regional domain that is closed and in which neither inflow nor outflow is permitted. Hence, the ultimate intent of imposing the normal stress is to provide dynamical behaviour that is similar to that which would occur if the models were performed in a much larger domain. Moreover, the reference pressure can also be used as an initial guess for the pressure when solving linear or nonlinear system flow problems with iterative methods.
The common approach to compute a reference pressure P is to define a set of depth columns and integrate the rock density ρ(x) over each column to obtain the pressure at depth. Thus, to compute the pressure P at some point of x^{′}, we evaluate the 1D integral as follows:
where ${\mathit{x}}_{\mathrm{s}}^{\prime}$ is the projection of x^{′} onto the surface of the Earth in the direction opposite to the gravity vector g and P_{s} is the reference pressure at the surface ${\mathit{x}}_{\mathrm{s}}^{\prime}$. For the case of a constant density and gravity, this expression reduces to $P\left({\mathit{x}}^{\prime}\right)={P}_{\mathrm{s}}+\mathit{\rho}gD$, where D is the distance (depth) given by $D=\Vert {\mathit{x}}_{\mathrm{s}}^{\prime}{\mathit{x}}^{\prime}\left\right$ and $g=\left\right\mathit{g}\left\right$. When the density is a function of space and gravity is constant, the 1D integral is decomposed into different segments D_{i} and a suitable quadrature rule is applied over each segment. For example using a onepoint Gauss quadrature rule we have
where ρ_{i} is the density at the centroid of the segment D_{i}. For the case of a uniform mesh with cell edges aligned with the gravity vector, all the cell edges and vertices are located along straight lines that are parallel to the direction of gravity. Hence, Eq. (2) can be simply evaluated by traversing along a column associated with a set of cells (or vertices). In this special case, the subdivision of the integral is naturally defined by mesh cells. If the column sweep is performed from the surface to depth, then only a single pass over each cell in a column is required to compute the pressure at any depth within that column by accumulating values from cells at shallower depths. Therefore, if we traverse from segments i = 0, 1, 2, …, N, where the segments are ordered such that D_{i+1} is located at greater depth than D_{i}, then we have the following sequence ${P}_{\mathrm{0}}={P}_{\mathrm{s}}+{\mathit{\rho}}_{\mathrm{0}}g{D}_{\mathrm{0}}$, ${P}_{\mathrm{1}}={P}_{\mathrm{s}}+{\mathit{\rho}}_{\mathrm{0}}g{D}_{\mathrm{0}}+{\mathit{\rho}}_{\mathrm{1}}g{D}_{\mathrm{1}}={P}_{\mathrm{0}}+{\mathit{\rho}}_{\mathrm{1}}g{D}_{\mathrm{1}}$, …, ${P}_{N}={P}_{\mathrm{s}}+{\sum}_{i}^{N}{\mathit{\rho}}_{i}g{D}_{i}={P}_{N\mathrm{1}}+{\mathit{\rho}}_{N}g{D}_{N}$.
Although evaluating Eq. (2) may appear simple, its implementation may be inefficient or too algorithmically complex for general use. Below we outline some common use cases that render the columnwise integration difficult (or expensive):

a mesh with cell edges (2D) or faces (3D) that are not aligned with the gravity vector (Fig. 1a),

an unstructured mesh (Fig. 1b),

a density structure (or gravity vector) that is spatially varying,

a parallel decomposition of the mesh (Fig. 1c),

time dependence in the density or mesh coordinates that requires continual reevaluation of the reference pressure.
To compute P(x^{′}) we first have to define the location ${\mathit{x}}_{\mathrm{s}}^{\prime}$. In general this is nontrivial for use cases (1) and (2). If both the density and gravity are constant, then the only complexity associated with meshes identified in points (1) and (2) relate to computing ${\mathit{x}}_{\mathrm{s}}^{\prime}$. Due to the fact that the path of the integral (i.e. the “column”) does not coincide in general with a set of mesh cells or vertices, the line integral must be performed for each point x^{′} in the mesh, meaning that the single pass approach used in the gravityaligned mesh is not possible. If the density (or gravity) vary in space throughout the domain, the integral must be approximated via a suitable subdivision in space and/or a quadrature rule. Assuming that the density is a piecewise constant over each cell, the simplest approximation would be to determine the intersection between the line segment $[{\mathit{x}}^{\prime},{\mathit{x}}_{\mathrm{s}}^{\prime}]$ and each cell and apply a onepoint quadrature rule over the intersecting segment times. In Fig. 1a and b we depict the complexity of this procedure for a noncoordinatealigned and unstructured mesh. When performing simulations in parallel where the mesh is distributed across multiple Message Passing Interface (MPI) ranks, even for the case of a uniform mesh aligned with the gravity vector, the columnwise integration approach is somewhat complicated. Individual MPI ranks may compute their local contribution to the sum of accumulated pressures; however, the final pressure requires a partial sum to be performed over mesh subdomains that intersect the 1D line integral. The global reduction (with the chosen MPI ranks overlapping with each 1D line integral) is complicated to define for mesh types identified in points (1) and (2). Lastly, if the reference pressure associated with some density structure is to be used as a boundary condition in a mechanical model, time dependence of that density structure (or mesh) will require one to recompute the reference pressure at each time step. Hence, the efficiency of the implementation used to compute the pressure is important.
Moreover, when the density structure evolves with time as deformation occurs, the pressure gradients may no longer be aligned with the gravitational acceleration vector. In these nonhydrostatic cases, this pressure is not lithostatic. However, to be able to provide an approximation for the total pressure or to use stress boundary conditions, it is still important to approximate the total pressure in the best possible way.
For these reasons, we propose an efficient mesh and numerical method (finite elements, finite differences, finite volumes, etc.) to compute a reference pressure associated with the density structure of a domain in hydrostatic cases or an approximation of the total pressure for nonhydrostatic cases for all scenarios given above by solving a partial differential equation (PDE) derived from the conservation of the noninertial momentum equation for an incompressible fluid. We also present thermomechanical numerical models and static numerical models applied to Earth sciences and geodynamics to show the usefulness of this approach.
For an incompressible fluid in a domain Ω, the noninertial form of the conservation of momentum is given by the Stokes equation
where v is the velocity of the fluid, P is the total pressure, $\mathit{\tau}=\mathrm{2}\mathit{\eta}\dot{\mathit{\epsilon}}\left(\mathit{v}\right)$ is the deviatoric stress tensor with η the viscosity, $\dot{\mathit{\epsilon}}\left(\mathit{v}\right)$ is the strain rate tensor, $\mathit{\rho}:=\mathit{\rho}(\mathit{x},t)$ is the density, $\mathit{g}:=\mathit{g}(\mathit{x},t)$ is the gravity vector, and x and t denote the space and time, respectively. The incompressibility constraint is given as follows:
In the context of our problems we will decompose the boundary of the domain into two nonoverlapping segments: ∂Ω_{surf} that we will regard as the free surface and prescribe that tangential and normal stresses are zero, i.e. $\mathit{\tau}P\mathbb{I}=\mathbf{0}$, and ∂Ω_{i}, which denotes the interior parts of the boundary along which we may impose any valid combination of velocity or stress in the normal and tangential directions. Furthermore, $\partial \mathrm{\Omega}=\partial {\mathrm{\Omega}}_{\text{i}}\cup \partial {\mathrm{\Omega}}_{\mathrm{surf}}$ and $\partial {\mathrm{\Omega}}_{\text{i}}\cap \partial {\mathrm{\Omega}}_{\mathrm{surf}}=\mathrm{\varnothing}$. The outwardpointing unit normal vector to ∂Ω will be denoted via $\widehat{\mathit{n}}$.
To define the pressure associated with the density structure we make the “ansatz” that v=0; hence, Eq. (4) is trivially satisfied, and Eq. (3) reduces to the usual hydrostatic equilibrium problem
For spatial dimensions ${n}_{\mathrm{d}}=\mathrm{2},\mathrm{3}$, Eq. (5) is overdetermined as there are more equations (n_{d}) than unknowns:
As such, there is no unique solution to Eq. (5). To obtain a unique solution to Eq. (5), we require a single equation for P. This can be achieved by taking the divergence of Eq. (5):
Eq. (7) will be referred to as the pressure Poisson equation (PPE).
Equation (7) can be obtained in an alternative manner with less restrictive assumptions. First, we assume that we have a solution ((v),p) for Eqs. (3), (4). Then we take the divergence of the momentum equation (as before) and integrate over Ω
Then, applying the divergence theorem to the first term we obtain
If we assume the boundary term on the lefthand side is small, we have
which must be true for any arbitrary domain, thus resulting in Eq. (7).
Assuming that $\widehat{\mathit{n}}$ is constant or slowly varying and that τ is symmetric yields
Hence, dropping the first term in Eq. (9) is equivalent to saying that $\mathit{\tau}\widehat{\mathit{n}}\approx \mathbf{0}$ for all $\mathit{x}\in \partial \mathrm{\Omega}$. Alternatively, the term is zero if $\mathrm{\nabla}\cdot \mathit{\tau}=\mathbf{0}$ for all $\mathit{x}\in \partial \mathrm{\Omega}$. This condition is satisfied if the fluid experiences either rigid body translations or rigid body rotation along ∂Ω.
2.1 Boundary conditions
A unique solution to Eq. (7) requires boundary conditions to be specified on P. Our choice of boundary conditions for Eq. (7) is motivated by Earthlike bodies. The boundary conditions will be specified in the usual manner, i.e. in terms of a Dirichlet constraint in which we impose P and Neumann constraints in which we impose the behaviour of $\mathrm{\nabla}P\cdot \widehat{\mathit{n}}$.
Along the surface of the domain, which represents the free surface of the Earth, we impose
Equation (12) is a Dirichlet constraint and specifies that the reference (or datum) pressure should be zero on the surface of our geological body. This is consistent with the observation that the mean pressure on all points on the surface (above sea level) are approximately equal. This Dirichlet boundary condition is a natural extension of the free surface boundary condition used for the flow problem in Eqs. (3), (4), namely $\widehat{\mathit{n}}\cdot (\mathit{\tau}P\mathbb{I})\widehat{\mathit{n}}=\widehat{\mathit{n}}\cdot (\mathit{\tau}P\mathbb{I})\widehat{\mathit{t}}=\mathrm{0}$, which reduces to $\widehat{\mathit{n}}\cdot \mathit{\tau}\widehat{\mathit{n}}P=\mathrm{0}$ and $\widehat{\mathit{n}}\cdot \mathit{\tau}\widehat{\mathit{t}}=\mathrm{0}$ with $\widehat{\mathit{t}}$ a tangent unit vector to the boundary such that $\widehat{\mathit{n}}\cdot \widehat{\mathit{t}}=\mathrm{0}$. Equation (12) is consistent with a fluid at rest since τ=0. In the nonhydrostatic case, we require $\mathit{\tau}\widehat{\mathit{n}}\approx \mathrm{0}$ to arrive at Eq. (12). We also note that $\widehat{\mathit{n}}\cdot \mathit{\tau}\widehat{\mathit{n}}$ is proportional to the mean curvature κ of the boundary ∂Ω (Barth and Carey, 2007). Hence, if there is zero topography, κ=0 and P=0 on ∂Ω_{surf}, and if the change in topography is small, then κ≈0 and P≈0.
Two different boundary conditions are introduced to constrain $\mathrm{\nabla}P\cdot \widehat{\mathit{n}}$. These are defined as a direct extension of the 1D hydrostatic assumptions to 2D and 3D domains. We first introduce some additional quantities that will aid the definition of the Neumann boundary condition. First, we split ∂Ω_{i} into two parts, such that $\partial {\mathrm{\Omega}}_{i}=\partial {\mathrm{\Omega}}_{\u27c2}\cup \partial {\mathrm{\Omega}}_{\parallel}$. Next, we define the gravity unit vector $\widehat{\mathit{g}}$ such that
and the unit vector ${\widehat{\mathit{g}}}_{\u27c2}$, which is perpendicular to $\widehat{\mathit{g}}$, i.e.
The first constraint states that P should increase only in the direction of gravity. Hence, from Eq. (5) we have
The second constraint states that P should not change along directions perpendicular to the gravity; hence,
Since the unit vector $\widehat{\mathit{n}}$ normal to the boundary ∂Ω can be decomposed according to
we have
Hence, the two Neumann boundary conditions may be stated as
and
Equations (19) and (20) may appear peculiar since both the lefthand side and righthand side involve the gradient of pressure. In principle, to obtain a unique solution to Eq. (7) one can constrain the gradient in any direction, independent of the boundary normal $\widehat{\mathit{n}}$, and our 1D inspired gradients do exactly that.
We note for domains with boundaries parallel to either $\widehat{\mathit{g}}$ or ${\widehat{\mathit{g}}}_{\u27c2}$ that the Neumann conditions (19) and (20) simplify (and in some cases do not provide any) the information to constrain $\mathrm{\nabla}P\cdot \widehat{\mathit{n}}$. For example, consider a 2D Cartesian domain's right and left boundaries Γ_{r,l} with normal $\widehat{\mathit{n}}=(\pm \mathrm{1},\mathrm{0})$ and a bottom boundary Γ_{b} with normal $\widehat{\mathit{n}}=(\mathrm{0},\mathrm{1})$ and $\widehat{\mathit{g}}=g(\mathrm{0},\mathrm{1})$. Invoking Eq. (19), we obtain
Invoking Eq. (20) we obtain
Clearly conditions (21a) and (22b) are redundant. As such, the usage of Eqs. (19) and (20) cannot be used arbitrarily.
One may also consider employing both Eqs. (19) and (20) simultaneously. In this way we would obtain
which is identical to the result obtained by computing the dot product of Eq. (5) with $\widehat{\mathit{n}}$. Equations (23) certainly avoids the potential issue of using Eqs. (19) and (20). If we again consider the 2D Cartesian domain example, imposing Eq. (23) on all of ∂Ω_{i} is equivalent to imposing Eqs. (21b) and (22a).
Nevertheless, for an arbitrarily shaped domain, using boundary conditions (15), (16) or (23) does not yield the same result (see Sect. 3.3). For general use (i.e. when considering arbitrarily shaped domains), we suggest employing Eqs. (15) and (16) as they are a direct extension of the 1D hydrostatic assumptions to 2D and 3D domains.
2.2 Weak formulation
To define the weak formulation of the PPE we will use functions that are square integrable in the sense of Lebesgue, i.e.
and functions from the H_{1}(Ω) Sobolev space
Finally we will require the space of functions in H_{1}(Ω) that vanish on the Dirichlet boundary ∂Ω_{surf}:
Given a test function $q\in {H}_{\mathrm{1}}^{\mathrm{d}}\left(\mathrm{\Omega}\right)$, the weak form of the PPE is obtained by multiplying Eq. (7) by q and integrating both sides over Ω
Applying integration by parts to the left and righthand sides yields
Note that the boundary ∂Ω_{surf} does not appear in Eq. (25) since the test function q vanishes along the Dirichlet boundary. We also note that Eq. (25) only requires ρg∈L_{2}(Ω), and thus the formulation is valid for cases when the density ρ is discontinuous.
Splitting the surface integrals over the two segments $\partial {\mathrm{\Omega}}_{i}=\partial {\mathrm{\Omega}}_{\u27c2}\cup \partial {\mathrm{\Omega}}_{\parallel}$ and using Eqs. (19), (20) we have
Noting that the second term of the lefthand side and part of the last term on the righthand side exactly cancel each other yields the following
From Eq. (25), the weak formulation obtained if using Eq. (23) applied over all of ∂Ω_{i} is simply
2.3 Implementation
The strong (Eq. 7) and weak (Eq. 25) formulations of the pressure Poisson problem can be solved using the standard spatial discretization techniques, e.g. finite differences or finite elements. Moreover, since the equation is of Poisson type, it is readily amenable to being solved using standard iterative multigrid and/or direct solvers. Lastly, because the formulation is expressed in terms of a PDE, it is also straightforward to compute the pressure on parallel computing architecture as we can reuse existing discretization implementations that support domain decomposition.
In this section we provide several numerical models to show the following items:

the accuracy and consistency of the method for hydrostatic cases,

the accuracy of the approximation of the total pressure in nonhydrostatic cases,

the effect of using the depthintegrated approach and the pressure Poisson problem approach to impose boundary conditions on the momentum equation,

the usefulness of the method for 3D geodynamic thermomechanical modelling.
3.1 Hydrostatic pressure
To compare the numerical solution of Eq. (7) with the analytical solution obtained with Eq. (2), we designed four hydrostatic models (in 1D and 2D) for which the analytical solution is easily obtained (Figs. 2 and 3).
3.1.1 Box domain
We define the domain $\mathrm{\Omega}=x\in [\mathrm{0},\mathrm{1}]\times y\in [\mathrm{0},\mathrm{1}]$ and assume that $\mathit{g}=(\mathrm{0},\mathrm{1})$. We consider three depthdependent density structures ρ=ρ(y) that thus admit a hydrostatic pressure solution, i.e. satisfy $\partial P/\partial x=\mathrm{0}$, $\partial P/\partial y=\mathit{\rho}\left(y\right)g$.
 Case 1.

This case assumes a constant density, ρ(y)=1 (Fig. 2a, b). The analytic pressure solution is given by
$$\begin{array}{}\text{(29)}& \begin{array}{rl}P\left(y\right)& =\underset{y}{\overset{{y}_{\mathrm{s}}=\mathrm{1}}{\int}}\mathit{\rho}g\phantom{\rule{0.25em}{0ex}}\mathrm{d}y={\left[\mathit{\rho}gy\right]}_{y}^{{y}_{\mathrm{s}}}=\mathit{\rho}g({y}_{\mathrm{s}}y)\\ & ={y}_{\mathrm{s}}y=\mathrm{1}y.\end{array}\end{array}$$  Case 2.

This case assumes a continuous depthvarying density $\mathit{\rho}\left(y\right)=\mathrm{2}y$ (Fig. 2c, d). The analytic pressure solution is given by
$$\begin{array}{}\text{(30)}& \begin{array}{rl}P\left(y\right)& =\underset{y}{\overset{{y}_{\mathrm{s}}=\mathrm{1}}{\int}}\mathit{\rho}g\phantom{\rule{0.25em}{0ex}}\mathrm{d}y=\underset{y}{\overset{{y}_{\mathrm{s}}=\mathrm{1}}{\int}}g(\mathrm{2}y)\phantom{\rule{0.25em}{0ex}}\mathrm{d}y\\ & ={\left[g\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{y}^{\mathrm{2}}+\mathrm{2}y\right)\right]}_{y}^{{y}_{\mathrm{s}}}\\ & =g\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({y}^{\mathrm{2}}{y}_{\mathrm{s}}\right)+\mathrm{2}({y}_{\mathrm{s}}y)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{y}^{\mathrm{2}}\mathrm{2}y+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}.\end{array}\end{array}$$  Case 3.

This case assumes a discontinuous density such that ρ(y)=1 for $y\in [\mathrm{0.5},\mathrm{1}]$ and ρ(y)=2 for $y\in [\mathrm{0},\mathrm{0.5})$ (Fig. 2e, f). The analytic pressure solution is given by
$$\begin{array}{}\text{(31)}& P\left(y\right)=\left\{\begin{array}{ll}\mathit{\rho}\left(y\right)g(\mathrm{1}y)=(\mathrm{1}y)& \forall y\mathit{\u2a7e}\mathrm{0.5}\\ & \\ \mathit{\rho}\left(y\mathit{\u2a7e}\mathrm{0.5}\right)gD& \forall y\mathrm{0.5}\\ \phantom{\rule{1em}{0ex}}+\mathit{\rho}\left(y\right)g(\mathrm{0.5}y)& \\ \phantom{\rule{1em}{0ex}}=D+\mathrm{1}\mathrm{2}y,& \end{array}\right)\end{array}$$with D the distance between the surface and the y coordinate at which the pressure is computed.
A finiteelement (FE) method employing an unstructured triangular mesh with a P_{2} function space was used to obtain the numerical solution for each case. The FE method was applied to Eq. (27) using boundary conditions described by Eq. (15) at the base and Eq. (16) on the lateral sides. Along the upper surface we impose the Dirichlet constraint P=0. Unless otherwise stated, when solving the PPE with FEs the Dirichlet constraints are imposed strongly (i.e. pointwise), whereas Neumann constraints are imposed weakly via surface integrals defined on facets of the FE cells that live on the boundary of the domain. Accordingly, all points living on ∂Ω_{surf} (including corner points) will be associated with the Dirichlet constraint. Figure 2a–f shows the 1D and 2D solution of these three models. On the 1D models both the numerical and analytical solutions of Eqs. (7) and (2) are shown, whereas on the 2D models only the numerical solution is provided. Since the P_{2} FE approximation contains the monomials 1, y, and y^{2}, the FE solution exactly reproduces the analytic solution for case 1 and case 2 independent of the number of finite elements used in the domain (e.g. subdividing the box into two triangles would be sufficient to obtain an exact solution). For case 3, the analytic pressure solution is piecewise linear, and provided that the density discontinuity is exactly resolved by the faces of the triangular FE mesh (which was the case here), the FE method exactly reproduces the analytic solution.
3.1.2 Halfannulus domain
The 2D halfannulus model aims to show the efficiency of the method when computing the lithostatic pressure in a body with a radial gravity vector and concentric density structure (Fig. 3a). This model represents a domain extending from
where θ is the polar angle, r is the radius in polar coordinates, R_{E} is the approximative Earth radius (6371 km), and R_{E}−2891 km is the approximative coremantle boundary. Mapped into Cartesian coordinates this gives x=rsin (θ) and y=rcos (θ). The gravity vector pointing to the centre ($\mathit{x}=(\mathrm{0},\mathrm{0})$) is defined as $\mathit{g}=\mathrm{9.8}\left(\frac{x}{\sqrt{{x}^{\mathrm{2}}+{y}^{\mathrm{2}}}},\phantom{\rule{0.25em}{0ex}}\frac{y}{\sqrt{{x}^{\mathrm{2}}+{y}^{\mathrm{2}}}}\right)$ m s^{−2}, and the density is defined as five concentric layers with a constant density in each (Fig. 3a). The pressure was computed by solving Eq. (27) with the boundary condition of Eq. (15) on the core–mantle boundary and the conditions of Eq. (16) on the sides parallel to $\widehat{\mathit{g}}$. At the surface of the domain we impose the Dirichlet constraint P=0. As in Sect. 3.1 we use a finiteelement discretization employing an unstructured mesh of triangles employing a P_{2} function space. Pressure values vary from 0 to 180 GPa with a concentric distribution following the density distribution and the gravity vector orientation (Fig. 3b).
The numerical solution extracted at x=0 along a line parallel to the gravity vector field reproduces the analytical solution computed for a 1D profile using Eq. (2) for the density distribution displayed in the halfannulus model (Fig. 3c). The difference between the pressure obtained using the depthintegrated approach and the PPE is very small (Fig. 3d). Unlike the analytic solution for the box models, the analytic solution for P here is nonpolynomial (in Cartesian coordinates); hence, the FE solution (which was discretized in Cartesian coordinates) cannot exactly reproduce the analytic solution.
The four hydrostatic models clearly illustrate that the solutions obtained using the depthintegrated approach and the pressure Poisson equation (with one set of boundary constraints) are equivalent for scenarios that admit a hydrostatic solution.
3.2 Nonhydrostatic pressure
Here, we show the differences and accuracy of the depthintegrated equation (Eq. 2) and the pressure Poisson equation (Eq. 7) approaches to approximate the total pressure. First, we compute and compare the pressure from the different methods in a large domain (referred to as the “global” model) containing a topographic perturbation (Fig. 4). Following this, we compute the pressure in a smaller domain (referred to as the “regional” model) and show the accuracy of the different methods to approximate the total pressure from the large domain (Figs. 5a, c and 6). Finally, we show the velocity field resulting from applying these approximated pressures as boundary conditions to solve the conservation of momentum in the small domain (Fig. 5b, d). The pressure Poisson problem was discretized and solved using the same FE method described in Sect. 3.1. The flow field was computed using the same underlying FE mesh and a mixed P_{2}P_{1} function space for velocity and pressure, respectively.
3.2.1 “Global” model
We define a large domain ${\mathrm{\Omega}}_{\mathrm{G}}=x\in [\mathrm{10},\mathrm{10}]\times y\in [\mathrm{0},\mathrm{1}]$ representing a global model, i.e. 20 times larger than the domain of interest, in order to avoid boundary condition influence (Fig. 4). In this domain, we introduce a topographic perturbation through a slope defined as ${y}_{\mathrm{s}}\left(x\right)=\frac{\mathrm{1}}{\mathrm{4}}x$, where y_{s} is the surface between x=0 and x=1. For demonstration purposes this topography is highly exaggerated with respect to the depth of the domain compared with the actual regional geodynamic models. Moreover, we use a constant density ρ=1 and a vertical gravity vector $\mathit{g}=(\mathrm{0},\mathrm{1})$.
We solve the flow problem described by Eqs. (3), (4) using a constant viscosity and density, along with the following boundary conditions: noslip at base, freeslip on the right and left faces, and a free surface along the top of the model domain. Due to the topography, a nonvertical pressure gradient that will drive flow is generated below that perturbation (Fig. 4a). The generated flow shows a velocity field characteristic of a gravitational collapse. Figure 4b shows the difference between the total pressure solution from Eq. (3) with the pressure P_{d} obtained by solving Eq. (7) using the boundary conditions described by Eq. (16) on the vertical sides, Eq. (15) on the bottom boundary, and P_{d}=0 on the upper surface. Figure 4c shows the difference with the pressure P_{a} obtained with Eq. (2).
The difference between the total pressure and the approximated pressure P_{d} is negligible in the nonperturbed domain and increases with the topography perturbation. It shows that as the system tends towards a hydrostatic state the total pressure and the approximated pressure P_{d} tends to the same value, i.e. the hydrostatic pressure. However, in the vicinity of the topographic perturbation the differences between the total pressure and the approximated pressure can be extremely small (Fig. 4b). In contrast, the difference between the total pressure and depthintegrated approximated pressure P_{a} is larger below the topographic slope, particularly below the points at which the slope begins and ends (Fig. 4c). In general, the pressure P_{a} obtained by applying the 1D solution is less accurate than P_{d} within both the interior and along the left, right, and bottom boundaries.
3.2.2 “Regional” model
Nevertheless, modelling a domain 20 times larger than the domain of interest can hardly be achieved in practice, mainly due to the numerical cost it represents. Thus, the boundary conditions are a firstorder component of regional models in order to best capture the global behaviour and interactions in a region without having to model the whole Earth. Therefore, we define a smaller domain ${\mathrm{\Omega}}_{R}=x\in [\mathrm{0},\mathrm{1}]\times y\in [\mathrm{0},\mathrm{1}]$ representing a regional model (Fig. 5). This domain represents the portion of the large domain in which the topographic slope is defined.
In this case we aim to apply a normal stress on the boundaries of our regional model that will generate a flow similar (or close) to the flow generated in the global model. To achieve this we present two models using Eqs. (7) and (2), respectively, to compute an approximated pressure that is then used as a boundary condition for the momentum equation (Eq. 3) on vertical boundaries as follows:
where ${P}_{\mathit{\alpha}}=\mathit{\{}{P}_{\mathrm{d}},{P}_{\mathrm{a}}\mathit{\}}$ denote the pressure computed using the pressure Poisson equation and the 1D depthintegrated approach, respectively. The viscous flow problem for the regional domain setting employs the following additional boundary conditions: a free surface on the top boundary and a noslip condition at the base. To solve for the pressure P_{d} using Eq. (7), we impose the following boundary conditions: P_{d}=0 on the surface, $\mathrm{\nabla}{P}_{\mathrm{d}}\cdot \widehat{\mathit{n}}=(\widehat{\mathit{n}}\cdot \widehat{\mathit{g}})\mathrm{\nabla}{P}_{\mathrm{d}}\cdot \widehat{\mathit{g}}$ on vertical sides, and $\mathrm{\nabla}{P}_{\mathrm{d}}\cdot \widehat{\mathit{n}}=(\widehat{\mathit{n}}\cdot \widehat{\mathit{g}})\mathit{\rho}g+(\widehat{\mathit{n}}\cdot \widehat{{\mathit{g}}_{\u27c2}})\mathrm{\nabla}{P}_{\mathrm{d}}\cdot \widehat{{\mathit{g}}_{\u27c2}}$ at the bottom.
Figure 5a shows the pressure field P_{d} computed on the regional domain and the total pressure P extracted from the global model. The approximated pressure P_{d} highlights differences with the total pressure from the global model, especially those regarding depth along the boundaries (blue lines in Fig. 6b, d). These differences are mainly due to the size of the domain that defines only the perturbed region without providing information about the domain in which it is enclosed. The boundary condition described by Eq. (16) used to solve Eq. (7) enforces the idea that the pressure gradient must be colinear with the gravity vector on the boundary. Therefore, given the definition of g, ∇P is enforced to be vertical on the vertical sides, whereas the deflection of the pressure field due to the topographic perturbation should occur in spatial offset from the topographic perturbation, as shown by the total pressure in the global model (Fig. 4a).
In Fig. 5b we show the velocity field resulting from solving the Stokes equation (Eqs. 3, 4) using P_{d} in Eq. (33). The flow field highlights velocities of around $\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}$ (velocity unit) at the surface pointing toward the right side of the domain, i.e. the bottom of the slope. Velocities progressively decrease at depth. While the orientation of the velocity field is more laminar than in the global model, its magnitude is very close, with the highest velocities being only 1.5× higher than in the global model.
Figure 5c displays the pressure field P_{a} computed from the depthintegrated approach (Eq. 2) within the regional domain and the total pressure P extracted from the global model. Because of the 1D behaviour of that method, the differences between the approximated pressure P_{a} and the total pressure P (red lines in Figs. 5c and 6b, d) are independent of the domain size and boundary conditions, and thus they are exactly the same as in the global model (Fig. 4c). The velocity field (Fig. 5d) resulting from using these approximated pressure values as a boundary condition to solve the momentum equation is approximatively 4 times higher than velocities obtained while using the pressure P_{d} (Fig. 5b) and approximatively 6 times higher than in the global model (Fig. 4a). As for the velocity field orientation, the vectors at the top of the slope show that the material uplifts, whereas in the global model velocities at the same location show a gravitational sliding. This orientation results from a stress value that is too high being imposed in the boundary conditions by the value of the pressure P_{a} at the boundaries (Fig. 6).
These simple tests demonstrate that the approximated pressure P_{d} computed from the Eq. (7) is more accurate for approximating the total pressure P than the pressure P_{a} computed from the Eq. (2). The only area where this is not true is located in the bottomleft corner of the regional model, and this is again due to the boundary condition not capturing the deflection of the pressure due to the size of the domain. Thus, the pressure computed with the pressure Poisson problem should be preferred for use as a boundary condition for the momentum equation. Moreover, as the domain size increases the error with the total pressure decreases, which is not the case with the depthintegrated approach due to its 1D nature.
3.3 Influence of the boundary condition type
The boundary conditions used to solve the pressure Poisson problem are stated in Eqs. (15), (16), and (23). To show the influence of the boundary conditions choices on the resulting pressure field, we define an irregular quadrilateral domain with the same topographic perturbation as the one described in Sect. 3.2, with a constant density ρ=1, a gravity vector $\mathit{g}=(\mathrm{0},\mathrm{1})$, and a Dirichlet boundary condition P_{d}=0 at the top surface (Γ_{s}). The irregular domain (shown in Fig. 7) is constructed such that none of the three boundary segments defining $\partial {\mathrm{\Omega}}_{i}={\mathrm{\Gamma}}_{\mathrm{1}}\cup {\mathrm{\Gamma}}_{\mathrm{2}}\cup {\mathrm{\Gamma}}_{\mathrm{3}}$ are parallel or perpendicular to g. The PPE was solved using the same FE method described in Sect. 3.1.
In Fig. 7 we show different pressure fields P_{d} obtained using several different boundary condition configurations. These results show that the boundary conditions described by Eq. (16) (or Eq. 20) force the pressure gradient to be parallel to g, i.e. ${\partial}_{x}{P}_{\mathrm{d}}=\mathrm{0}$ (Fig. 7b segment Γ_{3}, Fig. 7c segments Γ_{1,3}, and Fig. 7d segments Γ_{1,2}). We also confirm that boundary conditions described by Eq. (15) (or Eq. 19) constrain ∇P_{d} to be equivalent to the 1D solution of Eq. (2) on the boundary (Fig. 7a segments ${\mathrm{\Gamma}}_{\mathrm{1},\mathrm{2},\mathrm{3}}$, Fig. 7b segments Γ_{1,2}, Fig. 7c segment Γ_{2}, and Fig. 7d segment Γ_{3}). Figure 7a highlights that the solution of the PPE can be identical to that obtained using the depthintegrated approach with a specific choice of boundary conditions.
Moreover, Fig. 8a shows the pressure field P_{d} obtained using the boundary condition described by Eq. (23). The pressure solution is relatively similar to the solution obtained in Fig. 7c. However, as shown in Fig. 8b, the solution along the boundary (and therefore also in the interior) differs since Eq. (23) does not enforce the pressure gradient to be vertical on the boundary compared with the condition (16).
As previously noted, discretizations employing the weak form with Eq. (23) are certainly simpler to implement than imposing the Eqs. (19) and (20) as all the surface integrals cancel (see Eq. 28), and thus no surface integrals appear in either the linear form or bilinear form. We also reiterate that when the domain has boundaries that are aligned with g, using the boundary conditions (15) on the boundaries orthogonal to g and the condition (16) on the boundaries that are parallel to g results in an identical formulation.
3.4 Thermomechanical model
3.4.1 Physical model
To simulate the longterm evolution of the deformation of the lithosphere, we solve the stationary, noninertial form of the conservation of momentum described by Eq. (3) with the incompressible constraint (Eq. 4). Moreover, to consider the temperature variations in the domain, the following timedependent conservation of energy is solved:
where T is the temperature, t is the time, k is the thermal conductivity, H is the heat source, ρ_{0} is the reference density, and C_{p} the thermal heat capacity.
The numerical solution of Eqs. (3) and (4) is obtained using a mixed finiteelement method that independently discretizes the velocity and pressure fields. Hence, the numerical velocity and pressure obtained are solutions of the weak form of the Stokes problem given by
where w∈H_{1}(Ω) and q∈L_{2}(Ω) are test functions for the velocity and pressure, respectively, Γ_{N} denotes the Neumann boundary, T denotes the traction vector given by $\mathit{T}(\mathit{v},p)=\left(\mathit{\tau}\right(\mathit{u})p\mathbb{I})\widehat{\mathit{n}}$, and $\widehat{\mathit{n}}$ is the outward pointing normal vector from the boundary. The bilinear forms for the Stokes problem are given by (Elman et al., 2014):
Both the Stokes and thermal problems were solved using the parallel finiteelement code pTatin3D
(May et al., 2014, 2015),
which employs a mixed Q_{2}P_{1} discretization for velocity and pressure.
3.4.2 Initial conditions and rheology
To model the strain localization we use nonlinear viscoplastic rheologies expressed in term of viscosity. The ductile parts of the domain are simulated using an Arrhenius flow law for dislocation creep
where A, Q, and n are materialdefined parameters (see Table 1); R is the universal gas constant; V is the activation volume; and ${\dot{\mathit{\epsilon}}}^{\mathrm{II}}$ is the square root of the strain rate second invariant, computed as
The brittle parts of the domain are simulated using a Drucker–Prager yield criterion adapted to continuum mechanics, which is given by
where C is the cohesion of the material and ϕ is the friction angle.
The modelled domain contains four initial flat layers representing the upper continental crust, the lower continental crust, the lithosphere mantle, and the asthenosphere mantle, respectively (Fig. 9a). The upper crust extends from the surface of the domain (y=0 km) to $y=\mathrm{25}$ km and is modelled with a dislocation creep quartz rheology (Ranalli, 1997). The lower crust extends from $y<\mathrm{25}$ km to $y=\mathrm{35}$ km and is modelled with a dislocation creep anorthite rheology (Rybacki and Dresen, 2000). The lithosphere mantle extends from $y<\mathrm{35}$ km to $y=\mathrm{120}$ km, whereas the asthenosphere mantle extends from $y<\mathrm{120}$ km to $y=\mathrm{450}$ km. They are both modelled using a dislocation creep olivine flow law (Hirth and Kohlstedt, 2003).
The initial density distribution follows the lithologies and is reported in Table 1. In addition, the density varies with pressure and temperature following the Boussinesq approximation
where ρ_{0} is the reference density at T_{0} and P_{0}, P is the total pressure computed from the conservation of momentum (Eq. 3) and continuity equation (Eq. 4), α is the thermal expansion and β the compressibility. The Boussinesq approximation states that perturbations of density, if sufficiently small, can only be considered in the buoyancy term and neglected elsewhere regardless of the origin of the perturbation.
Moreover, the initial temperature field is computed as a steadystate solution of the heat equation
using a surface temperature of T=0 ^{∘}C at y=0 km and T=1450 ^{∘}C at $y=\mathrm{450}$ km. Moreover, to simulate an adiabatic thermal gradient in the asthenosphere due to thermal convection, the initial temperature field is solved with a conductivity of k=70 W m^{−1} K^{−1} in the asthenospheric mantle. However, for the actual model run we used a more realistic conductivity of k=3.3 W m^{−1} K^{−1} to solve Eq. (34). Other thermal parameters are reported in Table 1.
3.4.3 Boundary conditions
To show the influence of the normal stress boundary condition, we compare two rift models. In the reference model, an extension velocity of v_{x}=1 cm yr^{−1} is applied on the whole faces of normal x, whereas on faces of normal z a freeslip boundary condition is applied (Fig. 9c). To ensure mass conservation we impose an inflow velocity on the bottom face of normal y to balance any outflow that occurs due to the imposed extension. Along the surface of the model we use a free surface (zero normal stress, zero tangential stress) boundary condition.
The second rift model (Fig. 9b) uses the same Dirichlet boundary conditions on faces of normal x. On the faces of normal z we impose a Neumann boundary condition as
where P_{d} is the pressure computed with Eq. (7) and $\widehat{\mathit{n}}$ is the normal vector pointing outward from the domain.
To account for the density evolution through time due to the deformation and material advection, Eq. (7) is solved at every nonlinear iteration for each time step, and the Neumann boundary condition described by Eq. (40) is evaluated at every nonlinear iteration. Using ρ(P,T) computed from Eq. (38) to evaluate the pressure P_{d} going into the boundary condition described by Eq. (40) adds a new nonlinearity to the system.
The bottom of the domain is prescribed as an inflow condition balancing the outflow, and the surface of the domain is a free surface where the mesh deforms according to the computed velocity field. These Neumann boundary conditions allow material to flow both in and out through the boundary depending only on the Dirichlet boundary conditions and deformation that occurs inside the modelled domain.
3.4.4 Pressure Poisson problem in the 3D geodynamic model
In the context of our finiteelement forward model, we also solve the pressure Poisson problem using finite elements. As such, to compute P_{d} we employ the weak formulation given by Eq. (27) using boundary conditions from Eqs. (15) and (16) on bottom boundary and vertical boundaries, respectively. In our particular implementation, we employ Q_{1} for P_{d}, and these Q_{1} elements overlap the Q_{2} elements used to approximate the velocity.
As a demonstration of the computed P_{d} using this approach, in Fig. 10c and d we show the approximated pressure in our rift model at 8.7 Myr after large deformations that led to mantle exhumation and differential thinning of the continental crust, causing a variable topography. In this model, P_{d} was evaluated on a mesh consisting of $\mathrm{256}\times \mathrm{64}\times \mathrm{128}{Q}_{\mathrm{1}}$ finite elements on 1024 MPI ranks. The discrete pressure Poisson system was solved using geometric multigrid. As a rough estimate, solving for P_{d} required ∼0.2 % of the time required to solve the nonlinear viscous flow problem. Obviously this value is strongly dependent on both the physical model (linear viscous versus nonlinear viscous) and the implementation details, as well as the efficiency of how the discrete flow problem is solved. However, when considering even the simplest flow problem imaginable (i.e. linear isoviscous flow laws), it remains true that solving the Poisson problem will be far less expensive than solving either the linear or nonlinear viscous flow problem.
3.4.5 Tectonics evolution
The model using freeslip boundary conditions displays a cylindrical deformation pattern that could be reduced to a twodimensional model. As shown by the shear zone orientation and strain regime, the deformation is only extensional and perpendicular to the extension direction (Figs. 11a–d and 12). This strain localization is directly due to the freeslip boundary condition stating that any flow perpendicular to the boundary is prohibited.
In contrast, the model using the P_{d} pressure as a boundary condition displays a noncylindrical deformation. While extensional shear zones perpendicular to the extension direction develop in the central part of the domain, the edges of the rift experience oblique and strikeslip deformation (Fig. 11e to h). As the extension goes on, the extensional deformation localizes along a spreading centre, causing an increasing inflow on the boundaries of the domain with the normal stress boundary condition (Fig. 13c). As a result, near these boundaries the velocity field introduces noncylindrical features that are accommodated by strikeslip faults (Fig. 11g, h). These strikeslip faults delimit a triangular region terminating on a triple junction between two strikeslip faults and a ridge (ridge–fault–fault, RFF, triple junction). Along these strikeslip faults, the deformation is partitioned between purely vertical strikeslip shear zones and shallowdipping normal shear zones rooting into the strikeslip shear zones (Fig. 12).
4.1 Alternative PDEbased approaches
Recall that the starting point of defining the PPE was purely algebraic, with the sole intention of removing the nonuniqueness associated with Eq. (5). Here we discuss two alternative PDEbased approaches constructed with a similar rationale.
Rather than enforcing Eq. (15) only along the boundary, suppose we wished to enforce it everywhere throughout the domain,
This constraint can be interpreted as the steadystate solution of the following scalar hyperbolic PDE:
where τ plays the role of a timelike parameter having units of length and $\widehat{\mathit{g}}$ plays the role of a velocitylike quantity. Along the “inflow segments” ${\mathrm{\Gamma}}_{\mathrm{in}}=\mathit{\{}\mathit{x}\in \partial \mathrm{\Omega}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\widehat{\mathit{g}}\cdot \widehat{\mathit{n}}<\mathrm{0}\mathit{\}}$ we will impose P=0^{1}. On the outflow segments Γ_{out} = ∂Ω \ Γ_{in}; due to the hyperbolic nature of the PDE, no boundary constraint are required. Since we seek the steadystate solution of Eq. (42), no initial condition is required, but for completeness we chose $P(\mathit{\tau}=\mathrm{0})=\mathrm{0}$. Hence, the solution of Eq. (42) is equivalent to solving
along the family of characteristics given by
with $P(\mathit{\tau}=\mathrm{0})=\mathrm{0}$ and $\mathit{x}(\mathit{\tau}=\mathrm{0})\in {\mathrm{\Gamma}}_{\mathrm{in}}$.
Compared to the PPE, the hyperbolic formulation has several disadvantages.

We have less freedom to specify how ∇P varies along the boundary. The choice of boundary conditions (BCs) is largely dictated by the “inflow” and “outflow” segments. Along outflow segments, the only constraint available is Eq. (41) evaluated on ∂Ω. If inflow occurs on any part of ∂Ω not contained in ∂Ω_{surf}, we have to choose a flux BC, as using P=0 does not make physical sense. Consistency may require the use of a constraint that is independent of the PDE, e.g. Eq. (16).

The formulation may place restrictions on the shape of ∂Ω. If we wish to avoid the definition of new flux BCs (described above), the domain must be defined such that for every ${\mathit{x}}_{b}\in \partial \mathrm{\Omega}\backslash \partial {\mathrm{\Omega}}_{\mathrm{surf}}$ there exists a characteristic that intersects both x_{b} and ∂Ω_{surf}.

The lack of flexibility in controlling the boundary behaviour of ∇P will in general result in solutions of Eq. (42) being identical to a family of 1D solutions to Eq. (2) applied in directions parallel to the direction of gravity; see Fig. 14b, d.

The spatial discretization required for the accurate solution of Eq. (42) is arguably more complicated to implement (on unstructured meshes) compared with discretizations for the Poisson equation. Scalable multilevel solvers for the steadystate hyperbolic problem are much more challenging to develop in comparison to the Poisson problem.
From a linear algebra perspective, the nonuniqueness of Eq. (5) can alternatively be addressed by (i) discretizing Eq. (5) in space, yielding Gp=F, and then (ii) solving the normal equations
This approach obtains a unique solution that minimizes the following objective function:
In some senses, this approach is like the discrete counterpart of the PPE, in so much as G^{T}G are the discrete Laplacian defined on the approximation space used to represent pressure. The similarity between the pressure solution from the PPE and $\stackrel{\mathrm{\u0303}}{\mathbf{p}}$ are shown in Fig. 14a, c. Similar to the hyperbolic formulation, obtaining the pressure via the normal equations is restrictive from a modelling perspective as the formulation does not allow for control of the behaviour of ∇P on the boundary, and as such only the Dirichlet data on ∂Ω_{surf} can be specified. In Fig. 14a, c we provide snapshots of the pressure obtained from solving Eq. (45) in two different domains. Interestingly, the solutions indicate that $\mathrm{\nabla}P\cdot \widehat{\mathit{g}}$ appears to be approximately zero along the boundary despite the lack of a constraint enforcing this.
4.2 Flexibility of the PPE with respect to boundary conditions
Since the PPE is a secondorder PDE, the formulation permits a range of possible boundary constraints on ∇P to be imposed. Specific choices allow one to define whether (i) the equivalent of 1D pressure profiles as would be obtained by applying Eq. (2) along a boundary face or (ii) an approximation of the pressure field on the boundary would be obtained if a complete flow field was computed in a much larger global domain.
Experiments showed that the PPE approach better approximates the total pressure computed from the momentum equation. Therefore, its use as a boundary condition (or as an initial guess) for the pressure field to solve the momentum equation is preferred over hydrostatic solutions associated with Eq. (2). Moreover, as the domain size increases, the PPE formulation gives a more accurate approximation of the total pressure than the 1D depthintegrated approach.
4.3 Implications for lithosphere deformation
In the geodynamic rift model, using the pressure computed with Eq. (7) as a boundary condition produces a velocity field perpendicular to the extension direction in the rift axis (Fig. 13). This velocity field introduces noncylindrical deformation accommodated for by oblique and strikeslip structures (Fig. 11). The results of this study are very similar to previous studies directly applying an inflow perpendicular to the extension direction (Le Pourhiet et al., 2018; Jourdon et al., 2020). At the tip of the rift, a triangular region delimited by strikeslip faults or very oblique rift develops to accommodate the oblique velocity field. The similarity of these results shows that using open boundary conditions instead of kinematic boundary conditions in 3D may reveal firstorder implications for the lithosphere strain localization. In the case of geodynamic systems presenting the characteristics of a propagating rift (or ridge) with oblique and strikeslip deformation at its tip, considering the forces applied by the surrounding material weight could be the firstorder process at the origin of noncylindrical deformation.
4.4 Linear and nonlinear density
Two approaches can be considered to compute the pressure using Eq. (7). The first approach (also the simplest) is to consider that the density ρ is defined as a reference density ρ_{0} that only depends on rock type for a reference state, e.g. T_{0} and P_{0}. In that case, Eq. (7) is linear and only depends on the reference density structure. The second approach is to consider that the density ρ can vary with respect to other parameters. In our geodynamic example, we considered that the density can vary with pressure and temperature according to Eq. (38). In that case the Eq. (7) becomes nonlinear. This kind of nonlinearity is not rare in geodynamics, and the formulation of Eq. (38) may appear in the pressure computation when solving for an incompressible Boussinesq approximation or a compressible Stokes problem (e.g. Dannberg and Heister, 2016; King et al., 2010; Tackley, 2008).
In this study we presented a method to compute a reference pressure associated with the density structure of a domain in which we cast the problem in terms of a partial differential equation (PDE). From a practical standpoint, the PDE approach is generic (it is applicable to all spatial discretization and on any type of computational grid), efficient, and applicable in parallel computing environments. From the modelling perspective, the PDE approach has specific advantages, for example in models with a variable density structure (stationary or timedependent) and models that employ a reference pressure as a boundary condition of the flow problem (stationary or timedependent problems). Reevaluating that pressure in timedependent problems is not problematic (even if the mesh deforms) since solving the Poisson problem can be performed using optimal preconditioners (e.g. geometric or algebraic multigrid preconditioners). Importantly, the time to solve the pressure Poisson problem is a small fraction of the time required to solve the linear (or nonlinear) incompressible viscous flow problem. Moreover, we also demonstrate that the PDE formulation results in a better approximation of the total pressure than the 1D depthintegrated approach in nonhydrostatic cases.
Lastly, we showed in the context of 3D geodynamic models of continental rifting that using a reference pressure as a boundary condition within the flow problem resulted in noncylindrical velocity fields. These 3D velocity fields produced strain localization in the lithosphere along largescale strikeslip shear zones and the formation and evolution of triple junctions.
The code pTatin3D
used in this study to produce the 3D thermomechanical models is an opensource free software licensed under GPL3. The Supplement contains the version of the code used to produce the models presented in this study. To run the same models, users should use the driver named test_driver_checkpoint_fv.app
and the options files (.opts
) provided in the Supplement.
We also provide Firedrake
code (Firedrake team, 2022; Balay et al., 2019, 1997; Dalcin et al., 2011; Rathgeber et al., 2016) to compute the pressure Poisson problem in a halfannulus domain, in a deformed domain, and in the large and small domains with a topography perturbation used in this study. The version of Firedrake
used is 0.13.0+4944.g22178416
and is freely available.
We also provide FEniCS
code (FEniCS team, 2022; Alnaes et al., 2015; Logg et al., 2012b; Logg and Wells, 2010; Logg et al., 2012a) to reproduce the models solving for the normal and hyperbolic equations. The version of FEniCS
used is 2016.1.0
and is freely available.
The supplement related to this article is available online at: https://doi.org/10.5194/se1311072022supplement.
AJ and DAM conceptualized the project and made the mathematical and software development. AJ and DAM conceived, designed, and performed the experiments. AJ analysed the results. AJ and DAM drafted and corrected the manuscript. DAM procured funding and supervised the project.
The contact author has declared that neither they nor their coauthor has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gausscentre.eu, last access: 27 June 2022) for providing computing time on the GCS Supercomputer SuperMUCNG at the Leibniz Supercomputing Centre (https://www.lrz.de, last access: 27 June 2022) through project pr63qo. The authors also thank Cedric Thieulot and Rene Gassmoeller for their constructive reviews.
This project was supported by NSF Award EAR2121666. Anthony Jourdon has been supported by the European Union's Horizon 2020 research and innovation programme (TEAR ERC Starting, grant no. 852992).
This paper was edited by Taras Gerya and reviewed by Rene Gassmoeller and Cedric Thieulot.
Alnaes, M. S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M. E., and Wells, G. N.: The FEniCS Project Version 1.5, Archive of Numerical Software, 3, 9–23, https://doi.org/10.11588/ans.2015.100.20553, 2015. a
Baes, M., Sobolev, S. V., and Quinteros, J.: Subduction initiation in midocean induced by mantle suction flow, Geophys. J. Int., 215, 1515–1522, https://doi.org/10.1093/gji/ggy335, 2018. a
Balay, S., Gropp, W. D., McInnes, L. C., and Smith, B. F.: Efficient Management of Parallelism in Object Oriented Numerical Software Libraries, in: Modern Software Tools in Scientific Computing, edited by: Arge, E., Bruaset, A. M., and Langtangen, H. P., 163–202, Birkhäuser Press, https://doi.org/10.1007/9781461219866_8, 1997. a
Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D., Knepley, M. G., May, D. A., McInnes, L. C., Mills, R. T., Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H.: PETSc Users Manual, Tech. Rep. ANL95/11 – Revision 3.11, Argonne National Laboratory, 2019. a
Barth, W. L. and Carey, G. F.: On a boundary condition for pressuredriven laminar flow of incompressible fluids, Int. J. Numer. Meth. Fl., 54, 1313–1325, 2007. a
Brune, S.: Evolution of stress and fault patterns in oblique rift systems: 3D numerical lithosphericscale experiments from rift to breakup, Geochem. Geophy. Geosy., 15, 3392–3415, https://doi.org/10.1002/2014GC005446, 2014. a
Brune, S., Popov, A. A., and Sobolev, S. V.: Modeling suggests that oblique extension facilitates rifting and continental breakup, J. Geophys. Res., 117, 1–16, https://doi.org/10.1029/2011JB008860, 2012. a
Brune, S., Heine, C., PérezGussinyé, M., and Sobolev, S. V.: Rift migration explains continental margin asymmetry and crustal hyperextension, Nat. Commun., 5, 1–9, https://doi.org/10.1038/ncomms5014, 2014. a
Brune, S., Heine, C., Clift, P. D., and PérezGussinyé, M.: Rifted margin architecture and crustal rheology: Reviewing IberiaNewfoundland, Central South Atlantic, and South China Sea, Mar. Petrol. Geol., 79, 257–281, https://doi.org/10.1016/j.marpetgeo.2016.10.018, 2017. a
Chertova, M. V., Geenen, T., van den Berg, A., and Spakman, W.: Using open sidewalls for modelling selfconsistent lithosphere subduction dynamics, Solid Earth, 3, 313–326, https://doi.org/10.5194/se33132012, 2012. a
Chertova, M. V., Spakman, W., Geenen, T., Van Den Berg, A., and van Hinsbergen, D. J. J.: Underpinning tectonic reconstructions of the western Mediterranean region with dynamic slab evolution from 3D numerical modeling, J. Geophys. Res.Sol. Ea., 119, 5876–5902, https://doi.org/10.1002/2014JB011150, 2014. a
Dalcin, L. D., Paz, R. R., Kler, P. A., and Cosimo, A.: Parallel distributed computing using Python, new Computational Methods and Software Tools, Adv. Water Resour., 34, 1124–1139, https://doi.org/10.1016/j.advwatres.2011.04.013, 2011. a
Dannberg, J. and Heister, T.: Compressible magma/mantle dynamics: 3D, adaptive simulations in ASPECT, Geophys. J. Int., 207, 1343–1366, https://doi.org/10.1093/gji/ggw329, 2016. a
Elman, H. C., Silvester, D. J., and Wathen, A. J.: Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics, Numer. Math. Sci. Comp., Oxford Scholarship Online, https://doi.org/10.1093/acprof:oso/9780199678792.001.0001, 2014. a
FEniCS team: FEniCS, https://fenicsproject.org/download/, last access: 27 June 2022. a
Firedrake team: Firedrake, https://www.firedrakeproject.org/download.html, last access: 27 June 2022. a
Glerum, A., Thieulot, C., Fraters, M., Blom, C., and Spakman, W.: Nonlinear viscoplasticity in ASPECT: benchmarking and applications to subduction, Solid Earth, 9, 267–294, https://doi.org/10.5194/se92672018, 2018. a
Hirth, G. and Kohlstedt, D. L.: Rheology of the Upper Mantle and the Mantle Wedge: A View from the Experimentalists, Geophys. Monogr., 138, 83–105, 2003. a
IsmailZadeh, A., Honda, S., and Tsepelev, I.: Linking mantle upwelling with the lithosphere decent and the Japan Sea evolution: A hypothesis, Sci. Rep.UK, 3, 1–7, https://doi.org/10.1038/srep01137, 2013. a
Jourdon, A., Le Pourhiet, L., Mouthereau, F., and May, D.: Modes of Propagation of Continental Breakup and Associated Oblique Rift Structures, J. Geophys. Res.Sol. Ea., 125, 1–27, https://doi.org/10.1029/2020JB019906, 2020. a
King, S. D., Lee, C., Van Keken, P. E., Leng, W., Zhong, S., Tan, E., Tosi, N., and Kameyama, M. C.: A community benchmark for 2D Cartesian compressible convection in the Earth's mantle, Geophys. J. Int., 180, 73–87, https://doi.org/10.1111/j.1365246X.2009.04413.x, 2010. a
Le Pourhiet, L., ChamotRooke, N., Delescluse, M., May, D. A., Watremez, L., and Pubellier, M.: Continental breakup of the South China Sea stalled by farfield compression, Nat. Geosci., 11, 605–609 https://doi.org/10.1038/s4156101801785, 2018. a
Logg, A. and Wells, G. N.: DOLFIN: Automated Finite Element Computing, ACM T. Math. Software, 37, 1–8, https://doi.org/10.1145/1731022.1731030, 2010. a
Logg, A., Mardal, K., and Wells, G. N. (Eds.): Automated Solution of Differential Equations by the Finite Element Method, Springer, https://doi.org/10.1007/9783642230998, 2012a. a
Logg, A., Wells, G. N., and Hake, J.: DOLFIN: a C++/Python Finite Element Library, in: Automated Solution of Differential Equations by the Finite Element Method, edited by: Logg, A., Mardal, K., and Wells, G. N., Vol. 84 of Lecture Notes in Computational Science and Engineering, chap. 10, Springer, https://doi.org/10.1007/9783642230998, 2012b. a
May, D. A., Brown, J., and Le Pourhiet, L.: pTatin3D : HighPerformance Methods for LongTerm Lithospheric Dynamics, Proceeding SC′14 Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, IEEE (Institute of Electrical and Electronics Engineers, 274–284, https://doi.org/10.1109/SC.2014.28, 2014. a
May, D. A., Brown, J., and Le Pourhiet, L.: A scalable, matrixfree multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow, Comput. Method. Appl. M., 290, 496–523, https://doi.org/10.1016/j.cma.2015.03.014, 2015. a
Popov, A. A. and Sobolev, S. V.: SLIM3D: A tool for threedimensional thermomechanical modeling of lithospheric deformation with elastoviscoplastic rheology, Phys. Earth Planet. In., 171, 55–75, https://doi.org/10.1016/j.pepi.2008.03.007, 2008. a
Quinteros, J., Sobolev, S. V., and Popov, A. A.: Viscosity in transition zone and lower mantle: Implications for slab penetration, Geophys. Res. Lett., 37, L09307, https://doi.org/10.1029/2010gl043140, 2010. a
Ranalli, G.: Rheology of the lithosphere in space and time, Geol. Soc. Lond. Spec. Publ., 121, 19–37, https://doi.org/10.1144/GSL.SP.1997.121.01.02, 1997. a
Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., McRae, A. T. T., Bercea, G.T., Markall, G. R., and Kelly, P. H. J.: Firedrake: automating the finite element method by composing abstractions, ACM T. Math. Softw., 43, 24:1–24:27, https://doi.org/10.1145/2998441, 2016. a
Rybacki, E. and Dresen, G.: Dislocation and diffusion creep of synthetic anorthite aggregates, J. Geophys. Res.Sol. Ea., 105, 26017–26036, https://doi.org/10.1029/2000JB900223, 2000. a
Tackley, P. J.: Modelling compressible mantle convection with large viscosity contrasts in a threedimensional spherical shell using the yinyang grid, Phys. Earth Planet. In., 171, 7–18, https://doi.org/10.1016/j.pepi.2008.08.005, 2008. a
Yamato, P., Burov, E., Agard, P., Le Pourhiet, L., and Jolivet, L.: HPUHP exhumation during slow continental subduction: Selfconsistent thermodynamically and thermomechanically coupled model with application to the Western Alps, Earth Planet. Sc. Lett., 271, 63–74, https://doi.org/10.1016/j.epsl.2008.03.049, 2008. a