Coupling Heat Transfer with Subsurface Porous Media Flow

Phillip Oberdorfer April 24, 2014
Share this on Facebook Share this on Twitter Share this on Google+ Share this on LinkedIn

In the second part of our Geothermal Energy series, we focus on the coupled heat transport and subsurface flow processes that determine the thermal development of the subsurface due to geothermal heat production. The described processes are demonstrated in an example model of a hydrothermal doublet system.

Deep Geothermal Energy: The Big Uncertain Potential

One of the greatest challenges in geothermal energy production is minimizing the prospecting risk. How can you be sure that the desired production site is appropriate for, let’s say, 30 years of heat extraction? Usually, only very little information is available about the local subsurface properties and it is typically afflicted with large uncertainties.

Over the last decades, numerical models became an important tool to estimate risks by performing parametric studies within reasonable ranges of uncertainty. Today, I will give a brief introduction to the mathematical description of the coupled subsurface flow and heat transport problem that needs to be solved in many geothermal applications. I will also show you how to use COMSOL software as an appropriate tool for studying and forecasting the performance of (hydro-) geothermal systems.

Governing Equations in Hydrothermal Systems

The heat transport in the subsurface is described by the heat transport equation:


(\rho C_p)_{eq} \frac{\partial T}{\partial t} + \rho C_p {\bf u } \cdot \nabla T = \nabla \cdot (k_{eq} \nabla T ) + Q + Q_{geo}

Heat is balanced by conduction and convection processes and can be generated or lost through defining this in the source term, Q. A special feature of the Heat Transfer in Porous Media interface is the implemented Geothermal Heating feature, represented as a domain condition: Q_{geo}.

There is also another feature that makes the life of a geothermal energy modeler a little easier. It’s possible to implement an averaged representation of the thermal parameters, composed from the rock matrix and the groundwater using the matrix volume fraction, \theta, as a weighting factor. You may choose between volume and power law averaging for several immobile solids and fluids.

In the case of volume averaging, the volumetric heat capacity in the heat transport equation becomes:


(\rho C_p )_{eq} = \sum_{i} ( \theta_{pi}\rho_{pi}C_{p,pi})+(1-\sum_{i}\theta_{pi})\rho C_p

and the thermal conductivity becomes:


k_{eq}=\sum_{i} \theta_{pi} k_{pi} + ( 1-\sum_{i} \theta_{pi} ) \rho C_p

Solving the heat transport properly requires incorporating the flow field. Generally, there can be various situations in the subsurface requiring different approaches to describe the flow mathematically. If the focus is on the micro scale and you want to resolve the flow in the pore space, you need to solve the creeping flow or Stokes flow equations. In partially saturated zones, you would solve Richards’ equation, as it is often done in studies concerning environmental pollution (see our past Simulating Pesticide Runoff, the Effects of Aldicarb blog post, for instance).

However, the fully-saturated and mainly pressure-driven flows in deep geothermal strata are sufficiently described by Darcy’s law:


{\mathbf u} = -\frac{\kappa}{\mu} \nabla p

where the velocity field, \mathbf{u}, depends on the permeability, \kappa, the fluid’s dynamic viscosity, \mu, and is driven by a pressure gradient, p. Darcy’s law is then combined with the continuity equation:


\frac{\partial}{\partial t} (\rho \epsilon_p) + \nabla \cdot ( \rho {\bf u} ) = Q_m

If your scenario concerns long geothermal time scales, the time dependence due to storage effects in the flow is negligible. Therefore, the first term on the left-hand side of the equation above vanishes because the density, \rho, and the porosity, \epsilon_p, can be assumed to be constant. Usually, the temperature dependencies of the hydraulic properties are negligible. Thus, the (stationary) flow equations are independent of the (time-dependent) heat transfer equations. In some cases, especially if the number of degrees of freedom is large, it can make sense to utilize the independence by splitting the problem into one stationary and one time-dependent study step.

Fracture Flow and Poroelasticity

Fracture flow may locally dominate the flow regime in geothermal systems, such as in karst aquifer systems. The Subsurface Flow Module offers the Fracture Flow interface for a 2D representation of the Darcy flow field in fractures and cracks.

