Tutorial 9: Hydraulic Fracturing

Context

In this tutorial, we use a fully coupled hydrofracture solver from GEOSX to solve for the propagation of a single fracture within a reservoir with hetrogeneous in-situ properties. Advanced xml features will be used throughout the example.

Objectives

At the end of this tutorial you will know:

• how to use multiple solvers for hydraulic fracturing problems,
• how to specify pre-existing fractures and where new fractures can develop,
• how to construct a mesh with bias,
• how to specify heterogeneous in-situ properties and initial conditions,
• how to use parameters, symbolic math, and units in xml files.

Input files

This tutorial uses a set of input files and table files located at:

examples/hydraulicFracturing/heterogeneousInSituProperties

Note: because these files use the advanced xml features, they must be preprocessed using the geosx_xml_tools package. To install geosx_xml_tools, see Advanced XML Features (geosx_xml_tools)

Description of the case

Here, our goal is to demonstrate how hydraulic fractures are modeled in a typical environment. The in-situ properties and initial conditions are based upon a randomly generated, fractal, 1D layer-cake model.

Preparing the input files

The inputs for this case are contained inside a case-specific (heterogeneousInSitu_singleFracture.xml) and base (heterogeneousInSitu_base.xml) XML file. The tables directory contains the pre-constructed geologic model. This tutorial will first focus on the case-specific input file, which contains the key parameter definitions, then consider the base-file.

Included: including external xml files

At the head of the case-specific xml file is a block that will instruct GEOSX to include an external file. In our case, this points to the base hydraulic fracturing input file.

  <Included>
    <File
      name="./heterogeneousInSitu_base.xml"/>
  </Included>

Parameters: defining variables to be used throughout the file

The Parameters block defines a series of variables that can be used throughout the input file. These variables allow a given input file to be easily understood and/or modified for a specific environment, even by non-expert users. Parameters are specified in pairs of names and values. The names should only contain alphanumeric characters and underlines. The values can contain any type (strings, doubles, etc.).

Parameters can be used throughout the input file (or an included input file) by placing them in-between dollar signs. Barring any circular-definition errors, parameters can be used within other parameters. For example, see the parameter mu_upscaled. The value of this parameter is a symbolic expression, which is denoted by the surrounding back-ticks, and is dependent upon two other parameters. During pre-processing, geosx_xml_tools will subsititue in the parameter definitions, and evaluate the symbolic expression using a python-derived syntax.

A number of the input parameters include optional unit definitions, which are denoted by the square brackets follwing a value. For example, the parameter t_max is used to set the maximum time for the simulation to 20 minutes.

  <Parameters>
    <!-- Use the swarm upscaling law -->
    <Parameter
      name="Nperf"
      value="5"/>

    <Parameter
      name="Nswarm"
      value="8"/>

    <Parameter
      name="mu_init"
      value="0.005"/>

    <Parameter
      name="K_init"
      value="1e6"/>

    <Parameter
      name="mu_upscaled"
      value="`$mu_init$*($Nswarm$**2)`"/>

    <Parameter
      name="K_upscaled"
      value="`$K_init$*($Nswarm$**0.5)`"/>

    <Parameter
      name="ContactStiffness"
      value="1e10"/>

    <!-- Event timing -->
    <Parameter
      name="t_init_a"
      value="1 [min]"/>

    <Parameter
      name="dt_max_a"
      value="10 [s]"/>

    <Parameter
      name="t_init_b"
      value="2 [min]"/>

    <Parameter
      name="dt_max_b"
      value="2 [s]"/>

    <Parameter
      name="t_max"
      value="20 [min]"/>

    <Parameter
      name="dt_max_c"
      value="4 [s]"/>

    <!-- Etc. -->
    <Parameter
      name="table_root"
      value="./tables"/>

    <Parameter
      name="t_allocation"
      value="28 [min]"/>
  </Parameters>

Mesh: building a mesh with biased boundaries

The mesh block for this example uses a biased mesh along the simulation boundaries to reduce the size of the problem, while maintaining the desired spatial extents. For a given element region with bias, the left-hand element will be x% smaller and the right-hand element will be x% larger than the average element size. Along the x-axis of the mesh, we select a value of zero for the first region to indicate that we want a uniform-sized mesh, and we select a bias of -0.6 for the second region to indicate that we want the element size to smoothly increase in size as it moves in the +x direction. The other dimensions of the mesh follow a similar pattern.

  <Mesh>
    <InternalMesh
      name="mesh1"
      xCoords="{ 0, 200, 250 }"
      yCoords="{ -100, 0, 100 }"
      zCoords="{ -150, -100, 0, 100, 150 }"
      nx="{ 50, 5 }"
      ny="{ 10, 10 }"
      nz="{ 5, 25, 25, 5 }"
      xBias="{ 0, -0.6 }"
      yBias="{ 0.6, -0.6 }"
      zBias="{ 0.6, 0, 0, -0.6 }"
      cellBlockNames="{ cb1 }"
      elementTypes="{ C3D8 }"/>
  </Mesh>

