Technical note: Finite element formulations to map discrete fracture elements in three-dimensional groundwater models

Typically, in finite element groundwater models, fractures are represented by two-dimensional triangular or quadrilateral elements. When embedded in a three-dimensional space, the Jacobian matrix governing the transformation from the global three-dimensional space to the local two-dimensional space is rectangular and thus not invertible. There exist different approaches to obtain a unique mapping from local to global space even though the Jacobian matrix is not invertible. These approaches are discussed in this study. It is illustrated that all approaches yield the same result and may be applied to curved elements. The mapping of anisotropic hydraulic conductivity tensors for possibly curved fracture elements is also discussed.

1 Introduction

The finite element method is well suited for accommodating fractures in groundwater models. Typically, fractures are represented by discrete two-dimensional elements, and these fracture elements can be embedded within a three-dimensional continuum consisting of three-dimensional elements. For example, within a tetrahedral mesh, fractures can be embedded by using triangular elements such that each triangle corresponds to a face shared by two adjacent tetrahedral elements. Similarly, quadrilaterals can be embedded within a hexahedral mesh. Indeed, such discrete-continuum models with embedded fractures are routinely applied (Blessent et al., 2009, 2011; Li et al., 2020; Watanabe, 2011).

A key component in the finite element method is the mapping of the gradient matrix from local to global space, where the global space is typically defined by a standard orthogonal coordinate system. The local space within a finite element can be curvilinear and has the same dimension as the element itself. If the global space has the same dimension as the local space, then the mapping is defined by the inverse of the Jacobian matrix. However, in the case of two-dimensional fractures embedded in a global three-dimensional space, the Jacobian matrix is non-square and thus not invertible (Juanes et al., 2002; Perrochet, 1995). A couple of different techniques enable a mapping from two-dimensional local to three-dimensional global space.

A first approach is based on using contravariant base vectors and the contravariant metric tensor (Kiraly, 1985; Cornaton et al., 2004; Juanes et al., 2002; Perrochet, 1995). This approach requires some understanding of tensor calculus, and the few studies that describe this approach refer to mathematical textbooks for more details. Nonetheless, this approach yields a rather simple expression for the mapping and is directly applicable to curved elements.

A second approach uses the right Penrose–Moore inverse of the Jacobian matrix. As shown in this study, the derivation of this pseudo-inverse is relatively straightforward. Within the field of finite elements, the left Penrose–Moore inverse has been applied for the reverse mapping from a three-dimensional global space to a two-dimensional local space (Rognes et al., 2013). One study mentions the pseudo-inverse for mapping finite elements to higher dimensions (Reichenberger et al., 2004) but only within the context of non-curved elements and without much further detail.

A third approach is to introduce an intermediate mapping to an orthonormal two-dimensional space tangent to the fracture space. The Jacobian of such a mapping is invertible. A matrix of directional cosines is used for a subsequent mapping to the global space. This approach is widely used, and the available literature is quite detailed (Diersch et al., 2005; Watanabe, 2011; Kolditz, 2002). However, the approach as discussed in the available literature is only applicable to non-curved finite elements.

The existence of multiple approaches, which are quite different from a mathematical point of view, makes it difficult to navigate the literature for those in need of implementing the mapping of a gradient matrix to higher dimensions. This study provides a comprehensive discussion of the three approaches. It is shown that all approaches yield the exact same result. It is illustrated that the third approach can be applied to curved elements by a minor adjustment. Although rarely discussed, this study highlights that the right Penrose–Moore inverse is an elegant alternative approach to find the gradient matrix in global coordinates.

The mapping of locally defined hydraulic conductivity tensors to the global space is also discussed. Although this mapping is discussed in the existing literature for non-curved elements (Kolditz, 2002), here a more general mapping is presented that is also applicable to curved fracture elements. This is useful, as such a mapping for curved elements is not discussed in the existing literature.

2 Preliminary on the geometry of a fracture finite element

Figure 1 illustrates a curved quadrilateral fracture finite element. The orientation of the fracture element can be defined by the normal, strike and dip directions. The local space within the curved quadrilateral is defined by local coordinates s^k with $-\mathrm{1}\le s^k\le \mathrm{1}$. To describe this curved space, some differential geometry of surfaces is needed (Farrashkhalvat and Miles, 2003; Nguyen-Schäfer and Schmidt, 2014; Lebedev et al., 2010; Itskov, 2007). The covariant base vectors are tangent to the local coordinate axes and are given by

$$\begin{array}{}\text{(1)}& {\mathit{a}}_k={\displaystyle \frac{\partial x^j}{\partial s^k}}{\mathit{e}}_j.\end{array}$$

The contravariant base vectors (a^k) are perpendicular to planes along which s^k varies and are given by

$$\begin{array}{}\text{(2)}& {\mathit{a}}^k={\displaystyle \frac{\partial s^k}{\partial x^i}}{\mathit{e}}_i,\end{array}$$

such that

$$\begin{array}{}\text{(3)}& {\mathit{a}}^j\cdot {\mathit{a}}_i={\mathit{\delta }}_j^i,\end{array}$$

where ${\mathit{\delta }}_j^i$ is the Kronecker delta symbol. The contravariant base vectors and the covariant base vectors are related by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle }{\mathit{a}}_i& {\displaystyle }=G_{ij}{\mathit{a}}^j,\\ {\displaystyle }{\mathit{a}}^i& {\displaystyle }=H^{ij}{\mathit{a}}_j,\end{array}\end{array}$$

