Articles | Volume 13, issue 12
Model description paper
23 Dec 2020
Model description paper |  | 23 Dec 2020

HydrothermalFoam v1.0: a 3-D hydrothermal transport model for natural submarine hydrothermal systems

Zhikui Guo, Lars Rüpke, and Chunhui Tao

Herein, we introduce HydrothermalFoam, a three-dimensional hydro-thermo-transport model designed to resolve fluid flow within submarine hydrothermal circulation systems. HydrothermalFoam has been developed on the OpenFOAM platform, which is a finite-volume-based C++ toolbox for fluid-dynamic simulations and for developing customized numerical models that provides access to state-of-the-art parallelized solvers and to a wide range of pre- and post-processing tools. We have implemented a porous media Darcy flow model with associated boundary conditions designed to facilitate numerical simulations of submarine hydrothermal systems. The current implementation is valid for single-phase fluid states and uses a pure-water equation of state (IAPWS-97). We here present the model formulation; OpenFOAM implementation details; and a sequence of 1-D, 2-D, and 3-D benchmark tests. The source code repository further includes a number of tutorials that can be used as starting points for building specialized hydrothermal flow models. The model is published under the GNU General Public License v3.0.

1 Introduction

High-temperature hydrothermal circulation through the ocean floor plays a key role in the exchange of mass and energy between the solid earth and the global ocean (German and Seyfried2014; Elderfield and Schultz1996). It influences the thermal evolution of young oceanic plates (Stein and Stein1994; Theissen-Krah et al.2016), modulates global ocean biogeochemical cycles (German et al.2016; Tagliabue et al.2010), and is associated with massive sulfide ore deposits that form around vent sites (Hannington et al.2011). Hydrothermal convection occurs over large spatial and temporal scales. At fast-spreading ridges, convection cells may either be confined to the upper extrusive crust above the axial melt lens (Faak et al.2015; Coumou et al.2008; Fontaine et al.2009) or extend all the way down to the crust–mantle boundary at approx. 6 km depth (Hasenclever et al.2014; Dunn et al.2000; Cathles1993). At slow-spreading ridges, fluid circulation may extend much deeper (up to 35 km) into the ultramafic mantle (Schlindwein and Schmid2016), although the maximum extent of the brittle layer remains debated and may be confined to the upper 15 km (Grevemeyer et al.2019). Such deep-reaching fluid flow can sometimes be channelized along deep detachments at fault-controlled systems such as Trans-Atlantic Geotraverse (TAG) and Lonqi (Tao et al.2020; deMartin et al.2007) and/or may propagate to greater depths via thermal cracking (Olive and Crone2018; Lister1974). Temperatures can also vary over large ranges, with high-temperature systems typically being driven by a magmatic heat source of 1000 C or more and porosity and permeability staying open up to 600–800 C (Lister1974). Finally, hydrothermal systems can evolve over long timescales of up to 50–100 kyr (Jamieson et al.2014) but also respond to shorter events like glacial sea level changes (Middleton et al.2016) or magmatic and seismic events (Germanovich et al.2000; Wilcock2004; Singh and Lowell2015) and even tidal pressure changes (Crone et al.2011; Barreyre et al.2018). These spatial and temporal scales, in combination with the extreme pressure (P) and temperature (T) conditions (up to 300 MPa and 1000 C) of submarine hydrothermal systems, make direct and long-term observations challenging and pose a problem for laboratory work (Ingebritsen et al.2010). Hence, numerical simulations have become indispensable tools for understanding and characterizing fluid flow and for relating seafloor observations to physicochemical processes at depth.

In the last few decades, significant progress has been made in hydrothermal flow modeling both theoretically and numerically (Lowell1991; Ingebritsen et al.2010). Due to the high complexity of the heterogeneous sub-seafloor, a continuum porous medium approach is typically used, in which the conservation equations are written for control volumes with effective properties such as Darcy velocity, permeability, and porosity. Such approaches have been successfully used to make fundamental progress in our understanding of the nature and mechanisms of hydrothermal transport, including a thermodynamic explanation of black smoker temperatures (Jupp and Schultz2000), the three-dimensional structure of hydrothermal circulation cells at mid-ocean ridges (Coumou et al.2009; Hasenclever et al.2014), and phase separation phenomena and salinity variations of hydrothermal fluids (Lewis and Lowell2009a, b; Coumou et al.2009; Weis et al.2014).

Current numerical simulators of hydrothermal flow can be divided into two families: (1) multiphase codes that thrive towards resolving saltwater convection and associated phase separation phenomena and (2) single-phase hydrothermal codes that focus on sub-critical low-temperature fluid flow and/or super-critical high-temperature flow of pure water, i.e., codes that only “work” within single-phase fluid states. Multiphase saltwater codes are at the forefront of what is currently feasible in numerical simulations, as accounting for the complexity of the equation-of-state (EOS) of seawater (Driesner and Heinrich2007; Driesner2007) in combination with multiphase transport is a challenge (Ingebritsen et al.2010). Existing codes of this type include CSMP++, which is capable of treating salt water up to magmatic temperatures on unstructured finite-element–finite-volume (FEFV) meshes (Weis et al.2014). The hydrothermal multiphase version of CSMP++ is currently 2-D and it is a closed-source project. FISHES is a 2-D open-source academic code, which uses the finite-volume method to solve thermohaline convection on structured meshes (Lewis and Lowell2009a, b) but has some restrictions on the phase states that can be resolved. Currently there is no 3-D model that can resolve multiphase saltwater convection but there are some developments efforts on the way. In addition, there are a number of geothermal modeling codes that can handle two-phase behavior that have not (yet) been adapted to handle the complex EOS of saltwater over sufficiently large pressure and temperature ranges. HYDROTHERM (Kipp et al.2008), FEHM (Zyvoloski et al.1997), HT2_NR (Vehling et al.2018), and TOUGH2 (Pruess et al.1999) are examples of such codes. The second type of code family refers to somewhat simpler models that circumvent the numerical challenges of multiphase phenomena by staying in P-T regions, where the simulated fluid is in a single phase. A popular approach is to use a pure water instead of a saltwater EOS at pressures beyond the critical end-point (22 MPa). These models, despite making simplifying assumptions, continue to be widely used in the submarine hydrothermal system community and have been successfully applied to solve a wide range of problems. Examples include Jupp and Schultz (2000, 2004), who showed that hydrothermal systems operate close to optimal efficiency with their maximum vent temperatures set by the thermodynamic properties of water, studies that revealed the complex 3-D structure of recharge and discharge flow in mid-ocean-ridge hydrothermal systems (Hasenclever et al.2014; Coumou et al.2008; Fontaine et al.2014), dedicated case studies for individual vent systems (Tao et al.2020; Andersen et al.2015; Lowell et al.2012), and models exploring tidal forcing of hydrothermal circulation (Crone and Wilcock2005; Barreyre et al.2018). This list is nowhere near complete and there are many more examples. The bottom line is that single-phase circulation models continue to be widely used and highly useful “workhorses” in the hydrothermal community. Somewhere in the hopefully not so distant future, 2-D and 3-D multiphase models will be the new standard but for now robust and tested single-phase codes continue to be useful tools for a variety of applications.

