Crustal Conditions Favoring Convective Downward Migration of Fractures in Deep Hydrothermal Systems

Cooling magma plutons and intrusions are the heat sources for hydrothermal systems in volcanic settings. To explain system longevity and observed heat transfer at rates higher than those explained by pure conduction, the concept of fluid convection in fractures that deepen because of thermal rock contraction has been proposed as a heat‐source mechanism. While recent numerical studies have supported this half a century old hypothesis, understanding of the various regimes where convective downward migration of fractures can be an effective mechanism for heat transfer is lacking. Using a numerical model for fluid flow and fracture propagation in thermo‐poroelastic media, we investigate scenarios for which convective downward migration of fractures may occur. Our results support convective downward migration of fractures as a possible mechanism for development of hydrothermal systems, both for settings within active zones of volcanism and spreading and, under favorable conditions, in older crust away from such zones.

• Numerical modeling supports convective downward migration of fractures as a source mechanism for hydrothermal systems • Fluid flow, fracture opening and propagation in a thermo-poroelastic rock mass are simulated in different geological settings in the crust • Crustal stresses are key to understanding whether a hydrothermal system can evolve in regions away from active zones of volcanism

Supporting Information:
Supporting Information may be found in the online version of this article. 10.1029/2023GL105380 2 of 11 Weir, 2001).For hydrothermal systems in Iceland, Bodvarsson (1982) argued that the intrusive intensity would be low, and, hence, hypothesized CDM to be an important mechanism for heat transfer.Bodvarsson (1982Bodvarsson ( , 1983) ) and Axelsson (1985) suggested that CDM could be effective as a source mechanism also in systems with elevated heat flux away from central magmatic heat sources and, hence, may account for some of the long-lasting low-enthalpy hydrothermal activity in the Icelandic crust.For such settings, they both proposed that the onset of the process would depend on local stress conditions.Such settings are not only limited to Iceland but can also be found in other areas (Jolie et al., 2021;Limberger et al., 2018).A better understanding of the settings controlling the highly coupled processes governing CDM may therefore also shed light on the existence of hydrothermal systems in other parts of the world.
While a comprehensive understanding of the geological settings where CDM can occur is still lacking, two recent numerical modeling studies support the hypothesis.Patterson and Driesner (2021) present a model of large-scale natural convection in a downward propagating fracture zone (of dimensions: H = 3 km, L = 8 km, W = 1 m) in a thermo-elastic medium, considering what we in the present context will denote as a low geothermal gradient of 55°C/km.They investigate the effect of thermoelastic rock stresses and fracture fluid pressure on fracture-zone transmissivity by use of a Barton-Bandis relationship between fracture transmissivity and effective normal stress (Bandis et al., 1983;Barton et al., 1985).Expanding on the conceptual model by Lister (1974) and Axelsson (1985), Stefansson, Keilegavlen, et al. (2021) present a fully coupled numerical approach for fracture propagation and deformation that allow for multiple fractures in a thermo-poroelastic medium.The effect of thermo-poroelastic stresses on fracture transmissivity is incorporated through a fracture contact mechanics model.In the work of Stefansson, Keilegavlen, et al. (2021), the downward migration of fractures due to convective cooling is considered for a test case with a set of smaller fractures (H = 200 m, L = 200 m, W = 2 mm) at the bottom of a geothermal system.The study applied parameters which are representative for a geothermal system in a geological setting with a high geothermal gradient of 150°C/km.
In this paper we use numerical modeling to investigate CDM as a source mechanism for hydrothermal activity.We consider both young crust in active zones of volcanism and spreading as well as older crust away from such zones.The numerical approach builds on the methodology by Stefansson, Keilegavlen, et al. (2021), accounting for flow and fracture propagation in thermo-poroelastic media.This enables us to study effects of the stress regime and thermal gradients, as well as effects of important rock parameters such as the permeability in the medium surrounding the fractures.The results by Stefansson, Keilegavlen, et al. (2021) indicate that CDM can be a plausible source mechanism in systems away from central heat sources if the thermal gradient is sufficiently high.Furthermore, it is clear that the local stress setting (e.g., in spreading systems) can be favorable for CDM.This leads to the following hypotheses: 1.No central magmatic heat sources are needed for CDM; a high geothermal gradient is sufficient to maintain the process.2. With lower geothermal gradients, crustal stress conditions in a range of geological settings are still favorable for CDM as a mechanism for heat transfer.
Using simulation models, we test these hypotheses with different thermal gradients, and show how, for lower thermal gradients found away from volcanic belts, the local stress setting is a dominating factor for the onset of CDM.
As a basis for the simulations, geological conditions found in Iceland are chosen.Iceland is famous for its hydrothermal systems providing its nation with both electricity and space heating through the utilization of geothermal fluids.A zone of spreading and volcanism cuts through the center of the country along the Mid-Atlantic Ridge, explaining the elevated heat content in its crust.The volcanic belts mainly consist of very long fracture-zones dominated by spreading and injection of extremely long dikes at depth.Several central volcano complexes are interspersed in the zone, but these take up a very small part of its total area.For the study of CDM two regional settings are chosen, (a) within the spreading zone but away from any central volcanoes and (b) away from the spreading zone.Those two settings provide scenarios with elevated heat flux, that we in the present context will denote by a high geothermal gradient on one hand and low geothermal gradient on the other hand.As discussed in Section 2, similar settings can be found in other regions of the world.