where G_ij and H^ij are the covariant and contravariant metric tensors, respectively. These tensors are given by

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle }{G}_{ij}& {\displaystyle }={\mathit{a}}_i\cdot {\mathit{a}}_j,\\ {\displaystyle }{H}^{ij}& {\displaystyle }={\mathit{a}}^i\cdot {\mathit{a}}^j=G_{ij}^{-\mathrm{1}}.\end{array}\end{array}$$

The unit normal vector is simply defined by the cross product of the covariant base vectors:

$$\begin{array}{}\text{(6)}& \mathit{n}={\displaystyle \frac{{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}}{\left|{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}\right|}}.\end{array}$$

Making use of Lagrange's identity, the area spanned by the covariant base vectors can be shown to equal the square root of the determinant of G:

$$\begin{array}{}\text{(7)}& \left|{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}\right|=\sqrt{({\mathit{a}}_{\mathrm{1}}\cdot {\mathit{a}}_{\mathrm{1}})({\mathit{a}}_{\mathrm{2}}\cdot {\mathit{a}}_{\mathrm{2}})-({\mathit{a}}_{\mathrm{1}}\cdot {\mathit{a}}_{\mathrm{2}}{)}^{\mathrm{2}}}=\sqrt{\text{det}\mathbf{G}}.\end{array}$$

The local two-dimensional space can be expanded to a local three-dimensional space with the following base vectors all normal to the fracture surface:

$$\begin{array}{}\text{(8)}& {\mathit{a}}_{\mathrm{3}}={\mathit{a}}^{\mathrm{3}}=\mathit{n}.\end{array}$$

Then Eq. (3) implies that the contravariant base vectors can also be expressed as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle }{\mathit{a}}^{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\sqrt{g}}}({\mathit{a}}_{\mathrm{2}}\times {\mathit{a}}_{\mathrm{3}}),\\ {\displaystyle }{\mathit{a}}^{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\sqrt{g}}}({\mathit{a}}_{\mathrm{3}}\times {\mathit{a}}_{\mathrm{1}}),\\ {\displaystyle }{\mathit{a}}^{\mathrm{3}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\sqrt{g}}}({\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}),\end{array}\end{array}$$

where g=det G. It is noted that covariant and contravariant base vectors as well as metric tensors can similarly be defined for triangular finite elements.

Figure 1 Geometry of a curved fracture element.

3 The basic mapping problem

The finite element formulations for groundwater flow result in element matrices that require the element shape functions and their partial derivatives with respect to global Cartesian coordinates. These matrices also involve an integration over the finite element domain Ω_e. For the objective of this study, it suffices to consider the element conductance matrix for saturated groundwater flow:

$$\begin{array}{}\text{(10)}& \mathbf{G}={\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}\mathrm{\nabla }\mathbf{NK}\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}d{\mathrm{\Omega }}_{\mathrm{e}},\end{array}$$

where K is the hydraulic conductivity tensor defined with respect to a global Cartesian coordinate system xⁱ, and ∇N is the gradient matrix often denoted by B (Perrochet, 1995):

$$\begin{array}{}\text{(11)}& B_{ni}={\displaystyle \frac{\partial N_n}{\partial x^i}},\end{array}$$

where N_n is the nth nodal shape function. Typically, however, the shape functions are provided with respect to a local coordinate system s^k. To find the partial derivatives of the shape functions with respect to global coordinates, the standard approach is to use the Jacobian matrix of the coordinate transformation between local and global space. Following the chain rule,

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial N_n}{\partial s^k}}={\displaystyle \frac{\partial N_n}{\partial x^i}}{\displaystyle \frac{\partial x^i}{\partial s^k}}.\end{array}$$

The Jacobian is defined as follows:

$$\begin{array}{}\text{(13)}& J_{ki}={\displaystyle \frac{\partial x^i}{\partial s^k}}.\end{array}$$

The components of the Jacobian are computed using the derivatives of the shape functions with respect to local coordinates and the nodal coordinates:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial x^i}{\partial s^k}}={\displaystyle \frac{\partial N_n}{\partial s^k}}{x}_n^i.\end{array}$$

It can be observed that the Jacobian contains the covariant base vectors per row, and Eq. (14) illustrates how to compute these vectors from local shape functions and nodal coordinates. If the Jacobian is invertible, then the derivatives with respect to global coordinates can be computed as follows:

$$\begin{array}{}\text{(15)}& {\mathbf{B}}^{\mathrm{T}}={\mathbf{J}}^{-\mathrm{1}}{\mathbf{B}}^{\ast \mathrm{T}},\end{array}$$

where B^∗ denoted the gradient matrix with respect to local coordinates:

$$\begin{array}{}\text{(16)}& B_{nk}^{\ast }={\displaystyle \frac{\partial N_n}{\partial s^k}}.\end{array}$$

Once B^T has been computed, the matrix B can be computed easily by taking the transpose of B^T. Introducing the coordinate matrix X containing the nodal coordinates per row, it follows from Eq. (14) that the Jacobian can be computed using

$$\begin{array}{}\text{(17)}& \mathbf{J}={\mathbf{B}}^{\ast \mathrm{T}}\mathbf{X}.\end{array}$$