Interestingly, even single-phase models are not that easily accessible to the hydrothermal community. Many research groups maintain 2-D research codes that resolve hydrothermal flow but single-phase 3-D models continue to be rare. To our knowledge there are basically three single-phase code families that are routinely used in 3-D studies (Coumou et al.2008; Hasenclever et al.2014; Fontaine et al.2014) and none of them are open source. There are some major open-source initiatives that provide 3-D porous flow simulators or libraries such as RichardsFoam2 (Orgogozo2015), porousMultiphaseFoam (Horgue et al.2015), Dumux (Flemisch et al.2011), MRST (Lie2019), and OpenGeoSys (Kolditz et al.2012) that can be used to simulate hydrothermal flow but none of them has been adapted and documented for simulating submarine hydrothermal systems. In this paper, we present a toolbox, named HydrothermalFoam, to simulate 2-D and 3-D hydrothermal circulation in a single-phase regime for seafloor hydrothermal systems. The toolbox is build upon the open-source platform OpenFOAM® (Jasak1996; Weller et al.1998), which is not only a widely used simulator for solving Navier–Stokes-type problems but also a general toolbox for solving partial differential equations. OpenFOAM is based on the cell-centroid finite-volume method (FVM) and is written in C++. It provides high-level interfaces to field operations and includes a series of features such as support for flexible meshes (e.g., structured meshes, unstructured meshes, and mixed meshes), utilities of pre- and post-processing, and parallel computing in 2-D and 3-D (Moukalled et al.2016). Based on this established framework, we present a toolbox to simulate flow in submarine hydrothermal systems. We solve the porous flow problem using a continuum porous medium approach in which the fluid velocity is expressed by Darcy's law and the pressure equation is constructed from Darcy's law and the mass conservation equation. All the partial differential equations are solved implicitly in the framework of OpenFOAM and in a sequential scheme, and the thermal–physical models are developed using a pure-water EOS. HydrothermalFoam inherits all kinds of basic features of OpenFOAM, including boundary conditions. In addition, we have also customized several special boundary conditions for seafloor hydrothermal system modeling. The purpose of this toolbox is to provide the hydrothermal community with a state-of-the-art, yet easy-to-use and well-documented, simulator for resolving hydrothermal flow in submarine hydrothermal systems.

The paper is organized as follows. In Sect. 2, we present the mathematical model, its implementation in OpenFOAM, information on initial and boundary conditions, and the thermal–physical model selection. In Sect. 3 we describe the different toolbox components, in Sect. 4 we describe the installation options and procedures, and in Sect. 5 we validate the toolbox using several published benchmark tests.

Table 1Definitions and values of variables used in this study.

Download Print Version | Download XLSX

2 Model development

2.1 Mathematical model of hydrothermal flow

We use a continuum porous media approach and describe creeping flow in hydrothermal circulation systems using Darcy's law, where the Darcy velocity of the fluid is given by

(1) U = - k μ f ( p - ρ g ) ,

in which k denotes permeability, μf the fluid's dynamic viscosity, p total fluid pressure, and g gravitational acceleration. All variables and symbols are listed in Table 1. Considering a compressible fluid in a porous medium with given porosity structure, the mass balance is expressed by Eq. (2) (Theissen-Krah et al.2011),

(2) ε ρ f t + ( U ρ f ) = 0 ,

where ε is the porosity of the rock. Note that we assume the matrix to be incompressible, and thus the porosity is outside the time derivative. The equation for pressure can be derived by substituting Darcy's law (Eq. 1) into the continuity Eq. (2) and treating the fluid's density as a function of temperature T and pressure p:

(3) ε ρ f β f p t - α f T t = ρ f k μ f ( p - ρ f g ) ,

with αf and βf being the fluid's thermal expansivity and compressibility, respectively. Again there is no rock compressibility, as we consider the incompressible matrix case. Energy conservation of a single-phase fluid can be expressed using a temperature formulation (Hasenclever et al.2014)

(4) ( ε ρ f C pf + ( 1 - ε ) ρ r C pr ) T t = ( k r T ) - ρ f C pf U T + μ f k U 2 - ln ρ f ln T p ε p t + U p ,

where Cp is heat capacity and kr is the bulk thermal conductivity of porous rock. Subscripts of r and f refer to rock and fluid, respectively. As the matrix is incompressible, ρr and Cpr are constant in time. Fluid and rock are assumed to be in local thermal equilibrium (T=Tr=Tf) so that the mixture appears on the left-hand side of equation (4). Note that the assumption of thermal equilibrium is valid for most practical applications in submarine hydrothermal system modeling, as the equilibration timescale is short being related to grain size in a porous medium. However, this assumption should be carefully reviewed when simulating special geometries like fluid-filled cracks (Schmeling et al.2018) or very short timescales like the response of a hydrothermal circulation cell to a seismic event (Wilcock2004). For such specialized cases OpenFOAM offers support for multi-physics models that can resolve different physics in differing parts of the modeling domain, and thus heat transfer between solid and fluid can be explicitly resolved. Under the equilibrium assumption, changes in temperature depend on conductive heat transport, advective heat transport by fluid flow, heat generation by internal friction of the fluid, and pressure–volume work. Note that the third term on the right-hand side of Eq. (4) describes viscous dissipation, while the last describes pressure–volume work (and the pressure dependence of enthalpy, which shows up in the temperature formulation of the energy equation). Whether these terms are important depends on the application (Garg and Pritchett1977); our tests have shown that they matter, for example, in the pure vapor state, when simulating large vertical extents, and in supercritical states close to the critical end point of water. All fluid properties are functions of both pressure and temperature and are calculated using the IAPWS-IF97 formulation of water and steam properties as implemented in the freesteam project (Pye2010). Further details on the derivation of the governing equations can be found in the appendix of Hasenclever et al. (2014).

2.2 Implemented formulation

We solve for pressure (Eq. 3), velocity (Eq. 1), and temperature (Eq. 4) separately. Based on the finite-volume method implemented in OpenFOAM, the primary variables (p and T) related equations are discretized on an cell-centroid computational grid. The transient temperature term in Eq. (3) is evaluated explicitly and is treated numerically as a source term resulting in a Poisson-type equation

(5) ε ρ f β f p t - ρ f k μ f p = ε ρ f α f T t - ρ f k μ f ρ f g ,

where the left-hand side terms are pressure transient term and Laplacian term (or diffusion term), respectively. The first term on the right-hand side is evaluated explicitly using a known temperature field, and the second term on the right-hand side is a divergence term of gravity-related flux (ϕg), which is defined on each face of the computational grid.

To apply the finite-volume method, the advection term (the second term on the right-hand side) in the temperature Eq. (4) should be reformulated as a divergence term

(6) ρ f C pf U T = ( ρ f C pf U T ) - T ( ρ f C pf U ) .

Then substituting Eq. (6) in Eq. (4), the temperature equation can be rearranged as

