PDF
Geothermal Doublet
Introduction
The use of geothermal energy has become an important topic worldwide, not only sparked by the debate on climate change. Particularly in regions with lower geothermal energy, the requirements for geothermal plants are high in order to ensure a continuous and long-term energy supply.
This example of a hydrothermal doublet studies the coupled porous media flow and heat transfer problem using the Darcy’s Law and Heat Transfer in Porous Media interfaces. It shows how to set up a system of a layered subsurface with different geological and thermal properties including fractures and the injection and production side of a doublet.
Model Definition
The model geometry (Figure 1) consists of three geological layers ranged by a fault zone. The layers, their elevation, and the fault zone are interpolation functions from an artificial dataset. Different hydraulic and thermal properties are defined for the layers. The evolution of the flow and temperature field over 10 years is studied.
Figure 1: Geometry of the model including three geologic layers, fracture, and the doublet’s production and injection side (edges).
Fluid Flow
The flow in the fracture is in general much faster than in the surrounding porous matrix. The cubic law is a common correlation for modeling fracture flow. It defines the permeability κf (m2) in the fracture according to
where df (m) is the fracture’s aperture and ff the roughness factor. With this definition, the name cubic law results from the definition of the fracture’s transmissivity
.
A hydraulic gradient in x direction is applied as boundary condition on the vertical boundaries. Top and bottom boundary are defined as impermeable. Two edges represent the geothermal doublet injection and production side. Injection and pumping is modeled using the Well feature. Water injection is defined using the mass flow rate option with M0 = ρW·120 l/s.
Heat Transfer
An initial geothermal gradient of 0.03 K/m is applied. The vertical boundaries are defined as open boundaries. This means that a temperature condition using the same geothermal gradient is active if n · u < 0 (inflow) or otherwise −n · q = 0 (outflow). Analogously to the Fracture Flow boundary condition in the Darcy’s Law interface, the heat transfer in the fracture is modeled using the fracture boundary condition of the heat transfer interface.
To model the heat source term caused by the injection well, a line heat source feature is applied according to
where Cp (J/(kg·K)) is the water heat capacity, l the well length, Tinj = 278 K the injection temperature, and T the current temperature.
Results and Discussion
The pressure field is shown in Figure 2.
Figure 2: Pressure distribution after 10 years.
Figure 3 shows the temperature distribution in the fracture together with streamlines that visualize the flow field, showing the flow from the injection to the production side of the doublet.
Figure 3: Temperature distribution in the fracture and streamlines for the Darcy velocity field after 10 years.
It is interesting to evaluate the production temperature over time. As shown in Figure 4, the production temperature decreases about 20 K over 10 years. This indicates that the doublet at the operating conditions will not provide a stable long term energy supply and that a different configurations should be tested. This can be done by varying one or more of the parameters in the model to see which setup is more appropriate.
Figure 4: Production temperature over 10 years.
Figure 5 shows the pressure distribution and streamlines of the Darcy velocity field along the fracture.
Figure 5: Pressure distribution and velocity streamlines just along the fracture after 10 years.
Application Library path: Subsurface_Flow_Module/Heat_Transfer/geothermal_doublet
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Fluid Flow > Porous Media and Subsurface Flow > Darcy’s Law (dl).
3
Click Add.
4
In the Select Physics tree, select Heat Transfer > Porous Media > Heat Transfer in Porous Media (ht).
5
Click Add.
6
Click  Done.
Root
Start with loading the parameterized geometry sequence into the model.
Geometry 1
1
In the Geometry toolbar, click Insert Sequence and choose Insert Sequence.
2
Browse to the model’s Application Libraries folder and double-click the file geothermal_doublet_geom_sequence.mph.
3
In the Geometry toolbar, click  Build All.
The geometry parameters that are used were added to the Parameters list automatically. Add a few more parameters that are used to set up the physics interfaces.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
r_bore
1[m]
1 m
Borehole skin zone radius
pump
120[l/s]
0.12 m³/s
Pumping rate
deltaH
1[mm/m]
0.001
Hydraulic head gradient
T_top
283[K]
283 K
Surface temperature
T_inj
278[K]
278 K
Injection temperature
delta_Tz
0.03[K/m]
0.03 K/m
Geothermal gradient
d_f
0.2[cm]
0.002 m
Fracture thickness
f_f
1.6
1.6
Fracture roughness factor
Selections help to improve the whole modeling process. Define them now and use them later where needed.
Definitions
Injection well
1
In the Model Builder window, expand the Component 1 (comp1) > Definitions node.
2
Right-click Definitions and choose Selections > Explicit.
3
In the Settings window for Explicit, locate the Input Entities section.
4
From the Geometric entity level list, choose Edge.
5
Click the  Wireframe Rendering button in the Graphics toolbar.
6
Select Edge 19 only.
7
In the Label text field, type Injection well.
Production well
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, locate the Input Entities section.
3
From the Geometric entity level list, choose Edge.
4
Select Edge 37 only.
5
In the Label text field, type Production well.
Top layer
1
In the Definitions toolbar, click  Explicit.
2
Select Domains 3 and 6 only.
3
In the Settings window for Explicit, type Top layer in the Label text field.
Middle layer
1
In the Definitions toolbar, click  Explicit.
2
Select Domains 2 and 5 only.
3
In the Settings window for Explicit, type Middle layer in the Label text field.
Bottom layer
1
In the Definitions toolbar, click  Explicit.
2
Select Domains 1 and 4 only.
3
In the Settings window for Explicit, type Bottom layer in the Label text field.
Fracture
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, locate the Input Entities section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 21 only.
5
Select the Group by continuous tangent checkbox. This automatically adds all boundaries that belong to the large fault zone.
6
In the Label text field, type Fracture.
Outer boundaries
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, type Outer boundaries in the Label text field.
3
Locate the Input Entities section. From the Geometric entity level list, choose Boundary.
4
Select the Group by continuous tangent checkbox.
5
Select any face on each side. Automatically all faces of this side will be added to the selection.
Global Definitions
Add materials to the materials node, but do not define their properties at this point. After setting up the physics interface COMSOL Multiphysics automatically detects which material properties are needed.
Add Material
1
In the Materials toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in > Water, liquid.
4
Click the Add to Global Materials button in the window toolbar.
5
In the Materials toolbar, click  Add Material to close the Add Material window.
Materials
Upper Aquitard
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, type Upper Aquitard in the Label text field.
Aquifer
1
Right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, type Aquifer in the Label text field.
3
Locate the Geometric Entity Selection section. From the Selection list, choose Middle layer.
Lower aquitard
1
Right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, locate the Geometric Entity Selection section.
3
From the Selection list, choose Bottom layer.
4
In the Label text field, type Lower aquitard.
Fracture
1
Right-click Materials and choose More Materials > Material Link.
2
In the Settings window for Material Link, type Fracture in the Label text field.
3
Locate the Geometric Entity Selection section. From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose Fracture.
Darcy’s Law (dl)
Now, set up the physics interfaces. This determines which material properties are required.
Fracture 1
1
In the Physics toolbar, click  Boundaries and choose Fracture.
2
In the Settings window for Fracture, locate the Boundary Selection section.
3
From the Selection list, choose Fracture.
4
Locate the Aperture section. In the df text field, type d_f.
Fracture Material 1
1
In the Model Builder window, click Fracture Material 1.
2
In the Settings window for Fracture Material, locate the Fracture Material Properties section.
3
From the εp list, choose User defined. In the associated text field, type 0.6.
4
From the Permeability model list, choose Cubic law.
5
In the ff text field, type f_f.
Hydraulic Head 1
1
In the Physics toolbar, click  Boundaries and choose Hydraulic Head.
2
In the Settings window for Hydraulic Head, locate the Boundary Selection section.
3
From the Selection list, choose Outer boundaries.
4
Locate the Hydraulic Head section. In the H0 text field, type -deltaH*x.
Use the Well feature to define the injection and production side of the doublet.
Well 1
1
In the Physics toolbar, click  Edges and choose Well.
2
In the Settings window for Well, locate the Edge Selection section.
3
From the Selection list, choose Injection well.
4
Locate the Well section. In the dw text field, type 2*r_bore.
5
From the Specify list, choose Mass flow.
6
Locate the Mass Flow section. In the M0 text field, type pump*dl.rho.
The expression dl.rho refers to the water density that is defined by Darcy’s Law interface.
Well 2
1
In the Physics toolbar, click  Edges and choose Well.
2
In the Settings window for Well, locate the Edge Selection section.
3
From the Selection list, choose Production well.
4
Locate the Well section. In the dw text field, type 2*r_bore.
5
From the Well type list, choose Production.
6
From the Specify list, choose Mass flow.
7
Locate the Mass Flow section. In the M0 text field, type pump*dl.rho.
Heat Transfer in Porous Media (ht)
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Heat Convection section.
3
From the u list, choose Total Darcy velocity field (dl).
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the Define list, choose Solid phase properties.
Initial Values 1
Define the geothermal gradient as initial temperature distribution.
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T text field, type T_top-0.03[K/m]*z.
Open Boundary 1
1
In the Physics toolbar, click  Boundaries and choose Open Boundary.
2
In the Settings window for Open Boundary, locate the Boundary Selection section.
3
From the Selection list, choose Outer boundaries.
4
Locate the Upstream Properties section. In the Tustr text field, type T_top-delta_Tz*z.
Fracture 1
1
In the Physics toolbar, click  Boundaries and choose Fracture.
2
In the Settings window for Fracture, locate the Boundary Selection section.
3
From the Selection list, choose Fracture.
4
Locate the Shell Properties section. From the Shell type list, choose Nonlayered shell. In the Lth text field, type d_f.
5
Locate the Fluid Material section. From the list, choose Water, liquid (mat1).
6
Locate the Heat Convection section. From the u list, choose Total Darcy velocity field (dl).
7
Locate the Porous Material section. From the θp list, choose Volume fraction (dl/frac1).
8
Locate the Heat Conduction, Porous Matrix section. From the kp list, choose User defined. In the associated text field, type 3.
9
Locate the Thermodynamics, Porous Matrix section. From the ρp list, choose User defined. In the associated text field, type 1.2e3.
10
From the Cp,p list, choose User defined. In the associated text field, type 800.
Line Heat Source 1
1
In the Physics toolbar, click  Edges and choose Line Heat Source.
2
In the Settings window for Line Heat Source, locate the Edge Selection section.
3
From the Selection list, choose Injection well.
4
Locate the Line Heat Source section. In the Ql text field, type mat1.def.Cp*dl.well1.Ml*(T_inj-T).
5
Locate the Heat Source Radius section. Select the Specify heat source radius checkbox.
6
In the R text field, type r_bore.
Materials
Now you can define the material properties.
Upper Aquitard (pmat1)
1
In the Model Builder window, under Component 1 (comp1) > Materials click Upper Aquitard (pmat1).
2
In the Settings window for Porous Material, locate the Homogenized Properties section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Permeability
kappa_iso ; kappaii = kappa_iso, kappaij = 0
1e-10[cm^2]
m²
Basic
4
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Solid 1 (pmat1.solid1)
1
In the Model Builder window, click Solid 1 (pmat1.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 0.9.
4
From the Material list, choose Locally defined.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
1300
kg/m³
Basic
Heat capacity at constant pressure
Cp
900
J/(kg·K)
Basic
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
2
W/(m·K)
Basic
Aquifer (pmat2)
1
In the Model Builder window, under Component 1 (comp1) > Materials click Aquifer (pmat2).
2
In the Settings window for Porous Material, locate the Homogenized Properties section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Permeability
kappa_iso ; kappaii = kappa_iso, kappaij = 0
1e-6[cm^2]
m²
Basic
4
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Solid 1 (pmat2.solid1)
1
In the Model Builder window, click Solid 1 (pmat2.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
From the Material list, choose Locally defined.
4
In the θs text field, type 0.6.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
1900
kg/m³
Basic
Heat capacity at constant pressure
Cp
850
J/(kg·K)
Basic
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
3
W/(m·K)
Basic
Lower aquitard (pmat3)
1
In the Model Builder window, under Component 1 (comp1) > Materials click Lower aquitard (pmat3).
2
In the Settings window for Porous Material, locate the Homogenized Properties section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Permeability
kappa_iso ; kappaii = kappa_iso, kappaij = 0
1e-7[cm^2]
m²
Basic
4
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Solid 1 (pmat3.solid1)
1
In the Model Builder window, click Solid 1 (pmat3.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 0.7.
4
From the Material list, choose Locally defined.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
2300
kg/m³
Basic
Heat capacity at constant pressure
Cp
850
J/(kg·K)
Basic
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
3.5
W/(m·K)
Basic
Mesh 1
Adjust the default mesh settings slightly to make sure, that the injection and production well are resolved properly.
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Fine.
4
Locate the Sequence Type section. From the list, choose User-controlled mesh.
Distribution 1
1
In the Model Builder window, right-click Free Tetrahedral 1 and choose Distribution.
2
Select Edges 19 and 37 only.
3
In the Settings window for Distribution, click  Build All.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Stationary.
4
Find the Physics interfaces in study subsection. In the table, clear the Solve checkbox for Heat Transfer in Porous Media (ht).
5
Click the Add Study button in the window toolbar.
6
In the Home toolbar, click  Add Study to close the Add Study window.
Study 1
Step 1: Stationary
1
In the Settings window for Stationary, locate the Physics and Variables Selection section.
2
Select the Modify model configuration for study step checkbox.
3
In the tree, select Component 1 (comp1) > Darcy’s Law (dl) > Well 1 and Component 1 (comp1) > Darcy’s Law (dl) > Well 2.
4
Click  Disable.
This first stationary step results in a groundwater flow field which serves as a starting point for the subsequent time dependent analysis of the performance of the geothermal doublet.
Step 2: Time Dependent
1
In the Study toolbar, click  Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
From the Time unit list, choose a.
4
In the Output times text field, type range(0,0.2,10).
5
In the Model Builder window, click Study 1.
6
In the Settings window for Study, locate the Study Settings section.
7
Clear the Generate default plots checkbox.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
4
From the Maximum step constraint list, choose Constant.
Constraining the maximum time step ensures that the transient behavior of the whole system is resolved properly.
5
In the Study toolbar, click  Compute.
Result Templates
Add the pressure plot from Result Templates, and modify it to create Figure 2.
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Study 1/Solution 1 (sol1) > Darcy’s Law > Pressure (dl).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Selection 1
1
In the Model Builder window, expand the Pressure (dl) node.
2
Right-click Surface and choose Selection.
3
In the Settings window for Selection, locate the Selection section.
4
From the Selection list, choose All boundaries, and remove the top and front boundaries, which corresponds to:
5
Select Boundaries 3, 6, 9, 11–16, 18, 20, 21, 23, and 25–29 only.
6
In the Pressure (dl) toolbar, click  Plot.
Temperature and Flow Field
Create the plot in Figure 3 as follows:
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Temperature and Flow Field in the Label text field.
Surface 1
1
Right-click Temperature and Flow Field and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type T.
4
Locate the Coloring and Style section. From the Color table list, choose HeatCameraLight.
Selection 1
1
Right-click Surface 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Outer boundaries.
4
Select Boundaries 3, 11–13, 15, 16, and 25–29 only.
Surface 2
1
In the Model Builder window, under Results > Temperature and Flow Field right-click Surface 1 and choose Duplicate.
2
In the Settings window for Surface, click to expand the Inherit Style section.
3
From the Plot list, choose Surface 1.
Selection 1
1
In the Model Builder window, expand the Surface 2 node, then click Selection 1.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Fracture.
Volume 1
1
In the Model Builder window, right-click Temperature and Flow Field and choose Volume.
2
In the Settings window for Volume, locate the Expression section.
3
In the Expression text field, type T.
4
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
Selection 1
1
Right-click Volume 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Bottom layer.
Temperature and Flow Field
In the Model Builder window, under Results click Temperature and Flow Field.
Streamline Multislice 1
1
In the Temperature and Flow Field toolbar, click  More Plots and choose Streamline Multislice.
2
In the Settings window for Streamline Multislice, locate the Multiplane Data section.
3
Find the X-planes subsection. In the Planes text field, type 2.
4
Find the Y-planes subsection. In the Planes text field, type 0.
5
Find the Z-planes subsection. From the Entry method list, choose Coordinates.
6
In the Coordinates text field, type depth_w.
7
Locate the Streamline Positioning section. From the Positioning list, choose Uniform density.
8
Locate the Coloring and Style section. Find the Point style subsection. From the Type list, choose Arrow.
9
Click to expand the Quality section. From the Evaluation settings list, choose Manual.
10
From the Recover list, choose Everywhere.
11
In the Temperature and Flow Field toolbar, click  Plot.
This recovers the accurate derivatives for the streamlines.
Color Expression 1
1
Right-click Streamline Multislice 1 and choose Color Expression.
2
In the Settings window for Color Expression, locate the Coloring and Style section.
3
From the Color table list, choose GrayScale.
Streamline Multislice 1
1
In the Model Builder window, click Streamline Multislice 1.
2
In the Settings window for Streamline Multislice, locate the Coloring and Style section.
3
Find the Point style subsection.
4
Select the Scale factor checkbox. In the associated text field, type 4000000.
Temperature and Flow Field
1
In the Model Builder window, click Temperature and Flow Field.
2
In the Settings window for 3D Plot Group, click to expand the Title section.
3
From the Title type list, choose Manual.
4
In the Title text area, type Temperature distribution and Darcy's velocity field.
5
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
6
Locate the Color Legend section. Select the Show units checkbox.
7
In the Temperature and Flow Field toolbar, click  Plot.
Evaluate the production temperature and compare with Figure 4.
Line Average 1
1
In the Results toolbar, click  More Derived Values and choose Average > Line Average.
2
In the Settings window for Line Average, click Replace Expression in the upper-right corner of the Expressions section. From the menu, choose Component 1 (comp1) > Heat Transfer in Porous Media > Temperature > T - Temperature - K.
3
Select Edge 37 only.
4
Click  Evaluate.
Table 1
1
Go to the Table 1 window.
2
Click the Table Graph button in the window toolbar.
Results
Production Temperature
1
In the Model Builder window, under Results click 1D Plot Group 3.
2
In the Settings window for 1D Plot Group, type Production Temperature in the Label text field.
Temperature and Flow Field
To get a plot showing the pressure and velocity distribution along the fracture, add a plot from the Result Templates window and slightly adjust it to look like Figure 5.
Result Templates
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Study 1/Solution 1 (sol1) > Darcy’s Law > Pressure Along Fracture (dl).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Streamline Surface 1
1
In the Pressure Along Fracture (dl) toolbar, click  More Plots and choose Streamline Surface.
2
In the Settings window for Streamline Surface, locate the Surface Selection section.
3
From the Selection list, choose Fracture.
4
Locate the Streamline Positioning section. From the Positioning list, choose Magnitude controlled.
5
In the Minimum density level text field, type 0.
6
In the Maximum density level text field, type 7.
7
Locate the Coloring and Style section. Find the Point style subsection. From the Color list, choose White.
8
In the Pressure Along Fracture (dl) toolbar, click  Plot.