Typically, the element matrices are computed using Gaussian quadrature, although for a limited number of element types the integration can be carried out analytically (Diersch, 2013). The advantage of numerical integration is that it can be applied to any element type, including curved elements. To perform Gaussian quadrature, the integration limits need to be defined with respect to the local domain dΩ^∗. If the Jacobian is invertible, then (Perrochet, 1995)

$$\begin{array}{}\text{(18)}& d\mathrm{\Omega }=\text{det}\left(\mathbf{J}\right)d{\mathrm{\Omega }}^{\ast }.\end{array}$$

However, if the Jacobian is not a square matrix, then the Jacobian matrix it is not invertible, and Eqs. (15) and (18) cannot be used for the finite element computations. This occurs when the local space has a lower dimension than the global space. Thus, for two-dimensional fracture elements embedded within a three-dimensional model space, the problem is that the Jacobian is not a square matrix.

In Eq. (10), the hydraulic conductivity tensor for fractures is to be defined with respect to the global Cartesian space. In general, it is more convenient to start with tensors which are defined with respect to the strike and dip directions along a fracture. The strike, dip and normal directions provide a locally orthogonal coordinate system. On curvilinear elements, this local coordinate system varies from point to point.

4 Gradient mapping using contravariant and covariant bases

Similar to Eq. (12), it follows from the chain rule that

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial N_n}{\partial x^i}}={\displaystyle \frac{\partial N_n}{\partial s^k}}{\displaystyle \frac{\partial s^k}{\partial x^i}}.\end{array}$$

This indicates that the gradient matrix with respect to global coordinates can be obtained using the contravariant base vectors. Introducing a matrix D in which the columns contain the contravariant base vectors,

$$\begin{array}{}\text{(20)}& D_{ik}={\displaystyle \frac{\partial s^k}{\partial x^i}},\end{array}$$

and it follows that

$$\begin{array}{}\text{(21)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}=\mathbf{D}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}.\end{array}$$

The components in matrix D can be rewritten in terms of covariant vectors using Eq. (5):

$$\begin{array}{}\text{(22)}& D_{ij}=({\mathit{a}}^j{)}_i=(H^{jk}{\mathit{a}}_k{)}_i.\end{array}$$

Since the Jacobian J contains the covariant vectors per row, this can be written as

$$\begin{array}{}\text{(23)}& \mathbf{D}={\mathbf{J}}^{\mathrm{T}}\mathbf{H}.\end{array}$$

The contravariant metric tensor H can also be written in terms of the Jacobian matrices:

$$\begin{array}{}\text{(24)}& \mathbf{H}={\mathbf{G}}^{-\mathrm{1}}=({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}},\end{array}$$

where it is noted that JJ^T is an invertible square matrix. Thus, the gradient matrix in global coordinates is given by

$$\begin{array}{}\text{(25)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}.\end{array}$$

Once the Jacobian is available, Eq. (25) provides a straightforward solution for the gradient matrix in global coordinates.

The differential volume follows from Eq. (7):

$$\begin{array}{}\text{(26)}& d\mathrm{\Omega }=\sqrt{\text{det}\left(\mathbf{G}\right)}d{\mathrm{\Omega }}^{\ast }=\sqrt{\text{det}\left({\mathbf{JJ}}^{\mathrm{T}}\right)}d{\mathrm{\Omega }}^{\ast }\end{array}$$

It is interesting to observe that Eq. (9) permits one to write the matrix D as

$$\begin{array}{}\text{(27)}& \mathbf{D}={\displaystyle \frac{\mathrm{1}}{\sqrt{g}}}\left[\begin{array}{cc}({\mathit{a}}_{\mathrm{2}}\times {\mathit{a}}_{\mathrm{3}}{)}_{\mathrm{1}}& ({\mathit{a}}_{\mathrm{3}}\times {\mathit{a}}_{\mathrm{1}}{)}_{\mathrm{1}}\\ ({\mathit{a}}_{\mathrm{2}}\times {\mathit{a}}_{\mathrm{3}}{)}_{\mathrm{2}}& ({\mathit{a}}_{\mathrm{3}}\times {\mathit{a}}_{\mathrm{1}}{)}_{\mathrm{2}}\\ ({\mathit{a}}_{\mathrm{2}}\times {\mathit{a}}_{\mathrm{3}}{)}_{\mathrm{3}}& ({\mathit{a}}_{\mathrm{3}}\times {\mathit{a}}_{\mathrm{1}}{)}_{\mathrm{3}}\end{array}\right].\end{array}$$

Using the vector triple product, it can be shown that

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\displaystyle }{\mathit{a}}_{\mathrm{2}}\times {\mathit{a}}_{\mathrm{3}}& {\displaystyle }={\mathit{a}}_{\mathrm{2}}\times {\displaystyle \frac{{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}}{\left|{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}\right|}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\sqrt{g}}}\left(\right({\mathit{a}}_{\mathrm{2}}\cdot {\mathit{a}}_{\mathrm{2}}){\mathit{a}}_{\mathrm{1}}-({\mathit{a}}_{\mathrm{2}}\cdot {\mathit{a}}_{\mathrm{1}}\left){\mathit{a}}_{\mathrm{2}}\right),\\ {\displaystyle }{\mathit{a}}_{\mathrm{3}}\times {\mathit{a}}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}}{\left|{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}\right|}}\times {\mathit{a}}_{\mathrm{1}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\sqrt{g}}}\left(\right({\mathit{a}}_{\mathrm{1}}\cdot {\mathit{a}}_{\mathrm{1}}){\mathit{a}}_{\mathrm{2}}-({\mathit{a}}_{\mathrm{1}}\cdot {\mathit{a}}_{\mathrm{2}}\left){\mathit{a}}_{\mathrm{1}}\right).\end{array}\end{array}$$