(7) ( ε ρ f C pf + ( 1 - ε ) ρ r C pr ) T t + ( ρ f C pf U T ) = ( k r T ) + T ( ρ f C pf U ) + μ f k U 2 + T α f ε p t + U p ,

where on the left-hand side, the first term is temperature transient term and the second one is the advection term. On the right-hand side, the first two terms represent temperature diffusion and the source term resulting from the reformulation given in Eq. (6), respectively. The last two source terms are calculated explicitly.

Figure 1Phase diagram and density of pure water in temperature–pressure space. Red dot denotes supercritical point. The supercritical fluid region is outlined by critical pressure and critical temperature lines and is divided by critical density isoline into liquid-like fluid and vapor-like fluid. Three solid curves in different colors represent pressure-temperature path of 1-D benchmark examples shown in Fig. 4a, c, and e, respectively.


2.3 Time step limitations

To determine the time step, we adopt the limitation related to the Courant number Co, which is defined for compressible fluid as

(8) C o = max cell F = 1 m ϕ F 2 ρ f V cell Δ t ,

with m being the number of neighbor faces to a specific cell, ϕF=ρfUSF being mass flow through the cell faces, Vcell being volume of the cell, and Sf being area of face F. Then the coefficient for time step change is written as follows:

(9) C Δ t = C o max C o .

To avoid changes of the time step that are too large, which could lead to numerical instabilities, the time step is defined as follows (Jasak1996)

(10) Δ t = min ( min ( C Δ t , 1 + 0.1 C Δ t ) , 1.2 ) Δ t last ,

where Δtlast is the time step length of the previous time step. Implementation details can be found in the OpenFOAM documentation and the OpenFOAM source files included by the main source code file HydrothermalSinglePhaseDarcyFoam.C.

2.4 Boundary conditions

To solve the pressure and temperature equations, we have to impose suitable boundary conditions for T and p. The “typical” boundary conditions, e.g., fixed value, fixed gradient, and the mixing of both, are directly inherited from the basic boundary conditions of OpenFOAM. In submarine hydrothermal system modeling,some special adaptations of these basic boundary conditions can also be useful.

The hydrothermal heat flux (qh) boundary condition is a fixed gradient boundary condition that is often used to approximate heat input from a crustal magma chamber and is commonly used for simulations of mid-ocean ridge hydrothermal systems (e.g., Coumou et al.2009; Weis et al.2014). This Neumann boundary condition is called hydrothermalHeatFlux in the toolbox and can be used for the temperature field. Using it, the imposed gradient of temperature can be expressed as

(11) n T = - n q h k r ,

where n denotes the normal vector of the face boundary. In addition, we implement two options for heat flux distribution. Using the keyword shape allows for modifying the functional form of the heat flux boundary. The available options are fixed, gaussian2d, gaussian3d. The default option is fixed, and if gaussian2d or gaussian3d is specified, the Gaussian shape (Eq. 12) related parameters (qmin, qmax, c, x0 and/or z0) have to be specified (see Sect. 5.3).

(12) q h ( x , z ) = q min + ( q max - q min ) e - ( x - x 0 ) 2 + ( y - z 0 ) 2 2 c 2

A similar boundary condition called hydrothermalMassFluxPressure is defined for the pressure field to prescribe a mass influx into the modeling domain (ϕm=ρfU). The corresponding gradient of the pressure field can be derived from Darcy's law (Eq. 1)

(13) n p = n ρ f g - μ f k U = μ f ρ f k n ( ϕ g - ϕ m ) .

Where ϕm and ϕg denote mass flux and gravity-related flux, respectively. Further, we define another Neumann boundary condition (named noFlux) for the pressure field for impermeable boundaries, which is a special case of HydrothermalMassFluxPressure when ϕm=0. In addition, a Dirichlet boundary condition (named submarinePressure) for the pressure field is defined to describe hydrostatic pressure at the seafloor boundary due to bathymetric relief. Another commonly used boundary condition on hydrothermal venting boundary (e.g., seafloor) is OpenFOAM's inletOutlet boundary condition, which allows for the setting of a constant temperature for inflow and zero heat flux for outflow nodes – this type of boundary condition is often used to mimic free venting at the seafloor.

2.5 Fluid properties and equation-of-state

Numerical solutions of hydrothermal flow are known to strongly depend on the used thermodynamic properties of the simulated fluid. A series of studies using realistic thermodynamic properties of pure and salt water, rather than making a Boussinesq approximation or using linearized properties, have shown that realistic results depend critically on using a realistic EOS (Jupp and Schultz2000; Hasenclever et al.2014; Driesner2010; Carpio and Braack2012). Note that we here do not address any issues related to using pure versus saltwater EOSs, as outlined in the introduction. We use an EOS for pure water based on the IAPWS-IF97 parameterization and have created a corresponding OpenFOAM thermo-physical model in a single phase regime. The phase diagram is shown in Fig. 1.

Figure 2Schematic of the sequential algorithm.


2.6 Solution algorithm

The governing equations of pressure and temperature are solved in a sequential approach. The primary variables (pressure and temperature) and transport properties (such as permeability, porosity, etc.) have to be initialized before the time loop, and then the initial Darcy velocity and thermodynamic properties of fluid can be updated according to the temperature and pressure fields. The main computational sequence for a single time step is described below and sketched in Fig. 2.

Listing 1Implementation of temperature Eq. (7) with OpenFOAM (in EEqn.H).


Listing 2Implementation of pressure Eq. (5) with OpenFOAM (in pEqn.H).


  1. The time step size Δtn+1 is calculated from the condition related to the Courant number (Eq. 8).

  2. Temperature field Tn+1 is implicitly computed by solving the energy conservation equation (Eq. 7). The syntax of a partial differential equation (PDE) in OpenFOAM is very close to mathematical formulation, and a code snippet of the temperature Eq. (7) implementation is shown in Listing 1. All variable symbols, names, and OpenFOAM types (classes) are shown in Table 1. The transient term T/t can be implicitly discretized using the OpenFOAM operator fvm::ddt(T) with the discretization scheme (e.g., Euler scheme) being specified under the keyword ddtSchemes in the system/fvSchemes dictionary file (see Listing A1). The divergence term, Laplace term, and source term are implicitly discretized using fvm::div, fvm::laplacian and fvm::Sp operators. The last term on the right-hand side is explicitly calculated using known field values from the current or previous time step; the corresponding time derivative and gradient can be programmed using fvc::ddt and fvc::grad, respectively.

  3. The pressure field pn+1 is implicitly computed by solving the pressure equation (Eq. 5), the code snippet is shown in Listing 2. The temperature temporal term and divergence of ϕg on the right-hand side are evaluated explicitly by using fvc::ddt(T) and fvc::div(phig) (see line 7 in Listing 2). Although pressure boundary conditions are customized by flux directly (see Sect. 2.4), in order to specify pressure boundary conditions through velocity boundary conditions, e.g., OpenFOAM's fixedFluxPressure boundary condition, the OpenFOAM function of constrainPressure has to be called before solving the pressure equation (see line 3 in Listing 2). For non-orthogonal mesh, a non-orthogonal correction algorithm (line 4 in Listing 2) is commonly adopted to improve accuracy for gradient computation. The number of non-orthogonal corrections is specified by the nNonOrthogonalCorrectors key in the PIMPLE sub-dictionary in the system/fvSolution file.

  4. The velocity field is calculated explicitly using the latest pressure field based on Darcy's law (Eq. 1). Instead of calculating the velocity directly, we implement an indirect approach based on OpenFOAM's function fvc::reconstruct to reconstruct the velocity field from the computed mass flux (see Listing 3), which performs under higher numerical stability and benefits from the flux conservation characteristics of the finite-volume method. In addition, the boundary conditions of the velocity field have to be updated (line 3 in Listing 3) if OpenFOAM's fixedFluxPressure boundary condition is applied for the pressure field.

  5. Thermodynamic properties of fluid are updated by the thermo-physical model after solving the temperature and pressure field. The implementation code snippet is shown in Listing 4, in which thermo.correct() is used to update the temperature and pressure value for all the calculating nodes. Then the thermodynamic properties of the fluid, for example density (ρ) at each node, are calculated based on IAPWS-IF97 (see line 2–6 in Listing 4).