Geological Settings and Conceptual Model
Areas of elevated heat flow and heat content that favor the development of hydrothermal systems are located at divergent plate boundaries, where the spreading of the lithosphere leads to ascending magma and intrusion into 10.1029/2023GL105380 3 of 11 the crust.Regional extension tectonics create both regional and local structures that furthermore affect permeability.Some examples of divergent plate boundaries include the Mid-Atlantic Ridge, Red Sea Rift, Baikal Rift Zone, East African Rift (Great Rift Valley), East Pacific Rise, Gakkel Ridge, Galapagos Rise, Explorer Ridge, Juan de Fuca Ridge, Pacific-Antarctic Ridge, and West Antarctic Rift System.Many hydrothermal systems are located in those settings, for example, in Iceland (mid-Atlantic Ridge), Eritrea, Djibouti, Ethiopia, Uganda, Kenya, Tanzania and Malawi (western and eastern branches of the East Africa Rift) (Hochstein, 2005).Hydrothermal vents along the seabed are known to form along divergent mid-ocean ridges, such as the East Pacific Rise and the Mid-Atlantic Ridge (Petersen et al., 2018).
Divergent plate boundaries are examples of extension-controlled tectonics.Extension deformation that can enhance vertical permeability can also be found in locations away from convergent or transform plate boundaries, such as in back-arc basins (Lau Basin in the East Pacific), in continental extension zones (Rio Grande Rift Zone, East-African Rift Zone, Western Turkey (Bozkurt & Mittwede, 2005), and in releasing bends along strike-slip faults and in zones of thickened crust (gravitational spreading).Known areas of hydrothermal activity associated with continental rifting are, for example, located in the Great Basin in Western USA, the countries in the western and eastern branches of the East African Rift, in west Turkey and in Cyprus (Bettison-Varga et al., 1992;Murat Özler, 2000).Areas of hydrothermal activity associated with back-arc activity are for example, located in the Taupo Volcanic Zone New Zealand (Kissling & Weir, 2005), and the Okinawa Trough Japan (Halbach et al., 1989;Yang et al., 2020).
Extension deformation due to regional tectonics alone can act to enhance vertical permeability and therefore give rise to the circulation of fluids at intermediate depths in the crust.This, however, does not explain increased hydrothermal activity in tectonically less active areas.Pálmason (1981) speculated on the effect of cooling of the lithosphere as it moves away from the rift axis and suggests that thermal stresses consequently induce enhanced vertical permeability for geothermal fluids at depths on the flanks of the active rifting.Hence, we consider the effects of fluid circulation in existing fractures and how increased thermal stress may cause fracture propagation, thus increasing vertical permeability.It is expected that both temperature settings in the crust and regional stresses influence the initiation and maintenance of CDM.
We study a model of vertical fractures in the roots of geothermal systems, first proposed by Lister (1974) and Bodvarsson (1982).Convection of geothermal fluid in the fractures cools down the surrounding rock, causing horizontal tensile stresses in the rock, which causes (a) the rock to contract and (b) the fractures to propagate downward (Figure S1a in the Supporting Information S1).By this self-sustained process, the convection of thermal fluid extends downward constantly reaching fresh hot rock, enhancing the heat flux to the geothermal system above (Bjornsson & Stefansson, 1987;Björnsson et al., 1982;White, 1968).The process, combined with heat conduction, could sustain the hydrothermal systems in accordance with observed heat output.

Mathematical and Numerical Modeling of CDM
We use a mathematical model based on a discrete fracture matrix representation.The medium is 3D, consisting of the rock matrix and embedded fractures modeled as 2D planes.The model describes energy transport and fluid flow in rock matrix and fractures, thermo-poroelastic deformation of the rock, and fracture deformation and propagation.We assume single phase fluid conditions with the reservoir rock fully saturated, and with local thermal equilibrium between the fluid and the solid.We impose a balance of momentum, mass and energy in the rock matrix and a balance of mass and energy in the fractures, along with kinematic constraints and consti tutive laws for fracture deformation (Stefansson, Berre, & Keilegavlen, 2021; see also Barton, 1976).In the following, key components of the model are reviewed.A full description of the mathematical model, including coupling between variables in the rock matrix and the fracture (Jaffré et al., 2011;Martin et al., 2005;Stefansson, Keilegavlen, et al., 2021), is provided in the Supporting Information S1 of this paper.
Dimension reduction of the balance equations which is necessary to derive the flow of mass and energy in the fracture is detailed in Stefansson, Berre, and Keilegavlen (2021) and Keilegavlen, Berge, et al. (2021).
The aperture of a dimensionally reduced fracture is given by Here, a 0 is the residual hydraulic aperture in the undeformed state, and ⟦u⟧ n the normal component of a displacement-jump over the fracture, defined as the difference in the displacement, u, computed on the opposing fracture walls.The fracture aperture, a, will be affected by fluid pressure as well as thermo-poromechanical 10.1029/2023GL105380 4 of 11 forces in the matrix.The flow in the fracture is described by Darcy's law with a cubic law for the permeability, giving a strongly non-linear relation between the aperture and the fluid flow.
We consider propagation due to tensile forces, modeled by the stress intensity factor (Stefansson, Keilegavlen, et al., 2021; se also Nejati et al., 2015), where R d is the distance between the point where the displacement jump is evaluated and the fracture tip, μ is the shear modulus of the rock, and κ is the Kolosov constant for plain strain described as a function of the shear and the bulk moduli of the rock (see Supporting Information S1).The fracture tip propagates when SI I exceeds a critical value, which can be viewed as the rock toughness.
The mathematical model is implemented in the open-source simulation tool PorePy, which is tailored for representing complex multiphysics processes in fractured porous media (Keilegavlen, Berge, et al., 2021).The fractures are explicitly represented in the computational grid which allows for direct modeling of processes in the fracture and on the fracture walls.In the computational grid, pressure and temperature are represented in both rock matrix and fractures, the displacement is confined to the rock matrix and on the fracture walls, and contact tractions are represented on the fractures.Fracture propagation is represented by extending the fracture grid, with minimal adjustments needed to the rest of the computational model.

Simulations of CDM in Different Geological Settings
Two regional settings are considered to investigate the process of CDM as a source mechanism for hydrothermal activity.They are both considered representative of temperature conditions in the crust of Iceland: (a) within the Icelandic active zone of volcanism and spreading and (b) in older crust away from the active zone of spreading.
The first represents regions rich with geothermal resources while the second represents regions which are considered to include fewer geothermal systems (Axelsson, et al., 2005).The heat source is thermal energy within the Earth's crust, accumulated over time by heat flow from the lithosphere and intrusive activity.The background thermal gradient is 100-150°C/km (Flóvenz & Saemundsson, 1993) in young crust within the rifting zone, and varies between 50° and 100°C/km in older crust further away from the divergent ridge axis.These settings are distinct from known volcano complexes with associated high-enthalpy hydrothermal activity and heat sources of magmatic origin.
The depth at which fissures are assumed to be open varies between those regions: Close to the divergent axis, within the active zone, we assume fissures to be open down to 2 km depth, whereas away from the axis, in colder crust, we assume fissures to be open down to 1 km depth.Isothermal temperature with depth profiles from geothermal fields in Iceland strengthen this assumption (Arnórsson et al., 2008;Axelsson et al., 2000;Björnsson et al., 2000;Xianghui, 2012).We assume that the depth of open fissures reflects the depth of the geothermal systems in the initial state of the simulations.Since we are investigating the longevity of the systems and the CDM as a key mechanism to maintain the systems, we consider the opening and propagation of fractures beyond those depths.

Simulation Domain and Setup
For the numerical investigation of different geological settings, we choose a cube-shaped domain with side-lengths of 400 m.Five fractures are initially located in the top half of the domain, parallel to the yz-plane and evenly spaced, each spanning 200 m in length and depth (see Figure 1 top right for the initial fracture geometry).
The formation depth and thermal gradients are set in accordance with the two regional settings, and geothermal gradient is set accordingly to (a) 130°C/km and (b) 80°C/km, representing temperature conditions within and away from the active zone, respectively.The top of the simulation domain is located at (a) a depth of 2,000 m and (b) a depth of 1,000 m, respectively, thus the temperature at the top of the domain is 260°C and 80°C, respectively, for the two different cases.We denote settings corresponding to (a) as high-temperature (HT) regimes and settings corresponding to (b) as low-temperature (LT) regimes, respectively.
The background stress field is aligned with respect to the fractures, with vertical stress (S V ) equal to the weight of the overburden, the maximum horizontal stress (S H ) in direction of the fractures (along the y-axis), and the minimum horizontal stress (S h ) perpendicular to the fractures (x-axis).This background stress implies that fractures are favorably oriented for opening and propagation.10.1029/2023GL105380 6 of 11 see Table 1).The choice of the background stress, which corresponds to either a strike-slip or a normal fault regime, is based on stress settings observed in Iceland (Ziegler et al., 2016), away from the central volcanism and therefore also away from the complex stress conditions associated with volcanic activity.
The parameters representing the regional settings are given in Table 1).Motivated by the assumption of a low rock porosity of 5%, the bulk density of the overburden is set equal to that of the rock.The matrix permeability and the residual aperture are set in the lower range of what is suggested by measurements in active systems (Keilegavlen, Duboeuf, et al., 2021;Massiot et al., 2017;Sigurdsson et al., 2000).The linear thermal expansion of the rock is chosen according to the temperature conditions at the two different settings, with the thermal expansion coefficient for the high temperature case set five times larger than for the second case (Yin et al., 2021).The thermal expansion coefficient of water and water viscosity is chosen according to the overall temperature, while other parameters are constant, across the two different regional settings (1 and 2).
Dirichlet boundary conditions are chosen for the temperature and pressure.Temperature gradient and hydrostatic pressure are fixed on the sides of the numerical domain.At the top of the domain the boundary conditions are imposed on both the rock matrix and the fractures.Neumann boundary conditions are chosen for the displacement by imposing the anisotropic tractions defined in Table 1, assuming zero displacement in the center of the bottom boundary.Considering the onset of natural convection, these boundary conditions represent a conservative choice given the prescribed linear temperature profile on the vertical boundaries.
Uniform Cartesian grids are used for the matrix and the fracture domains.For the matrix, the grid has 36, 40, and 40 cells along the x-, y-, and z -axis respectively.The fracture grids conform to the matrix grid; hence, there are 400 cells in each fracture initially.As the fracture propagates, 2D cells are added between 3D grid cell edges along the predefined, vertical fracture path (Figures S1b and S1c in Supporting Information S1).