Eventually, after expanding the cross products in Eq. (27) using the vector triple products in Eq. (28), it can be shown that this eventually yields the same result $\mathbf{D}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}$.

5 Gradient mapping using the right Penrose–Moore inverse

Equation (15) can be written as

$$\begin{array}{}\text{(29)}& \mathbf{J}\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}.\end{array}$$

Since the Jacobian is rectangular, Eq. (29) represents an underdetermined system with infinite solutions. However, the particular solution that represents the desired mapping needs to be a solution that lies in the row space of J, namely, the row space of J contains the covariant base vectors spanning the local fracture space. To reflect this condition, Eq. (29) is written as

$$\begin{array}{}\text{(30)}& {\mathbf{JJ}}^{\mathrm{T}}\mathbf{M}={\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}},\end{array}$$

where the matrix ∇N^T=J^TM now lies within the row space of J. Equation (30) has a unique solution:

$$\begin{array}{}\text{(31)}& \mathbf{M}=({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}.\end{array}$$

Thus, the same result as in Eq. (25) is obtained:

$$\begin{array}{}\text{(32)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}.\end{array}$$

This can also be written as

$$\begin{array}{}\text{(33)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{J}}^{\mathit{†}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}},\end{array}$$

where J^† is the so-called right Penrose–Moore inverse given by

$$\begin{array}{}\text{(34)}& {\mathbf{J}}^{\mathit{†}}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}.\end{array}$$

A more in-depth background on the Penrose–Moore inverse is provided in Appendix A.

6 Gradient mapping using directional cosines

For each point on a possibly curved two-dimensional discrete element, it is possible to construct a two-dimensional orthonormal coordinate system tangent to the fracture defined by unit vectors ${\widehat{\mathit{e}}}_{\mathrm{1}}$ and ${\widehat{\mathit{e}}}_{\mathrm{2}}$. There are several possibilities, but here the procedure starts with taking the vector ${\widehat{\mathit{e}}}_{\mathrm{1}}$ parallel to the first covariant basis a₁:

$$\begin{array}{}\text{(35)}& {\widehat{\mathit{e}}}_{\mathrm{1}}={\displaystyle \frac{{\mathit{a}}_{\mathrm{1}}}{\left|{\mathit{a}}_{\mathrm{1}}\right|}}.\end{array}$$

The vector ${\widehat{\mathit{e}}}_{\mathrm{2}}$ can be easily obtained making use of the normal n.

$$\begin{array}{}\text{(36)}& {\widehat{\mathit{e}}}_{\mathrm{2}}={\widehat{\mathit{e}}}_{\mathrm{1}}\times \mathit{n}\end{array}$$

This two-dimensional orthonormal coordinate system can be expanded into three dimensions by adding a third unit vector:

$$\begin{array}{}\text{(37)}& {\widehat{\mathit{e}}}_{\mathrm{3}}=\mathit{n}.\end{array}$$

The differential volume simply follows from the covariant base vectors:

$$\begin{array}{}\text{(38)}& d\mathrm{\Omega }=\left|{\mathit{a}}_{\mathrm{1}}\times {\mathit{a}}_{\mathrm{2}}\right|d{\mathrm{\Omega }}^{\ast }.\end{array}$$

The transformation from the global coordinate system to the new coordinate system ${\widehat{x}}^i$ is given by a 2×3 matrix of directional cosines:

$$\begin{array}{}\text{(39)}& T_{ij}^{\mathrm{2}\times \mathrm{3}}={\displaystyle \frac{\partial {\widehat{x}}^i}{\partial x^j}}={\widehat{\mathit{e}}}_i\cdot {\mathit{e}}_j=\cos ({\widehat{x}}^i,x^j).\end{array}$$

The matrix of directional cosines can also be expressed as

$$\begin{array}{}\text{(40)}& \mathbf{T}=\left[\begin{array}{c}{\widehat{\mathit{e}}}^{\mathrm{1}}\\ {\widehat{\mathit{e}}}^{\mathrm{2}}\end{array}\right],\end{array}$$

where it is noted that the inverse of T is T^T.

The gradient matrix with respect to the new two-dimensional orthonormal coordinate system ∇^∧N^T is given by

$$\begin{array}{}\text{(41)}& {\mathrm{\nabla }}^{\wedge }{\mathbf{N}}^{\mathrm{T}}={\widehat{\mathbf{J}}}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}},\end{array}$$

where the Jacobian matrix $\widehat{\mathbf{J}}$ is an invertible 2×2 matrix:

$$\begin{array}{}\text{(42)}& {\widehat{J}}_{ki}={\displaystyle \frac{\partial {\widehat{x}}^i}{\partial s^k}}={\displaystyle \frac{\partial N_n}{\partial s^k}}{\widehat{x}}_n^i.\end{array}$$

Using the global coordinate matrix, this Jacobian is computed with