Listing 3Implementation of Darcy velocity calculation with OpenFOAM (in pEqn.H).


Listing 4Update of fluid thermodynamic properties (in updateProps.H).


2.7 Numerical schemes

Since the numerical evaluation of the divergence and gradient terms in the governing equations has great influence on heat and mass transfer, a suitable solution strategy regarding discretization and linear solver schemes needs to be chosen to ensure accuracy, robustness, and stability. In the presented solver HydrothermalSinglePhaseDarcyFoam, the discretization and interpolation scheme of the primary fields (T,p) can be defined in the simulation configuration files. In the following benchmark tests (Sect. 5), the advective discretization scheme is set to upwind to ensure consistency with HYDROTHERM. It should be noted that all of the basic numerical schemes of OpenFOAM are also valid for HydrothermalSinglePhaseDarcyFoam solver.

3 Description of toolbox components

The organization of the HydrothermalFoam toolbox is shown in Fig. 3. The toolbox contains five parts: a HydrothermalFoam solver, thermo-physical models, boundary conditions, cookbooks, and a manual.

  • HydrothermalSinglePhaseDarcyFoam. This block compiles the solver (an executable file) that solves the seafloor hydrothermal convection equations described in Sect. 2.1. It can be used to simulate single-phase hydrothermal circulation in an isotropic porous medium.

  • ThermoModels. This block compiles the libHydroThermoPhysicalModels library, containing the EOS of pure water, which is used to formulate the used thermo-physical model; see Sect. 2.5.

  • BoundaryConditions. This block compiles libHydrothermalBoundaryConditions library, containing four customized boundary conditions explained in Sect. 2.4. The example usage of each boundary condition can be found in the cookbooks and manual in the following GitLab repository (, last access: 21 December 2020).

  • benchmarks. This block contains the input files of all the benchmark tests (see Sect. 5) presented in this paper.

  • cookbooks. This block contains some example cases of parallel computing, user-defined boundary conditions, and post-processing.

  • manual. The manual presents detailed usage and all the numerical examples in cookbooks.

Figure 3Structure and components of the HydrothermalFoam toolbox.


4 Installation

We provided two options for installation: one is building from source and the other is using a pre-compiled docker image.

4.1 Building from source code

4.1.1 OpenFOAM

The HydrothermalFoam v1.0 is developed based on, which can be installed according to the installation instructions (, last access: 21 December 2020) given by the development team at (last access: 21 December 2020), (last access: 21 December 2020), (last access: 21 December 2020), and (last access: 21 December 2020) for each platform.


Once OpenFOAM is built successfully, the source code of HydrothermalFoam can be downloaded from (Guo and Rüpke2020). The directory structure and components of HydrothermalFoam are shown in Fig. 3 and the components can be built following the three steps given below.

  1. Build a freesteam-2.1 library. The freesteam project is constructed by (last access: 21 December 2020), which is an open-source software construction tool dependent on (last access: 21 December 2020), and based on (last access: 21 December 2020) (GNU Scientific Library, GSL). Therefore, Python 2, SCons, and GSL have to be installed first, followed by changing the directory to freesteam-2.1 in the HydrothermalFoam source code and type command of scons install to compile the freesteam library.

  2. Build libraries of customized boundary conditions and thermo-physical model. Change directory to libraries and type command of ./Allmake to compile the libraries.

  3. Build solver. Change directory to HydrothermalSinglePhaseDarcyFoam and type command of wmake to compile the solver.

All the library files and executable application (solver) file will be generated in directories defined by OpenFOAM's path variables of FOAM_USER_LIBBIN and FOAM_USER_APPBIN, respectively.

4.2 Pre-compiled docker image

In order to use all the tools directly without any compiling and development skills, we have published a pre-compiled Docker® image in a repository at (last access: 21 December 2020). The docker image can be used on any operation system (e.g., Windows, Mac OS, and Linux) to run HydrothermalFoam cases following the five steps below.

  1. Install Docker, then open Docker and keep it running.

  2. Pull the docker image by using command of docker pull zguo/hydrothermalfoam in the shell terminal, e.g., bash shell in Mac OS or PowerShell in Windows system.

  3. Install a container from the docker image by running a shell script, e.g., Unix shell script, shown in Listing 5. The directory named HydrothermalFoam_runs is a shared folder between the container and host machine.

  4. Start the container by running the command docker start hydrothermalfoam.

  5. Attach the container by running the command docker attach hydrothermalfoam. The user is now in a Ubuntu Linux environment with pre-compiled HydrothermalFoam tools that are located at the directory  /HydrothermalFoam. We recommend that the user runs HydrothermalFoam cases in the directory HydrothermalFoam_runs in the container, then the results will be synchronized in the shared directory in the host and thus can be visualized by ParaView®, Tecplot®, or other software.

Following above five steps, one can run HydrothermalFoam tools in the current shell terminal. The following links can be found in the Docker Hub repository: (last access: 21 December 2020) and (last access: 21 December 2020) (5 min tutorial video).

Listing 5Script for installing a container from the HydrothermalFoam docker image.


4.3 Run the first case of HydrothermalFoam

The basic directory structure for a HydrothermalFoam case that contains the mandatory files to run an application is shown in Fig. 4. There is a bash script file named in every HydrothermalFoam case provided by this paper. A HydrothermalFoam case can be run by executing ./ In addition, we provide a 5 min quick start tutorial video (, last access: 21 December 2020) to run the first case of HydrothermalFoam in Docker.

Figure 4Case directory structure of the HydrothermalFoam toolbox.


4.3.1 Mesh generation

The mesh information containing boundary patches definitions, cell face indices, and connections is located in the polyMesh subdirectory in the constant directory in a specific case folder (Fig. 4). All the OpenFOAM mesh generation approaches can be applied to HydrothermalFoam as well. For example, blockMesh generates a simple mesh defined by the blockMeshDict dictionary file in the system directory, and gmshToFoam transforms a mesh file generated by (last access: 21 December 2020) (Geuzaine and Remacle2009) to polyMesh.

