COMSOL 6.4 - Biot Poroelasticity

Poroelastic models describe the linked interaction between fluids and deformation in porous media. The fluids in a reservoir absorb stress, which registers as fluid pressure or equally hydraulic head. For example, if pumping significantly reduces pore fluid pressures, sediments could shift due to the increased load. Because the reduction in the pore space brings about more fluid movement, the reservoir could compact further. It follows that lateral stretching must compensate for the vertical compaction.

Model Definition

This example analyzes fluid and solid behavior within three sedimentary layers overlying an impermeable bedrock in a basin. The bedrock is faulted, which creates a step near a mountain front. The sediment stack totals 420 m at the centerline of the basin (x = 0 m) and thins to 120 m above the step (x > 4000 m). The top two layers are each 20 m thick.

Figure 1: Model geometry showing boundary segments (from Leake and Hsieh, Ref. 1).

Pumping from the lower aquifer reduces hydraulic head down the centerline of the basin by 6 m per year. The head drop moves fluid away from the step. The middle layer is relatively impermeable. The pumping does not diminish the supply of fluids in the unpumped reservoir above it. The flow field is initially at steady state. The period of interest is 10 years.

Governing Equations

This example uses the Poroelasticity multiphysics node, which includes the time rate of change in strain from the Solid Mechanics interface in the Darcy’s Law interface and adds the fluid pressure gradient as stress contribution in the Solid Mechanics interface. The following sections describe the built-in equations when using the Poroelastic Material feature as well as the initial and boundary conditions.

Fluid Flow

Use the Darcy’s Law interface to estimate the flow field in the poroelastic model with the pressure head formulation

(1)

where ρf is the fluid density, H is the pressure head (m), K is the hydraulic conductivity (m/s), εvol is the volumetric strain of the porous matrix, and αB is the Biot–Willis coefficient; in this model the Biot–Willis coefficient equals 1, which is a common value for “soft soils.” You can interpret the term on the right-hand side as the time rate of expansion of the porous matrix. The volume fraction available for liquid increases and thereby gives rise to liquid sink, which is why the sign is reversed in the source term. Leake and Hsieh (Ref. 1) defined Sα using coefficients from the solids equation, the Young’s modulus, E, and Poisson’s ratio, ν. Debate over poroelastic storage coefficients is heated (Ref. 2), and the subscript α here denotes that conventional storage terms might need redefinition for poroelasticity models.

The main differences between the Terzaghi compaction model and this poroelastic analysis lie in material coefficients and sources; the boundary conditions are identical: H is the offset in hydraulic head since an initial steady-state distribution of H0. This subtle twist simplifies describing the boundary conditions: the value at the outlet boundary becomes the decline in hydraulic head with time, and at the upper aquifer, the hydraulic head is fixed at H0. All other boundaries have no-flow conditions. The boundary and initial conditions for the Darcy’s law analysis are

where n is the normal to the boundary. The letters A through E come from Leake and Hsieh (Ref. 1), each letter denoting a specific boundary (see Figure 1).

Porous matrix Deformation

The governing equation for the poroelastic material model is

(2)

here, σ is the stress tensor, ρ is the total density, and g is acceleration of gravity. The poroelastic material model uses Equation 2 to describe changes in the stress tensor σ and porous matrix displacement u due to boundary conditions and changes in pore pressure. This focus on changes in displacement is standard poroelasticity and greatly simplifies specifying loads, boundaries, and initial conditions. Equation 2 defines a state of static equilibrium because the changes in the solid equilibrate quickly, unlike vibrations or waves, that is, there are no time-dependent terms. Still the time rate of change in strain ∂εvol/∂t appears as a coupling term in Darcy’s Law because the solids equation becomes quasi-static when solved simultaneously with a time-dependent flow model.

The governing equation for the bedrock step problem describes the deformation state of plane strain, which is the norm for 2D poroelasticity problems (Ref. 1 and Ref. 2), so the strain normal to the xy-plane equals zero. COMSOL Multiphysics solves Equation 2 using the Solid Mechanics interface.