$$\begin{array}{}\text{(43)}& \widehat{\mathbf{J}}={\mathbf{B}}^{\ast T}{\mathbf{XT}}^{\mathrm{T}}.\end{array}$$

The gradient matrix with respect to global coordinates follows from using the chain rule:

$$\begin{array}{}\text{(44)}& {\displaystyle \frac{\partial N_n}{\partial x^i}}={\displaystyle \frac{\partial N_n}{\partial {\widehat{x}}^k}}{\displaystyle \frac{\partial {\widehat{x}}^k}{\partial x^i}}\end{array}$$

which can be expressed using the transformation matrix:

$$\begin{array}{}\text{(45)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{T}}^{\mathrm{T}}{\widehat{\mathbf{J}}}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}\end{array}$$

This expression looks quite different compared to the expressions obtained using the first and second approach. However, it can be illustrated that the result is identical. From the chain rule:

$$\begin{array}{}\text{(46)}& {\displaystyle \frac{\partial s^k}{d{\widehat{x}}^i}}={\displaystyle \frac{\partial x^j}{d{\widehat{x}}^i}}{\displaystyle \frac{\partial s^k}{d{x}^j}},\end{array}$$

and it follows that

$$\begin{array}{}\text{(47)}& {\widehat{\mathbf{J}}}^{-\mathrm{1}}=\mathbf{TD}.\end{array}$$

Therefore, Eq. (45) is identical to Eq. (21):

$$\begin{array}{}\text{(48)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{T}}^{\mathrm{T}}{\widehat{\mathbf{J}}}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}={\mathbf{T}}^{\mathrm{T}}\mathbf{TD}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}=\mathbf{D}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}\end{array}$$

Since the covariant bases are used to construct a two-dimensional orthonormal coordinate system, the approach as discussed here is applicable to curved fracture elements. In the existing literature (Diersch, 2013; Kolditz, 2002; Watanabe, 2011), the two-dimensional orthonormal space is often constructed using the edges of non-curved fracture elements. That is, the unit normal is constructed from two element edges: the first unit vector is taken parallel to the first edge and finally a cross product of the unit normal and the first unit vector is used to compute the second unit vector. Such an approach assumes that the two-dimensional orthonormal space is constant across the fracture element, which is only valid for non-curved fracture elements.

7 Coordinate transformations for the hydraulic conductivity tensor

Here, it is assumed that a hydraulic tensor is initially provided with respect to the local strike and dip directions for each fracture element. On curvilinear elements, the strike and dip directions vary from point to point. Given the normal n, which also varies from point to point within a curved fracture element and a vertical unit vector v, the unit vector in the strike direction is given by

$$\begin{array}{}\text{(49)}& {\stackrel{\mathrm{‾}}{\mathit{e}}}_{\mathrm{1}}=\mathit{n}\times \mathit{v}.\end{array}$$

The unit vector in the dip direction follows directly from the following cross product:

$$\begin{array}{}\text{(50)}& {\stackrel{\mathrm{‾}}{\mathit{e}}}_{\mathrm{2}}=\mathit{n}\times {\stackrel{\mathrm{‾}}{\mathit{e}}}_{\mathrm{1}}.\end{array}$$

Finally, the unit vector normal to the fracture is given by

$$\begin{array}{}\text{(51)}& {\stackrel{\mathrm{‾}}{\mathit{e}}}_{\mathrm{3}}=\mathit{n}.\end{array}$$

The transformation from the orthonormal local coordinate system aligned with the strike and dip directions to the global coordinate system is defined by the following 3×2 matrix:

$$\begin{array}{}\text{(52)}& \mathbf{Q}=\left[\begin{array}{cc}{\widehat{\mathit{e}}}^{\mathrm{1}}& {\widehat{\mathit{e}}}^{\mathrm{2}}\end{array}\right]\end{array}$$

Denoting the two-dimensional hydraulic conductivity tensor in local coordinates by $\stackrel{\mathrm{‾}}{\mathbf{K}}$, the hydraulic conductivity tensor in global coordinates is given by

$$\begin{array}{}\text{(53)}& \mathbf{K}=\mathbf{Q}\stackrel{\mathrm{‾}}{\mathbf{K}}{\mathbf{Q}}^{\mathrm{T}}.\end{array}$$

For curved elements, the normal is to be computed from the covariant vectors using Eq. (6). For non-curved elements, the normal is constant across the element and can be computed by taking the cross product between two element edges.

8 Example

To illustrate how the different gradient mappings are applied in practice, a curved quadratic triangular element is considered. Three Gauss points, each at the midpoint on an edge, are used for numerical integration. The nodal shape functions are defined (Oñate, 2010):

$$\begin{array}{}\text{(54)}& \begin{array}{rl}{\displaystyle }{N}_{\mathrm{1}}& {\displaystyle }=(\mathrm{1}-r-s)(\mathrm{1}-\mathrm{2}r-\mathrm{2}s),\\ {\displaystyle }{N}_{\mathrm{2}}& {\displaystyle }=r(\mathrm{2}r-\mathrm{1}),\\ {\displaystyle }{N}_{\mathrm{3}}& {\displaystyle }=s(\mathrm{2}s-\mathrm{1}),\\ {\displaystyle }{N}_{\mathrm{4}}& {\displaystyle }=\mathrm{4}r(\mathrm{1}-r-s),\\ {\displaystyle }{N}_{\mathrm{5}}& {\displaystyle }=\mathrm{4}rs,\\ {\displaystyle }{N}_{\mathrm{6}}& {\displaystyle }=\mathrm{4}s(\mathrm{1}-r-s),\end{array}\end{array}$$