4.3.2 Input field data

Much of the input–output data in HydrothermalFoam are fields, e.g., temperature or pressure data, that are read from and written into the time directories. For example, the initial time directory is commonly named 0 (see Fig. 4). HydorthermalFoam writes field data as dictionary files using keyword entries described in Table 2.

Table 2Main keywords used in field dictionary files.

Download Print Version | Download XLSX

The required input field data of HydrothermalFoam are temperature, pressure, and permeability. Each kind of input field data begins with an entry for its dimensions, which is expressed by a vector of seven basic SI units (International System of Units) in the following order: kg, m, s, K, mol, A, and cd. The file names are the same as the variable names shown in Table 1. Following this is the internalField, described in one of the following ways.

  1. Uniform field. A single value is assigned to all elements within the field, taking the form:
    internalField uniform <entry>;

  2. Nonuniform field. Each field element is assigned a unique value from a list, taking the form:
    internalField nonuniform <List>;
    The nonuniform list of an internal field, e.g., permeability, can be specified in different mesh regions by using setFields transform setFieldsDict in the system directory (see the (last access: 21 December 2020) benchmark example described in Sect. 5.3.2).

Listing 6Example dictionary file for temperature field T.


The boundaryField is a dictionary containing a set of entries whose names are listed in the polyMesh/boundary file for each boundary patch. The compulsory entry, type, describes the boundary condition specified for the field. Aside from the OpenFOAM internal basic boundary condition types of fixedValue, zeroGradient, codedFixedValue, etc., the customized boundary condition types, e.g., noFlux for pressure and HydrothermalHeatFlux for temperature, are described in Sect. 2.4. An example set of field dictionary entries for temperature T are shown in Listing 6. While a boundary condition for permeability is not mathematically required, its internal OpenFOAM data type requires a corresponding boundaryField dictionary with specified boundary conditions. Therefore, we suggest that the type entry of all boundary patches for permeability field should be always set to zeroGradient.

4.3.3 Thermo-physical model

The thermophysicalProperties is a compulsory file in the constant directory, it contains keywords (Listing 7) of the newly defined thermo-physical model for water, which is described in Sect. 2.5 and contains constant properties of rock, e.g., porosity and density. It should be noted that porosity is implemented as volScalarField type just like permeability. It means that user can define a porosity field in the modeling domain. If the porosity file does not exist in the start time folder, e.g., the 0 folder, the solver will initialize the porosity field by using the constant porosity value specified in the porousMedia sub-dictionary in the constant/thermophysicalProperties file.

Listing 7thermophysicalProperties file.


4.3.4 Discretized schemes and solution control

The discretization schemes for primary variables in the PDEs (partial difference equations) and solver for linear equations are specified in fvSchemes and fvSolution files in the system directory, respectively. According to the implementation of temperature equation (Listing 1) and pressure equation (Listing 2), we have to specify discretization schemes for the transient terms, Laplace terms, and gradient and divergence terms, which are shown in Listing A1. A example of solver, preconditioner, and tolerance settings for linear equations of temperature and pressure fields are shown in Listing A2. We recommend keeping these two files the same for different cases unless one attempts to try different options available in OpenFOAM.

5 Benchmark tests

We have conducted a number of one-dimensional (1-D), two-dimensional (2-D), and three-dimensional benchmark tests and compared the results to other established software packages to validate HydrothermalFoam and to highlight some of its advantages. The reference software we used is version 3.1 of HYDROTHERM, a simulation tool developed and maintained by the US Geological Survey (USGS), which can be downloaded from the internet (, last access: 21 December 2020) for free. All the parameters used in the 1-D and 2-D examples are taken from Weis et al. (2014), who presented a sequence of well-defined and highly useful benchmarks designed to test code performance within different key thermodynamic fluid states. In those benchmarks the transport properties and rock properties are constant and uniform (values are also listed in Table 1); an isotropic permeability of k=10-15 m2, a porosity of ε=0.1, a heat capacity of Cpr=880 J/(kg C), a thermal conductivity of kr=2 W/(m C), and a rock density of ρr=2700 kg/m3 are used in all simulations below.

5.1 One-dimensional simulations

We conducted six 1-D simulations to test the code performance along the three pT paths in the phase diagram of pure water shown in Fig. 1. These runs are designed with constant pressure and temperature conditions on both ends of a domain with 2 km length and 10 m grid spacing. The boundary conditions and initial conditions of each 1-D test are listed in Table 2. For comparison, we use the same parameters for the HYDROTHERM simulations. The computational domain is oriented horizontally (model index are A, C, E) without gravity and vertically (model index are B, D, F) with gravity to evaluate gravitational effects on fluid flow. All input files can be found in the (last access: 21 December 2020) and (last access: 21 December 2020) directories.

Table 3Model parameters of 1-D benchmark examples.

Download Print Version | Download XLSX

Simulation results of the six 1-D examples are shown in Fig. 5. The example runs A and B describe invasion of a hot fluid into an initially colder domain and the fluids stay in single-phase liquid state. The thermal front moves from the start point (left and bottom for the horizontal and vertical example, respectively) towards the end point. As the vertical flow is opposing gravity, the flow is about 3 times slower. The fluids remain at pressures beyond the critical end point of pure water and therefore in the single-phase regime. Results calculated by HYDROTHERM and HydrothermalFoam are almost identical. In Examples A (horizontal) and B (vertical), the fluid remains in a liquid-like state, while in the examples C (horizontal) and D (vertical) the fluid flows along the pressure gradient from a cold liquid-like state to a hot vapor-like state. In C and D the fluids moves about 2 times faster than in A and B, resulting in a sharper thermal front. The sharpness of this front is, unfortunately, often affected by the numerical scheme (e.g., mesh geometry, upwinding scheme, and advection scheme). In fact, OpenFOAM seems to cope a bit better with resolving the sharpness of the front despite also using an upwind advection scheme. Benchmarks E and F explore sub-critical vapor flow. The results of horizontal and vertical flow look very similar because the density of single-phase vapor fluid is very low and thus the gravitational effects are relatively small with respect to the liquid cases. The results of HydrothermalFoam have a good agreement with that of HYDROTHERM for single-phase vapor flow.

Figure 5Snapshots of one-dimensional benchmark examples in horizontal (a, c, e) and vertical (b, d, f) orientation. Results (temperature in red and pressure in blue) of HydrothermalFoam and HYDROTHERM are plotted as solid lines and dashed lines, respectively.


5.2 Two-dimensional simulations

The two-dimensional models are performed on a rectangular domain with a length of 9 km in the x direction and 3 km in the y direction (Fig. 6), loosely representative of a vertical section through the upper oceanic crust with uniform permeability. The top boundary represents the seafloor and is kept at a constant pressure of 30 MPa, which is equivalent to about 3 km water depth, and a constant temperature of 5 C. At the bottom, a constant heat flux of Qb=0.05 W/m2 is applied. Further, we assume a magmatic heat source with constant heat flux of Qm=5 W/m2 extending 1 km wide along the x direction located at the center of the bottom boundary (shown in red line in Fig. 6).