Results
In total, eight simulations were run, representing eight different geological settings, based on the two regional settings (1-2) and the four background stress regimes (A-D).For the HT regimes, onset of propagation occurs for all the four background stress conditions shown in Table 1.For the LT regimes, onset of propagation only occurs for cases B and C.
The modeled fracture evolution can be seen in Figure 1, for setting 1A (HT, strike-slip) and for setting 2C (LT, normal stress) in the left and right columns respectively.The figure also shows the modeled fracture temperature and the temperature change compared to the initial state in the matrix after 10 years (top row) and 60 years (center row) simulation time.Note that in the HT case the fractures have already begun to propagate after 10 years but in the LT case onset of propagation is in the following time-step.The fracture propagation and associated convection causes significant cooling in the rock matrix, more than 10°C for the LT case and more than 20°C for the HT case after 60 years.The bottom row of Figure 1 shows the temperature in the middle fracture at the end of the simulation period; the fracture has propagated 170 and 160 m downward, respectively, in the HT and LT case.
In Figure 2 the Darcy velocity in the middle fracture is shown for the LT case after 50 and 60 years.As the figure shows, fluid circulation in the fracture forms convection cells, however the flow pattern and number of cells change as the fracture propagates.This is observed in both HT and LT cases and is due to changes in the fracture geometry as well as coupling to the surrounding rock matrix.
The onset of propagation and the propagation speed of the five fractures is shown in Figure 3 for the six simulations where propagation occurs.The onset of propagation in HT settings (1) is after 7-8 years (1A and 1D), and after 3-4 years (1B and 1C), with the thermal expansion coefficient set according to the predicted thermal conditions at 2 km depth within the crust.However, with thermal expansion coefficient set according to the predicted thermal conditions at 1 km depth within the crust in the LT setting (2), there is no propagation for stress regimes A and D, with a background minimum horizontal stress at 60% of the vertical stress.When the background minimum horizontal stress is 40% of the vertical stress (2B and 2C), the cooling-induced thermal stress in the vicinity of the fracture is sufficient to overcome the background stress, with propagation starting after 12-13 Based on additional simulation studies, the transition to a regime of propagation is estimated to occur when the minimum horizontal stress is between 45% and 40% of the vertical stress.