Hydrothermal heat extraction systems usually consist of one or more injection and production wells. Those are in many cases realized as separate boreholes, but the modern approach is to create one (or more) multilateral wells. There are even tactics that consist of single boreholes with separate injection and production zones.

Note that artificial pressure changes due to water injection and extraction can influence the structure of the porous medium and produce hydraulic fracturing. To take these effects into account, you can perform poroelastic analyses, but we will not consider these here.

COMSOL Model of a Hydrothermal Application: A Geothermal Doublet

It is easy to set up a COMSOL Multiphysics model that features long time predictions for a hydro-geothermal application.

The model region contains three geologic layers with different thermal and hydraulic properties in a box with a volume V≈500 [m³]. The box represents a section of a geothermal production site that is ranged by a large fault zone. The layer elevations are interpolation functions from an external data set. The concerned aquifer is fully saturated and confined on top and bottom by aquitards (impermeable beds). The temperature distribution is generally a factor of uncertainty, but a good guess is to assume a geothermal gradient of 0.03 [°C/m], leading to an initial temperature distribution T0(z)=10 [°C] – z·0.03 [°C/m].

Geometry of a hydrothermal doublet system in a layered subsurface domain created with COMSOL Multiphysics
Hydrothermal doublet system in a layered subsurface domain, ranged by a fault zone. The edge is about 500 meters long. The left drilling is the injection well, the production well is on the right. The lateral distance between the wells is about 120 meters.

COMSOL Multiphysics creates a mesh that is perfectly fine for this approach, except for one detail — the mesh on the wells is refined to resolve the expected high gradients in that area.

Screenshot of a mesh of the Hydrothermal doublet system

Now, let’s crank the heat up! Geothermal groundwater is pumped (produced) through the production well on the right at a rate of 50 [l/s]. The well is implemented as a cylinder that was cut out of the geometry to allow inlet and outlet boundary conditions for the flow. The extracted water is, after using it for heat or power generation, re-injected by the left well at the same rate, but with a lower temperature (in this case 5 [°C]).

The resulting flow field and temperature distribution after 30 years of heat production are displayed below:

Plot showing the results of the Hydrothermal doublet system after 30 years heat production
Result after 30 years of heat production: Hydraulic connection between the production and injection zones and temperature distribution along the flow paths. Note that only the injection and production zones of the boreholes are considered. The rest of the boreholes are not implemented, in order to reduce the meshing effort.

The model is a suitable tool for estimating the development of a geothermal site under varied conditions. For example, how is the production temperature affected by the lateral distance of the wells? Is it worthwhile to reach a large spread or is a moderate distance sufficient?

This can be studied by performing a parametric study by varying the well distance:

Graph showing the production temperature after stationary conditions have been reached
Flow paths and temperature distribution between the wells for different lateral distances. The graph shows the production temperature after reaching stationary conditions as a function of the lateral distance.

With this model, different borehole systems can easily be realized just by changing the positions of the injection/production cylinders. For example, here are the results of a single-borehole system:

Plot showing the results of a single-borehole approach after 30 years of heat production
Results of a single-borehole approach after 30 years of heat production. The vertical distance between the injection (up) and production (down) zones is 130 meters.

So far, we have only looked at aquifers without ambient groundwater movement. What happens if there is a hydraulic gradient that leads to groundwater flow?

The following figure shows the same situation as the figure above, except that now there is a hydraulic head gradient of \nablaH=0.01 [m/m], leading to a superposed flow field:

Plot showing the results of a single-borehole approach results after 30 years of overlapping groundwater flow and heat production
Single borehole after 30 years of heat production and overlapping groundwater flow due to a horizontal pressure gradient.