Figure 6Mesh and model configuration of the two-dimensional examples. (a) The grid for HYDROTHERM has to be set to a uniform regular mesh with resolution of 100 m and a total number of 2700 elements. (b) HydrothermalFoam uses an unstructured mesh with 7074 triangular elements and an edge length in the range of 60–142 m. Red lines on the bottom boundaries denotes heat flux boundary condition for temperature.


The simulation results are shown in Fig. 7. Fluid pressure and temperature calculated by HydrothermalFoam and HYDROTHERM agree very well. The evolution of the hydrothermal plume is summarized in Fig. 7a–c. A hot plume is forming at the base and rises upwards after 5 kyr (Fig. 7a) reaching the seafloor at 15 kyr (Fig. 7b). After 50 kyr, a quasi-steady hydrothermal circulation cell is established (Fig. 7c). The results of HydrothermalFoam and HYDROTHERM are again very similar. All input files can be found in the following directories: (last access: 21 December 2020) and (last access: 21 December 2020).

Figure 7Snapshots of two-dimensional examples. Contour lines for fluid pressure (blue) and temperature (red) calculated from HydrothermalFoam are plotted as solid lines and those from HYDROTHERM are plotted as dotted lines (black). The visualization window is shown by the black box in Fig. 6


5.3 Three-dimensional simulation

5.3.1 Homogeneous model

Similar to the two-dimensional model in Sect. 5.2, the three-dimensional models are performed on a cubic domain with a length of 9 km in both the x direction and the z direction and 3 km in y direction (vertical direction). Note that the vertical coordinate in HydrothermalFoam or OpenFOAM is y rather than z. It can be imagined as representing a three-dimensional section of oceanic crust with uniform permeability. The top boundary is the seafloor and is kept at a constant pressure of 30 MPa and a constant temperature of 5 C. At the bottom boundary, a zero mass flux is applied for pressure, and a constant Gaussian-shaped heat flux (see Eq. 12) is applied for temperature. The corresponding parameters in Eq. (12) are x0=0,z0=0,qmax=5 W/m2, qmin=0.05 W/m2, and c=500. All input files can be found in the 3d/Homogeneous directory in benchmarks/HydrothermalFoam and benchmarks/USGS_HYDROTHERMAL directories, respectively.

The simulation results at 50 kyr are shown in Fig. 8. Vertical slices at x=0 km and x=0.5 km, and horizontal slices at z=0.5 km and z=2.5 km are shown in Fig. 8a–d. Fluid pressure (blue contours) and the temperature field calculated by HydrothermalFoam and HYDROTHERM (dashed contours) agree very well, and the results of the central vertical slice are very close to two-dimensional model results shown in Fig. 7c. Three-dimensional flow path and isothermal surfaces of 300, 200, and 100 C are shown in Fig. 8e.

Figure 8Snapshots of three-dimensional model results. Contour lines for fluid pressure (blue) and temperature calculated from HydrothermalFoam are plotted as solid lines, and those from HYDROTHERM are plotted as dotted black lines. Panels (a) and (b) show pressure and temperature contours on vertical slices, and panels (c) and (d) show temperature contours on horizontal slices. Panel (e) shows three-dimensional isothermal surfaces and stream lines.


5.3.2 Heterogeneous model

The heterogeneous model with two-layer permeability structure is modified from the homogeneous model described in Sect. 5.3.1. Permeability values of two layers are k1=10-14 m2 and k2=10-15 m2, respectively (see Fig. 9a). The thickness of first layer is 1.1 km, and the other parameters are the same as the homogeneous model. All input files can be found in the 3d/Heterogeneous directory in benchmarks/HydrothermalFoam and benchmarks/USGS_HYDROTHERMAL directories, respectively. The simulation results at 50 kyr are shown in Fig. 9. Vertical slices at x=0 km and x=0.5 km and horizontal slices at z=0.5 km and z=2.5 km are shown in Fig. 9a–d. The higher permeability in the upper layer results in mixing with colder ambient fluids and focusing of the upflow zone. Three-dimensional isothermal surfaces of 300, 200, and 100 C are shown in Fig. 9e.

Figure 9Results of three-dimensional model with two-layer permeability structure. Two layers are separated by the dashed green line. Contour lines for fluid pressure (blue) and temperature calculated from HydrothermalFoam are plotted as solid lines, and that from HYDROTHERM are plotted as dotted black lines. Panels (a) and (b) show pressure and temperature contours on vertical slices, panels (c) and (d) show temperature contours on horizontal slices. Panel (e) shows three-dimensional isothermal surfaces and stream lines.


5.4 Cookbooks

In addition to the presented benchmarks, we have added a number of cookbooks to the code repositories that can be used as starting points for more complex models. These include simple 2-D and 3-D box models, 2-D single-pass loop models, and time-dependent permeability models. They also include examples of how to use more complex meshes generated by Gmsh (Geuzaine and Remacle2009). We intend to add additional cookbooks in the future and hope to receive contributions from users of HydrothermalFoam.

6 Conclusions

We have presented a toolbox for simulating flow in submarine hydrothermal circulation systems. Being based on the widely used fluid-dynamic simulation platform OpenFOAM, the toolbox provides the user with robust parallelized 3-D solvers and a whole suite of pre- and post-processing tools. The toolbox is meant to provide the interdisciplinary submarine hydrothermal systems community with an accessible and easy-to-use open-source platform for testing ideas on how hydrothermal systems “work”. The benchmark tests have shown that model matches previously published models and that the cookbooks provide the user with starting points for building more sophisticated models. By following an open-source approach and by providing extensive code documentation, we hope that the presented model will facilitate integrative studies that combine models with data to better assess the role of submarine hydrothermalism in the Earth system.

Appendix A: Key code snippets

Listing A1Example of fvSchemes file.


Listing A2Example of fvSolution file.


Code availability

The supplement related to this article is available online at:

Author contributions

ZG, LR, and CT designed the project. ZG developed the source code, ran simulations, and wrote the paper. LR provided suggestions for the benchmarks and co-wrote the paper and manual. CT provided further suggestions for the manuscript and the model. All authors discussed and contributed to the final paper.

Competing interests

The authors declare that they have no conflict of interest.


We would like to thank the reviewers Cyprien Soulaine, Matteo Cerminara, and Rene Gassmoeller for their helpful and constructive comments. Special thanks also go to the topical editor Thomas Poulet for his professional handling of the manuscript.

Financial support

This research has been supported by the National Key R&D Program of China (grant nos. 2018YFC0309901, 2017YFC0306603, 2017YFC0306803, and 2017YFC0306203), the COMRA Major Project (grant nos. DY135-S1-01-01 and DY135-S1-01-06), and the German Science Foundation (DFG) (grant no. 428603082).

The article processing charges for this open-access
publication were covered by a Research
Centre of the Helmholtz Association.

Review statement

This paper was edited by Thomas Poulet and reviewed by Cyprien Soulaine, Matteo Cerminara, and Rene Gassmoeller.