with r and s being the local coordinates ($\mathrm{0}\le r\le \mathrm{1}$ and $\mathrm{0}\le s\le \mathrm{1}$). The gradient matrix with respect to local coordinates, which varies within the triangle, is given by

$$\begin{array}{}\text{(55)}& {\mathbf{B}}^{\ast }=\left[\begin{array}{cc}-\mathrm{3}+\mathrm{4}r+\mathrm{4}s& -\mathrm{3}+\mathrm{4}r+\mathrm{4}s\\ \mathrm{4}r-\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{4}s-\mathrm{1}\\ \mathrm{4}-\mathrm{8}r-\mathrm{4}s& -\mathrm{4}r\\ \mathrm{4}s& \mathrm{4}r\\ -\mathrm{4}s& \mathrm{4}-\mathrm{4}r-\mathrm{8}s\end{array}\right].\end{array}$$

In this example, the coordinates of the element nodes are given by the following coordinate matrix:

$$\begin{array}{}\text{(56)}& \mathbf{X}=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{2}}& -\frac{\mathrm{1}}{\mathrm{8}}& \frac{\mathrm{1}}{\mathrm{2}}\\ \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{8}}& \frac{\mathrm{1}}{\mathrm{2}}\end{array}\right].\end{array}$$

Figure 2 illustrates the triangle. The conductivity tensor in this example is defined as

$$\begin{array}{}\text{(57)}& \stackrel{\mathrm{‾}}{\mathbf{K}}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right].\end{array}$$

The Jacobian for each Gauss point is computed using $\mathbf{J}={\mathbf{B}}^{\ast \mathrm{T}}\mathbf{X}$. For the first Gauss point given by $(r,s)=(\mathrm{0},\mathrm{1}/\mathrm{2})$, the Jacobian is

$$\begin{array}{}\text{(58)}& \mathbf{J}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{2}}& \mathrm{1}\end{array}\right].\end{array}$$

This Jacobian provides the covariant base vectors per row. Once the Jacobian is computed, the gradient matrix in global coordinates is easily computed using $\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}$. This is the expression that results from either using contravariant and covariant bases or a right Penrose–Moore inverse. For the first Gauss point, the gradient matrix is given by

$$\begin{array}{}\text{(59)}& \mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}=\left[\begin{array}{cccccc}-\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{\mathrm{2}}{\mathrm{5}}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{5}}& \frac{\mathrm{4}}{\mathrm{5}}& -\frac{\mathrm{4}}{\mathrm{5}}& -\frac{\mathrm{4}}{\mathrm{5}}\\ -\frac{\mathrm{4}}{\mathrm{5}}& \mathrm{0}& -\frac{\mathrm{4}}{\mathrm{5}}& -\frac{\mathrm{8}}{\mathrm{5}}& \frac{\mathrm{8}}{\mathrm{5}}& \frac{\mathrm{8}}{\mathrm{5}}\end{array}\right]\end{array}$$

and

$$\begin{array}{}\text{(60)}& \partial \mathrm{\Omega }=\sqrt{\text{det}\left({\mathbf{JJ}}^{\mathrm{T}}\right)}={\displaystyle \frac{\sqrt{\mathrm{5}}}{\mathrm{2}}}.\end{array}$$

The transformation matrix for the hydraulic conductivity tensor is obtained from the vectors ${\stackrel{\mathrm{‾}}{\mathbf{e}}}_{\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{\mathbf{e}}}_{\mathrm{1}}$. For the first Gauss point,

$$\begin{array}{}\text{(61)}& \mathbf{Q}=\left[\begin{array}{cc}\mathrm{0}& -\mathrm{1}\\ -\frac{\sqrt{\mathrm{5}}}{\mathrm{5}}& \mathrm{0}\\ \frac{\mathrm{2}\sqrt{\mathrm{5}}}{\mathrm{5}}& \mathrm{0}\end{array}\right].\end{array}$$

Applying Eq. (53) gives

$$\begin{array}{}\text{(62)}& \mathbf{K}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{5}}& \frac{-\mathrm{2}}{\mathrm{5}}\\ \mathrm{0}& \frac{-\mathrm{2}}{\mathrm{5}}& \frac{\mathrm{4}}{\mathrm{5}}\end{array}\right].\end{array}$$

The contribution to G from each Gauss point i is computed from

$$\begin{array}{}\text{(63)}& w_i(\mathrm{\nabla }\mathbf{NK}\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}\partial \mathrm{\Omega }{)}_i,\end{array}$$

where w_i is the Gauss weight for each Gauss point i. For the triangular element considered here, this weight is 1/6 for all Gauss points (Diersch, 2013; Oñate, 2010). For the first Gauss point, this contribution equals