Geometry: defining a fracture nodeset

For this example, we want to propagate a single hydraulic fracture along the y=0 plane. To acheive this, we need to define three nodesets:

• source_a: The location where we want to inject fluid. Typically, we want this to be a single face in the x-z plane.
• perf_a: This is the initial fracture for the simulation. This nodeset needs to be at least two-faces wide in the x-z plane (to represent the fracture at least one internal node needs to be open).
• fracturable_a: This is the set of faces where we will allow the fracture to grow. For a problem where we expect the fracture to curve out of the y=0 plane, this could be replaced.

  <Geometry>
    <Box
      name="source_a"
      xMin="-0.1, -0.1, -0.1"
      xMax="4.1, 0.1, 4.1"/>

    <Box
      name="perf_a"
      xMin="-4.1, -0.1, -4.1"
      xMax="4.1, 0.1, 4.1"/>

    <ThickPlane
      name="fracturable_a"
      normal="0, 1, 0"
      origin="0, 0, 0"
      thickness="0.1"/>
  </Geometry>

Boundary Conditions: defining boundary conditions

The boundary conditions for this problem are defined in the case-specific and the base xml files. The case specific block includes four instructions:

• frac: this marks the initial perforation.
• separableFace: this marks the set of faces that are allowed to break during the simualtion.
• waterDensity: this initializes the fluid in the perforation.
• sourceTerm: this instructs the code to inject fluid into the source_a nodeset. Note the usage of the symbolic expression and parameters in the scale. This boundary condition is also driven by a function, which we will define later.

  <FieldSpecifications>
    <!-- Fracture-related nodesets -->
    <FieldSpecification
      name="frac"
      fieldName="ruptureState"
      initialCondition="1"
      objectPath="faceManager"
      scale="1"
      setNames="{ perf_a }"/>

    <FieldSpecification
      name="separableFace"
      fieldName="isFaceSeparable"
      initialCondition="1"
      objectPath="faceManager"
      scale="1"
      setNames="{ fracturable_a }"/>

    <!-- Fluid Initial Conditions -->
    <FieldSpecification
      name="waterDensity"
      fieldName="water_density"
      initialCondition="1"
      objectPath="ElementRegions"
      scale="1000"
      setNames="{ perf_a }"/>

    <!-- Fluid BC's -->
    <!-- Note: the units of the source flux BC are in kg/s not m3/s -->
    <SourceFlux
      name="sourceTerm"
      objectPath="ElementRegions/Fracture"
      scale="`-1000.0/$Nperf$`"
      functionName="flow_rate"
      setNames="{ source_a }"/>
  </FieldSpecifications>

The base block includes instructions to set the initial in-situ properties and stresses. It is also used to specify the external mechanical boundaries on the system. In this example, we are using roller-boundary conditions (zero normal-displacement). Depending upon how close they are to the fracture, they can significantly effect its growth. Therefore, it is important to test whether the size of the model is large enough to avoid this.

  <FieldSpecifications>
    <FieldSpecification
      name="bulk_modulus"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions"
      fieldName="rock_BulkModulus"
      functionName="bulk_modulus"
      scale="1.0"/>

    <FieldSpecification
      name="shear_modulus"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions"
      fieldName="rock_ShearModulus"
      functionName="shear_modulus"
      scale="1.0"/>

    <FieldSpecification
      name="sigma_xx"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions"
      fieldName="rock_stress"
      component="0"
      functionName="sigma_xx"
      scale="1.0"/>

    <FieldSpecification
      name="sigma_yy"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions"
      fieldName="rock_stress"
      component="1"
      functionName="sigma_yy"
      scale="1.0"/>

    <FieldSpecification
      name="sigma_zz"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions"
      fieldName="rock_stress"
      component="2"
      functionName="sigma_zz"
      scale="1.0"/>

    <!-- Mechanical BC's -->
    <FieldSpecification
      name="x_constraint"
      component="0"
      fieldName="TotalDisplacement"
      objectPath="nodeManager"
      scale="0.0"
      setNames="{ xneg, xpos }"/>

    <FieldSpecification
      name="y_constraint"
      component="1"
      fieldName="TotalDisplacement"
      objectPath="nodeManager"
      scale="0.0"
      setNames="{ yneg, ypos }"/>

    <FieldSpecification
      name="z_constraint"
      component="2"
      fieldName="TotalDisplacement"
      objectPath="nodeManager"
      scale="0.0"
      setNames="{ zneg, zpos }"/>
  </FieldSpecifications>