Further Reading



  1. Bastiano Deschmann May 30, 2015   8:45 am

    Where can I find the model file of this geothermal doublet?

  2. Bastiano Deschmann May 30, 2015   8:46 am

    Where can I find the model file of this example?

  3. Phillip Oberdorfer June 2, 2015   5:58 am

    Dear Bastiano,

    thank you for your interest! I can send you the model file. Please send a request to and mention this blog (and my name).

    Kind regards

  4. Musa Aliyu July 9, 2015   1:09 pm

    Hoping is well with you? I’ve seen one of your models in COMSOL Blog with title “Coupling Heat Transfer with Subsurface Porous Media Flow” more especially “COMSOL Model of a Hydrothermal Application: A Geothermal Doublet”. I was wondering if you can help with the steps involve in the creation of the model or the COMSOL files; because I’m going to attend one of your free workshop in London next month. So, I’ll try it with them and if it works out it’ll be a very good step in convincing my university to get the product for our research team. Thank you and waiting for your response.

  5. GyeongHun Kim April 14, 2016   12:22 am

    Where can i find the model file? Can i get it here?

  6. Phillip Oberdorfer June 2, 2016   4:53 am

    The model is available in the meantime. Just follow this link:

  7. Martin Hasal October 4, 2016   3:04 am

    I was trying to find a coupling of Richardson equation (unsaturated flow) + Heat transfer in COMSOL, but I was not succesfull. Does anybody know about or better can share .mph file with this problematic?

  8. Caty Fairclough October 19, 2016   10:34 am

    Hi Martin,

    Thank you for your interest in this blog post! For your question, I suggest reaching out to our support team.

    Online support center:


  9. Harish Puppala October 24, 2017   12:46 pm

    Hii all, Can some one suggest me any blog or a post to learn how to couple 1D and 2D elements.
    Considering well bore as a 1D element, how can we couple to a 2D porous element.

  10. Phillip Oberdorfer October 25, 2017   2:42 am

    Hi Harish,

    I suggest to use linear extrusion operator for your approach. There is a blog post from my colleague Temesgen who describes the use of that operator in detail:


  11. Harish Puppala October 31, 2017   1:40 pm

    Hi phillip, I have gone through your geothermal doublet model. What if I want to know the variation of temperature within the extraction well. For example if the extraction well is till a depth of 2.5 km, temperature at 2.5 km might be 200 degC. however after extracted fluid reaches the surface, some fraction of heat might be lost. How to capture this?

  12. Phillip Oberdorfer November 1, 2017   7:12 am

    Hi Harish,

    you can take into account and elevate the heat loss through the borehole using the features of the Pipe Flow Module, check out

  13. chenchen lee January 29, 2018   5:02 pm

    Dear Prof. Oberdorfer, I am very interested in this blog of yours. In fact, I am working on flow including both free surface flow and porous flow, I note that you talked about stokes equation for free surface flow and Darcy’s law for porous flow, but there is one order difference between these two laws. How did you couple them together? And do you have some models files to describe similar problems? I have been confused by this problem for a long time. Thanks a lot!

    Best regards

  14. Harish Puppala February 27, 2018   12:55 pm

    Dear oberdorfer, if we maintain very low injection velocity, there may be a chance that the production well may not get hydraulic support from the injection well. How to understand this concept.

  15. Molika Linkolnli March 7, 2018   7:49 am

    Is it possible to use a dual permeability model here?

    That is, fluid flow for both matrix and fracture, they have defined two separate equations and connected to each other with a transfer function and shape factor as a fracture matrix interaction.

    If yes, which modules should be reset?

    Please advise.

  16. Phillip Oberdorfer April 3, 2018   9:03 am

    Dear Chenchen Lee,

    it is easily possible to couple free surface flow and porous media flow in COMSOL Multiphysics. The great thing is that you do not have to worry about the order differences because that is done internally by the software. You can couple the flow fields manually or use the interface Free and Porous Media Flow (fp) which takes care of that. There is a nice and simple tutorial model in our Application Library:

  17. Phillip Oberdorfer April 3, 2018   9:04 am

    Dear Molika,

    we do actually use a dual permeability approach here. As described in the section Fracture flow and Poroelasticity, there is a Fracture Flow boundary condition with re-definable fluid and matrix properties that is applied to the fracture in the model. I think that this is pretty much what you were asking for.

    All fracture flow functionality is available with the Subsurface Flow Module. Please note that the corresponding representation of heat transport in fractures is part of the Heat Transfer Module.

Loading Comments...