Verification of Induced Stresses Along a Fault

Context

In this example, we evaluate the induced stresses in a pressurized reservoir displaced by a normal fault (permeable or impermeable). This problem is solved using the poroelastic solver in GEOS to obtain the stress perturbations along the fault plane, which are verified against the corresponding analytical solution (Wu et al., 2020).

Input file

The xml input files for the test case with impermeable fault are located at:

inputFiles/poromechanics/faultPoroelastic_base.xml
inputFiles/poromechanics/impermeableFault_benchmark.xml

The xml input files for the test case with permeable fault are located at:

inputFiles/poromechanics/faultPoroelastic_base.xml
inputFiles/poromechanics/permeableFault_benchmark.xml

A mesh file and a python script for post-processing the simulation results are also provided:

inputFiles/poromechanics/faultMesh.vtu

src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/faultVerification/faultVerificationFigure.py

Description of the case

We simulate induced stresses along a normal fault in a pressurized reservoir and compare our results against an analytical solution. In conformity to the analytical set-up, the reservoir is divided into two parts by an inclined fault. The fault crosses the entire domain, extending into the overburden and the underburden. The domain is horizontal, infinite, homogeneous, isotropic, and elastic. The reservoir is pressurized uniformely upon injection, and we neglect the transient effect of fluid flow. A pressure buildup is applied to: (i) the whole reservoir in the case of a permeable fault; (ii) the left compartment in the case of an impermeable fault. The overburden and underburden are impermeable (no pressure changes). Due to poromechanical effects, pore pressure changes in the reservoir cause a mechanical deformation of the entire domain. This deformation leads to a stress perturbation on the fault plane that could potentially trigger sliding of the fault. Here, the fault serves only as a flow boundary, and the mechanical separation of the fault plane (either by shear slippage or normal opening) is prohibited, like in the analytical example. For verification purposes, a plane strain deformation is considered in the numerical model.

../../../../../../_images/geometry.PNG — Fig. 24 Sketch of the problem

In this example, we set up and solve a poroelastic model to obtain the spatial solutions of displacement and stress fields across the domain upon pressurization. Changes of total stresses along the fault plane are evaluated and compared with the corresponding published work (Wu et al., 2020).

For this example, we focus on the Mesh, the Constitutive, and the FieldSpecifications tags.

Mesh

The following figure shows the mesh used in this problem.

Here, we load the mesh with VTKMesh. The syntax to import external meshes is simple: in the XML file, the mesh file faultMesh.vtu is included with its relative or absolute path to the location of the GEOS XML file and a user-specified label (here FaultModel) is given to the mesh object. This mesh contains quadrilateral elements and local refinement to conform with the fault geometry, and two reservoir compartments displaced by the fault. The size of the reservoir should be large enough to avoid boundary effects.

  <Mesh>
     <VTKMesh         
      name="FaultModel"      
      file="faultMesh.vtu"
      regionAttribute="CellEntityIds"/>
  </Mesh> 

Solid mechanics solver

GEOS is a multi-physics platform. Different combinations of physics solvers available in the code can be applied in different regions of the domain and be functional at different stages of the simulation. The Solvers tag in the XML file is used to list and parameterize these solvers.

To specify a coupling between two different solvers, we define and characterize each single-physics solver separately. Then, we customize a coupling solver between these single-physics solvers as an additional solver. This approach allows for generality and flexibility in constructing multi-physics solvers. The order in which solvers are specified is not important in GEOS. Note that end-users should give each single-physics solver a meaningful and distinct name, as GEOS will recognize these single-physics solvers based on their customized names to create the expected couplings.

As demonstrated in this example, to setup a poromechanical coupling, we need to define three different solvers in the XML file:

the mechanics solver, a solver of type SolidMechanics_LagrangianFEM called here mechanicsSolver (more information here: Solid Mechanics Solver),

    <SolidMechanics_LagrangianFEM 
      name="mechanicsSolver" 
      timeIntegrationOption="QuasiStatic"
      logLevel="1"
      discretization="FE1"
      targetRegions="{ Domain }">
      <NonlinearSolverParameters
        newtonTol = "1.0e-5"
        newtonMaxIter = "15"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-10"/> 
    </SolidMechanics_LagrangianFEM>

the single-phase flow solver, a solver of type SinglePhaseFVM called here singlePhaseFlowSolver (more information on these solvers at Singlephase Flow Solver),

    <SinglePhaseFVM 
      name="singlePhaseFlowSolver"
      logLevel="1"
      discretization="singlePhaseTPFA"
      targetRegions="{ Domain }"> 
      <NonlinearSolverParameters
        newtonTol = "1.0e-6"
        newtonMaxIter = "8"
      />
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-12"/>   
    </SinglePhaseFVM>
  </Solvers>