Solvers: setting up the coupled hydraulic fracturing solver

The Solvers block is located in the base xml file. Note that the gravityVector attribute indicates that we are applying gravity in the -z direction in this problem.

Similar to other coupled physics solver, the Hydrofracture solver is specified in three parts:

• Hydrofracture: this is the primary solver, which will be called by the event manager. Two of its key attributes are the names of the dependent solid and fluid solvers.
• SolidMechanicsLagrangianSSLE: this is the solid mechancis solver.
• SinglePhaseFVM: this is the fluid solver.

The final solver present in this example is the SurfaceGenerator, which manages how faces in the model break.

  <Solvers
    gravityVector="0.0, 0.0, -9.81">
    <Hydrofracture
      name="hydrofracture"
      solidSolverName="lagsolve"
      fluidSolverName="SinglePhaseFlow"
      couplingTypeOption="FIM"
      logLevel="1"
      discretization="FE1"
      targetRegions="{ Region2, Fracture }"
      contactRelationName="fractureContact">
      <NonlinearSolverParameters
        newtonTol="1.0e-5"
        newtonMaxIter="50"
        lineSearchMaxCuts="10"/>
      <LinearSolverParameters
        logLevel="0"
        solverType="gmres"
        preconditionerType="amg"/>
    </Hydrofracture>

    <SolidMechanicsLagrangianSSLE
      name="lagsolve"
      timeIntegrationOption="QuasiStatic"
      logLevel="0"
      discretization="FE1"
      targetRegions="{ Region2 }"
      solidMaterialNames="{ rock }"
      contactRelationName="fractureContact">
      <NonlinearSolverParameters
        newtonTol="1.0e-6"
        newtonMaxIter="5"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-10"
        logLevel="0"/>
    </SolidMechanicsLagrangianSSLE>

    <SinglePhaseFVM
      name="SinglePhaseFlow"
      logLevel="0"
      discretization="singlePhaseTPFA"
      targetRegions="{ Fracture }"
      fluidNames="{ water }"
      solidNames="{ rock }"
      meanPermCoeff="0.8">
      <NonlinearSolverParameters
        newtonTol="1.0e-5"
        newtonMaxIter="10"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-12"
        logLevel="0"/>
    </SinglePhaseFVM>

    <SurfaceGenerator
      name="SurfaceGen"
      logLevel="0"
      fractureRegion="Fracture"
      targetRegions="{ Region2 }"
      solidMaterialNames="{ rock }"
      rockToughness="$K_upscaled$"
      nodeBasedSIF="1"
      mpiCommOrder="1"/>
  </Solvers>

Events: setting up flexible events

Rather than explicitly specify the desired timestep behavior, this example uses a flexible approach for timestepping. The hydrofracture solver is applied in three segments, where maxEventDt indicates the maximum allowable timestep:

• solverApplications_a: this corresponds to the problem initialization, where we request $dt_max_a$=10.
• solverApplications_b: this corresponds to the period of initial fluid injection, where we request $dt_max_b$=2.
• solverApplications_c: this corresponds to rest of the problem where we request $dt_max_c$=4.

Depending upon how well the solution converges, the timestep may be smaller than the maximum requested value. Other key events in this problem include:

• preFracture: this calls the surface generator at the beginning of the problem and helps to initialize the fracture.
• outputs: this produces output silo files.
• restarts: this is a HaltEvent, which tracks the external clock. When the runtime exceeds the spefified value (here $t_allocation$=28 minutes), the code will call the target (which writes a restart file) and instruct the code to exit.

  <Events
    maxTime="$t_max$">
    <!-- Generate the initial fractures -->
    <SoloEvent
      name="preFracture"
      target="/Solvers/SurfaceGen"/>

    <!-- Primary outputs -->
    <PeriodicEvent
      name="outputs"
      target="/Outputs/siloOutput"
      timeFrequency="1 [min]"/>

    <!-- Apply the hydrofracture solver, limiting the timesteps during certain time intervals -->
    <PeriodicEvent
      name="solverApplications_a"
      maxEventDt="$dt_max_a$"
      endTime="$t_init_a$"
      target="/Solvers/hydrofracture"/>

    <PeriodicEvent
      name="solverApplications_b"
      maxEventDt="$dt_max_b$"
      beginTime="$t_init_a$"
      endTime="$t_init_b$"
      target="/Solvers/hydrofracture"/>

    <PeriodicEvent
      name="solverApplications_c"
      maxEventDt="$dt_max_c$"
      beginTime="$t_init_b$"
      target="/Solvers/hydrofracture"/>

    <!-- Watch the wall-clock, write a restart, and exit gracefully if necessary -->
    <HaltEvent
      name="restarts"
      maxRuntime="$t_allocation$"
      target="/Outputs/restartOutput"/>
  </Events>