$$\begin{array}{}\text{(64)}& \left[\begin{array}{cccccc}\frac{\mathrm{3}}{\mathrm{20}}& \frac{-\mathrm{1}}{\mathrm{12}}& \frac{\mathrm{1}}{\mathrm{15}}& \frac{\mathrm{2}}{\mathrm{15}}& \frac{-\mathrm{2}}{\mathrm{15}}& \frac{-\mathrm{2}}{\mathrm{15}}\\ \frac{-\mathrm{1}}{\mathrm{12}}& \frac{\mathrm{1}}{\mathrm{12}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{15}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{15}}& \frac{\mathrm{2}}{\mathrm{15}}& \frac{-\mathrm{2}}{\mathrm{15}}& \frac{-\mathrm{2}}{\mathrm{15}}\\ \frac{\mathrm{2}}{\mathrm{15}}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{15}}& \frac{\mathrm{4}}{\mathrm{15}}& \frac{-\mathrm{4}}{\mathrm{15}}& \frac{-\mathrm{4}}{\mathrm{15}}\\ \frac{-\mathrm{2}}{\mathrm{15}}& \mathrm{0}& \frac{-\mathrm{2}}{\mathrm{15}}& \frac{-\mathrm{4}}{\mathrm{15}}& \frac{\mathrm{4}}{\mathrm{15}}& \frac{\mathrm{4}}{\mathrm{15}}\\ \frac{-\mathrm{2}}{\mathrm{15}}& \mathrm{0}& \frac{-\mathrm{2}}{\mathrm{15}}& \frac{-\mathrm{4}}{\mathrm{15}}& \frac{\mathrm{4}}{\mathrm{15}}& \frac{\mathrm{4}}{\mathrm{15}}\end{array}\right]\sqrt{\mathrm{5}}.\end{array}$$

For brevity, the calculations for the other two Gauss points are not presented but are computed similarly.

Figure 2 The curved quadratic triangle as used in the example. Nodes are represented as red circles.

Using the gradient mapping based on directional cosines requires more work, as illustrated below. After obtaining the Jacobian, the covariant base vectors are used to compute ${\widehat{\mathbf{e}}}_{\mathrm{1}}$ and ${\widehat{\mathbf{e}}}_{\mathrm{2}}$. For the first Gauss point, this results in the following 2×3 transformation matrix:

$$\begin{array}{}\text{(65)}& \mathbf{T}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\sqrt{\mathrm{5}}}{\mathrm{5}}& \frac{-\mathrm{2}\sqrt{\mathrm{5}}}{\mathrm{5}}\end{array}\right].\end{array}$$

Using $\widehat{\mathbf{J}}={\mathbf{B}}^{\ast T}{\mathbf{XT}}^{\mathrm{T}}$ (Eq. 43) and $\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{T}}^{\mathrm{T}}{\widehat{\mathbf{J}}}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}$ (Eq. 45) then yields the same result as before (Eq. 60).

9 Application

The application considers steady-state flow towards a pipe penetrating a curved fracture using quadratic triangles. Although relatively simple, this model application could be used, for example, to simulate the water inrush during tunnel construction when a fracture is penetrated, provided the geometry and hydraulic properties of the fracture are known. The triangular mesh is generated with Triangle (Shewchuk, 1996) on a flat surface with a length of a 100 m and a height of 50 m. Subsequently, the mesh is wrapped on a curved surface. The pipe is represented by a hole approximating a circular pipe that intersects the fracture horizontally. A constant pressure head of 0 m is prescribed at the nodes around this hole. At the top and bottom of the fracture, constant hydraulic heads of 10 m are prescribed. The hydraulic conductivity tensor is given by

$$\begin{array}{}\text{(66)}& \stackrel{\mathrm{‾}}{\mathbf{K}}=\left[\begin{array}{cc}\mathrm{0.75}& \mathrm{0.25}\\ \mathrm{0.25}& \mathrm{0.75}\end{array}\right]{\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}.\end{array}$$

Thus, the principal axes of this tensor make an angle of 45° with respect to local strike and dip directions.

The simulation is carried out with a finite element groundwater flow model developed by the author. Figure 3 illustrates the solution in terms of hydraulic heads inside the fracture. It is observed that the hydraulic head contours in the vicinity of the pipe follow an ellipsoid as a result of the anisotropic hydraulic conductivity tensor. The pipe penetrating the fracture is represented by a hole instead of a single boundary node for two reasons. Firstly, by using a hole, the radius of the pipe is accounted for (approximately). Secondly, since the hole is represented by relatively short edges, relatively small triangles are generated around the pipe, such that the model can better resolve the steep hydraulic gradients around this part of the model domain.

Figure 3 Hydraulic heads simulated in a curved fracture with an anisotropic hydraulic conductivity tensor.

It is noticed that because the fracture does not lie in a flat plane, the fracture is embedded in three-dimensional space. Therefore, gradient matrices must be mapped from two-dimensional local to three-dimensional global space. In the finite element code applied here, this mapping is preferably based the expression $\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}$ as the routines for this mapping are most concise. The alternative gradient mapping based on directional cosines provides exactly the same result. The hydraulic conductivity tensor must also be mapped to global space.

10 Discussion and conclusion

A key component in the finite element method is the mapping of the gradient matrix from local to global space. If the global space has the same dimension as the local space, then the mapping is defined by the inverse of the Jacobian matrix. However, in the case of two-dimensional fractures embedded in a global three-dimensional space, the Jacobian matrix is non-square and thus not invertible (Juanes et al., 2002; Perrochet, 1995). A couple of different techniques enable a mapping from two-dimensional local to three-dimensional global space.

