An efﬁcient partial-differential-equation-based method to compute pressure boundary conditions in regional geodynamic models

. 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 ﬂuid by using a hydrostatic scenario or a uniformly mov-ing ﬂuid 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 non-hydrostatic cases, we show that this PDE-based approach better approximates the total pressure than the classical 1D depth-integrated approach. To illustrate the performance of this PDE-based formulation we present several hydrostatic and non-hydrostatic 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 time-dependent boundary condition to simulate far-ﬁeld normal stresses. This model shows a high degree of non-cylindrical deformation, resulting from the stress boundary condition, that is accommodated by strike-slip shear zones. We compare the result of this numerical model with a traditional rift model employing free-slip boundary conditions to demonstrate the ﬁrst-order implications of considering “open” boundary conditions in 3D thermo-mechanical rift models.


Introduction
In Earth sciences and geodynamic modelling, computing the pressure can be essential. Specifically, numerous regional thermo-mechanical 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., 2012Brune et al., , 2014Brune et al., , 2017Chertova et al., 2012Chertova et al., , 2014Glerum et al., 2018;Ismail-Zadeh 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 non-linear 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 x s 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 x s . For the case of a constant density and gravity, this expression reduces to P (x ) = P s + ρgD, where D is the distance (depth) given by D = x s − x || and g = ||g||. When the density is a function Published by Copernicus Publications on behalf of the European Geosciences Union.

1108
A. Jourdon and D. A. May: An efficient PDE-based method to compute pressure boundary conditions 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 sub-division 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 0 = P s +ρ 0 gD 0 , P 1 = P s +ρ 0 gD 0 +ρ 1 gD 1 = P 0 + ρ 1 gD 1 , . . . , P N = P s + N i ρ i gD i = P N−1 + ρ N gD 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 column-wise integration difficult (or expensive): 1. a mesh with cell edges (2D) or faces (3D) that are not aligned with the gravity vector ( Fig. 1a), 2. an unstructured mesh (Fig. 1b), 3. a density structure (or gravity vector) that is spatially varying, 4. a parallel decomposition of the mesh (Fig. 1c), 5. time dependence in the density or mesh coordinates that requires continual re-evaluation of the reference pressure.
To compute P (x ) we first have to define the location x s . In general this is non-trivial 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 x s . 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 gravity-aligned mesh is not possible. If the density (or gravity) vary in space throughout the domain, the integral must be approximated via a suitable sub-division 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 [x , x s ] and each cell and apply a one-point quadrature rule over the intersecting segment times. In Fig. 1a and b we depict the complexity of this procedure for a noncoordinate-aligned 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 column-wise 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 re-compute 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 non-hydrostatic 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 non-hydrostatic 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 thermo-mechanical numerical models and static numerical models applied to Earth sciences and geodynamics to show the usefulness of this approach.

PDE-based pressure formulation
For an incompressible fluid in a domain , the non-inertial 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, τ = 2ηε(v) is the deviatoric stress tensor with η the viscosity,ε(v) is the strain rate tensor, ρ := ρ(x, t) is the density, g := g(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 non-overlapping segments: ∂ surf that we will regard as the free surface and prescribe that tangential and normal stresses are zero, i.e. τ −P I = 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, ∂ = ∂ i ∪ ∂ surf and ∂ i ∩ ∂ surf = ∅. The outward-pointing unit normal vector to ∂ will be denoted vian.
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 d = 2, 3, Eq. (5) is over-determined 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 left-hand side is small, we have which must be true for any arbitrary domain, thus resulting in Eq. (7).

1110
A. Jourdon and D. A. May: An efficient PDE-based method to compute pressure boundary conditions Assuming thatn is constant or slowly varying and that τ is symmetric yields Hence, dropping the first term in Eq. (9) is equivalent to saying that τn ≈ 0 for all x ∈ ∂ . Alternatively, the term is zero if ∇ · τ = 0 for all x ∈ ∂ . This condition is satisfied if the fluid experiences either rigid body translations or rigid body rotation along ∂ .

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 Earth-like 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 ∇P ·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), namelyn · (τ − P I)n =n · (τ − P I)t = 0, which reduces tô n·τn−P = 0 andn·τt = 0 witht a tangent unit vector to the boundary such thatn ·t = 0. Equation (12) is consistent with a fluid at rest since τ = 0. In the non-hydrostatic case, we require τn ≈ 0 to arrive at Eq. (12). We also note thatn · τ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 ∇P ·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 ∂ i = ∂ ⊥ ∪ ∂ . Next, we define the gravity unit vectorĝ such that and the unit vectorĝ ⊥ , which is perpendicular toĝ, i.e.
The first constraint states that P should increase only in the direction of gravity. Hence, from Eq. (5) we have ∇P ·ĝ = ρg ĝ 2 = ρg, for all x ∈ ∂ .
The second constraint states that P should not change along directions perpendicular to the gravity; hence, Since the unit vectorn normal to the boundary ∂ can be decomposed according tô we have ∇P ·n = (n ·ĝ)∇P ·ĝ + (n ·ĝ ⊥ )∇P ·ĝ ⊥ .
Hence, the two Neumann boundary conditions may be stated as and ∇P ·n = (n ·ĝ)∇P ·ĝ for all x ∈ ∂ ⊥ .
Equations (19) and (20) may appear peculiar since both the left-hand side and right-hand 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 normaln, and our 1D inspired gradients do exactly that.
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.

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 ∈ H d 1 ( ), 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 right-hand 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 ∂ i = ∂ ⊥ ∪ ∂ and using Eqs.
Noting that the second term of the left-hand side and part of the last term on the right-hand 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

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 re-use existing discretization implementations that support domain decomposition.