For an isotropic porous material under plane strain conditions, this simplifies to

(3)

where E is Young’s modulus (Pa) and ν is Poisson’s ratio of the drained porous matrix. The term αBp (Pa) amounts to the fluid pressure contribution and is often described as the fluid-to-structure coupling expression. This contribution is added by the Poroelasticity multiphysics node.

With small deformations for plane strain analyses, the normal strains εxx, εyy, εzz and shear strains εxy, εyz, εxz relate to the displacements u and v as

(4)

Inserting the relationships from Equation 3 and Equation 4 into Equation 2 gives the equation that COMSOL Multiphysics solves.

The boundary conditions on the porous matrix are a series of constraints on the displacement that allow for horizontal movement at the surface and throughout the basin. The base of the sediments is fixed, which means you constrain horizontal and vertical displacement to zero. The upper surface is free to vary in the horizontal and the vertical directions. The boundary on left of Figure 1 is a roller condition, so the sediments are free to move in the vertical direction, but there is no horizontal displacement. These conditions result in the following boundary expressions

where, once again, the letters A through E come from Leake and Hsieh, and denote the boundaries in Figure 1.

Model Data

The coefficients and parameters for the poroelasticity model are as follows:


Variable	Description	Value
g	Gravity	9.82 m/s2
ρf	Fluid density	1000 kg/m3
ρd	Drained density of porous layer	2750 kg/m3
Sα	Poroelastic storage coefficient, aquifer layers	1·10-6 m-1
Sα	Poroelastic storage coefficient, confining layer	1·10-5 m-1
K	Hydraulic conductivity, aquifer layers	25 m/day
K	Hydraulic conductivity, confining layer	0.01 m/day
αB	Biot–Willis coefficient	1
H0	Initial hydraulic head	0 m
H(t)	Declining head boundary	(6 m/year)·t
E	Drained Young’s modulus, aquifer layers	8·108 Pa
E	Drained Young’s modulus, confining layer	8·107 Pa
ν	Drained Poisson’s ratio, all regions	0.25

Results and Discussion

Figure 2, Figure 3, and Figure 4 are Year 2, Year 5, and Year 10 snapshots, respectively, from the COMSOL Multiphysics solution to the problem (Ref. 1) of linked fluid flow and solid deformation near a bedrock step in a sedimentary basin. The shading and arrows, respectively, represent the change in hydraulic head and velocities brought about by pumping from the basin interior at x = 0 m. The fluid moves from the surface toward the well screens in the lower aquifer, with an abrupt change in direction and velocity near the bedrock step. In this way, the flow solution here is similar to the results in Terzaghi compaction model.

The results for the solids displacement tell another story. In Figure 2, Figure 3, and Figure 4, contours and deformed shapes illustrate the total displacement. The plot sequence illustrates the evolution of lateral deformations that compensate for the changing surface elevation above the bedrock step.

Figure 2: Solution for a poroelasticity analysis of the bedrock step problem of Leake and Hsieh (Ref. 1): Hydraulic head (surface plot), displacement (contours and deformations) at Year 2. The arrows indicate the hydraulic head gradient. The vertical axis is expanded for clarity.

Figure 3: Hydraulic head (surface plot), displacement (contours and deformations), and hydraulic head gradient (arrows) at Year 5. The vertical axis and deformation are exaggerated for clarity.

Figure 4: Hydraulic head (surface plot), displacement (contours and deformations), and hydraulic head gradient (arrows) at Year 10. The vertical axis and deformation are exaggerated for clarity.

Figure 5, represents the coupling terms that link the fluid and structure equations. The shading gives the structure-to-fluid link, which is the negative of the time rate of change in strain (see Equation 1).

Figure 5: The solid-to-fluid coupling term evaluated at Year 10 with the poroelastic analysis. The squeeze of porous matrix results in an apparent mass source term.