It is shown in this work that applying the right Penrose–Moore inverse is an efficient, elegant and relatively simple alternative to find an expression to map the gradient matrix. This alternative avoids the use of tensor calculus or the use of cumbersome rotation matrices. Instead, it uses the concept of subspaces associated with matrices. It is also shown that the mapping approach based on an intermediate mapping to a two-dimensional orthonormal space and a subsequent mapping to the global space can be applied to curved elements. The approach based on the right Penrose–Moore inverse, the approach based on covariant and contravariant vectors, and the approach based on an intermediate mapping to a two-dimensional orthonormal space and a subsequent mapping to the global space all yield the same mapping result. If the Jacobian J is readily available, the expression $\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}={\mathbf{J}}^{\mathrm{T}}({\mathbf{JJ}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\mathrm{\nabla }}^{\ast }{\mathbf{N}}^{\mathrm{T}}$, as derived from the right Penrose–Moore inverse or from the approach based on covariant and contravariant vectors, is particularly straightforward to implement in computer code. While the approach based on an intermediate mapping to a two-dimensional orthonormal space and a subsequent mapping to the global space is easier to understand from a mathematical point of view, it involves extra steps that involve cumbersome rotation matrices. The Jacobian J is readily available if the finite element code uses numerical integration. It is noted that analytical integration avoids the need to define the local space, and as such the Jacobian J is not defined. Thus, if analytical integration is used, then there is no alternative for implementing the intermediate mapping to a two-dimensional orthonormal space and the subsequent rotation to the global space. In general, however, it can be argued that numerical integration is to be preferred, since it is far easier to implement (even without considering the mapping problem for fracture elements). Moreover, numerical integration is more general as it can be applied to all finite element types including curved elements.

Finally, this work provides a general approach, applicable to curved elements, to map hydraulic tensors as defined in a local orthonormal coordinate system aligned with the strike, dip and normal directions to the global coordinate system.

Appendix A

The Penrose–Moore inverse is widely used to solve overdetermined and underdetermined linear systems. By definition, the Penrose–Moore inverse satisfies the following conditions (Penrose, 1955):

$$\begin{array}{}\text{(A1)}& \begin{array}{rlrl}{\displaystyle }& {\displaystyle }{\mathbf{AA}}^{\mathit{†}}\mathbf{A}=\mathbf{A}& {\displaystyle }& {\displaystyle }\left(\text{I}\right)\\ {\displaystyle }& {\displaystyle }{\mathbf{A}}^{\mathit{†}}{\mathbf{AA}}^{\mathit{†}}=\mathbf{A}& {\displaystyle }& {\displaystyle }\left(\text{II}\right)\\ {\displaystyle }& {\displaystyle }{\mathbf{AA}}^{\mathit{†}}=({\mathbf{AA}}^{\mathit{†}}{)}^{\mathrm{T}}& {\displaystyle }& {\displaystyle }\left(\text{III}\right)\\ {\displaystyle }& {\displaystyle }{\mathbf{A}}^{\mathit{†}}\mathbf{A}=({\mathbf{A}}^{\mathit{†}}\mathbf{A}{)}^{\text{T}}& {\displaystyle }& {\displaystyle }\left(\text{IV}\right)\end{array}\end{array}$$

Condition (I) implies that A^†A is idempotent (A^†AA^†A=A^†A), and condition (IV) implies that A^†A is Hermitian. Therefore, A^†A is an orthogonal projection matrix. Using the right Penrose–Moore inverse ${\mathbf{A}}^{\mathit{†}}={\mathbf{A}}^{\mathrm{T}}({\mathbf{AA}}^{\mathrm{T}}{)}^{-\mathrm{1}}$ as used for an underdetermined system, A^†A is expressed as

$$\begin{array}{}\text{(A2)}& {\mathbf{A}}^{\mathit{†}}\mathbf{A}={\mathbf{A}}^{\mathrm{T}}({\mathbf{AA}}^{\mathrm{T}}{)}^{-\mathrm{1}}\mathbf{A}.\end{array}$$

This last expression illustrates that A^†A is the orthogonal projection matrix ${\mathbf{P}}_{\mathcal{R}\left(A^{\mathrm{T}}\right)}$ onto the column space or range of A^T (Strang, 2022) which equals the row space or range of A, and as such I−AA^† is the orthogonal projection matrix P_𝒩(A) onto the null-space of A:

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\displaystyle }{\mathbf{P}}_{\mathcal{R}\left(A^{\mathrm{T}}\right)}& {\displaystyle }={\mathbf{A}}^{\mathit{†}}\mathbf{A},\\ {\displaystyle }{\mathbf{P}}_{\mathcal{N}\left(A\right)}& {\displaystyle }=\mathbf{I}-{\mathbf{A}}^{\mathit{†}}\mathbf{A}.\end{array}\end{array}$$

Using these orthogonal projection matrices, a solution to an underdetermined system Ax=b can thus be expressed as

$$\begin{array}{}\text{(A4)}& \mathit{x}=\left({\mathbf{A}}^{\mathit{†}}\mathbf{A}\right)\mathit{x}+(\mathbf{I}-{\mathbf{A}}^{\mathit{†}}\mathbf{A})\mathit{x}={\mathbf{A}}^{\mathit{†}}\mathit{b}+(\mathbf{I}-{\mathbf{A}}^{\mathit{†}}\mathbf{A})\mathit{x}.\end{array}$$

This illustrates that the right Penrose–Moore inverse provides the solution x=A^†b that lies within the row space of A.