Andersen, C., Rüpke, L., Hasenclever, J., Grevemeyer, I., and Petersen, S.: Fault geometry and permeability contrast control vent temperatures at the Logatchev 1 hydrothermal field, Mid-Atlantic Ridge, Geology, 43, 51–54,, 2015. a

Barreyre, T., Olive, J. A., Crone, T. J., and Sohn, R. A.: Depth-Dependent Permeability and Heat Output at Basalt-Hosted Hydrothermal Systems Across Mid-Ocean Ridge Spreading Rates, Geochem. Geophy. Geosy., 19, 1259–1281,, 2018. a, b

Carpio, J. and Braack, M.: The effect of numerical methods on the simulation of mid-ocean ridge hydrothermal models, Theor. Comp. Fluid Dyn., 26, 225–243,, 2012. a

Cathles, L. M.: A Capless 350°C Flow Zone Model to Explain Megaplumes, Salinity Variations, and High-Temperature Veins in Ridge Axis Hydrothermal Systems, Econ. Geol., 88, 1977–1988, 1993. a

Coumou, D., Driesner, T., and Heinrich, C. A.: The structure and dynamics of mid-ocean ridge hydrothermal systems, Science, 321, 1825–1828,, 2008. a, b, c

Coumou, D., Driesner, T., Weis, P., and Heinrich, C. A.: Phase separation, brine formation, and salinity variation at Black Smoker hydrothermal systems, J. Geophys. Res.-Sol. Ea., 114, B03212,, 2009. a, b, c

Crone, T. J. and Wilcock, W. S. D.: Modeling the effects of tidal loading on mid-ocean ridge hydrothermal systems, Geochem. Geophy. Geosy., 6, 1–25,, 2005. a

Crone, T. J., Tolstoy, M., and Stroup, D. F.: Permeability structure of young ocean crust from poroelastically triggered earthquakes, Geophys. Res. Lett., 38, 1–5,, 2011. a

deMartin, B. J., Canales, R. A. R., Canales, J. P., and Humphris, S. E.: Kinematics and geometry of active detachment faulting beneath the Trans-Atlantic Geotraverse (TAG) hydrothermal field on the Mid-Atlantic Ridge, Geology, 35, 711–714,, 2007. a

Driesner, T.: The system H2O–NaCl. Part II: Correlations for molar volume, enthalpy, and isobaric heat capacity from 0 to 1000 C, 1 to 5000 bar, and 0 to 1 XNaCl, Geochim. Cosmochim. Ac., 71, 4902–4919,, 2007. a

Driesner, T.: The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems, Geofluids, 10, 132–141,, 2010. a

Driesner, T. and Heinrich, C. A.: The system H2O–NaCl. Part I: Correlation formulae for phase relations in temperature–pressure–composition space from 0 to 1000 C, 0 to 5000 bar, and 0 to 1 XNaCl, Geochim. Cosmochim. Ac., 71, 4880–4901,, 2007. a

Dunn, R. A., Toomey, D. R., and Solomon, S. C.: Three-dimensional seismic structure and physical properties of the crust and shallow mantle beneath the East Pacific Rise at 9 degrees 30'N, J. Geophys. Res., 105, 23537–23555, 2000. a

Elderfield, H. and Schultz, A.: Mid-ocean ridge hydrothermal fluxes and the chemical composition of the ocean, Annu. Rev. Earth Pl. Sc., 24, 191–224,, 1996. a

Faak, K., Coogan, L. A., and Chakraborty, S.: Near conductive cooling rates in the upper-plutonic section of crust formed at the East Pacific Rise, Earth Planet. Sc. Lett., 423, 36–47,, 2015. a

Flemisch, B., Darcis, M., Erbertseder, K., Faigle, B., Lauser, A., Mosthaf, K., Müthing, S., Nuske, P., Tatomir, A., Wolff, M., and Helmig, R.: DuMux: DUNE for multi-phase, component, scale, physics,... flow and transport in porous media, Adv. Water Resour., 34, 1102–1112,, 2011. a

Fontaine, F. J., Wilcock, W. S. D., Foustoukos, D. E., and Butterfield, D. A.: A Si-Cl geothermobarometer for the reaction zone of high-temperature, basaltic-hosted mid-ocean ridge hydrothermal systems, Geochem. Geophy. Geosy., 10, 5,, 2009. a

Fontaine, F. J., Cannat, M., Escartin, J., and Crawford, W. C.: Along-axis hydrothermal flow at the axis of slow spreading Mid-Ocean Ridges: Insights from numerical models of the Lucky Strike vent field (MAR), Geochem. Geophy. Geosy., 15, 2918–2931,, 2014. a, b

Garg, S. and Pritchett, J.: On pressure-work, viscous dissipation and the energy balance relation for geothermal reservoirs, Adv. Water Resour., 1, 41–47,, 1977. a

German, C. and Seyfried, W.: Hydrothermal Processes, in: Treatise on Geochemistry, Elsevier, 191–233,, 2014. a

German, C. R., Casciotti, K. A., Dutay, J.-C., Heimbürger, L. E., Jenkins, W. J., Measures, C. I., Mills, R. A., Obata, H., Schlitzer, R., Tagliabue, A., Turner, D. R., and Whitby, H.: Hydrothermal impacts on trace element and isotope ocean biogeochemistry, Philos. T. Roy. Soc. A, 374, 20160035,, 2016. a

Germanovich, L. N., Lowell, R. P., and Astakhov, D. K.: Stress-dependent permeability and the formation of seafloor event plumes, J. Geophys. Res.-Sol. Ea., 105, 8341–8354,, 2000. a

Geuzaine, C. and Remacle, J.-F.: Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities, Int. J. Nume. Meth. Eng., 79, 1309–1331,, 2009. a, b

Grevemeyer, I., Hayman, N. W., Lange, D., Peirce, C., Papenberg, C., Van Avendonk, H. J. A., Schmid, F., de La Peña, L. G., and Dannowski, A.: Constraining the maximum depth of brittle deformation at slow- and ultraslow-spreading ridges using microseismicity, Geology, 47, 1069–1073,, 2019. a

Guo, Z. and Rüpke, L.: HydrothermalFoam v1.0, Zenodo,, 2020. a, b, c, d, e

Hannington, M., Jamieson, J., Monecke, T., Petersen, S., and Beaulieu, S.: The abundance of seafloor massive sulfide deposits, Geology, 39, 1155–1158,, 2011. a

Hasenclever, J., Theissen-Krah, S., Rüpke, L. H., Morgan, J. P., Iyer, K., Petersen, S., and Devey, C. W.: Hybrid shallow on-axis and deep off-axis hydrothermal circulation at fast-spreading ridges, Nature, 508, 508–512,, 2014. a, b, c, d, e, f, g

Horgue, P., Soulaine, C., Franc, J., Guibert, R., and Debenest, G.: An open-source toolbox for multiphase flow in porous media, Comput. Phys. Commun., 187, 217–226,, 2015. a

Ingebritsen, S., Geiger, S., Hurwitz, S., and Driesner, T.: Numerical simulation of magmatic hydrothermal systems, Rev. Geophys., 48, RG1002,, 2010. a, b, c