the coupling solver (SinglePhasePoromechanics) that will bind the two single-physics solvers above, named poromechanicsSolver (more information at Poromechanics Solver).

  <Solvers gravityVector="{0.0, 0.0, 0.0}">
    <SinglePhasePoromechanics 
      name="poromechanicsSolver" 
      solidSolverName="mechanicsSolver"
      flowSolverName="singlePhaseFlowSolver"
      logLevel="1"
      targetRegions="{ Domain }">
      <LinearSolverParameters 
        solverType="gmres"
        preconditionerType="mgr"
        logLevel="1"
        krylovAdaptiveTol="1"
      />
      <NonlinearSolverParameters
        newtonMaxIter = "40"
      />
    </SinglePhasePoromechanics>

The two single-physics solvers are parameterized as explained in their corresponding documents.

In this example, let us focus on the coupling solver. This solver (poromechanicsSolver) uses a set of attributes that specifically describe the coupling process within a poromechanical framework. For instance, we must point this solver to the designated fluid solver (here: singlePhaseFlowSolver) and solid solver (here: mechanicsSolver). These solvers are forced to interact with all the constitutive models in the target regions (here, we only have one, Domain). More parameters are required to characterize a coupling procedure (more information at Poromechanics Solver). This way, the two single-physics solvers will be simultaneously called and executed for solving the problem.

Discretization methods for multiphysics solvers

Numerical methods in multiphysics settings are similar to single physics numerical methods. In this problem, we use finite volume for flow and finite elements for solid mechanics. All necessary parameters for these methods are defined in the NumericalMethods section.

As mentioned before, the coupling solver and the solid mechanics solver require the specification of a discretization method called FE1. In GEOS, this discretization method represents a finite element method using linear basis functions and Gaussian quadrature rules. For more information on defining finite elements numerical schemes, please see the dedicated Finite Element Discretization section.

The finite volume method requires the specification of a discretization scheme. Here, we use a two-point flux approximation scheme (singlePhaseTPFA), as described in the dedicated documentation (found here: Finite Volume Discretization).

  <NumericalMethods>
    <FiniteElements>
      <FiniteElementSpace
        name="FE1"
        order="1"/>
    </FiniteElements>

    <FiniteVolume>
      <TwoPointFluxApproximation
        name="singlePhaseTPFA"
        />
    </FiniteVolume>
  </NumericalMethods>

Constitutive laws

For this problem, a homogeneous and isotropic domain with one solid material is assumed for both the reservoir and its surroundings. The solid and fluid materials are named as rock and water respectively, and their mechanical properties are specified in the Constitutive section. PorousElasticIsotropic model is used to describe the linear elastic isotropic response of rock when subjected to fluid injection. And the single-phase fluid model CompressibleSinglePhaseFluid is selected to simulate the flow of water.

  <Constitutive>
    <PorousElasticIsotropic
      name="porousRock"
      solidModelName="rock"
      porosityModelName="rockPorosity"
      permeabilityModelName="rockPerm"
    />

    <ElasticIsotropic
      name="rock"
      defaultDensity="2700"
      defaultYoungModulus="14.95e9"
      defaultPoissonRatio="0.15"
    />

    <CompressibleSinglePhaseFluid 
      name="water"
      defaultDensity="1000"
      defaultViscosity="0.001"
      referencePressure="0e6"
      referenceDensity="1000"
      compressibility="2.09028227021e-10"
      referenceViscosity="0.001"
      viscosibility="0.0"
    />

    <BiotPorosity
      name="rockPorosity"
      defaultGrainBulkModulus="7.12e10"
      defaultReferencePorosity="0.3"
    />  
  
    <ConstantPermeability
      name="rockPerm"
      permeabilityComponents="{1.0e-18, 1.0e-18, 1.0e-18}"
    /> 
  </Constitutive>

All constitutive parameters such as density, viscosity, and Young’s modulus are specified in the International System of Units.

Initial and boundary conditions

The next step is to specify fields, including:

The initial value (the in-situ stresses and pore pressure have to be initialized),

The boundary conditions (pressure buildup within the reservoir and constraints of the outer boundaries have to be set).