Functions: building functions to set in-situ properties

The function definitions are in the base xml file, and rely upon the files in the tables directory. The functions in this example include the flow rate over time, the in-situ principal stress, and the bulk/shear moduli of the rock. Note the use of the table_root parameter, which contains the root path to the table files.

The flow_rate TableFunction is an example of a 1D function. It has a single input, which is time. The table is defined using a single coordinateFile:

0.00000e+00
6.00000e+01
1.20000e+02
3.72000e+03
3.78000e+03
1.00000e+09

And a single voxelFile:

0.00000e+00
0.00000e+00
5.00000e-02
5.00000e-02
0.00000e+00
0.00000e+00

Given the specified linear interpoation method, these values define a simple trapezoidal-function. Note: since this is a 1D table, these values could alternately be given within the xml file using the coordinates and values attributes.

The sigma_xx TableFunction is an example of a 3D function. It uses elementCenter as its input, which is a vector. It is specified using a set of three coordinate files (one for each axis), and a single voxel file. The geologic model in this example is a layer-cake, which was randomly generated, so the size of the x and y axes are 1. The interpolation method used here is upper, so the values in the table indicate those at the top of each layer.

  <Functions>
    <!-- Pumping Schedule -->
    <TableFunction
      name="flow_rate"
      inputVarNames="{ time }"
      coordinateFiles="{ $table_root$/flowRate_time.csv }"
      voxelFile="$table_root$/flowRate.csv"
      interpolation="linear"/>

    <!-- Geologic Model -->
    <TableFunction
      name="sigma_xx"
      inputVarNames="{ elementCenter }"
      coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
      voxelFile="$table_root$/sigma_xx.csv"
      interpolation="upper"/>

    <TableFunction
      name="sigma_yy"
      inputVarNames="{ elementCenter }"
      coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
      voxelFile="$table_root$/sigma_yy.csv"
      interpolation="upper"/>

    <TableFunction
      name="sigma_zz"
      inputVarNames="{ elementCenter }"
      coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
      voxelFile="$table_root$/sigma_zz.csv"
      interpolation="upper"/>

    <TableFunction
      name="init_pressure"
      inputVarNames="{ elementCenter }"
      coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
      voxelFile="$table_root$/porePressure.csv"
      interpolation="upper"/>

    <TableFunction
      name="bulk_modulus"
      inputVarNames="{ elementCenter }"
      coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
      voxelFile="$table_root$/bulkModulus.csv"
      interpolation="upper"/>

    <TableFunction
      name="shear_modulus"
      inputVarNames="{ elementCenter }"
      coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
      voxelFile="$table_root$/shearModulus.csv"
      interpolation="upper"/>
  </Functions>

Running the case and inspecting the results

Preprocessing the input file

Because we are using advanced xml features in this example, the input file must be pre-processed using geosx_xml_tools. To build the final input file hydrofracture_processed.xml, run the following:

geosx_bin_dir/preprocess_xml examples/hydraulicFracturing/heterogeneousInSituProperties/heterogeneousInSitu_singleFracture.xml -o hydrofracture_processed.xml

Running the case

This is a moderate-sized example, so it is recommended to run this problem in parallel. For example, this will run the code on the debug partition using a total of 36 cores:

srun -n 36 -ppdebug geosx_bin_dir/geosx -i hydrofracture_processed.xml -x 6 -y 2 -z 3

Inspecting results

In the above example, we requested silo-format output files every minute. We can therefore import these into VisIt or python and visualize the outcome. The following figure shows the extents of the generated fracture over time:

Because we did not explicitly specify any fracture barriers in this example, the fracture dimensions are controlled by the in-situ stresses. During the first couple of minutes of growth, the fracture quickly reaches its maximum/minimum height, which corresponds to a region of low in-situ minimum stress.

The following figures show the aperture and pressure of the hydraulic fracture (near the source) over time:

To go further

Feedback on this tutorial

This concludes the hydraulic fracturing tutorial. For any feedback on this tutorial, please submit a GitHub issue on the project’s GitHub page.

For more details

• More on advanced xml features, please see Advanced XML Features (geosx_xml_tools).
• More on functions, please see Functions.
• More on biased meshes, plase see Mesh Bias.