Jamieson, J. W., Clague, D. A., and Hannington, M. D.: Hydrothermal sulfide accumulation along the Endeavour Segment, Juan de Fuca Ridge, Earth Planet. Sc. Lett., 395, 136–148,, 2014. a

Jasak, H.: Error analysis and estimation for the finite volume method with applications to fluid flows., PhD thesis, Imperial College of Science, Technology and Medicine, available at: (last access: 21 December 2020), 1996. a, b

Jupp, T. and Schultz, A.: A thermodynamic explanation for black smoker temperatures, Nature, 403, 880–883,, 2000. a, b, c

Jupp, T. E. and Schultz, A.: Physical balances in subseafloor hydrothermal convection cells, J. Geophys. Res.-Sol. Ea., 109, B05101,, 2004. a

Kipp, K. L., Hsieh, P. A., and Charlton, S. R.: Guide to the Revised Ground-Water Flow and Heat Transport Simulator: HYDROTHERM–Version 3, US Department of the Interior, US Geological Survey, available at: (last access: 21 December 2020), 2008. a

Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J.-O., Fischer, T., Görke, U. J., Kalbacher, T., Kosakowski, G., McDermott, C., et al.: OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media, Environ. Earth Sci., 67, 589–599,, 2012. a

Lewis, K. and Lowell, R.: Numerical modeling of two-phase flow in the NaCl-H2O system: Introduction of a numerical method and benchmarking, J. Geophys. Res.-Sol. Ea., 114, B05202,, 2009a. a, b

Lewis, K. and Lowell, R.: Numerical modeling of two-phase flow in the NaCl-H2O system: 2. Examples, J. Geophys. Res.-Sol. Ea., 114, B08204,, 2009b. a, b

Lie, K.-A.: An introduction to reservoir simulation using MATLAB/GNU Octave: user guide for the MATLAB Reservoir Simulation Toolbox (MRST), Cambridge University Press,, 2019. a

Lister, C. R. B.: Penetration of water into hot rock, Geophys. J. Roy. Astr. S., 39, 465–509,, 1974. a, b

Lowell, R. P.: Modeling continental and submarine hydrothermal systems, Rev. Geophys., 29, 457–476,, 1991. a

Lowell, R. P., Farough, A., Germanovich, L. N., Hebert, L. B., and Horne, R.: A Vent-Field-Scale Model of the East Pacific Rise 9 degrees 50' N Magma-Hydrothermal System, Oceanography, 25, 158–167, 2012. a

Middleton, J. L., Langmuir, C. H., Mukhopadhyay, S., McManus, J. F., and Mitrovica, J. X.: Hydrothermal iron flux variability following rapid sea level changes, Geophys. Res. Lett., 43, 3848–3856,, 2016. a

Moukalled, F., Mangani, L., and Darwish, M.: The finite volume method in computational fluid dynamics, Springer, vol. 113,, 2016. a

Olive, J. A. and Crone, T. J.: Smoke Without Fire: How Long Can Thermal Cracking Sustain Hydrothermal Circulation in the Absence of Magmatic Heat?, J. Geophys. Res.-Sol. Ea., 123, 4561–4581,, 2018. a

Orgogozo, L.: RichardsFoam2: A new version of RichardsFoam devoted to the modelling of the vadose zone, Comput. Phys. Commun., 196, 619–620,, 2015. a

Pruess, K., Oldenburg, C., and Moridis, G.: TOUGH2 user's guide, version 2.1. Earth sciences division, Lawrence Berkeley National Laboratory University of California, Berkeley, California (USA), available at: (last access: 21 December 2020), 1999. a

Pye, J.: Freesteam project, available at: (last access: 21 December 2020), 2010. a

Schlindwein, V. and Schmid, F.: Mid-ocean-ridge seismicity reveals extreme types of ocean lithosphere, Nature, 535, 276–279,, 2016. a

Schmeling, H., Marquart, G., and Grebe, M.: A porous flow approach to model thermal non-equilibrium applicable to melt migration, Geophys. J. Int., 212, 119–138,, 2018. a

Singh, S. and Lowell, R. P.: Thermal response of mid-ocean ridge hydrothermal systems to perturbations, Deep-Sea Res. Pt. II, 121, 41–52,, 2015.  a

Stein, A. and Stein, S.: Constraints on hydrothermal heat flux through the oceanic lithosphere from global heat flow, J. Geophys. Res., 99, 3081–3095, 1994. a

Tagliabue, A., Bopp, L., Dutay, J. C., Bowie, A. R., Chever, F., Jean-Baptiste, P., Bucciarelli, E., Lannuzel, D., Remenyi, T., Sarthou, G., Aumont, O., Gehlen, M., and Jeandel, C.: Hydrothermal contribution to the oceanic dissolved iron inventory, Nat. Geosci., 3, 252–256,, 2010. a

Tao, C., Seyfried, W. E., J., Lowell, R. P., Liu, Y., Liang, J., Guo, Z., Ding, K., Zhang, H., Liu, J., Qiu, L., Egorov, I., Liao, S., Zhao, M., Zhou, J., Deng, X., Li, H., Wang, H., Cai, W., Zhang, G., Zhou, H., Lin, J., and Li, W.: Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge, Nat. Commun., 11, 1300,, 2020. a, b

Theissen-Krah, S., Iyer, K., Rüpke, L. H., and Morgan, J. P.: Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to fast spreading ridges, Earth Planet. Sc. Lett., 311, 275–286, 2011. a

Theissen-Krah, S., Rüpke, L. H., and Hasenclever, J.: Modes of crustal accretion and their implications for hydrothermal circulation, Geophys. Res. Lett., 43, 1124–1131,, 2016. a

Vehling, F., Hasenclever, J., and Rüpke, L.: Implementation Strategies for Accurate and Efficient Control Volume-Based Two-Phase Hydrothermal Flow Solutions, Transport Porous Med., 121, 233–261,, 2018. a

Weis, P., Driesner, T., Coumou, D., and Geiger, S.: Hydrothermal, multiphase convection of H2O-NaCl fluids from ambient to magmatic temperatures: a new numerical scheme and benchmarks for code comparison, Geofluids, 14, 347–371,, 2014. a, b, c, d

Weller, H. G., Tabor, G., Jasak, H., and Fureby, C.: A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys., 12, 620–631,, 1998. a

Wilcock, W. S. D.: Physical response of mid-ocean ridge hydrothermal systems to local earthquakes, Geochem. Geophy. Geosy., 5, 11,, 2004. a, b

Zyvoloski, G. A., Robinson, B. A., Dash, Z. V., and Trease, L. L.: Summary of the models and methods for the FEHM application-a finite-element heat-and mass-transfer code, Tech. rep., Los Alamos National Lab., NM (US), available at: (last access: 21 December 2020), 1997. a

Short summary
We present the 3-D hydro-thermo-transport model HydrothermalFoam v1.0, which we designed to provide the marine geosciences community with an easy-to-use and state-of-the-art tool for simulating mass and energy transport in submarine hydrothermal systems. HydrothermalFoam is based on the popular open-source platform OpenFOAM, comes with a number of tutorials, and is published under the GNU General Public License v3.0.