Conclusions
Our results contribute to ongoing discussion on (a) the convective downward migration of fractures (CDM) in the roots of geothermal systems, and (b) CDM's importance in explaining the origin and sustainability of hydrothermal activity.Many have contributed to the development of the mathematical model and description of the phenomenon; however, the complexity of the coupled processes of fluid flow, heat transfer and fracture and rock deformation, have put limitations on its understanding.Based on numerical simulations, we show that this process is plausible in different geological settings known to host hydrothermal systems.
The results of the numerical study show that the proposed CDM is highly relevant in understanding the nature of hydrothermal systems and the origin of hydrothermal activity.Notably, the study shows that CDM is possible in settings away from central heat sources, such as magma or cooling intrusions.As proposed by Bodvarsson (1982Bodvarsson ( , 1983) ) and Axelsson (1985), a relatively low geothermal gradient is sufficient to initiate the process.This is also supported by a resent numerical study (Patterson & Driesner, 2021).We have further shown that, in the absence of a high thermal gradient (e.g., in old crust away from active zones of volcanism and spreading), the local stress setting is important.The present numerical simulations show that CDM is possible in both strike-slip and normal stress regimes.In the case of lower thermal gradient, the stress perpendicular to the fracture must be low compared to the vertical stress-as already pointed out by previous studies, but now strengthened by results of the present numerical simulation.Therefore, our results indicate that crustal stresses are a controlling facture as to whether a hydrothermal system can evolve in regions away from volcano-tectonic activity.In the search for hidden hydrothermal activity, this could be key.