Numerical examples
In this section we provide several numerical models to show the following items: 1. the accuracy and consistency of the method for hydrostatic cases, 2. the accuracy of the approximation of the total pressure in non-hydrostatic cases, 3. the effect of using the depth-integrated approach and the pressure Poisson problem approach to impose boundary conditions on the momentum equation, 4. the usefulness of the method for 3D geodynamic thermo-mechanical modelling.

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).
Case 1. This case assumes a constant density, ρ(y) = 1 (Fig. 2a, b). The analytic pressure solution is given by Case 2. This case assumes a continuous depth-varying density ρ(y) = 2 − y (Fig. 2c, d). The analytic pressure solution is given by Case 3. This case assumes a discontinuous density such that ρ(y) = 1 for y ∈ [0.5, 1] and ρ(y) = 2 for y ∈ [0, 0.5) (Fig. 2e, f). The analytic pressure solution is given by with D the distance between the surface and the y coordinate at which the pressure is computed.
A finite-element (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. point-wise), 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. sub-dividing 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.

Half-annulus domain
The 2D half-annulus 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 = r sin(θ ) and y = r cos(θ ). The gravity vector pointing to the centre (x = (0, 0)) is defined as g = −9.8 x √ x 2 +y 2 , y √ x 2 +y 2 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ĝ. At the surface of the domain we impose the Dirichlet constraint P = 0. As in Sect. 3.1 we use a finite-element 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 half-annulus model (Fig. 3c). The difference between the pressure obtained using the depth-integrated approach and the PPE is very small (Fig. 3d). Unlike the analytic solution for the box models, the analytic solution for P here is non-polynomial (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 depth-integrated approach and the pressure Poisson equation (with one set of boundary constraints) are equivalent for scenarios that admit a hydrostatic solution.

Non-hydrostatic 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.

"Global" model
We define a large domain G = x ∈ [−10, 10] × y ∈ [0, 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 s (x) = − 1 4 x, where y s is the surface between x = 0 and x = 1. For demonstra- tion 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 g = (0, −1).
We solve the flow problem described by Eqs.
(3), (4) using a constant viscosity and density, along with the following boundary conditions: no-slip at base, free-slip on the right and left faces, and a free surface along the top of the model domain. Due to the topography, a non-vertical 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 non-perturbed 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 ex-tremely small (Fig. 4b). In contrast, the difference between the total pressure and depth-integrated 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.

"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 first-order 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 R = x ∈ [0, 1] × y ∈ [0, 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 α = {P d , P a } denote the pressure computed using the pressure Poisson equation and the 1D depth-integrated 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 no-slip 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, ∇P d ·n = (n ·ĝ)∇P d ·ĝ on vertical sides, and ∇P d ·n = (n ·ĝ)ρg + (n ·ĝ ⊥ )∇P d ·ĝ ⊥ 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 co-linear 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 5 × 10 −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 depth-integrated 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 bottom-left 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 depth-integrated approach due to its 1D nature. Figure 6. Curves showing the pressure from the regional and global models along the profiles (a) x = 0 and (c) x = 1. Absolute value of the difference between the total pressure P in the global model and the different pressures in the global and regional models along the profiles (b) x = 0 and (d) x = 1.

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 g = (0, −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 ∂ i = 1 ∪ 2 ∪ 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. ∂ x P d = 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 1,2,3 , Fig. 7b segments   1,2 , Fig. 7c segment 2 , and Fig. 7d segment 3 ).  highlights that the solution of the PPE can be identical to that obtained using the depth-integrated 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.

Physical model
To simulate the long-term evolution of the deformation of the lithosphere, we solve the stationary, non-inertial 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 finite-element 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 T (v, p) = (τ (u) − pI)n, andn 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 finite-element code pTatin3D (May et al., 2014(May et al., , 2015, which employs a mixed Q 2 -P 1 discretization for velocity and pressure.

Initial conditions and rheology
To model the strain localization we use non-linear 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 material-defined parameters (see Table 1); R is the universal gas constant; V is the activation volume; andε II is the square root of the strain rate second invariant, computed aṡ 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 = −25 km and is modelled with a dislocation creep quartz rheology (Ranalli, 1997). The lower crust extends from y < −25 km to y = −35 km and is modelled with a dislocation creep anorthite rheology (Rybacki and Dresen, 2000). The lithosphere mantle extends from y < −35 km to y = −120 km, whereas the asthenosphere mantle extends from y < −120 km to y = −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 steady-state solution of the heat equation using a surface temperature of T = 0 • C at y = 0 km and T = 1450 • C at y = −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.

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 free-slip 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) andn 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 non-linear iteration for each time step, and the Neumann boundary condition described by Eq. (40) is evaluated at every non-linear 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 non-linearity 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.

Pressure Poisson problem in the 3D geodynamic model
In the context of our finite-element 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.
Pa −1 10 −11 10 −11 10 −11 10 −11 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 256 × 64 × 128Q 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 non-linear viscous flow problem. Obviously this value is strongly dependent on both the physical model (linear viscous versus non-linear 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 non-linear viscous flow problem.

Tectonics evolution
The model using free-slip boundary conditions displays a cylindrical deformation pattern that could be reduced to a two-dimensional model. As shown by the shear zone orientation and strain regime, the deformation is only extensional and perpendicular to the extension direction . This strain localization is directly due to the free-slip 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 non-cylindrical 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 strike-slip 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 non-cylindrical features that are accommodated by strike-slip faults (Fig. 11g, h). These strike-slip faults delimit a triangular region terminating on a triple junction between two strike-slip faults and a ridge (ridge-fault-fault, RFF, triple junction). Along these strike-slip faults, the deformation is partitioned between purely vertical strike-slip shear zones and shallowdipping normal shear zones rooting into the strike-slip shear zones (Fig. 12).

Alternative PDE-based 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 PDE-based 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 steady-state solution of the following scalar hyperbolic PDE: where τ plays the role of a time-like parameter having units of length andĝ plays the role of a velocity-like quantity. Along the "inflow segments" in = {x ∈ ∂ :ĝ ·n < 0} we will impose P = 0 1 . On the outflow segments out = along the family of characteristics given by with P (τ = 0) = 0 and x(τ = 0) ∈ in . Compared to the PPE, the hyperbolic formulation has several disadvantages.
1. 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).
2. 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 x b ∈ ∂ \∂ surf there exists a characteristic that intersects both x b and ∂ surf .
3. 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.
4. 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 steady-state hyperbolic problem are much more challenging to develop in comparison to the Poisson problem.
From a linear algebra perspective, the non-uniqueness 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 andp 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 pres-sure obtained from solving Eq. (45) in two different domains. Interestingly, the solutions indicate that ∇P ·ĝ appears to be approximately zero along the boundary despite the lack of a constraint enforcing this.

Flexibility of the PPE with respect to boundary conditions
Since the PPE is a second-order 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 depth-integrated approach.

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 non-cylindrical deformation accommodated for by oblique and strike-slip structures (Fig. 11). The results of this study are very similar to previous studies directly applying an inflow perpen-dicular to the extension direction (Le Pourhiet et al., 2018;Jourdon et al., 2020). At the tip of the rift, a triangular region delimited by strike-slip 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 first-order 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 first-order process at the origin of non-cylindrical deformation.

Linear and non-linear 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 non-linear. This kind of non-linearity 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).

Conclusions
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 time-dependent) and models that employ a reference pressure as a boundary condition of the flow problem (stationary or time-dependent problems). Re-evaluating that pressure in time-dependent 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 non-hydrostatic 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 large-scale strikeslip shear zones and the formation and evolution of triple junctions.
Code availability. The code pTatin3D used in this study to produce the 3D thermo-mechanical models is an open-source 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 Figure 14. Pressure in non-dimensional models with constant density ρ = 1. Pressure computed for a hydrostatic case in an irregularly shaped domain using (a) the normal equations (Eq. 45) and (b) the hyperbolic equation (Eq. 42). Pressure computed for a nonhydrostatic case due to topography using (c) the normal equation (Eq. 45) and (d) the hyperbolic equation (Eq. 42). The contour lines show contours of iso-pressure values every 0.1. 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., 2019Balay et al., , 1997Dalcin et al., 2011;Rathgeber et al., 2016) to compute the pressure Poisson problem in a half-annulus 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.
Author contributions. 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.
Competing interests. The contact author has declared that neither they nor their co-author has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.