In this example, we need to specify isotropic horizontal total stress ( $\sigma_h$ = -60.0 MPa and $\sigma_H$ = -60.0 MPa), vertical total stress ( $\sigma_v$ = -70.0 MPa), and initial reservoir pressure ( $P_0$ = 35.0 MPa). When initializing the model, a normal traction (name="NormalTraction") of -70.0 MPa is imposed on the upper boundary (setNames="{ 91 }") to reach mechanical equilibrium. The lateral and lower boundaries are subjected to roller constraints. These boundary conditions are set up through the FieldSpecifications section.

  <FieldSpecifications>
    <FieldSpecification 
      name="initialPressure"
      initialCondition="1"
      setNames="{all}"
      objectPath="ElementRegions/Domain"
      fieldName="pressure"
      scale="35.0e6"
    />
    	     
    <FieldSpecification 
      name="stressXX"
      initialCondition="1"
      setNames="{all}"
      objectPath="ElementRegions/Domain"
      fieldName="rock_stress"
      component="0"
      scale="-28.499545e6"
    />

    <FieldSpecification 
      name="stressYY"
      initialCondition="1"
      setNames="{all}"
      objectPath="ElementRegions/Domain"
      fieldName="rock_stress"
      component="1"
      scale="-38.499545e6"
    />

    <FieldSpecification 
      name="stressZZ"
      initialCondition="1"
      setNames="{all}"
      objectPath="ElementRegions/Domain"
      fieldName="rock_stress"
      component="2"
      scale="-28.499545e6"
    />

    <FieldSpecification
      name="xconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="0"
      scale="0.0"
      setNames="{ 89, 88 }"/>

    <FieldSpecification
      name="yconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="1"
      scale="0.0"
      setNames="{ 90 }"/>

    <FieldSpecification
      name="zconstraintFront"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="2"
      scale="0.0"
      setNames="{ 92, 93 }"/>

    <Traction
      name="NormalTraction"
      objectPath="faceManager"
      tractionType="normal"
      scale="-70.0e6"      
      setNames="{ 91 }"/>  	  
  </FieldSpecifications>

In this example, the only difference between the impermeable fault and permeable fault cases is how to apply pressure buildup. For the impermeable fault case, a constant pressure buildup is imposed to the left compartment of the reservoir (objectPath="ElementRegions/Domain/97_hexahedra"):

  <FieldSpecifications>
    <FieldSpecification 
      name="injection"
      initialCondition="0"
      setNames="{all}"
      objectPath="ElementRegions/Domain/97_hexahedra"
      fieldName="pressure"      
      scale="55.0e6"/>   	  
  </FieldSpecifications>

For the permeable fault case, a constant pressure buildup is imposed to both compartments of the reservoir: (objectPath="ElementRegions/Domain/97_hexahedra" and objectPath="ElementRegions/Domain/96_hexahedra"):

  <FieldSpecifications>
   <FieldSpecification 
      name="injection"
      initialCondition="0"
      setNames="{all}"
      objectPath="ElementRegions/Domain/97_hexahedra"
      fieldName="pressure"      
      scale="55.0e6"/>

    <FieldSpecification 
      name="injection2"
      initialCondition="0"
      setNames="{all}"
      objectPath="ElementRegions/Domain/96_hexahedra"
      fieldName="pressure"      
      scale="55.0e6"/>	  	  
  </FieldSpecifications>

The parameters used in the simulation are summarized in the following table, which are specified in the Constitutive and FieldSpecifications sections. Note that stresses and traction have negative values, due to the negative sign convention for compressive stresses in GEOS.

Symbol	Parameter	Unit	Value
$E$	Young’s Modulus	[GPa]	14.95
$\nu$	Poisson’s Ratio	[-]	0.15
$\sigma_h$	Min Horizontal Stress	[MPa]	-60.0
$\sigma_H$	Max Horizontal Stress	[MPa]	-60.0
$\sigma_v$	Vertical Stress	[MPa]	-70.0
$p_0$	Initial Reservoir Pressure	[MPa]	35.0
${\Delta}p$	Pressure Buildup	[MPa]	20.0
$K_s$	Grain Bulk Modulus	[GPa]	71.2
$\theta$	Fault Dip	[Degree]	60.0
$\kappa$	Matrix Permeability	[m²]	1.0*10^-18
$\phi$	Porosity	[-]	0.3
$D_L$	Domain Length	[m]	4000.0
$D_W$	Domain Width	[m]	2000.0
$D_T$	Domain Thickness	[m]	1000.0
$Res_T$	Reservoir Thickness	[m]	300.0
$F_{off}$	Fault Vertical Offset	[m]	100.0

Inspecting results

We request VTK-format output files and use Paraview to visualize the results. The following figure shows the distribution of resulting shear stress ( $\sigma_{xy}$ ) in the computational domain for two different cases (a permeable vs. an impermeable fault). Numerical solutions for both cases are also compared with the corresponding analytical solutions.

../../../../../../_images/contour.PNG — Fig. 26 Simulation results of $\sigma_{xy}$ 

Fig. 26 Simulation results of $\sigma_{xy}$ 

The figure below compares the results from GEOS (marks) and the corresponding analytical solution (solid curves) for the change of total stresses ( $\sigma_{xx}$ , $\sigma_{yy}$ and $\sigma_{xy}$ ) along the fault plane. As shown, GEOS reliably captures the mechanical deformation of the faulted reservoir and shows excellent agreement with the analytical solutions for two different scenarios. Differences in the stress perturbations between the cases with permeable and impermeable fault are also noticeable, which suggests that fault permeability plays a crucial role in governing reservoir deformation for the problems with reservoir pressurization or depletion.

(Source code)

../../../../../../_images/faultVerificationFigure.png

To go further

Feedback on this example

For any feedback on this example, please submit a GitHub issue on the project’s GitHub page.