Figure 1 (
top right) shows the initial geometrical setup of the domain and the fractures.The background stress is defined such that ℎ = xx = 1zz,  = yy = 2zz,  = zz =  , (4) where ρ S is the bulk density of the overburden, g is the acceleration of gravity, z is the depth and b 1 and b 2 are positive constants.Four background stress regimes (A-D) are defined for b 1 = {0.4,0.6} and for b 2 = {0.8,1.2},

Figure 1 .
Figure 1.From left to right: Setting 1A (high temperature, strike-slip stress), and 2C (low temperature, normal stress).From top to bottom: Temperature distribution in fractures and temperature change in matrix after 10 years (top) and 60 years (center) simulation time and temperature in the middle fracture at the end of the simulation time (bottom).
Background-Stress for Regimes A to D and (b) Parameters for Regional Settings 1 and 2, Defining Conditions for the Eight Different Cases Modeled

Figure 2 .
Figure 2. Darcy velocity in the middle fracture in setting 2C (low temperature, normal stress) after 50 years (left) and 60years (right) simulation time.

Figure 3 .
Figure 3. Propagation of the five fractures depicted as growth in surface area of each fracture.Results for the high temperature regimes are presented in the upper two rows and for the low temperature regimes in the bottom row.The plots to the left are for a strike-slip background stress regime, and the plots to the right are for a normal background stress regime.