The Terzaghi and Biot solutions differ most when it comes to predicting the horizontal strain at the edge of the bedrock step. The Biot poroelasticity analysis predicts horizontal strain; the Terzaghi compaction analysis does not. The horizontal strains at the ground surface from the Biot poroelasticity approach appear in Figure 6. It depicts negative strain or compaction immediately on the basin side of the step; the positive strains correspond to tension or lengthening on the mountain side.

Figure 6: COMSOL Multiphysics estimates of horizontal strain from poroelastic analysis for the bedrock step problem of Leake and Hsieh: Year 2 (blue line), Year 5 (green line), and Year 10 (red line).

Failure criteria or expressions defining a critical threshold for stress, strain, or displacements facilitate evaluating whether the strain differential at the bedrock steps is big enough to produce fissures. Figure 7 plots von Mises stresses (surface plot), fluid velocities (streamlines), and total displacement (deformation). The von Mises stresses are postprocessing variables defined by COMSOL Multiphysics in all structural-deformation analyses. The von Mises stresses are variables in many failure expressions and you can also use them as yield functions in the elastoplastic materials dialogs.

Figure 7: COMSOL Multiphysics estimates of von Mises stresses (surface plot), fluid velocities (streamlines), and displacement (deformation) at Year 10. These results correspond to the poroelastic analysis from Leake and Hsieh (Ref. 1). The vertical axis and deformation are exaggerated for clarity.

Notes About the COMSOL Implementation

During the simulation an equilibrium between the gravity forces and the poroelastic forces is reached after about 0.2 years. After that, the simulation results change further only due to the changing hydraulic head at the left-hand side boundary. To reproduce this behavior correctly, it is important to start with a small time step in the beginning to reproduce the development of an equilibrium state correctly, whereas later on, a bigger time step is sufficient and should be chosen to save computational time.

References

1. S.A. Leake and P.A. Hsieh, Simulation of Deformation of Sediments from Decline of Ground-Water Levels in an Aquifer Underlain by a Bedrock Step, U.S. Geological Survey Open File Report 97-47, 1997.

2. H.F. Wang, Theory of Linear Poroelasticity with Application to Geomechanics and Hydrogeology, Princeton Univ. Press, 2000.

Subsurface_Flow_Module/Poromechanics/biot_poroelasticity

Modeling Instructions

Root

Start by opening the model containing the solution for the compaction according to Terzaghi theory and use it as entry point for the poroelasticity study according to Biot’s theory.

Application Libraries

From the menu, choose .

In the window, select > > in the tree.

Click

Add Physics

Add a interface and a multiphysics node.

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Find the subsection. In the table, clear the checkbox for .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Add Multiphysics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Find the subsection. In the table, clear the checkbox for .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Solid Mechanics (solid)

Gravity 1

In the toolbar, click

and choose .

Fixed Constraint 1

In the toolbar, click

and choose .

Select Boundaries 2, 8, and 9 only.

Roller 1

In the toolbar, click

and choose .

Select Boundaries 1, 3, and 5 only.

Multiphysics

Poroelasticity 1 (poro1)

In the window, under > click .

In the window for , locate the section.

From the αB list, choose . 1 is the default value for the Biot–Willis coefficient.

Materials

Fill the required fields in the material properties tables.

Porous Material 1 (pmat1)

In the window, expand the > node, then click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	2750	kg/m³	Basic
Young’s modulus	E	800[MPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.25	1	Young’s modulus and Poisson’s ratio

Porous Material 2 (pmat2)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	2750	kg/m³	Basic
Young’s modulus	E	80[MPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.25	1	Young’s modulus and Poisson’s ratio

Definitions (comp1)

Create a variable to calculate the solid-to-fluid coupling term.

Variables 3

In the window, expand the > node, then click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit
Q_biot	-poro1.alphaBd(solid.evol,TIME)dl.rho	kg/(m³·s)

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

Step 1: Time Dependent

Now set up the fully coupled study. As described in Notes About the COMSOL Implementation in the beginning, the time step has to be small until an equilibrium between gravity forces and poroelastic forces is reached. Studies have shown that this is the case after 0.2 years. After that, a bigger time step is sufficient.

In the window for , locate the section.

From the list, choose .

Click