Geothermal Doublet

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.

Subsurface_Flow_Module/Heat_Transfer/geothermal_doublet

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

In the tree, select .

Click .

Click

Root

Start with loading the parameterized geometry sequence into the model.

Geometry 1

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file geothermal_doublet_geom_sequence.mph.

In the toolbar, click

The geometry parameters that are used were added to the list automatically. Add a few more parameters that are used to set up the physics interfaces.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
r_bore	1[m]	1 m	Borehole skinzone 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

In the window, expand the node.

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Click the

button in the toolbar.

Select Edge 19 only.

In the text field, type Injection well.

Production well

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Edge 37 only.

In the text field, type Production well.

Top layer

In the toolbar, click

Select Domains 3 and 6 only.

In the window for , type Top layer in the text field.

Middle layer

In the toolbar, click

Select Domains 2 and 5 only.

In the window for , type Middle layer in the text field.

Bottom layer

In the toolbar, click

Select Domains 1 and 4 only.

In the window for , type Bottom layer in the text field.

Fracture

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 21 only.

Select the check box. This automatically adds all boundaries that belong to the large fault zone.

In the text field, type Fracture.

Outer boundaries

In the toolbar, click

In the window for , type Outer boundaries in the text field.

Locate the section. From the list, choose .

Select the check box.

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

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Materials

Upper Aquitard

In the window, under right-click and choose .

In the window for , type Upper Aquitard in the text field.

Aquifer

Right-click and choose .

In the window for , type Aquifer in the text field.

Locate the section. From the list, choose .

Lower aquitard

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Lower aquitard.

Fracture

Right-click and choose .

In the window for , type Fracture in the text field.

Locate the section. From the list, choose .

From the list, choose .

Darcy’s Law (dl)

Now, set up the physics interfaces. This determines which material properties are required.

Fracture Flow 1

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Fluid and Fracture Properties 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

In the ff text field, type f_f.

Aperture 1

In the window, click .

In the window for , locate the section.

In the df text field, type d_f.

Hydraulic Head 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the H0 text field, type -deltaH*x.

Use the feature to define the injection and production side of the doublet.

Well 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the dw text field, type 2*r_bore.

From the list, choose .

Locate the 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

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the dw text field, type 2*r_bore.

From the list, choose .

Locate the section. In the M0 text field, type pump*dl.rho.

Heat Transfer in Porous Media (ht)

Fluid 1

In the window, under click .

In the window for , locate the section.

From the u list, choose .

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Initial Values 1

Define the geothermal gradient as initial temperature distribution.

In the window, under click .

In the window for , locate the section.

In the T text field, type T_top-0.03[K/m]*z.

Open Boundary 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Tustr text field, type T_top-delta_Tz*z.

Fracture 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose . In the Lth text field, type d_f.

Locate the section. From the list, choose .

Locate the section. From the u list, choose .

Locate the section. From the θp list, choose .

Line Heat Source 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Ql text field, type mat1.def.Cp*dl.well1.Ml*(T_inj-T).

Locate the section. Select the check box.

In the R text field, type r_bore.

Materials

Now you can define the material properties.

Upper Aquitard (pmat1)

In the window, under click .

In the window for , locate the section.

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

Locate the section. Click

Solid 1 (pmat1.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 0.9.

From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1300	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	2	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	900	J/(kg·K)	Basic

Aquifer (pmat2)

In the window, under click .

In the window for , locate the section.

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

Locate the section. Click

Solid 1 (pmat2.solid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

In the θs text field, type 0.6.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1900	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	3	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	850	J/(kg·K)	Basic

Lower aquitard (pmat3)

In the window, under click .

In the window for , locate the section.

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

Locate the section. Click

Solid 1 (pmat3.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 0.7.

From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	2300	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	3.5	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	850	J/(kg·K)	Basic

Fracture (pmat4)

In the window, under click .

In the window for , locate the section.

In the εp text field, type 0.6.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	3	W/(m·K)	Basic
Density	rho	1.2e3	kg/m³	Basic
Heat capacity at constant pressure	Cp	800	J/(kg·K)	Basic