Hydromechanical embedded finite element for conductive and impermeable strong discontinuities in porous media

The pore pressure inside oil and gas reservoirs compartmentalized by sealing faults increases during injection processes. The rise in the pore pressure can induce fault reactivation, leading to hydraulic issues such as fluid leakage from the reservoir to other layers and seismicity. Therefore


Introduction
The use of deep geological formations to store carbon dioxide (CO 2 ) is agreed to be feasible and important to limit the increase in the average atmospheric temperature (IPCC, 2021;Metz et al., 2005).In these sites, pre-existing discontinuities, such as geological faults, can play a significant role as a trapping mechanism for storing CO 2 (Pereira et al., 2014;Silva et al., 2023;Metz et al., 2005).In the risk assessment and monitoring of the injection and storage of CO 2 , the risk of fault reactivation must be analyzed to guarantee that it will not escape to other rock strata and that no induced seismicity appears.Therefore, numerical simulations of such problems have been extensively studied in the literature (Cappa and Rutqvist, 2011;Chan and Zoback, 2007;Jha and Juanes, 2014;Pan et al., 2016;Pruess and Spycher, 2007;Quevedo et al., 2017;Rueda et al., 2014;Rutqvist et al., 2007Rutqvist et al., , 2013Rutqvist et al., , 2016;;Serajian et al., 2016).
The complexity of simulating a fault reactivation problem requires an accurate method that can capture the reactivation mechanisms.Different numerical approaches have been proposed to represent discontinuities (Berre et al., 2019;Jing, 2003).In the context of the Finite Element Method, the most accurate and computationally cumbersome is the Discrete Fracture Method (DFM), which uses zero-thickness interface elements to model the behavior of the discontinuities (Day and Potts, 1994;de-Pouplana and Oñate, 2018;Nguyen et al., 2017;Segura andCarol, 2008a, 2008b).Despite its computational cost, one of the main difficulties associated with DFMs is the requirement of discontinuity-conforming meshes (Fadakar et al., 2017;Lima et al., 2023).In realistic problems, DFM imposes a huge challenge in mesh generation.
To avoid the need to explicitly represent the discontinuity in the mesh, enriched finite element formulations have gained attention due to their capability to implicitly represent the discontinuity in the domain.
In the Extended Finite Element Method (XFEM), the primary fields of the problem, e.g., the displacement and the pore-pressure fields, are enriched with a discontinuous function using a global or local partition of unity (Belytschko et al., 2001;Crusat et al., 2019;De Borst et al., 2006;Khoei, 2015).Its application in geomechanics extends from fault reactivation problems (Hosseini et al., 2023;Liu, 2020;Prévost and Sukumar, 2016) to hydraulic fracturing problems (Cruz et al., 2018(Cruz et al., , 2019;;Hosseini and Khoei, 2020;Khoei et al., 2015Khoei et al., , 2016Khoei et al., , 2018;;Vahab et al., 2019).Similarly, the embedded finite element methods (EFEM) also implicitly incorporate a discontinuity by enriching the primary fields (Cazes et al., 2016;Damirchi et al., 2021;Jirásek, 2000).The difference between the two approaches is mainly the nature of the enrichment.The XFEM defines the enrichment degrees of freedom in a global sense at the nodes of the finite element mesh, meaning that the enrichment degrees of freedom at these nodes are shared by the elements connected to them.On the other hand, EFEM represents the discontinuity of the primary fields using the jumps that occur at the nodes of the intersection between the discontinuity and the finite element mesh (Cervera et al., 2022;Dias-Da-Costa et al., 2010;Wu, 2011).Therefore, the most straightforward advantage of EFEM over XFEM is the reduced number of enrichment degrees of freedom.In EFEM, if the jump degrees of freedom are taken as local variables, the enrichment degrees of freedom can be condensed by solving the governing equations associated with the discontinuity at an element level (Linder and Armero, 2007;Stanić et al., 2020).However, considering the jump degrees of freedom as global variables defined at the nodes of the discontinuity has the advantage of guaranteeing the continuity of the jump field along its extension (Bybordiani and Dias-da-Costa, 2021;Dias-da-Costa et al., 2009a, 2009b).
In the context of coupled hydromechanical problems, the existing embedded finite element formulations have limitations regarding their application in some practical reservoir simulations.For example, in a fault reactivation simulation, it is essential to model discontinuities that act as preferential paths or as barriers to the fluid flow.To meet the latter behavior, the pressure field needs to be enriched with a Heaviside function (Cavalcanti et al., 2024).Also, the enriched displacement field must include linear deformation modes to avoid stress locking (Linder and Armero, 2007).In the literature, none of the coupled hydromechanical embedded finite element formulations have such features (Armero and Callari, 1999;Benkemoun et al., 2015;Callari, 2004;Callari et al., 2010;Larsson et al., 1996;Lu et al., 2017).Recently, Damirchi et al. (2022) presented a hydromechanical embedded finite element formulation based on the strong discontinuity approach using a coupling finite element (Bitencourt et al., 2015) to assure compatibility between the continuum and discontinuity meshes.They applied the formulation to problems with discontinuities acting as a preferential path for fluid flow.
This paper proposes a fully coupled hydromechanical embedded finite element formulation, considering a strong discontinuity on both displacement and pore-pressure fields, for the transient behavior of deformable porous media.The present formulation allows the modeling of discontinuities that act as preferential paths or as barriers for the fluid flow in the porous media.In this formulation, no coupling finite element or coupling parameter is required.The paper is structured in the following sections.Section 2 presents the derivation of the weak form of the governing equations as well as the definition of the enrichment of the displacement and pore-pressure fields.Section 3 presents how the problem is solved with the proposed embedded finite element formulation.Section 4 presents two validation tests.Section 5 presents the application of the proposed formulation in solving a fault reactivation problem.Finally, Section 6 summarizes the main conclusions of the work.

Governing equations
A domain Ω crossed by discontinuity Γ d is adopted for the following derivations, as shown in Fig. 1.It separates the model into two subdomains (Ω + and Ω − ) based on the normal vector to the discontinuity surface, n d , such that Ω + is the sub-region at which the normal vector n d is pointing.
The boundary conditions associated with the displacement field, u(x, t), are: where σ is the Cauchy stress tensor, t is the traction vector on the boundary with normal n associated with the natural boundary conditions, u is the prescribed displacement, t d is the traction vector on the discontinuity boundary, with normal n d .
The boundary conditions associated with the pore-pressure field in the continuum porous media, p m (x, t) are: where v m is the fluid relative velocity vector, q m is the flux associated with the natural boundary conditions, q d is the flux through the discontinuity, p m is the prescribed pressure associated with the essential boundary conditions, and E F denotes the jump in the field along Γ d .

Approximation of primary fields
The problem is governed by three independent fields: the displacement, u(x,t), the pore-pressure of the continuum porous media, p m (x,t), and the pore-pressure along the discontinuity surface, p f (x,t).
To incorporate a discontinuity into the domain without having it explicitly represented requires that u(x,t) and p m (x,t) be defined as presented for the following generic field g(x) (Oliver, 1996): where g(x) is the regularized part of the field, g(x) is the enrichment part of the field, H Γd (x) is the Heaviside function, and φ(x) is an auxiliary function defined, respectively, as with ĝ(x) being the regular part of the field.Fig. 2 illustrates the composition of the generic field with a regularized enrichment.

Mechanical equilibrium problem
The differential equilibrium equation is: where b represents the body force per unit volume vector.Applying the Galerkin method followed by the Divergence theorem: The integral along the boundary of the domain is decomposed as and applying the boundary conditions, Eq. ( 1), and substituting Eq. ( 9) into Eq.( 8) Finally, adopting Biot's stress decomposition (Biot and Blot, 1941;Biot and Willis, 1957).
where σ' is the effective stress tensor, t d ʹ is the effective cohesive stresse vector, and α and α f are the Biot's coefficient of the porous media and the discontinuity, respectively.
Assuming that the continuum porous media is elastic, the effective stresses are where D is the elastic constitutive matrix.
A traction-separation elastic and an elastoplastic cohesive law governs the discontinuity behavior, For the elastic cohesive model, the effective traction vector is where R d is the rotation matrix from the global cartesian to the discontinuity local coordinate system, and T d is the constitutive tensor with k N and k S being the discontinuity normal and shear cohesive stiffness.To avoid the penetration of the two regions along the discontinuity, a penalty approach is used.The normal stiffness is defined based on the current opening of the discontinuity, w = EuF n + w 0 , where EuF n is the normal displacement jump and w 0 is the initial opening: where k 0 n is the initial normal stiffness, and p is the penalization coefficient, which is equal to 5 E/l e , with E being the Young modulus and l e the characteristic length (Hosseini et al., 2023).
For the elastoplastic cohesive model, an additive decomposition of the displacement jump vector is defined as where EuF e and EuF p are the elastic and plastic parts of the displacement jump, respectively.The yield criterion is a combination of the Mohr-Coulomb, f MC , and a tension cut-off, f TC , surface: where t S and t N are the shear and normal components of the cohesive traction vector, ϕ is the friction angle, c is the cohesion, and f t is the tensile strength.Fig. 3 presents the split of the cohesive stress space t N -t S based on the two yield surfaces and how an inadmissible trial elastic traction state is corrected considering these surfaces.
The plastic part of the displacement jump is computed using the flow rule combining the active surfaces.For instance, when both are active: where λ MC and λ TC are the plastic multipliers and g MC and g TC are the plastic potential surfaces of the Mohr-Coulomb and the tension cut-off yield criteria, respectively.The plastic potential surface for the tension cut-off equals its yield surface.For the Mohr-Coulomb, the plastic potential is where φ is the dilatancy angle.
For more details on the cohesive stress integration algorithm and the definition of the tangent constitutive matrix, see Abreu et al. (2022) and Potts and Zdravkovic (1999).

Fluid flow through the porous media
The continuity equation of the continuum porous media is, where M is Biot's modulus (Biot and Willis, 1957): where ϕ is the porosity, and K s and K f are the bulk moduli of the solid skeleton and the fluid, respectively.Applying the Galerkin method followed by the Divergence theorem gives: Decomposing the integral along the boundary and applying the boundary conditions, eq. ( 2) yields: Neglecting gravity effects, Darcy's law governs the relative fluid velocity in the porous medium where k m is the intrinsic permeability matrix and μ is the dynamic fluid viscosity.
The fluid flow exchange between the discontinuity and the porous medium, q d , can be decomposed into two parts: with where c T and c B are the leak-off parameters, which control the transversal flow through the discontinuity.

Fluid flow inside the discontinuity
The one-dimensional continuity equation of the discontinuity is where M f is Biot's modulus of the discontinuity.The divergence operator is defined at the discontinuity local system as The derivation of the weak form of the continuity equation for the continuum porous media follows the steps of the previous section.The final integral form is the following: The terms of Eq. ( 33) are assumed constant in the discontinuity normal direction, see Fig. 1, and, thus,can be simplified.The first term is rewritten as, Fig. 3. Illustration of the return mapping of cohesive stresses to an admissible state.This representation is valid for f t < c / tan ϕ.
The second term in Eq. ( 33) is The fourth term in Eq. ( 33) is: where k l is the discontinuity longitudinal permeability, which is defined by the cubic law as (Snow, 1965;Witherspoon et al., 1980) Finally, the weak form of the continuity equation results by substituting Eqs. ( 34)-(36) into Eq.( 33):

Embedded finite element formulation
This section presents the embedded finite element discretization as a reduction of the enrichment variable used in an XFEM discretization.In XFEM, the displacement and pore-pressure fields of the porous media are discretized as where, d and a are the standard and enrichment DOFs of the displacement field, p and b are the standard and enrichment DOFs of the pore pressure field in the porous media, N α , α = u, p, contains the standard shape functions of the finite element, and Ñα contains the enrichment shape functions defined as: where d is a diagonal matrix with the Heaviside function evaluated at the nodes associated with each degree of freedom of the finite element.
In the continuity equations, the discretization of the pore pressure field of the porous medium at a specific region of the domain, Ω + or Ω − , is defined as where Ñ(+) p (x) is the matrix with the enhanced shape functions at Ω + , with H Γd (x) = 1, and Ñ(− ) p (x) is the matrix with the enhanced shape functions at Ω − , with H Γd (x) = 0.
The pore pressure in the discontinuity mid-plane is: where p f is the pore pressure DOF vector at the discontinuity mid-plane and N f is a vector containing Lagrange shape functions defined in terms of the tangential coordinate of the discontinuity, s.
The displacement jump along the discontinuity is where Δd is the vector with the discontinuity nodal displacement jumps and N d is a matrix containing Lagrange shape functions defined in terms of the tangential coordinate of the discontinuity.The symmetric part of the displacement field gradient in the continuous part of the domain is where B u (x) = ∂N u (x) and Bu (x) = ∂ Ñu (x), with ∂ being the representation of the symmetric gradient operator in matrix form.The divergent of the displacement field rate is with m T = [1 1 0] T , for bidimensional problems.
The rate of the normal displacement jump is The mean rate of the derivative of the shear displacement jump in the tangential direction is 1 2 with r s = [cos 2 α sin 2 α sin α cos α] and The gradient of the pore-pressure field of the porous media is: where B p (x) = ∇N p (x) and Bp (x) = ∇ Ñp (x).

Mapping matrices
In the proposed EFEM, the enrichment part of displacement and pore-pressure field in the porous media are discretized using the displacement and pressure jumps at the intersection nodes between the discontinuity and the finite element borders, as illustrated in Fig. 4. In order to obtain such a discretization, we define a linear transformation that reduces the enrichment DOFs space to the jumps of each field, as: To obtain M u , the jump of the displacement field is assumed to have a linear variation along the tangential coordinate, s.Then, the transmission of the displacement jump to the continuum is assumed to be a combination of three different deformation modes: a constant translation, a small relative rotation around the discontinuity reference point, and a stretching mode, as illustrated in Fig. 5 (Dias-da-Costa et al., 2009).The linear deformation modes are important to avoid stresslocking behavior (Linder and Armero, 2007).
The displacement mapping matrix of a node, M u (x), can be written as a composition of two parts: the rigid-body deformation modes, M RB u (x), and the stretching modes, M S u (x). with The mapping matrix of the pressure jump of a node, M p (x), is defined assuming that the transmission to the continuum occurs entirely in the normal direction, i.e. (Damirchi et al., 2021): Finally, element mapping matrices are obtained by concatenating by rows the nodes mapping matrices evaluated at the nodes of the element, namely,

Linearization and time integration
Substituting Eqs. ( 39)-( 52) into the weak form of the equilibrium, Eq. ( 10), and both continuity equations, Eqs. ( 30) and (38), leads to the following residual vector, where U = [d T Δd T ] T is the set of displacement DOFs, P = [p T Δp T ] T is the set of the continuum pore-pressure DOFs.The sub-indices in the zeros in eq. ( 58) indicate the size of the vector/matrix, with n U , n P and n pf being the number of components of U, P, and P f , respectively.F ext U (t), F ext P (t) and F ext p f (t) are the nodal external loads and fluid discharge, vectors; F int U is the internal force vector, H m is the fluid flow matrix of the porous medium, H f is the fluid flow matrix of the discontinuity, L m , L mf and L f are the matrices that couple the flow between the porous medium and the discontinuity, S m is the compressibility matrix of the porous medium, S f is the compressibility matrix of the discontinuity, Q is the hydromechanical coupling matrix of the porous medium, Q f and E f are the hydromechanical coupling matrix of the discontinuity, where n d is the number of components of standard displacement DOFs d.
Using a fully implicit time integration scheme and then linearizing the system results in where K is the tangent stiffness matrix: The numerical integration of the matrices and vectors defined over the continuum porous media domain presented in eqs.( 59) -( 73) is performed by applying a sub-division of the element domain into triangular sub-domains, a common practice in XFEM (Khoei, 2015).This is required due to the presence of the Heaviside function in sub-matrices associated with the enrichment DOFs (Dias-da-Costa et al., 2009).The matrices and vectors defined over the discontinuity domain are integrated using a two-node trapezoidal rule (Dias- Da-Costa et al., 2010).
The enrichment DOFs Δd, Δp, and p f , are all considered global degrees of freedom.Therefore, they are not statically condensed.

Validation tests
The solution obtained with the embedded formulation here proposed is compared with that obtained using a discrete fracture model with triple-node interface elements implemented in the in-house computational framework GeMA (Mejia et al., 2020;Mendes et al., 2016;Rueda et al., 2020).
Despite the formulation being flexible and allowing different leak-off parameters at each side of the discontinuity, c T and c B , in the following examples, they are assumed to have the same value c, and are referred to as the leak-off parameter.

Consolidation of a block crossed by a discontinuity
Segura and Carol initially solved this consolidation problem (Segura and Carol, 2008) using fully coupled HM interface elements.The geometry and boundary conditions are summarized in Fig. 6.The displacement jumps in the discontinuity are fixed for comparison with the reference solution.The model properties are Young's modulus of 1.0 MPa, Poisson's ratio of 0.3, normal cohesive stiffness of 2.0e04 kPa/m, shear cohesive stiffness of 1.0e03 kPa/m, discontinuity opening of 2.4038e-04 m, hydraulic conductivity of 1.1574e-05 m/s, fluid dynamic viscosity of 1.0e-06 kPa s, and leak-off parameters of 1.0 m/(kPa s).
The domain was discretized with a regular finite element mesh of 11 x 11 linear quadrilaterals.A transient analysis was performed with a fixed time step of 0.5 s.Fig. 7 presents the EFEM results for the pore pressure profile at the discontinuity mid-plane at different times.It also includes the solutions proposed by Segura and Carol (Segura and Carol, 2008).
The numerical results show that the EFEM captures the influence of the discontinuity in the consolidation process and demonstrates good agreement with the solution obtained using the DFM, as illustrated in Fig. 7.

Consolidation of a block with an inclined internal discontinuity
This benchmark is inspired by the model proposed by Lamb et al. (2013) using XFEM and considering permeable and cohesionless discontinuities.Here, the problem is adapted considering a stressdependent discontinuity crossing the domain.The influence of the discontinuity on transversal flow is investigated.The discontinuity longitudinal permeability, Eq. ( 37), is updated during the analysis.
The geometry and boundary conditions are presented in Fig. 8(a).The nodal displacement jump in the x-direction at the lateral borders is zero.The finite element mesh of the DFM used for validation is presented in Fig. 8(b), which has 2703 nodes, 2500 linear quadrilateral finite elements, and 50 triple-node interface elements (Nguyen et al., 2017;Rueda et al., 2020).
The mechanical properties are Young's modulus of 40.0 MPa and Poisson ratio of 0.3.The discontinuity cohesive material is elastic with normal and shear cohesive stiffness of 5.0e03 kPa/m and an initial opening of 1.063e-03 m.The hydraulic conductivity is 1.0e− 08 m/s, the fluid dynamic viscosity is 1.0e− 06 kPa⋅s, the compressibility of both the fluid and the solid grains are neglected.
The domain was discretized with a structured mesh of 50 × 50 linear quadrilateral finite elements.A transient analysis was performed during 100 days, using a fixed time-step of 1 day.
Fig. 9 shows the pressure evolution over time at the node with coordinates (0.0,8.0) for EFEM and DFM models considering two leak-off parameter values, a high leak-off, c = 1.0 m/(kPa s), and a low leak-off, c = 1.0e− 15 m/(kPa s).The result for a model without discontinuity is   The results in Fig. 9 show good agreement between the EFEM and the DFM solutions, and how the discontinuity can affect the consolidation process.The discontinuity with a high leak-off parameter acts as a preferential path for fluid flow, allowing the pressure to decrease at a higher rate than the model without discontinuity.On the other hand, the discontinuity with a low leak-off parameter acts as a barrier for transversal fluid flow through the porous media.In this case, the pressure drop happens at an even higher rate.To better understand the influence of the discontinuity leak-off parameters over the pressure distribution, Fig. 10 shows the pressure field for the models without and with the discontinuity using EFEM, and Fig. 11 shows their pressure profile along the y-coordinate at x = 5 m.The values correspond to the final time step.
Fig. 10 shows how the discontinuity with different values of leak-off parameter influences the pressure distribution.In Fig. 10(b) a preferential path of fluid flow can be observed, following the direction of the discontinuity.In Fig. 10(c), since the discontinuity is acting as a barrier for the fluid flow, the pressure field on the bottom part of the domain does not dissipate.
Fig. 11 allows visualizing the pressure jump across the discontinuity for the low leak-off parameter scenario.
The influence of the discontinuity with different leak-off values on the displacement field is also analyzed, and the profiles along the ycoordinate at x = 5 m are shown in Fig. 12.
Fig. 12 shows that the discontinuity significantly impacts the displacement field components and that it is strongly affected by the leak-off parameter.The results showed an excellent agreement with the solutions from the DFM.
The cohesive stresses along the discontinuity are compared with solutions using a DFM (see Fig. 13).The results show spurious oscillations on the cohesive stresses along the discontinuity using the EFEM.The literature reports this behavior in association with XFEM, occurring when the discontinuity runs skewed to the mesh (Crusat et al., 2019;Moës et al., 2006;Sanders et al., 2009).Since the EFEM, as presented here, can be seen as a particular case of XFEM, it is assumed that the causes of the oscillations are the same.Despite these spurious oscillations in the cohesive stresses, the results of the displacement and pressure fields agree well with the solutions obtained with a DFM.
Finally, to analyze the effect of spatial discretization, the consolidation problem with a low leak-off parameter was solved using six different finite element meshes: two regular structured meshes with different levels of refinement, 50 × 50 (S50) and 100 × 100 (S100); one uniform unstructured quadrilateral mesh (USQ); one uniform unstructured triangular mesh (UST); one unstructured mesh refined at the region close to the discontinuity (USQR); and one non-regular 50 × 50 structured mesh that aligns to the discontinuity (SA50).Fig. 14 presents the non-uniform and unstructured meshes.
Fig. 15 presents the results for the displacements, cohesive stresses, and pressure for the six finite element meshes.
From the displacement (Fig. 15 (a) and (b)) and pressure (Fig. 15 (e)) profiles along x = 5 m and the pressure evolution over time (Fig. 15 (f)), the results show that they are not mesh-dependent.Meanwhile, the adopted mesh can strongly affect the cohesive stresses along the discontinuity.The unstructured mesh with triangular elements reported major spurious oscillations.This is expected since linear triangular finite elements have a constant strain and stress field and, hence, a poorer representation.The oscillations in the cohesive shear stresses tend to reduce with mesh refinement.However, the cohesive normal stresses for the unstructured quadrilateral mesh with a high refinement close to the discontinuity still presented high oscillations at some points.The points of higher oscillations occur when the discontinuity crosses the finite element very close to one of its vertices.Finally, the oscillations are eliminated with the use of a mesh that is aligned with the discontinuity.

Fault reactivation problem
The fault reactivation problem analyzed here is based on the one studied by Cappa and Rutqvist (2011).The geometry, boundary conditions and DFM mesh are presented in Fig. 16.It consists of a reservoir sealed by two caprock layers crossed by a fault.The fluid is injected with a volumetric discharge of 2.0e-5 m 3 /s.The nodal displacements and pressure jumps at the bottom edge of the model are zero.The origin of the cartesian system is located at the left bottom node of the model, as indicated in Fig. 16.The initial pressure at the top of the reservoir corresponds to a water level located at y = 2500 m, and an atmospheric pressure of 100 kPa acts at the ground surface.
All layers have the following common properties: Young's modulus of 10 GPa, Poisson's ratio of 0.25, fluid dynamic viscosity of 5.98e-7 kPa.s, fluid bulk modulus of 2.2e6 kPa, Biot coefficient of 1, saturated rock specific weight of 22.6 kPa/m, water specific weight of 10 kPa/m, and incompressible solid grains.Other parameters adopted for each layer are summarized in Table 1.
To avoid numerical instabilities due to the very low permeability of the caprock layers, the pressure degrees of freedom of the elements located internally in the caprock layers are fixed at the hydrostatic value.
The fault segment in the aquifers is taken as elastic, and in the reservoir and caprock layers as elastoplastic, following the Mohr-Coulomb criterion with tension cut-off.The normal and shear stiffnesses are 5.0e9 kPa/m, the initial opening is 1.0e-7 m, the friction angle is 25 • , the dilatancy angle is 20 • , and the cohesion and the tensile strength are zero.The pore-pressure field is initialized with a hydrostatic gradient of 10.0 kPa/m, and the vertical effective stress is initialized considering the self-weight with a resultant gradient of 12.6 kPa/m, as shown in Fig. 17 The horizontal effective stresses are initialized by σ ʹ x = K 0 σ ʹ y , with K 0 being the Earth-pressure coefficient.
The initial cohesive stresses along the discontinuity are computed by rotating the stress vector following the direction of the discontinuity:

{
Appendix A presents how these initial conditions are incorporated into the numerical model A simulation was performed for 365 days.The initial time step is 1000 s and is adapted during the analysis limited by a maximum value of 1 day.Two finite element meshes are tested, one structured aligned with the discontinuity and one regular structured, as shown in Fig. 18.
The problem is solved with two values for the leak-off parameter, 1.0e-9 and 1.0e-12 m/(kPa s).Fig. 19 shows the pressure field at the end of the analysis using EFEM with the aligned mesh, presented in Fig. 18 (a).First, results using EFEM with the aligned mesh, presented in Fig. 18 (a), are compared with those using the DFM.The displacement profiles along the y-coordinate at x = 1000 m are shown in Fig. 20 for two leakoff parameters, 1.0e-9 and 1.0e-12 m/(kPa s).The results show a good agreement between the EFEM and DFM solutions, and it evidence that the change in the leak-off parameter strongly affects the displacement field.Fig. 21 shows the fault internal pressure along its length in comparison with the pressure at the fault adjacent faces.The results show that EFEM captured a similar pressure build-up as the DFM, for both values of leak-off parameter.
Fig. 21 (a) shows that for the higher leak-off parameter, 1.0e-9 m/ (kPa s), the pore pressure inside and at the fault adjacent planes are approximately the same.However, for a lower leak-off parameter, 1.0e-12 m/(kPa s), the results show that there is a pressure jump across the fault and the internal pressure is an intermediate value concerning the pressure in the adjacent planes.Hosseini et al. (2023) and Jha and Juanes (2014) discuss that to be closer to safety in assessing fault reactivation, the pore pressure inside the discontinuity should be imposed as the maximum value regarding the adjacent planes.This could be accomplished through the present embedded formulation by imposing a constraint on the definition of the pressure field inside the discontinuity.
The increments in the cohesive stresses for the two leak-off parameters at the end of the analysis are presented in Fig. 22.The results show that the normal cohesive stress increment is higher for the lower leak-off parameter, which is related to a higher pore-pressure build-up observed in Fig. 21.Higher shear cohesive stresses at the layers outside the reservoir are also observed for a lower leak-off value.
To assess the fault reactivation, as performed by Quevedo et al. (2017), a slip tendency (ST) parameter, Eq. ( 75 Fig.23 shows that the leak-off parameter strongly influences the reactivation of the fault.For a lower leak-off parameter, since there is a higher pore-pressure build-up, the reactivation happens earlier.It also shows that the first point in the reservoir, in both cases, to reactivate is the one located at the edge with the upper caprock layer. Finally, the influence of the mesh in EFEM is evaluated by comparing the cohesive stress increments and the slip tendency parameter evolution along the fault for the two meshes shown in Fig. 18, for a leak-off parameter of 1.0e-9 m/(kPa s).The results in Fig. 25 show that the spurious oscillations in the cohesive stresses along the fault, observed in    Fig. 24 and discussed in the previous example, impact the fault reactivation analysis.Given that an elastoplastic cohesive model is used, for a same mesh discretization, the spurious oscillations in the cohesive stresses can indicate that the fault is pre-maturely reactivated.We have observed, by performing further tests, that using furtther refined meshes reduces this problem.

Conclusions
This paper presented a fully coupled embedded finite element method (EFEM) capable of modeling discontinuities that act as both preferential paths or barriers for fluid flow in permeable porous media.The coupling between the porous media and the discontinuity is achieved by deriving the weak form of the governing equations and the constitutive model.
The results obtained with the proposed EFEM were successfully verified using a discrete fracture finite element model with triple-node interface elements and were in good agreement with the reference solution.
This work also discussed the spurious oscillations of the cohesive stresses along the discontinuities associated with the mesh discretization when using the EFEM.These oscillations depend on the finite element mesh, happening when the discontinuity crosses the finite element close to one of its vertices.This issue is well reported in the XFEM literature but not that much in the EFEM one.We observed that these oscillations did not affect the global behavior of the model when the cohesive material is linear elastic.
Finally, the proposed EFEM was used to model fault reactivation problems, capturing the principal reactivation mechanisms.The results show that the spurious oscillations in the cohesive stresses can affect the evaluation of the fault reactivation since in this problem it was adopted an elastoplastic cohesive material model.The authors observed that this problem can be mitigated by using more refined meshes, but will further investigate ways to completely avoid it in future works.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
) dΓ with σ′ and t d ′ being the initial stress and cohesive stress tensors, respectively, and the pressures p and p f correspond to the initial hydrostatic state.
The vectors associated with the pressure DOFs are:

Fig. 2 .
Fig. 2. Behavior of a generic one-dimensional field with a strong embedded discontinuity: (a) regularized generic function and (b) discontinuous part of the generic field.
Fig. 4 presents the degrees of freedom (DOFs), based on both XFEM and EFEM, of a quadrilateral finite element crossed by a strong discontinuity with length l d , a tangential unit vector m d = [cos α sin α] T , and a normal unit vector n d = [− sin α cos α] T .

Fig. 6 .
Fig. 6.Geometry and boundary conditions of consolidation in discontinuous block.

Fig. 7 .
Fig. 7. Pore-pressure profile at the discontinuity mid-plane for different times.

Fig. 8 .
Fig. 8. Consolidation of a block with inclined discontinuity (a) Geometry and boundary conditions (Adapted from Lamb et al. (2013), (b) Finite element mesh used in the DFM.

Fig. 9 .
Fig. 9. Pressure evolution over time at the left bottom node.

Fig. 15 .
Fig. 15.Results of the consolidation problem using different meshes: (a) horizontal displacement profile at x = 5 m, (b) vertical displacement profile at x = 5 m, (c) shear cohesive stress along the discontinuity, (d) normal cohesive stress along the discontinuity, (e) pressure profile at x = 5 m, (f) pressure at the right bottom node over time. .

Fig. 21 .
Fig. 21.Fault internal pressure in comparison with the pore pressure at the adjacent faces of the fault for leak-off: (a) 1.0e− 9 and (b) 1.0e− 12 m/(kPa s).
), was evaluated at three points along the fault.The results are shown in Fig. 23.The points of the fault are labeled based on their y-coordinate.ST = |t s | t N tan ϕ + c (75)

Fig. 22 .
Fig. 22.Comparison of (a) the normal and (b) the shear stresses along the fault at the end of the simulation.

D
.Cavalcanti et al.

Table 1
Hydraulic properties of each layer.