Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Analyzing the influence of correlation length in permeability on convective systems in heterogeneous aquifers using entropy production

Analyzing the influence of correlation length in permeability on convective systems in... jan.niederau@erdw.ethz.ch Institute for Applied Hydrothermal convection in porous geothermal reservoir systems can be seen as Geophysics and Geothermal a double-edged sword. On the one hand, regions of upflow in convective systems Energy, E.ON Energy Research Center, RWTH can increase the geothermal energy potential of the reservoir; on the other hand, Aachen University, Aachen, convection introduces uncertainty, because it can be difficult to locate these regions Germany of upflow. Several predictive criteria, such as the Rayleigh number, exist to estimate Full list of author information is available at the end of the whether convection might occur under certain conditions. As such, it is of interest article which factors influence locations of upwelling regions and how these factors can be determined. We use the thermodynamic measure entropy production to describe the influence of spatially heterogeneous permeability on a hydrothermal convection pattern in a 2D model of a hot sedimentary aquifer system in the Perth Basin, Western Australia. To this end, we set up a Monte Carlo study with multiple ensembles. Each ensemble contains several hundred realizations of spatially heterogeneous permeabil- ity. The ensembles only differ in the horizontal spatial continuity (i.e., correlation length) of permeability. The entropy production of the simulated ensembles shows that the convection patterns in our models drastically change with the introduction and increase of a finite, lateral correlation length in permeability. An initial decrease of the average entropy production number with increasing lateral correlation length shows that fewer ensemble members show convection. When neglecting the purely conduc- tive ensembles in our analysis, no significant change in the number of convection cells is seen for lateral correlation lengths larger than 2000 m. The result suggests that the strength of convective heat transfer is not sensitive to changes in lateral correlation length beyond a specific factor. It does, however, change strongly compared to simula- tions with a homogeneous permeability field. As such, while the uncertainty in spatial continuity of permeability may not strongly influence the convective heat transfer, our findings show that it is important to consider spatial heterogeneity and continuity of permeability when simulating convective heat transfer in an aquifer. Keywords: Hydrothermal convection, Entropy production, Heterogeneous permeability, Perth Basin, Spatial continuity © The Author(s) 2019. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Niederau et al. Geotherm Energy (2019) 7:35 Page 2 of 27 Introduction Convective heat transfer in sedimentary reservoir systems can significantly increase their geothermal potential, as hot water is transported to shallower depths by zones of upflow (Barbier 2002). However, this increased potential is connected to a similarly increased uncertainty, as these zones can be difficult to locate due to uncertain reser - voir structure and properties (Rana et  al. 1979; Schilling et  al. 2013; Lipsey et  al. 2016; Niederau et al. 2017), or as they are influenced active exploitation (Grant et al. 1982). On the other hand, convection can increase the risk of success for geothermal projects, e.g., when boreholes are drilled in a zone of downflow. As such, knowledge about the state of the convection system and its pattern is crucial for assessing the uncertainty and associ- ated risk of an exploration study. Several parameters, such as aquifer thickness or perme- ability, control whether free convection can occur. However, those “global” parameters do not give unique insight into spatial distribution of up- and downwelling zones. It is, therefore, necessary to incorporate more detailed knowledge about these parameters for predicting the hydrothermal potential of (sedimentary) reservoir systems with convec- tion. While the spatial extent of a reservoir unit can be derived from geophysical meas- urements to a fine degree of detail, detailed information about the spatial distribution of permeability is often not available. It is usually also difficult to assess, as permeability of sedimentary reservoir systems correlates to the depositional paleo-environment. Depo- sitional structures may be inferrable from seismics to some extent, and seismic proper- ties of rocks can be related to their hydraulic properties (Rubin et al. 1992). But even if such methods are successful, an uncertainty in spatial permeability distribution remains. The Perth Basin, a sedimentary rift basin in Western Australia, was extensively studied for its geothermal potential with focus on one prominent reservoir unit, the Yarragadee Aquifer. While the hydrothermal conditions in this aquifer are largely unknown (Pujol et  al. 2015), multiple modeling studies suggest that hydrothermal convection occurs within the Yarragadee Aquifer (Sheldon et al. 2012; Schilling et al. 2013; Niederau et al. 2017). In particular, Sheldon et  al. (2012) show that temperature and salinity measure- ments in the Perth Basin indicate convective flow. For instance, they show that laterally varying temperature gradients which were deduced from corrected bottom hole tem- peratures, are indicative for regions of upflow or downflow, respectively. Furthermore, zones of uniform salinity were measured in the Yarragadee Aquifer. Such homogeniza- tion of salinity can be caused by convective flow. Horowitz et  al. (2008) present three arguments indicating convection in the Perth Basin: (i) anomalously low thermal gradi- ents in deep aquifers, (ii) thermal gradients varying from 28 to 68 K/km within a region of 55 km × 50 km , (iii) Rayleigh numbers above the critical one, even if conservative parameters are used for the Yarragadee Aquifer when calculating the Rayleigh number. Timms et  al. (2015) classified available core material from wells to several lithostrati - graphic units corresponding to different depositional paleoenvironments in the Perth Basin. While they found a general correlation between lithostratigraphy and average permeability, the spatial correlation of permeability in the Yarragadee Aquifer remains highly uncertain. However, for characterizing the convection pattern, knowledge about heterogeneity of permeability and its spatial continuity is crucial, as both influence the convection system in an aquifer (Nield and Simmons 2007; Nield et  al. 2009). For instance, Simmons et al. (2010) conclude that the scale of heterogeneity, i.e., the size of Niederau et al. Geotherm Energy (2019) 7:35 Page 3 of 27 standard deviation of the underlying permeability distribution, has a strong impact on the nucleation of a convection system. They observe further that a convection system is moderately sensitive to changes in the correlation length. With uncertain spatial heterogeneity of permeability, an associated convection pat- tern is also uncertain. For assessing the geothermal potential of the Perth Basin with a focus on the spatial variability of a convection pattern, it is, therefore, necessary to address heterogeneous permeability and its spatial correlation in the Yarragadee Aquifer. Based on stratigraphic forward modeling, Irvine et  al. (2014) created Rayleigh number maps of the Yarragadee Aquifer, addressing vertical heterogeneity of permeability in the Yarragadee Aquifer. They concluded that convection is likely in the Yarragadee Aquifer and that these maps are a suitable tool for locating free convection in the aquifer. Later- ally continuous permeability facilitates a characterization based on the Rayleigh num- ber (Rayleigh 1916). However, permeability in aquifers is often laterally discontinuous. Performing a standard Rayleigh number analysis in such heterogeneous porous media is difficult and not straightforward (Nield and Simmons 2007; Nield et al. 2009). In the presented work, a different, thermodynamically motivated measure is used for charac - terizing the state of a convection system: the so-called entropy production (Kleidon and Lorenz 2005; Siavashi et al. 2015; Börsing et al. 2017). While the porous Rayleigh num- ber is a sound and tested indicator for characterizing a convective system in an aquifer, it needs the two parameters of the porous medium: L (thickness of the porous medium) and k (the average permeability). These are not constant in real-world reservoir systems, and are usually averaged, i.e., simplified. By assessing the additional heat transport intro - duced by convection, the measure of entropy production bypasses these two simplifica - tions. This can be shown in a numerical experiment with a simple geometry. Figure  1 shows the heat and mass transport in a closed-box model similar to Börsing et al. (2017), with changes made only to the average vertical permeability, either in total or by intro- ducing layers of low permeability. In plot a, the average permeability is high enough to cause instabilities, which are indicated by both, Rayleigh number and entropy produc- tion number (Ra = Rayleigh number and N  = entropy production number in the figure, respectively). Both these numbers decrease, when heterogeneities in permeability (plot b) are introduced. While this is straightforward, it becomes interesting, when the (har- monic) average of (vertical) permeability in a model yields subcritical Rayleigh numbers, as in Fig. 1c, d. These illustrate that models with a similar Rayleigh number may behave completely differently due to a spatially heterogeneous permeability. A subcritical Ray - leigh number, as well as an entropy production number around 1, suggests that heat dif- fusion is the dominant transport process (Fig. 1c). A horizontal, barely permeable layer intersects the aquifer in Fig.  1d. The average vertical permeability in Model c and d is similar, and so are their Rayleigh numbers (1.43 and 1.54, respectively). The Rayleigh number does not suggest convection in model d, while the entropy production num- ber diagnoses convection, as it is greater than 1. Also, a Rayleigh number is a predictive measure of a system, while entropy production is a diagnostic one, thus corresponding more closely to the Nusselt number. In comparison to the Nusselt number, entropy pro- duction has the advantage of being a volumetric estimate. The non-dimensional Nusselt number, defined as the ratio of total heat flow to conductive heat flow over a defined domain, is estimated across a cross-sectional area. Furthermore, similar to the Rayleigh Niederau et al. Geotherm Energy (2019) 7:35 Page 4 of 27 Fig. 1 Four scenarios of temperature in a closed-box model similar to models in Börsing et al. (2017). a Model with homogeneous permeability. b Introduction of a layer with low permeability (black-outlined block) makes permeability heterogeneous. c, d Models with similar harmonic average vertical permeability. In c, the model has a homogeneous permeability field; in d, a barely permeable horizontal layer causes a heterogeneous permeability field. The average vertical permeability of the model, as well as the corresponding Rayleigh number and entropy production number is annotated for each scenario number, it is complicated to assess under non-idealized, i.e., more complex, model domains. Entropy production is mostly independent of the domain geometry as it is cal- culated from the spatial temperature distribution induced by heat transport processes. In this work, we estimate the influence of different correlation lengths on the asso - ciated convection pattern by systematically varying the lateral correlation length of porosity and permeability in a 2D cross section of the Perth Basin model by Niederau et al. (2017). For each correlation length, we generate a Monte Carlo ensemble com- prising 240 ensemble members. To retain comparability, all ensembles have the same random number seeds. This is illustrated in Fig.  2, which shows three exemplary per- meability fields with different horizontal correlation lengths. Due to the same random number seed, areas of low or high permeability are roughly at the same location. The presented manuscript is structured as follows: after reviewing the geological back - ground and introducing the concept of entropy production, we describe the geological model used in this study, as well as the corresponding stratigraphical and petrophysical model. The latter two are connected via the lithostratigraphy of the Yarragadee Aqui - fer described by Timms et al. (2015). We then present the structure of the Monte Carlo ensembles and the setup of the numerical model. The Monte Carlo ensembles are ana - lyzed in the results section and discussed afterwards. Niederau et al. Geotherm Energy (2019) 7:35 Page 5 of 27 Fig. 2 Three examples of permeability fields in the reservoir model used in this study. Permeability fields only differ in the prescribed horizontal correlation length (here, 500 m, 2500 m, and 5000 m). All models have the same random number seed, so the general spatial location of high- and low-permeable areas is preserved. Note that the model is strongly vertically exaggerated Methods Forward code For modeling heat and mass transport in a porous medium, we use the coupled balance equations for energy, mass, and momentum. These conservation equations are imple - mented in the code SHEMAT-Suite (Rath et al. 2006), which is developed from the SHE- MAT code presented in Clauser and Bartels (2003). SHEMAT-Suite models heat and mass transport in porous media in steady as well as in transient state. The transient flow equation used in our simulations is based on Darcy’s law and mass conservation: ∂h ρ g k S =∇ · (∇h + ρ ∇z) + W , (1) s r ∂t μ −3 −2 2 where ρ is fluid density ( kg m ), g gravity ( ms ), k permeability tensor ( m ), μ f f dynamic viscosity of the fluid (Pa  s), h hydraulic head (m), and z the vertical coordi - −1 nate space (m). W describes a volume source term ( s ); ρ is the relative density change ρ −ρ f 0 defined by ρ = , where ρ is density at reference conditions (i.e., atmospheric con- r 0 −1 −1 ditions). S = ρ g(α + φα ) ( m ) is the specific storage coefficient ( m ), ( α and α s f r f r f −1 Pa are compressibilities ( ) of the rock and the fluid phase, respectively). Equation  (1) is derived, using an extended Oberbeck–Boussinesq approximation (Kolditz et al. 1998; Diersch and Kolditz 2002). In our simulations, however, fluid properties are calculated as a function of temperature and pressure after the formulation given in Zyvoloski et al. (1996). The approximation made for Eq. (1) is valid for cases without very high-temper - ature gradients or brines with high salinity (Diersch and Kolditz 2002). These two condi - tions are not present in our simulations, and the approximation is, therefore, suitable for the conditions considered here. The heat transport equation is derived from Fourier’s law q = − ∇T and energy conservation: ∂T (ρ c ) = ∇ · [ ∇T ] + (ρ c ) v∇T + K , (2) p e e p f ∂t −2 ◦ where q is specific heat flow ( Wm ), T is temperature ( C ), and (ρ c ) and (ρ c ) are p e p the effective thermal capacities of the fluid-saturated porous medium and the fluid −3 −1 −1 −1 ( Jm K ), resp e ctively .  is the tensor of bulk thermal conductivity ( Wm K ), and e Niederau et al. Geotherm Energy (2019) 7:35 Page 6 of 27 ρ gk −3 f K is the radiogenic heat production rate ( Wm ); v = (∇h + ρ ∇z) is specific dis - 0 r −1 charge (or Darcy velocity) ( ms ). Parameter values for the bulk porous rock are derived from appropriate mixing laws. The arithmetic mean was used for estimating an average specific heat capacity and the geometric mean for thermal conductivity (Clauser 2011). A more detailed description of the transport equations is provided in Clauser and Bartels (2003). Entropy production The description of the concept of entropy production presented here follows the descrip - tion in Börsing et al. (2017). The second fundamental law of thermodynamics states that the entropy of a closed system is always increasing. The irreversible entropy produc - tion dS by thermal dissipation characterizes irreversible conductive heat transport dQ irr through a (porous) medium with different temperatures at its ends (from higher T to lower temperatures T ) (Tolman and Fine 1948). It can be described as: 1 1 dS = − dQ. irr (3) T T 2 1 Bejan (2013) formulates the production of entropy in a fluid-filled porous medium as the ′′ sum of contributions by thermal dissipation (i.e., conductive heat transport, ) and therm ′′ viscous dissipation (i.e., fluid friction S ): visc ′′ dS ′′ ′′ ′′ total ˙ ˙ ˙ (4) = S = S + S , total therm visc dt ′′ −3 −1 ′′ ′′ where S is in W m K and S the transient derivative. The notation indicates that the calculation is done across a grid face, i.e., between two cell-centered grid nodes. By the definition in Eq.  (3), the transient change of dS in a specified medium can be ′′ described by introducing a specific heat flow q (Tolman and Fine 1948). This vector describes the rate of heat transported per unit time and area, and can be expressed in terms of temperature gradient and thermal conductivity, according to Fourier’s law. The ′′ transient thermal dissipation term ( S in Eq. (4)) can, hence, be reformulated: therm ′′ 1 q  T + T b 1 2 ′′ ′′ 2 S = q ·∇ =− ·∇T = (∇T ) with T = . 0 (5) therm 2 2 T 2 T T 0 0 T and T are the lower and upper temperature boundaries, respectively (Börsing et  al. 2 1 2017; Siavashi et al. 2017), similar to Eq. (3). These temperatures need to be defined for assessing an average initial temperature T . In case of a conductive case (Fig.  3  left), it is reasonable to assume T as the arithmetic mean between upper and lower temper- ature boundaries. As all transient simulations presented in this manuscript start from the same conductive initial conditions, this assumption remains valid.  is the average thermal conductivity over the model domain. In a geothermal system dominated by heat conduction, the predominant heat transport direction is in the direction of the tempera- ture gradient, i.e., normally vertical (Fig.  3 left). If the dominant transport mechanism is advective heat transport by convection, lateral temperature gradients superpose the Niederau et al. Geotherm Energy (2019) 7:35 Page 7 of 27 Fig. 3 Cross sections showing a temperature field of a conductive solution (left) and a convective solution (right). Horizontal temperature profiles (dashed lines) show no lateral temperature gradient for the conductive state, but a significant lateral temperature gradient for the convective state underlying vertical temperature gradient. That is, in addition to a vertical heat trans - port, a lateral (conductive) heat transport occurs, from hotter areas (upwelling regions) towards colder areas (downwelling regions). This lateral heat flow can be described by Eq.  (5), and yields an increase in entropy in the colder areas, while entropy in the hot- ter areas decreases. It should be noted that this observation does not violate the second law of thermodynamics, as the areas described are open. Over the whole model domain, entropy cannot be reduced. In addition to thermal dissipation, entropy is also produced by viscous dissipation ′′ ( S in Eq. (4)), e.g., by internal friction of the fluid. Entropy production by viscous dis - visc sipation can be described as: μ μ f f ′′ 2 S = v + Φ, visc (6) T k −1 v ms k where is the specific discharge in ( ), and Φ the viscous dissipation function. is the average of the reduced permeability tensor ( k ), as we are working with a 2D model: k k xx xz k = . (7) k k zx zz −4 −1 −6 −1 Considering viscosity of water and realistic ranges (10 ms to 10 ms ) of the spe- cific discharge in aquifers, viscous dissipation will have a negligible contribution to the total entropy production. This was also shown by Costa (2006). Thus, in homogeneous porous media, viscous dissipation only has a small effect on the total entropy production in a system. However, if preferred flow paths of a convective system are blocked by zones of low permeability caused by spatially heterogeneous permeability, the importance of viscous dissipation effects can increase. In case of a two-dimensional model like in this work, Eq. (4) can be written as: Niederau et al. Geotherm Energy (2019) 7:35 Page 8 of 27 2 2 ∂T ∂T μ b f ′′ 2 2 dS = + + [v + v ]. (8) x z total T ∂x ∂z kT 0 0 By integrating over the model domain, one obtains the total entropy production S : total ′′ ˙ ˙ S = S dV . total (9) total Non-dimensionalising the total entropy production yields the entropy production Num- ber N (Siavashi et al. 2015), which is used in this study as a representative of the total entropy production in the system. Using a non-dimensionalized version of entropy pro- duction facilitates comparability of our results with the literature: S q total N = with Y = ˙ (10) VY  T , where q is the total specific heat flow (Siavashi et al. 2015). Several thermal rock proper - ties, such as thermal conductivity, are used for calculating the entropy production num- ber. Each grid cell has a bulk thermal conductivity value, which is the geometric mean of thermal conductivities of the rock matrix and pore-filling fluid. The temperature- dependent bulk thermal conductivity value  , used in Eq. (10), is the average of bulk thermal conductivities in the model grid. Temperature-dependent changes of matrix thermal conductivity are small. For water, however, thermal conductivity and viscosity vary strongly with temperature. Water thermal conductivity in the model cells is calcu- lated by the formulation in Zyvoloski et al. (1996), while its viscosity for calculating the entropy production number is calculated based on the IAPWS formulation (Wagner and Kretzschmar 2008). Reservoir model used in this study Different types of models of the Yarragadee Aquifer are presented in this study. Mod - els range from a geometrical geological model, over conceptual stratigraphic and petro- physical models to a numerical model of relevant partial differential equations presented in the methods section. The synthesis of information from all these models is a reservoir model of the Yarragadee Aquifer. In the following, we will briefly describe the geologi - cal model with focus on the main aquifer unit, the stratigraphic model underlying the spatial facies distribution of the main aquifer unit in this model, and the petrophysical model relating porosity and permeability in the aquifer. Finally, we will present the dis- cretized numerical model which combines the features of the previous in a rectilinear hexahedral grid for simulating the heat and mass transport in the Yarragadee Aquifer. Geological model The geological model in this manuscript is a 2D section extracted from a 3D geologi - cal model of the central Perth Basin presented in Niederau et  al. (2017). Its location is indicated in Fig. 4. The Perth Basin is the result of rifting between Australia and Greater India in the late Paleozoic to Mesozoic (Harris 1994) during the breakup of Gondwana. Sedimentary successions filling the Perth Basin can be up to 15 km thick. Thickness of Niederau et al. Geotherm Energy (2019) 7:35 Page 9 of 27 Fig. 4 Map showing well locations and model boundaries in the study area in the central Perth Basin (modified from Niederau et al. (2017)). Location of the 2D cross section is indicated by the dashed line the Yarragadee Aquifer in the central Perth Basin range from 500  m to about 3000  m. More detailed information about the geological and tectonic history of the Perth Basin can be found in Playford et  al. (1976), Backhouse (1984), Harris (1994) and Crostella and Backhouse (2000). Lithological units in the geological model are divided based on their hydraulic properties (see Niederau et al. 2017), i.e., divided into aquifer and aqui- tard units, based on geological and petrophysical data (Crostella and Backhouse 2000; Davidson and Yu 2008; Sheldon et al. 2012; Timms et al. 2015). This is often referred to as “zonation method” (Testerman 1962; Alpay 1972). The main aquifer unit of interest is the jurassic Yarragadee Aquifer, which comprises lithologies of three stratigraphical units, namely Cadda, Yarragadee, and Gage formations. Working groups of the Western Niederau et al. Geotherm Energy (2019) 7:35 Page 10 of 27 Australian Geothermal Center of Excellence (WAGCoE) published several detailed stud- ies covering the petrography of the Yarragadee Aquifer (see Timms et  al. (2015) and references therein), its estimated lateral extension and thickness, and estimated temper- atures in the aquifer (see Sheldon et  al. 2012 and references therein). Our 2D geologi- cal model comprises six geological units, (three aquifer units and three aquitard units), dissected by three major faults. One fault, the Darling Fault, displaces the sedimentary units of the Perth Basin against the Yilgarn Craton, which bounds the Perth Basin to the east. The Darling Fault is the east boundary of the 2D geological model. Two fur - ther faults displace the basement and the Yarragadee Aquifer, visible in the cross sec- tion in Fig.  4. The Yarragadee Aquifer is covered by Post-rift sediments, of which the South-Perth-Shale acts as a thermal blanket for the Yarragadee Aquifer. In Niederau et al. (2017), we showed that reservoir structure has a potentially stronger impact on a convection pattern than spatially heterogeneous permeability. There, we conclude that abrupt changes in reservoir thickness particularly impact a convective pattern. Conse- quently, when studying convection in the Perth Basin, fault displacements should be rep- resented in the model geometry. The two fault displacements considered in this study, are at x ≈ 30, 000 km , and, more pronounced, in the west, at x ≈ 14, 000 km in Fig. 4. However, Niederau et  al. (2017) kept the spatial correlation length of permeability constant; thus, their conclusion relies on a Monte Carlo study with a singular spatial cor- relation length. Next to differences in aquifer thickness due to tectonic activity, spatially heterogene - ous permeability caused by differential sedimentation and compaction was found to affect a convective pattern. In the following, we will present a stratigraphic model of the Yarragadee Aquifer, giving insight in the spatial distribution of lithofacies and cor- responding spatial differences in aquifer permeability. These spatial differences are the heterogeneity, whose spatial continuity is bound to the lithofacies. As their spatial distri- bution and continuity is uncertain, so is that of permeability. Stratigraphical model The Yarragadee Formation was deposited under a predominantly fluvial sedimentary environment, as first inferred by Playford et  al. (1976). Timms et  al. (2015) conducted a detailed study analyzing lithostratigraphy and petrography of multiple core samples from wells in the central Perth Basin (among others the Cockburn 1 well). They identi - fied a variable lithofacies distribution within the Yarragadee Aquifer, with predominant lithofacies corresponding to channel fill deposits and fluvial barforms. Furthermore, they state that formation porosity strongly correlates with lithofacies, which is sup- ported by the previous studies (Delle Piane et al. 2013). Delle Piane et al. (2013) divided core measurements into groups of different flow-zone indicators (FZI). This method is used for identifying classes of rocks which presumably belong to the same hydraulic unit based on porespace geometry and flow properties (Prasad 2003). As lithofacies are cor - related to porosity according to Timms et al. (2015), and FZI deduced from porosity data connect similar hydraulic units (Delle Piane et al. 2013), it may be concluded that litho- facies and hydraulic units correlate. As the depositional environment of the Yarragadee Aquifer varies spatially, types and changes of depositional environment correlate with different lithofacies, and porosity Niederau et al. Geotherm Energy (2019) 7:35 Page 11 of 27 and permeability significantly vary spatially. This may be obvious, considering the flu - vial nature of the Yarragadee Formation, and that hydraulic parameters in a fluviatile sedimentation environment can significantly vary in space. Studies using stratigraphic forward modeling to characterize spatial sediment distribution and distribution of hydraulic properties also show strong spatial variations of porosity and grain size (Cor- bel et  al. 2012a). Vertical variability in porosity can also be determined from wells and geophysical well-log analysis. In Niederau et  al. (2017), we deduced a vertical correla- tion length of around 150 m for porosity, based on porosity measurements on samples from the borehole Cockburn  1. In sedimentary depositional environments, horizontal correlation of lithofacies and connected properties is generally larger than their vertical correlation. Timms et al. (2015) defined facies based on metrics from Gamma ray logs, called GR facies. Those GR facies broadly correspond to lithofacies described in the Yar - ragadee Aquifer. They did not observe clear cross-well correlation of the GR facies and concluded that the corresponding lithofacies are not continuous over distances greater than 5 km. The horizontal correlation length in our study varies between 0.5  and 10 km and is, therefore, significantly larger than the conclusion of Timms et  al. (2015). We extend the horizontal correlation beyond 5 km to assess if larger correlation lengths still affect the occurrence of hydrothermal convection. In the case of the Perth Basin, how - ever, we focus our main analysis on horizontal correlation lengths between 0.5 and 5 km, considering the findings of Timms et al. (2015). Petrophysical model In Niederau et  al. (2017), we deduced porosity from the sonic log of Cockburn  1 after estimating clay content from the Gamma ray log. Using the Wyllie method (Wyllie et al. 1956), porosity was estimated from measured reciprocal interval velocity (“slowness”) Δt: Δt = φΔt + (1 − φ)(V Δt + (1 − V )Δt ) ma (11) f sh sh sh with φ being porosity, V being the shale volume deduced from the GR log, and t , t , sh f sh t the slowness of formation water, shale fraction, and remaining rock matrix, respec- ma tively. Equation (11) was calibrated using available porosity measurements on core samples from Cockburn  1. For relating porosity with permeability, an equation based on fractal theory (Pape et al. 1999) was previously calibrated to the Yarragadee Aquifer, using data from CSIRO (2007) (see Niederau 2014 and Niederau et  al. 2017 for more details): 2.8 8.01 −18 k =[166φ + 19904φ + 388(10φ) ]× 10 , (12) where k is permeability in m . We created a suite of porosity samples (Fig. 5 left) based on mean and standard deviation of the porosity log deduced from Eq. (11). A correlated permeability distribution was generated using Eq.  (12) (Fig.  5 right). We chose to gen- erate a porosity distribution based on mean and standard deviation of the porosity log over using the porosity log directly to have idealized distributions. Interpretation of the results would be significantly more difficult if the underlying distributions of spatial het - erogeneity are not normal, or log-normal distributions, respectively, but asymmetrical or skewed. Niederau et al. Geotherm Energy (2019) 7:35 Page 12 of 27 Fig. 5 Histograms of porosity and permeability samples used for populating the parameter space in the simulation ensembles. Mean and standard deviation of porosity was deduced by log analysis of the borehole Cockburn 1. The permeability distribution was calculated from the porosity distribution by applying Eq. (12) Monte Carlo ensembles Considering the stratigraphic and petrophysical models, we spatially vary porosity and permeability within the Yarragadee Aquifer using Sequential Gaussian Simulation (Deutsch and Journel 1998). This algorithm requires a distribution to draw samples from and a correlation length for defining the spatial variance. The previously described ver - tical correlation length cl of 150  m is kept constant, while we increase the horizontal correlation length cl from 500 to 10,000 m by 500 m steps. For each horizontal correla- tion length, we generate an ensemble comprising 240 members, all with different—but equally likely—porosity and permeability fields, based on the distributions presented in Fig.  5. Ensembles are simulated on the JURECA-HPC (Jülich Supercomputing Cen- tre 2016) cluster and each ensemble comprises up to 240 members to fully utilize the JURECA architecture, i.e., 24 cores per node with 10 nodes allocated for each ensem- ble. Therefore, the particular ensemble size is a result of the architecture of the JURECA cluster, where we ran our simulations on. For comparison, we assess a simulation with homogeneous porosity and permeability (i.e., infinite spatial correlation). Values for the homogeneous model are presented in Table 1. Porosity and permeability of the homoge- neous model resemble mean values of the corresponding distributions. To evaluate the influence of vertical correlation length, we simulate a series of ensembles with a vertical correlation length of 100  m. Again, we increase the horizontal correlation length from 500 to 5000 m by 500 m steps. Numerical model All presented studies on the impact of spatially heterogeneous permeability on the entropy production caused by convection are based on the same model geometry: a two- dimensional cross section taken from the model presented by Niederau et al. (2017). In the following, we describe the essential model parameters in more detail. Niederau et al. Geotherm Energy (2019) 7:35 Page 13 of 27 Fig. 6 Initial temperature (color coded) and boundary conditions for the 2D model of the Perth Basin. Black lines indicate interfaces between different model units Table 1 Geological units implemented in  the  2D model and  corresponding petrophysical properties which were used in the simulation Model unit Porosity Vertical Thermal conductivity Radiogenic −1 −1 permeability (matrix; W m K ) heat production −15 2 −3 (1 0 m ) (μ W m ) Leederville Aquifer 0.3 118.0 3.4 0.6 South Perth Shale 0.1 0.01 1.8 0.82 Yarragadee Aquifer Mean: 0.2 Mean: 11.8 4.3 0.52 Range: 0.14–0.33 Range: 0.3–4100 4.3 0.52 Cattamarra coal measures 0.1 0.018 4.1 0.45 −4 Basement/Yilgarn Craton 0.06 4 4 1.18 10 References 1 2 3 4,5 References: (1) CSIRO (2007); (2) Schilling et al. (2013); (3) HDRPL (2008); (4) Botman (2010); (5) Middleton (2013) Boundary conditions and initial values The lower ( T ) and upper (T ) boundary temperatures of the model domain are kept l u constant (Fig.  6). T equals 95 C , which is shown to be a representative temperature at a depth of 3.5 km in the central Perth Basin (Sheldon et al. 2012; Schilling et al. 2013), and 19 C , the average annual surface temperature around Perth. Lateral boundaries are set as zero Neumann boundary conditions, implying no heat or mass transport across the boundaries. The top boundary for the hydraulic potential was set as a function of topography and sea level respectively, varying around 3500 m, as the onshore part of the model does not have high topography. The initial conditions for temperature are results of a purely conductive steady-state solution, using model parameters stated in Table  1 and displayed in Fig.  6. For the hydraulic head, we applied a constant initial value of 3500 m, representing hydrostatic conditions in the entire model. Spatial discretization The calculated entropy production is affected by the grid discretization of a numerical model. Kleidon and Lorenz (2005) state that the value of entropy production decreases with smaller cell sizes until a threshold value is reached. However, a finer numerical grid, consequently, means more unknowns and longer simulation times, i.e., higher computa- tional cost. To find the optimal discretization with respect to a trade-off between com - putational cost and accuracy, we simulate the 2D model with different discretizations, ranging from 100 to 5  m in cell size. Figure  7 shows the resulting entropy production Niederau et al. Geotherm Energy (2019) 7:35 Page 14 of 27 Fig. 7 Entropy production number N plotted versus cell size for different discretizations of the 2D model shown in Fig. 6. N decreases with smaller cell sizes until it approaches a steady value of around 1.28 number N for different cell sizes. With equidistant sizes of ≤ 25 m , changes in N are ˙ ˙ S S minor compared to the exponential increase in unknowns. Therefore, an equidistant cell size ( dx = dy = dz ) of 25 m is used for our simulations, yielding a total of 196,000 cells. This represents a discretised domain of 35 km in length and 3.5 km in depth (1400 cells in x-direction times 140 cells in z-direction). Transient evolution Each transient simulation in this study starts from the same initial conditions, a purely conductive temperature field. From this conductive state, convection develops over a certain time interval. It is not unusual, that simulations with convection do not reach a steady state, i.e., the numerical solution does not change with time. A reason for this is often oscillation of convection cells (Niederau et  al. 2017), or their superposition with a regional groundwater flow (Schilling et  al. 2013). While a steady-state solution may not be found, simulations with convection can start to show oscillations, e.g., within a temperature interval. Methods exist to assess the timing when such oscillation within a “solution interval” starts. One method is to observe transient convection in a phase space and to see if the solution orbits around the so-called attractors (Lorenz 1963; Vin- cent and Yuen 1988). In Niederau et al. (2016), we showed that the concept of attractors is a viable method to determine the onset of transient oscillation in the Perth Basin. This onset, and the following oscillation, is called “dynamic steady state”. To assess the mini- mal simulation time needed for allowing a convective system in a dynamic steady state to develop, we study the entropy production over time for a small ensemble comprising 12 members with cl equal to 500  m. Generally, the smaller cl is, the longer simula- h h tions need to reach a dynamic steady state. This is why we used cl of 500  m as proxy for assessing the minimum simulation time for the Monte Carlo study. If no significant changes in the entropy production number with time can be observed, a dynamic steady state is reached. Using the entropy production as a measure based on thermodynam- ics, the dynamic steady state can be seen as a thermodynamic steady state. Figure  8 shows the transient evolution of the entropy production number of such an ensemble Niederau et al. Geotherm Energy (2019) 7:35 Page 15 of 27 Fig. 8 Change of entropy production number N with time of a test ensemble comprising 12 members with a horizontal correlation length of 500 m. The average entropy production number stabilizes after around 45,000 years, whereas the majority of ensemble members do not show a strong variation after 60,000 years with a horizontal correlation length of 500 m. The entropy production number increases steeply from the simulation start, suggesting that the whole ensemble is an supercriti- cal state, i.e., above the corresponding critical Rayleigh number for this system. Despite that all realizations show indications of free convection, peak entropy production num- bers differ as well as the transient development until 50,000–60,000 years. This suggests that spatially heterogeneous permeability not only influences whether convection cells develop, but also how long a convective system needs to stabilize. After about 60,000 years, the convective system stabilizes into a (thermo-) dynamic steady state. This behavior for the convective system in the Perth Basin is also described in Schilling et al. (2013) and Niederau et al. (2017). Based on these findings, we limit the simulated time of ensemble members to 60,000 years. Results In this section, we present the simulation results of the aforementioned Monte Carlo ensembles with focus on the influence of changing spatial correlation length, analyzed using the average entropy production number (Eq.  10). Per definition, higher entropy production numbers correspond to more vigorous free convection, i.e., stronger lateral heat transport. Following the hypothesis that increasing the finite horizontal correlation length of permeability coincides with a less likely development of a convective system, we assess the average entropy production number of each ensemble, and relate it to the scenario with homogeneous permeability. Average entropy production with correlation length Figure 9 shows the entropy production number N for ensembles with a horizontal cor- relation length up to 10 km. The average ensemble entropy production number ( N ) decreases with increasing horizontal correlation length ( cl ), suggesting that more ensemble members do not show convection at higher cl . The entropy production of ensembles with heterogeneous permeability is compared to a homogeneous case, where Niederau et al. Geotherm Energy (2019) 7:35 Page 16 of 27 Fig. 9 Variation of the ensemble-average entropy production number N with different horizontal correlation length cl (vertical correlation length cl = 150 m ). Left: including non-convecting ensemble members. Right: h v excluding non-convecting ensemble members aquifer permeability does not vary spatially and is equal to the mean of the permeabil- ity distribution used in the heterogeneous ensembles (Fig.  5 right). Although correla- tion length in a homogeneous case can theoretically be seen as infinite (i.e., cl >> model domain), we plot it at cl equal zero in Fig.  9. The reason is that infinite and finite cor - relation length are different. Both finite correlation lengths cl and cl would have to be h v magnitudes larger than the model extent to mimic a homogeneous case, which is practi- cally the same as not stating a finite correlation length at all. In that case, only the per - meability value, and not also its spatial change, would be a controlling factor whether a convection system forms, similar to a model with homogeneous permeability. However, due to the sequential nature of the SGSim algorithm, the first determined permeability value would be used for the whole domain in that case, and likely not represent the aver- age permeability of the reservoir. In a case with observable heterogeneous permeability, both, the permeability value, and its change in space, are important to consider, making a homogeneous approach not feasible. This means that in Fig.  9 and following visualiza- tions, a finite correlation length cl of zero means that no finite correlation length exists, i.e., a homogeneous case. The homogeneous case, i.e., no finite correlation length, is characterized by a higher entropy production number, associated with an ideal convection system, where cells are not distorted or diverted by lenses with low permeability. The entropy production of a homogeneous case is used as a reference point for emphasizing the importance of considering heterogeneous permeability in an aquifer (here specifically the Yarragadee Aquifer) when modeling hydrothermal convection. While the entropy production of the whole ensemble continuously decreases with increasing cl (Fig.  9 left), it stabilizes at around 2000 m horizontal correlation length if members without convection are ignored when calculating N (Fig.  9 right). No significant influence of cl on the aver- ˙ h age entropy production of an ensemble is observed for correlation lengths larger 4000 m (Fig.  9 right). This suggests that heat transport is not strongly influenced by changes in cl beyond this range. This can further mean that the number of convection cells is rela - tively stable for cl larger 4000 m. We will test this hypothesis later in this manuscript by assessing the relation between number of convection cells and entropy production. Niederau et al. Geotherm Energy (2019) 7:35 Page 17 of 27 Fig. 10 Boxplots of the entropy production number for each ensemble with different horizontal correlation length cl . At the beginning, average entropy production decreases with increasing cl , but stabilizes at larger h h cl Figure 10 shows boxplots for each ensemble, including only members showing convec- tion. While the mean entropy production decreases until cl = 2000 m , ensemble spread increases with increasing cl as indicated by the longer whiskers. As said before, the ensemble average of the entropy production number is constant for larger cl . Howe ver , analysis of all the ensembles shows oscillation of the ensemble spread for larger hori- zontal correlation lengths. Outliers (i.e., above 1.5 times the upper quartile or below 1.5 times the lower quartile) are mostly found at the lower end of the ensembles (diamonds in Fig. 10). This may be linked directly to the top and bottom Dirichlet boundary condi - tions for temperature, meaning that there is a fixed upper limit for heat transport in the model, thus reducing the likelihood of “positive” outliers. Simulation results suggest that the increase in horizontal correlation length does not have a strong influence once a convection system has formed. We do not observe a systematic connection between cl and ensemble variance, which is expressed by the varying whisker lengths for boxplots with increasing cl in Fig.  10. It should be noted, however, that variations in whisker lengths may also be due to the ensemble size, so con- clusions should be made cautiously. The ensemble median, however, is relatively con - stant over this range of cl , indicating no significant shift of the ensemble distribution, but a change in its width, i.e., its standard deviation. Influence of vertical correlation length In a sedimentary environment, lithofacies and associated petrophysical parameters do not only vary laterally, but also vertically. Such vertical variations due to layering have a stronger influence on the resulting vertical heat transport than lateral variations. This can be illustrated by comparison with a simplified circuit model. In case of variations along the transport path, i.e., vertical variations when looking at vertical heat transport, resistors in the electric circuit would be connected in series. That is, adding a resistor with a high resistivity would already change the resulting voltage significantly. For lat - eral variations, the analog circuit would be resistors connected in parallel, and the total Niederau et al. Geotherm Energy (2019) 7:35 Page 18 of 27 Fig. 11 Schematic showing the significant influence of vertical variations along the transport path compared to lateral variations, and analog electrical circuits as examples. Arrows represent the relative proportions of q passing through the material, not its magnitude resistance is lower than the smallest individual resistor in the system. Here, an additional resistor in series with a high resistivity would have a less significant influence. Figure  11 illustrates this analog. As a result, one would use the harmonic mean for averaging vertical variations along the transport path, and the arithmetic mean for averaging the lateral variations. Accordingly, changing the vertical correlation length of permeability should have a more significant influence on the associated convection system than changing the horizontal correlation length. To test this hypothesis, we performed the same Monte Carlo simulations with cl equal to 100  m. Results show that the vertical correlation length ( cl ) of permeability has a significant effect on the likelihood of convective heat transport. Again, increasing horizontal correlation length of permeability correlates with increasing number of realizations showing no convective heat transport for both vertical correlation lengths (Fig.  12). While over 80% of members in ensembles with cl = 150 m show convection, this fraction decreases drastically for ensembles with cl = 100 m down to < 50 % . For both vertical correlation lengths, the number of con- vecting ensemble members mainly decreases between 1500 and 3500 m horizontal correlation length. In ensembles with cl larger than 3500 m, the fraction of members showing convection stabilizes for cl = 100 m , while it continuously slightly decreases for cl = 150 m. Fig. 12 Variation of the fraction of realizations with free convection N with horizontal correlation convection length cl for two different cl . Ensembles with cl = 100 m show a significantly stronger decrease in h v v convecting realizations, compared to ensembles with cl = 150 m v Niederau et al. Geotherm Energy (2019) 7:35 Page 19 of 27 These results emphasize how important knowledge about the (vertical) spatial continuity of permeability in a modeled aquifer is for quantifying the probability of convective heat transport. Furthermore, we observe that cl has a stronger influence on the development of convection cells than cl , which does not influence the heat transfer, expressed by entropy production, over a certain range. This observation has important implications, as the vertical correlation length cl can usually be deduced from well logs. A shorter vertical correlation length means more frequent alterna- tion of different permeabilities, meaning more little permeable zones, which strongly influence the average permeability of the sedimentary column, as it is averaged by the harmonic mean. Number of convection cells and entropy production In the context of this work, one convection cell is understood as the connected pair of zones of upflow and downflow. To test the hypothesis that spatially heterogeneous permeability influences the convection pattern, i.e., the number of convection cells, we count the convection cells in each simulated model. For this purpose, we apply a peak count method to horizontal temperature profiles through the model. Each counted posi - tive peak can be associated with a zone of upflow. Before presenting the results of the cross section of the Perth Basin, we show, as an example, the influence of cl on the number of convective rolls in closed boxes of dimensions 2 km (x-direction) ×1 km (z-direction) in Fig.  13. In these simplified models without complicated reservoir structure, the change in the number of convec- tion cells is solely controlled by the change in horizontal correlation length of perme- ability. Figure 13 shows isotherms for three different cl . At cl = 50 m , three distinct h h convection cells develop. Increasing cl to 250 m and 500 m yields a reduction in the number of convection cells caused by wider zones of low permeability (gray-to-white fields in Fig.  13). Results of this analysis for the Perth Basin ensembles are presented as a joint plot of correlation length and number of convection cells (wave number) in Fig.  14. A joint plot shows the bivariate distribution between two variables, here cl and wave number, as well as their univariate distributions at the margin. The number of convection cells decreases at shorter horizontal correlation lengths, visualized by the darker colors in Fig. 14. Darker blue indicates that a wave number was counted in Fig. 13 Simulation of convection in a closed box with heterogeneous permeability. Horizontal correlation length cl of permeability from left to right: 50 m, 250 m, and 500 m. Permeability is displayed by gray fields, while isotherms are presented as colored lines Niederau et al. Geotherm Energy (2019) 7:35 Page 20 of 27 Fig. 14 Heatmap of wave number (number of convection cells) and horizontal correlation length for the ensembles. Color intensity represents the amount of ensemble members showing the according number of convection cells with lighter colors correlating to more models in a specific bin more ensemble members. At cl = 500 m , the majority of simulations have eight con- vection cells. At cl = 1000 m , the average wave number decreased to seven, and to six at cl = 4000 m . This observation agrees with results from simpler models, where no major changes in the number of convection rolls are observed for larger cl . Simil arly to the entropy production number (Fig. 10), the number of convection cells stabilizes at longer correlation lengths. While 12 peaks are counted in the model with homo- geneous permeability, the number of peaks undulates between 6 and 7 for ensembles with a correlation length greater than 1500 m. That is, at 1500 m, the average number of convection cells in the model is halved compared to a simulation with homogene- ous permeability. The similarity between the distribution of the entropy production number and the wave number suggests that the former can be related to the probable number of convection cells, given a closed model with a fixed size. Changes in the model size, i.e., its aspect ratio, also influence the number of convection cells (Beck 1972; Bejan 1979; Turcotte and Schubert 2014; Bejan 2013). Discussion The presented results are analyzed following the hypothesis that convective systems are less likely to develop when we increase the finite horizontal correlation length in an aquifer, and that this can be measured by the entropy production number. This hypothesis may need further explanation. We assume that all parameters relevant for convection, such as temperature difference across the porous medium, or its thick - ness, are kept constant. A model with homogeneous permeability can be treated as a case with an infinite spatial correlation length. In such a model, only the magni - tude of permeability determines if convection is likely or not. If we consider the same model with heterogeneous permeability, i.e., a finite correlation length, the magni - tude of permeability and its spatial variation would control the onset of free convec- tion. It may appear logical that longer correlation lengths might facilitate the onset of Niederau et al. Geotherm Energy (2019) 7:35 Page 21 of 27 convection, as permeable areas stretch over a longer distance. However, so do areas of low permeability. These work as barriers to the buoyancy-driven convective cur - rent, which is mainly vertical. At shorter finite horizontal correlation lengths, these currents just get diverted around lenses of low permeability, whereas a lateral elonga- tion of those lenses can completely inhibit the onset of convection. And such a lateral elongation corresponds to a longer horizontal correlation length. One controlling factor in a convective system is the temperature difference between top and bottom of the convecting layer. Consequently, the type of boundary conditions applied to a model influences the simulated convection. We apply Dirichlet boundary conditions at the top and bottom boundaries of the model. Thus, the applied boundary conditions supply a constant amount of energy to the system. It needs to be discussed whether these boundary conditions reflect a realistic state of the system. If tempera - tures at depth are known, applying Dirichlet boundary conditions is reasonable. How- ever, as this is usually not the case, thermal Neumann boundary conditions are applied as bottom boundary condition (Smith and Chapman 1983). In the Perth Basin region, several studies use a basal Dirichlet boundary condition for temperature and yield real- istic simulation results (Corbel et  al. 2012b; Schilling et  al. 2013; Niederau et  al. 2017). Their results suggest that, for the bottom boundary, a Dirichlet boundary condition for temperature is suitable for the Perth Basin. Still, it is worthy to discuss the implications which the choice of boundary condition has on the entropy production. Constant force boundary conditions (Dirichlet boundary conditions) yield a maximum of estimated entropy production and constant flux boundary conditions (Neumann boundary condi - tions) define a lower bound. We focused here on maximum conditions to obtain a com - parison to the previous work on cell structures (Börsing et al. 2017). Furthermore, under the scope of this paper, i.e., the possibility to use entropy production to obtain insight in the behavior of convection in heterogeneous systems, we generally expect similar results concerning the geometric interpretation of the entropy production. Additionally to permeability of a reservoir, its structure, i.e., abrupt changes in reser- voir thickness, has a strong influence on the convective system. Such changes are usually caused by fault zones which not only influence reservoir structure, but normally also the spatial permeability distribution, and can thus influence groundwater flow paths (Frey - mark et al. 2019; Inbar et al. 2019). In terms of their permeability relative to that of the surrounding rock, fault zones may be more permeable, similarly permeable, or less per- meable. For instance, fault-related processes, such as clay smear (Vrolijk et al. 2016), can drastically influence permeability of a fault, e.g., by acting as a barrier to fluid flow. The existence of clay smear alone, however, is not a defining indicator on the permeability of the fault, as potential breaching of the clay layer coincides with increased flow across the fault (Kettermann et  al. 2017). Furthermore, experimental data showed that the influ - ence of clay smear on permeability changes with fault activity, i.e., rock deformation (Takahashi 2003). We only consider fault displacement in our simulations which over- simplifies the impact of fault systems. Schilling et al. (2013) observe that zones of upflow occur in fault zones or in their vicinity, depending if the faults are highly permeable or impermeable respectively. While it is logical that upflow is focused within permeable fault zone, it can appear counterintuitive that zones of upflow occur next to fault zones if these are impermeable. However, aquifer thickness of the Yarragadee Aquifer is highest Niederau et al. Geotherm Energy (2019) 7:35 Page 22 of 27 in the hanging walls of normal faults. Thus, formation of convection cells in these areas is facilitated, especially when fluid flow across a fault is inhibited by its low permeability. However, fault zones are not homogeneously permeable over their whole extent, mean- ing that the results presented by Schilling et  al. (2013) reflect extreme conditions. But even if faults are modeled just as displacement, they have an influence on the convective pattern (Niederau et al. 2017). It is not known whether faults in the central Perth Basin are sealing or transmissive. Faults studied in the northern Perth Basin suggest that they may be sealing (Olierook et al. 2014), but this cannot be assumed to be true in the cen- tral Perth Basin. Considering the great uncertainty in fault permeability in general, and also in the Perth Basin, we decided to model only the fault displacement, as our study focuses on the influence of the correlation length of matrix permeability on the convec - tive system, and not on the influence of fault permeability. It would, however, be possible and important, to extend the analysis to simulations which consider fluid flow in and across fault zones, for example by representing faults as 2D elements in a finite-element discretization. The presented diagnostic equations for the entropy production would be equally applicable. Our simulation results suggest that varying the correlation length of permeability in the Yarragadee Aquifer changes the convection pattern in the aquifer. These results agree with the theoretical studies on the influence of the spatial continuity on the onset of convection (Nield et al. 2010; Simmons et al. 2010). For instance, Nield and Simmons (2007) state that weak vertical heterogeneity has presumably little or no effect on the onset of convection. In our simulations, reducing the vertical correlation length by 50 m caused a significant reduction of ensembles showing convection at higher horizon - tal correlation lengths. By reducing the vertical correlation length, we increase the ver- tical heterogeneity. However, the variance of the permeability distribution applied in our model (see Fig.  5 right) is too large to be considered weak heterogeneity after the definition provided by Nield and Simmons (2007) and Simmons et al. (2010). Thus, the influence of vertical heterogeneity is more significant in our study, but this does not con - tradict previous studies. In the case of the Perth Basin, it can further be concluded that a horizontal correla- tion length greater than 2000 m has a small influence on the convective pattern. Such a stabilization in convective systems was also found by Prasad and Simmons (2003) who observed a decrease of overall degree of instability with increasing horizontal correlation length. While the model which they studied is a heterogeneous Elder Problem (Elder 1967), their conclusions can be confirmed to some extent in this study with a geometri - cally more complex model. Furthermore, our findings coincide with those by Prasad and Simmons (2003) that reservoir models using representative averages or flow parameters will often significantly underestimate the importance of their spatial heterogeneities, e.g., their effect on the onset of convection. We use a geostatistical approach for varying horizontal correlation lengths, as we have no real information about the latter in the Yarragadee Aquifer. Cross correlation between wells is a usual approach for assessing a horizontal correlation length. How- ever, the distance between wells in the study area is between 40 and 100 km, thus too long to assess a meaningful horizontal correlation length. The vertical correlation length, however, as well as the value range of porosity and permeability used in our models are Niederau et al. Geotherm Energy (2019) 7:35 Page 23 of 27 representative for the Yarragadee Aquifer, as it is based on core- and log data from wells in the Perth Basin (Niederau et al. 2017). Another approach for modeling the spatial per- meability distribution would be stratigraphic forward modeling. Existing stratigraphic models for the Yarragadee Aquifer by Corbel et al. (2012a) provide a good basis for esti- mating a range of correlation lengths and a different approach for modeling a spatial permeability distribution. In their model, sediment influx in the Perth Basin was from the north and east. Despite yielding a probable permeability distribution, stratigraphic forward modeling in such a sedimentary environment is non-trivial and model results are just one possible realization of the sediment distribution in the subsurface. As such, it is comparable to geostatistical realizations of the same distribution. The key difference between the two approaches is that stratigraphic forward modeling is a process-based approach, while geostatistics is a data-based one. In contrast to analyzing one strati- graphic forward model, a Monte Carlo approach enables us to draw statistically more robust conclusions. Ideally, both approaches should be combined. That is, conditioned geostatistical ensembles of permeability distributions based on stratigraphic forward models could provide a model representative for the Yarragadee Aquifer. A classic approach in studying convecting system is assessing the systems’ Rayleigh number. While the Rayleigh number is very useful for estimating if a system is likely to show convection, it does not provide information about the likely convective pattern. The same Rayleigh number may characterize different convection patterns which have different entropy production numbers (Börsing et al. 2017). Also, assessing a meaning - ful Rayleigh number in a non-idealized setting (e.g., a convecting layer with non-parallel boundaries) is difficult, as the height of the convecting layer, or the temperature differ - ence across it is two input parameters in classical Rayleigh number formulations. Both parameters are not constant in a real-world example. For instance, estimating the Ray- leigh number of a simple model similar to the one shown in Fig.  1a is straightforward, as is the estimation of N . Introducing lenses of low permeability as in Fig.  1b changes the convective system. Such a change in the average permeability is represented by a change in the Rayleigh number and N of the model. If introduction of heterogeneity in permeability is not reflected by a change in the average permeability, a correct Ray - leigh number analysis gets problematic, especially if the harmonic mean of the vertical permeability of the studied aquifer is subcritical, as in Fig. 1c, d. Such a conclusion, that traditional Rayleigh analysis can be misleading, as convection can occur at subcritical Rayleigh numbers (due to subcritical average permeability), was also drawn by Prasad and Simmons (2003), who stated “the traditional Rayleigh stability criterion based upon an average permeability is inadequate for describing solute transport processes in het- erogeneous systems”. A traditional Rayleigh number analysis may only be valid when heterogeneity in permeability is small, i.e., a deduced log-permeability distribution has a small standard deviation. In a system with spatially heterogeneous permeability, a Rayleigh number analysis may lead to incorrect predictions, as low-permeable layers in a system significantly reduce the average vertical permeability (e.g., when using the harmonic mean) and thereby the estimated Rayleigh number. Thus, entropy production can be used as a thermodynam - ically-based, diagnostic metric which may vary in systems of equal Rayleigh number and thus represent different convection patterns. This is one advantage of the entropy Niederau et al. Geotherm Energy (2019) 7:35 Page 24 of 27 production number compared to a Rayleigh number assessment. Concerning the Nus- selt number, entropy production has the advantage that a volumetric estimation of con- vection is possible, due to the formulation by Bejan (2013). Furthermore, the entropy production number can be used for assessing a probable convection pattern in a res- ervoir of more complex geometry. In our simulations, the number of convection cells decreased with increasing horizontal correlation length of permeability, at least at lower correlation lengths. The average entropy production number similarly decreased at lower correlation lengths, suggesting that it correlates with the number of convec- tion cells. It should, however, be noted that entropy production provides the potential to observe whether major changes in a convective system occur. It does not directly help in locating zones of up flow or down flow in a realistic environment. In this sense, we can - not expect the measure of entropy production to indicate regions of upflow in an area. However, it provides an estimation about the expected number, and accuracy, of up flow and down flow regions in the entire simulation domain. Conclusion In this manuscript, we analyze how different correlation lengths in permeability affect a convection pattern by using a Monte Carlo approach. We use the measure of entropy production for quantifying this effect. Simulation results of the MC ensembles show that the number of ensemble members with dominant convective heat transport decreases with increasing correlation length, as expressed by a continuously decreasing entropy production number ( N ). Excluding non-convective ensemble members reveals that N ˙ ˙ S S stabilizes at a horizontal correlation length of around 2000 m. As expected, the num- ber of convection cells also decreases with increasing correlation length, similarly to the entropy production number, also stabilizing at a horizontal correlation length of around 2000 m. Vertical correlation length has a stronger influence on the possibility of convec - tion development in our model. This result has important implications for similar case studies, as deducing the vertical correlation length from well logs is usually feasible in reservoir studies. The measure presents itself as a powerful diagnostic tool when studying hydrothermal convection in real-world reservoir systems, and that in such settings, it can be more use- ful than a classical Rayleigh number assessment. It is not fully clear which conditions in our simulations facilitate the development of a convective system. It is likely that certain combinations of heterogeneous permeability and reservoir structure are an important factor for convection in our model of the Perth Basin. The evaluation of entropy produc - tion shows further that it is important to consider spatially heterogeneous permeability with an accurate correlation length to reliably simulate the probable convection present in the Perth Basin. We consider the use of entropy production as a diagnostic tool and focus here on 2-D analyses, as an extension and application to previous theoretical work (Börsing et  al. 2017). However, one strength of the measure itself is that it can be applied to any more complex setting,  both in the choice of boundary conditions, as well as in the dimen- sionality of the problem. Including faults as own hydraulic units in a finite-element dis - cretization, for instance, the influence of a constant flux boundary condition, or how entropy production is affected by the development of elongated convection rolls, rather Niederau et al. Geotherm Energy (2019) 7:35 Page 25 of 27 than cells, in three-dimensional problems. We already performed the first studies inves - tigating the development entropy production and its interpretability in three-dimen- sional settings (Huang 2017) and see this as an interesting path for future research. Acknowledgements The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium on the supercom- puter JURECA at Forschungszentrum Jülich for the JARA-HPC project 11061. The authors thank the reviewers and editor for their valuable comments which helped improving the manuscript. Authors’ contributions JN sets up the geological model, did the ensemble simulation and majority of analysis. He wrote the majority of the manuscript. JFW was a major contributor in structuring and revising the manuscript. NB significantly contributed to the manuscript, and provided the basis code for the entropy production analysis. All authors read and approved the final manuscript. Funding The conducted research was not part of a funded project. Computing time on the supercomputer JURECA were Granted by JARA-HPC under the project number 11061. Availability of data and materials The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request. Competing interests The authors declare that they have no competing interests. Author details Institute for Applied Geophysics and Geothermal Energy, E.ON Energy Research Center, RWTH Aachen University, 2 3 Aachen, Germany. Jülich-Aachen Research Alliance (JARA-HPC), RWTH Aachen University, Aachen, Germany. Present Address: Geothermal Energy and Geofluids Group, Department of Earth Sciences, ETH Zürich, Zurich, Switzerland. Com- putational Geoscience and Reservoir Engineering, RWTH Aachen University, Aachen, Germany. Institute of Geophysics, Department of Earth Sciences, ETH Zürich, Zurich, Switzerland. Received: 2 August 2019 Accepted: 8 November 2019 References Alpay OA, et al. A practical approach to defining reservoir heterogeneity. J Petrol Technol. 1972;24(07):841–8. Backhouse J. Revised late Jurassic and early Cretaceous stratigraphy in the Perth Basin. Geol Surv Western Aust Rep. 1984;12:1–6. Barbier E. Geothermal energy technology and current status: an overview. Renew Sustain Energy Rev. 2002;6(1):3–65. https ://doi.org/10.1016/S1364 -0321(02)00002 -3. Beck JL. Convection in a box of porous material saturated with fluid. Phys Fluids. 1972;15(8):1377–83. Bejan A. A study of entropy generation in fundamental convective heat transfer. J Heat Transf. 1979;101(4):718–25. Bejan A. Convection heat transfer. Hoboken: Wiley; 2013. Börsing N, Wellmann JF, Niederau J, Regenauer-Lieb K. Entropy production in a box: analysis of instabilities in confined hydrothermal systems. Water Resour Res. 2017;53(9):7716–39. https ://doi.org/10.1002/2017W R0204 27. Botman C. Determining the geothermal properties of the Perth Basin lithology to aid geothermal energy projects. Perth: Curtin University of Technology B.Sc (Geophysics) Honours thesis Unpublished. 2010. Clauser C. Thermal storage and transport properties of rocks, II: thermal conductivity and diffusivity. Dordrecht: Springer; 2011. p. 1431–48. https ://doi.org/10.1007/978-90-481-8702-7_67. Clauser C, Bartels J. Numerical simulation of reactive flow in hot aquifers: SHEMAT and processing SHEMAT. Berlin: Springer; 2003. Corbel S, Griffiths C, Dyt C, Ricard L. Hot sedimentary aquifer characterization using stratigraphic forward modelling, Perth Basin, Australia. In: Proceedings of the thirty-seventh workshop on geothermal reservoir engineering, Stanford University, Stanford, California. 2012a. Corbel S, Schilling O, Horowitz F, Reid L, Sheldon H, Timms N, Wilkes P. Identification and geothermal influence of faults in the Perth metropolitan area, Australia. In: Thirty-seventh workshop on geothermal reservoir engineering, Stanford, CA. 2012b. Costa V. On natural convection in enclosures filled with fluid-saturated porous media including viscous dissipation. Int J Heat Mass Transf. 2006;49(13–14):2215–26. Crostella A, Backhouse J. Geology and petroleum exploration of the central and southern Perth Basin, Western Australia. 57, Geological Survey of Western Australia. 2000. CSIRO Csiro Pressureplot(tm) software and Pressuredb (tm) database. 2007. www.press urepl ot.com. Accessed 08 Apr Davidson WA, Yu X. Perth regional aquifer modelling system (PRAMS) model development: Hydrogeology and ground- water modelling. Western Australia Department of Water. 2008. Delle Piane C, Esteban L, Timms L, Ramesh Israni S. Physical properties of mesozoic sedimentary rocks from the Perth Basin, Western Australia. Aust J Earth Sci. 2013;60(6–7):735–45. Niederau et al. Geotherm Energy (2019) 7:35 Page 26 of 27 Deutsch C, Journel A. GSLIB Geostatistical software library and user’s guide. New York: Oxford University Press; 1998. Diersch HJG, Kolditz O. Variable-density flow and transport in porous media: approaches and challenges. Adv Water Resour. 2002;25(8):899–944. https ://doi.org/10.1016/S0309 -1708(02)00063 -5. Elder J. Transient convection in a porous medium. J Fluid Mech. 1967;27(3):609–23. Freymark J, Bott J, Cacace M, Ziegler M, Scheck-Wenderoth M. Influence of the main border faults on the 3d hydraulic field of the central Upper Rhine Graben. Geofluids. 2019;. https ://doi.org/10.1155/2019/75207 14. Grant M, Donaldson I, Bixley P. Geothermal reservoir engineering. New York: Academic Press; 1982. https ://doi. org/10.1016/B978-0-122-95620 -1.X5001 -2. Harris L. Structural and tectonic synthesis for the Perth Basin, Western Australia. J Petrol Geol. 1994;17(2):129–56. HDRPL Geothermal energy potential in selected areas of Western Australia (Perth Basin). Hot Dry Rocks Pty Ltd. (HDRPL), South Yarra, Australia. 2008. Horowitz F, Regenauer-Lieb K, Wellmann J, Chua H, Wang X, Poulet T. Evidence for hydrothermal convection in the Perth Basin, Australia. In: Gurgenci H, Budd AR, editors. Proceedings of the Sir Mark Oliphant international frontiers of sci- ence and technology Australian geothermal energy conference. 2008;18. Huang P. Investigating different formulations for hydrothermal convection in geothermal systems. Master’s thesis, RWTH Aachen University. 2017. Inbar N, Rosenthal E, Magri F, Alraggad M, Möller P, Flexer A, Guttman J, Siebert C. Faulting patterns in the lower yarmouk gorge potentially influence groundwater flow paths. Hydrol Earth Syst Sci. 2019;23(2):763–71. Irvine DJ, Sheldon HA, Simmons CT, Werner AD, Griffiths CM. Investigating the influence of aquifer heterogeneity on the potential for thermal free convection in the Yarragadee Aquifer, Western Australia. Hydrogeol J. 2014;. https ://doi. org/10.1007/s1004 0-014-1194-1. Jureca Jülich Supercomputing Centre. General-purpose supercomputer at jülich supercomputing centre. J Large-Scale Res Facil. 2016;2:A62. https ://doi.org/10.17815 /jlsrf -2-121. Kettermann M, Urai JL, Vrolijk PJ. Evolution of structure and permeability of normal faults with clay smear: insights from water-saturated sandbox models and numerical simulations. J Geophys Res Solid Earth. 2017;122(3):1697–725. Kleidon A, Lorenz R. Entropy production by earth system processes. In: Kleidon A, Lorenz R, editors. Non-equilibrium thermodynamics and the production of entropy. Heidelberg: Springer; 2005. p. 1–20. Kolditz O, Ratke R, Diersch HJG, Zielke W. Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models. Adv Water Resour. 1998;21(1):27–46. Lipsey L, Pluymaekers M, Goldberg T, van Oversteeg K, Ghazaryan L, Cloetingh S, van Wees JD. Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands. Geothermics. 2016;64:135–51. Lorenz EN. Deterministic nonperiodic flow. J Atmos Sci. 1963;20(2):130–41. Middleton M. Radiogenic heat generation in the Darling Range, Western Australia. Explor Geophys. 2013;44:206–14. Niederau J. Calibration of a fractal model relating porosity to permeability and its use for modeling hydrothermal transport processes in the Perth Basin, Australia. Energy Procedia. 2014;59:293–300. https ://doi.org/10.1016/j.egypr o.2014.10.380. Niederau J, Börsing N, Wellmann F, Clauser C. Entropy production and attractors: Measures to quantify uncertainty and complexity introduced by convection. In: Proceedings European geothermal congress, 19–24 September, Stras- bourg, France. 2016. Niederau J, Ebigbo A, Marquart G, Arnold J, Clauser C. On the impact of spatially heterogenous permeability on free convection in the Perth Basin, Australia. Geothermics. 2017;66:119–33. Nield D, Simmons CT. A discussion on the effect of heterogeneity on the onset of convection in a porous medium. Transp Porous Media. 2007;68(3):413–21. Nield D, Kuznetsov A, Simmons CT. The effect of strong heterogeneity on the onset of convection in a porous medium. Transp Porous Media. 2009;77(2):169. Nield D, Kuznetsov A, Simmons CT. The effect of strong heterogeneity on the onset of convection in a porous medium: 2d/3d localization and spatially correlated random permeability fields. Transp Porous Media. 2010;83(3):465–77. Olierook HK, Timms NE, Hamilton PJ. Mechanisms for permeability modification in the damage zone of a normal fault, northern Perth Basin, Western Australia. Marine Petrol Geol. 2014;50:130–47. https ://doi.org/10.1016/j.marpe tgeo.2013.10.012. Pape H, Clauser C, Iffland J. Permeability prdiction based on fractal pore-space geometry. Geophysics. 1999;64(5):1447–60. Playford PE, Cockbain AE, Low G. Geology of the Perth Basin, Western Australia. Geological Survey of Western Australia, Perth, Australia. 1976. Prasad M. Velocity-permeability relations within hydraulic units. Geophysics. 2003;68(1):108–17. Prasad A, Simmons CT. Unstable density-driven flow in heterogeneous porous media: a stochastic study of the elder “short heater” problem. Water Resour Res. 2003;39(1):SBH 4–1–SBH 4–21. https ://doi.org/10.1029/2002W R0012 90. Pujol M, Ricard LP, Bolton G. 20 years of exploitation of the Yarragadee Aquifer in the Perth Basin of Western Australia for direct-use of geothermal heat. Geothermics. 2015;57:39–55. https ://doi.org/10.1016/j.geoth ermic s.2015.05.004. Rana R, Horne R, Cheng P. Natural convection in a multi-layered geothermal reservoir. J Heat Transf. 1979;101(3):411–6. Rath V, Wolf A, Bücker H. Joint three-dimensional inversion of coupled groundwater flow and heat transfer based on automatic differentiation: sensitivity calculation, verification, and synthetic examples. Geophys J Int. 2006;167(1):453–66. Rayleigh L. On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Lond Edinburgh Dublin Philos Mag J Sci. 1916;32(192):529–46. Rubin Y, Mavko G, Harris J. Mapping permeability in heterogeneous aquifers using hydrologic and seismic data. Water Resour Res. 1992;28(7):1809–16. Schilling O, Sheldon H, Reid L, Corbel S. Hydrothermal models of the perth metropolitan area, Western Australia: implica- tions for geothermal energy. Hydrogeol J. 2013;21(3):605–21. https ://doi.org/10.1007/s1004 0-012-0945-0. Sheldon HA, Florio B, Trefry MG, Reid LB, Ricard LP, Ghori KAR. The potential for convection and implications for geother- mal energy in the Perth Basin, Western Australia. Hydrogeol J. 2012;20(7):1251–68. Niederau et al. Geotherm Energy (2019) 7:35 Page 27 of 27 Siavashi M, Bahrami HRT, Saffari H. Numerical investigation of flow characteristics, heat transfer and entropy generation of nanofluid flow inside an annular pipe partially or completely filled with porous media using two-phase mixture model. Energy. 2015;93:2451–66. Siavashi M, Bordbar V, Rahnama P. Heat transfer and entropy generation study of non-darcy double-diffusive natural convection in inclined porous enclosures with different source configurations. Appl Therm Eng. 2017;110:1462–75. https ://doi.org/10.1016/j.applt herma leng.2016.09.060. Simmons CT, Kuznetsov AV, Nield DA. Eec ff t of strong heterogeneity on the onset of convection in a porous medium: importance of spatial dimensionality and geologic controls. Water Resour Res. 2010;. https ://doi.org/10.1029/2009W R0086 06. Smith L, Chapman DS. On the thermal effects of groundwater flow: 1. Regional scale systems. J Geophys Res Solid Earth (1978–2012). 1983;88(B1):593–608. Takahashi M. Permeability change during experimental fault smearing. J Geophys Res Solid Earth. 2003. https ://doi. org/10.1029/2002J B0019 84. Testerman J, et al. A statistical reservoir-zonation technique. J Petrol Technol. 1962;14(08):889–93. Timms NE, Olierook HK, Wilson ME, Delle Piane C, Hamilton PJ, Cope P, Stütenbecker L. Sedimentary facies analysis, mineralogy and diagenesis of the mesozoic aquifers of the central Perth Basin, Western Australia. Marine Petrol Geol. 2015;60:54–78. Tolman RC, Fine PC. On the irreversible production of entropy. Rev Mod Phys. 1948;20(1):51. Turcotte D, Schubert G. Geodynamics. Cambridge: Cambridge University Press; 2014. Vincent AP, Yuen DA. Thermal attractor in chaotic convection with high-prandtl-number fluids. Phys Rev A. 1988;38:328– 34. https ://doi.org/10.1103/PhysR evA.38.328. Vrolijk PJ, Urai JL, Kettermann M. Clay smear: review of mechanisms and applications. J Struct Geol. 2016;86:95–152. Wagner W, Kretzschmar HJ. Iapws industrial formulation 1997 for the thermodynamic properties of water and steam. In: International steam tables: properties of water and steam based on the industrial formulation IAPWS-IF97 2008. p. 7–150. Wyllie MRJ, Gregory AR, Gardner LW. Elastic wave velocities in heterogeneous and porous media. Geophysics. 1956;21(1):41–70. Zyvoloski G, Robinson B, Dash Z, Trease L. Models and methods summary for the FEHMN application. Tech. rep., Los Alamos National Lab., NM (United States). 1996. https ://www.osti.gov/servl ets/purl/36645 5. Accessed 12 Oct 2016. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geothermal Energy Springer Journals

Analyzing the influence of correlation length in permeability on convective systems in heterogeneous aquifers using entropy production

Geothermal Energy , Volume 7 (1) – Nov 21, 2019

Loading next page...
 
/lp/springer-journals/analyzing-the-influence-of-correlation-length-in-permeability-on-U4L0iKfJet
Publisher
Springer Journals
Copyright
Copyright © 2019 by The Author(s)
Subject
Earth Sciences; Geotechnical Engineering & Applied Earth Sciences; Renewable and Green Energy; Geoecology/Natural Processes
eISSN
2195-9706
DOI
10.1186/s40517-019-0151-6
Publisher site
See Article on Publisher Site

Abstract

jan.niederau@erdw.ethz.ch Institute for Applied Hydrothermal convection in porous geothermal reservoir systems can be seen as Geophysics and Geothermal a double-edged sword. On the one hand, regions of upflow in convective systems Energy, E.ON Energy Research Center, RWTH can increase the geothermal energy potential of the reservoir; on the other hand, Aachen University, Aachen, convection introduces uncertainty, because it can be difficult to locate these regions Germany of upflow. Several predictive criteria, such as the Rayleigh number, exist to estimate Full list of author information is available at the end of the whether convection might occur under certain conditions. As such, it is of interest article which factors influence locations of upwelling regions and how these factors can be determined. We use the thermodynamic measure entropy production to describe the influence of spatially heterogeneous permeability on a hydrothermal convection pattern in a 2D model of a hot sedimentary aquifer system in the Perth Basin, Western Australia. To this end, we set up a Monte Carlo study with multiple ensembles. Each ensemble contains several hundred realizations of spatially heterogeneous permeabil- ity. The ensembles only differ in the horizontal spatial continuity (i.e., correlation length) of permeability. The entropy production of the simulated ensembles shows that the convection patterns in our models drastically change with the introduction and increase of a finite, lateral correlation length in permeability. An initial decrease of the average entropy production number with increasing lateral correlation length shows that fewer ensemble members show convection. When neglecting the purely conduc- tive ensembles in our analysis, no significant change in the number of convection cells is seen for lateral correlation lengths larger than 2000 m. The result suggests that the strength of convective heat transfer is not sensitive to changes in lateral correlation length beyond a specific factor. It does, however, change strongly compared to simula- tions with a homogeneous permeability field. As such, while the uncertainty in spatial continuity of permeability may not strongly influence the convective heat transfer, our findings show that it is important to consider spatial heterogeneity and continuity of permeability when simulating convective heat transfer in an aquifer. Keywords: Hydrothermal convection, Entropy production, Heterogeneous permeability, Perth Basin, Spatial continuity © The Author(s) 2019. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Niederau et al. Geotherm Energy (2019) 7:35 Page 2 of 27 Introduction Convective heat transfer in sedimentary reservoir systems can significantly increase their geothermal potential, as hot water is transported to shallower depths by zones of upflow (Barbier 2002). However, this increased potential is connected to a similarly increased uncertainty, as these zones can be difficult to locate due to uncertain reser - voir structure and properties (Rana et  al. 1979; Schilling et  al. 2013; Lipsey et  al. 2016; Niederau et al. 2017), or as they are influenced active exploitation (Grant et al. 1982). On the other hand, convection can increase the risk of success for geothermal projects, e.g., when boreholes are drilled in a zone of downflow. As such, knowledge about the state of the convection system and its pattern is crucial for assessing the uncertainty and associ- ated risk of an exploration study. Several parameters, such as aquifer thickness or perme- ability, control whether free convection can occur. However, those “global” parameters do not give unique insight into spatial distribution of up- and downwelling zones. It is, therefore, necessary to incorporate more detailed knowledge about these parameters for predicting the hydrothermal potential of (sedimentary) reservoir systems with convec- tion. While the spatial extent of a reservoir unit can be derived from geophysical meas- urements to a fine degree of detail, detailed information about the spatial distribution of permeability is often not available. It is usually also difficult to assess, as permeability of sedimentary reservoir systems correlates to the depositional paleo-environment. Depo- sitional structures may be inferrable from seismics to some extent, and seismic proper- ties of rocks can be related to their hydraulic properties (Rubin et al. 1992). But even if such methods are successful, an uncertainty in spatial permeability distribution remains. The Perth Basin, a sedimentary rift basin in Western Australia, was extensively studied for its geothermal potential with focus on one prominent reservoir unit, the Yarragadee Aquifer. While the hydrothermal conditions in this aquifer are largely unknown (Pujol et  al. 2015), multiple modeling studies suggest that hydrothermal convection occurs within the Yarragadee Aquifer (Sheldon et al. 2012; Schilling et al. 2013; Niederau et al. 2017). In particular, Sheldon et  al. (2012) show that temperature and salinity measure- ments in the Perth Basin indicate convective flow. For instance, they show that laterally varying temperature gradients which were deduced from corrected bottom hole tem- peratures, are indicative for regions of upflow or downflow, respectively. Furthermore, zones of uniform salinity were measured in the Yarragadee Aquifer. Such homogeniza- tion of salinity can be caused by convective flow. Horowitz et  al. (2008) present three arguments indicating convection in the Perth Basin: (i) anomalously low thermal gradi- ents in deep aquifers, (ii) thermal gradients varying from 28 to 68 K/km within a region of 55 km × 50 km , (iii) Rayleigh numbers above the critical one, even if conservative parameters are used for the Yarragadee Aquifer when calculating the Rayleigh number. Timms et  al. (2015) classified available core material from wells to several lithostrati - graphic units corresponding to different depositional paleoenvironments in the Perth Basin. While they found a general correlation between lithostratigraphy and average permeability, the spatial correlation of permeability in the Yarragadee Aquifer remains highly uncertain. However, for characterizing the convection pattern, knowledge about heterogeneity of permeability and its spatial continuity is crucial, as both influence the convection system in an aquifer (Nield and Simmons 2007; Nield et  al. 2009). For instance, Simmons et al. (2010) conclude that the scale of heterogeneity, i.e., the size of Niederau et al. Geotherm Energy (2019) 7:35 Page 3 of 27 standard deviation of the underlying permeability distribution, has a strong impact on the nucleation of a convection system. They observe further that a convection system is moderately sensitive to changes in the correlation length. With uncertain spatial heterogeneity of permeability, an associated convection pat- tern is also uncertain. For assessing the geothermal potential of the Perth Basin with a focus on the spatial variability of a convection pattern, it is, therefore, necessary to address heterogeneous permeability and its spatial correlation in the Yarragadee Aquifer. Based on stratigraphic forward modeling, Irvine et  al. (2014) created Rayleigh number maps of the Yarragadee Aquifer, addressing vertical heterogeneity of permeability in the Yarragadee Aquifer. They concluded that convection is likely in the Yarragadee Aquifer and that these maps are a suitable tool for locating free convection in the aquifer. Later- ally continuous permeability facilitates a characterization based on the Rayleigh num- ber (Rayleigh 1916). However, permeability in aquifers is often laterally discontinuous. Performing a standard Rayleigh number analysis in such heterogeneous porous media is difficult and not straightforward (Nield and Simmons 2007; Nield et al. 2009). In the presented work, a different, thermodynamically motivated measure is used for charac - terizing the state of a convection system: the so-called entropy production (Kleidon and Lorenz 2005; Siavashi et al. 2015; Börsing et al. 2017). While the porous Rayleigh num- ber is a sound and tested indicator for characterizing a convective system in an aquifer, it needs the two parameters of the porous medium: L (thickness of the porous medium) and k (the average permeability). These are not constant in real-world reservoir systems, and are usually averaged, i.e., simplified. By assessing the additional heat transport intro - duced by convection, the measure of entropy production bypasses these two simplifica - tions. This can be shown in a numerical experiment with a simple geometry. Figure  1 shows the heat and mass transport in a closed-box model similar to Börsing et al. (2017), with changes made only to the average vertical permeability, either in total or by intro- ducing layers of low permeability. In plot a, the average permeability is high enough to cause instabilities, which are indicated by both, Rayleigh number and entropy produc- tion number (Ra = Rayleigh number and N  = entropy production number in the figure, respectively). Both these numbers decrease, when heterogeneities in permeability (plot b) are introduced. While this is straightforward, it becomes interesting, when the (har- monic) average of (vertical) permeability in a model yields subcritical Rayleigh numbers, as in Fig. 1c, d. These illustrate that models with a similar Rayleigh number may behave completely differently due to a spatially heterogeneous permeability. A subcritical Ray - leigh number, as well as an entropy production number around 1, suggests that heat dif- fusion is the dominant transport process (Fig. 1c). A horizontal, barely permeable layer intersects the aquifer in Fig.  1d. The average vertical permeability in Model c and d is similar, and so are their Rayleigh numbers (1.43 and 1.54, respectively). The Rayleigh number does not suggest convection in model d, while the entropy production num- ber diagnoses convection, as it is greater than 1. Also, a Rayleigh number is a predictive measure of a system, while entropy production is a diagnostic one, thus corresponding more closely to the Nusselt number. In comparison to the Nusselt number, entropy pro- duction has the advantage of being a volumetric estimate. The non-dimensional Nusselt number, defined as the ratio of total heat flow to conductive heat flow over a defined domain, is estimated across a cross-sectional area. Furthermore, similar to the Rayleigh Niederau et al. Geotherm Energy (2019) 7:35 Page 4 of 27 Fig. 1 Four scenarios of temperature in a closed-box model similar to models in Börsing et al. (2017). a Model with homogeneous permeability. b Introduction of a layer with low permeability (black-outlined block) makes permeability heterogeneous. c, d Models with similar harmonic average vertical permeability. In c, the model has a homogeneous permeability field; in d, a barely permeable horizontal layer causes a heterogeneous permeability field. The average vertical permeability of the model, as well as the corresponding Rayleigh number and entropy production number is annotated for each scenario number, it is complicated to assess under non-idealized, i.e., more complex, model domains. Entropy production is mostly independent of the domain geometry as it is cal- culated from the spatial temperature distribution induced by heat transport processes. In this work, we estimate the influence of different correlation lengths on the asso - ciated convection pattern by systematically varying the lateral correlation length of porosity and permeability in a 2D cross section of the Perth Basin model by Niederau et al. (2017). For each correlation length, we generate a Monte Carlo ensemble com- prising 240 ensemble members. To retain comparability, all ensembles have the same random number seeds. This is illustrated in Fig.  2, which shows three exemplary per- meability fields with different horizontal correlation lengths. Due to the same random number seed, areas of low or high permeability are roughly at the same location. The presented manuscript is structured as follows: after reviewing the geological back - ground and introducing the concept of entropy production, we describe the geological model used in this study, as well as the corresponding stratigraphical and petrophysical model. The latter two are connected via the lithostratigraphy of the Yarragadee Aqui - fer described by Timms et al. (2015). We then present the structure of the Monte Carlo ensembles and the setup of the numerical model. The Monte Carlo ensembles are ana - lyzed in the results section and discussed afterwards. Niederau et al. Geotherm Energy (2019) 7:35 Page 5 of 27 Fig. 2 Three examples of permeability fields in the reservoir model used in this study. Permeability fields only differ in the prescribed horizontal correlation length (here, 500 m, 2500 m, and 5000 m). All models have the same random number seed, so the general spatial location of high- and low-permeable areas is preserved. Note that the model is strongly vertically exaggerated Methods Forward code For modeling heat and mass transport in a porous medium, we use the coupled balance equations for energy, mass, and momentum. These conservation equations are imple - mented in the code SHEMAT-Suite (Rath et al. 2006), which is developed from the SHE- MAT code presented in Clauser and Bartels (2003). SHEMAT-Suite models heat and mass transport in porous media in steady as well as in transient state. The transient flow equation used in our simulations is based on Darcy’s law and mass conservation: ∂h ρ g k S =∇ · (∇h + ρ ∇z) + W , (1) s r ∂t μ −3 −2 2 where ρ is fluid density ( kg m ), g gravity ( ms ), k permeability tensor ( m ), μ f f dynamic viscosity of the fluid (Pa  s), h hydraulic head (m), and z the vertical coordi - −1 nate space (m). W describes a volume source term ( s ); ρ is the relative density change ρ −ρ f 0 defined by ρ = , where ρ is density at reference conditions (i.e., atmospheric con- r 0 −1 −1 ditions). S = ρ g(α + φα ) ( m ) is the specific storage coefficient ( m ), ( α and α s f r f r f −1 Pa are compressibilities ( ) of the rock and the fluid phase, respectively). Equation  (1) is derived, using an extended Oberbeck–Boussinesq approximation (Kolditz et al. 1998; Diersch and Kolditz 2002). In our simulations, however, fluid properties are calculated as a function of temperature and pressure after the formulation given in Zyvoloski et al. (1996). The approximation made for Eq. (1) is valid for cases without very high-temper - ature gradients or brines with high salinity (Diersch and Kolditz 2002). These two condi - tions are not present in our simulations, and the approximation is, therefore, suitable for the conditions considered here. The heat transport equation is derived from Fourier’s law q = − ∇T and energy conservation: ∂T (ρ c ) = ∇ · [ ∇T ] + (ρ c ) v∇T + K , (2) p e e p f ∂t −2 ◦ where q is specific heat flow ( Wm ), T is temperature ( C ), and (ρ c ) and (ρ c ) are p e p the effective thermal capacities of the fluid-saturated porous medium and the fluid −3 −1 −1 −1 ( Jm K ), resp e ctively .  is the tensor of bulk thermal conductivity ( Wm K ), and e Niederau et al. Geotherm Energy (2019) 7:35 Page 6 of 27 ρ gk −3 f K is the radiogenic heat production rate ( Wm ); v = (∇h + ρ ∇z) is specific dis - 0 r −1 charge (or Darcy velocity) ( ms ). Parameter values for the bulk porous rock are derived from appropriate mixing laws. The arithmetic mean was used for estimating an average specific heat capacity and the geometric mean for thermal conductivity (Clauser 2011). A more detailed description of the transport equations is provided in Clauser and Bartels (2003). Entropy production The description of the concept of entropy production presented here follows the descrip - tion in Börsing et al. (2017). The second fundamental law of thermodynamics states that the entropy of a closed system is always increasing. The irreversible entropy produc - tion dS by thermal dissipation characterizes irreversible conductive heat transport dQ irr through a (porous) medium with different temperatures at its ends (from higher T to lower temperatures T ) (Tolman and Fine 1948). It can be described as: 1 1 dS = − dQ. irr (3) T T 2 1 Bejan (2013) formulates the production of entropy in a fluid-filled porous medium as the ′′ sum of contributions by thermal dissipation (i.e., conductive heat transport, ) and therm ′′ viscous dissipation (i.e., fluid friction S ): visc ′′ dS ′′ ′′ ′′ total ˙ ˙ ˙ (4) = S = S + S , total therm visc dt ′′ −3 −1 ′′ ′′ where S is in W m K and S the transient derivative. The notation indicates that the calculation is done across a grid face, i.e., between two cell-centered grid nodes. By the definition in Eq.  (3), the transient change of dS in a specified medium can be ′′ described by introducing a specific heat flow q (Tolman and Fine 1948). This vector describes the rate of heat transported per unit time and area, and can be expressed in terms of temperature gradient and thermal conductivity, according to Fourier’s law. The ′′ transient thermal dissipation term ( S in Eq. (4)) can, hence, be reformulated: therm ′′ 1 q  T + T b 1 2 ′′ ′′ 2 S = q ·∇ =− ·∇T = (∇T ) with T = . 0 (5) therm 2 2 T 2 T T 0 0 T and T are the lower and upper temperature boundaries, respectively (Börsing et  al. 2 1 2017; Siavashi et al. 2017), similar to Eq. (3). These temperatures need to be defined for assessing an average initial temperature T . In case of a conductive case (Fig.  3  left), it is reasonable to assume T as the arithmetic mean between upper and lower temper- ature boundaries. As all transient simulations presented in this manuscript start from the same conductive initial conditions, this assumption remains valid.  is the average thermal conductivity over the model domain. In a geothermal system dominated by heat conduction, the predominant heat transport direction is in the direction of the tempera- ture gradient, i.e., normally vertical (Fig.  3 left). If the dominant transport mechanism is advective heat transport by convection, lateral temperature gradients superpose the Niederau et al. Geotherm Energy (2019) 7:35 Page 7 of 27 Fig. 3 Cross sections showing a temperature field of a conductive solution (left) and a convective solution (right). Horizontal temperature profiles (dashed lines) show no lateral temperature gradient for the conductive state, but a significant lateral temperature gradient for the convective state underlying vertical temperature gradient. That is, in addition to a vertical heat trans - port, a lateral (conductive) heat transport occurs, from hotter areas (upwelling regions) towards colder areas (downwelling regions). This lateral heat flow can be described by Eq.  (5), and yields an increase in entropy in the colder areas, while entropy in the hot- ter areas decreases. It should be noted that this observation does not violate the second law of thermodynamics, as the areas described are open. Over the whole model domain, entropy cannot be reduced. In addition to thermal dissipation, entropy is also produced by viscous dissipation ′′ ( S in Eq. (4)), e.g., by internal friction of the fluid. Entropy production by viscous dis - visc sipation can be described as: μ μ f f ′′ 2 S = v + Φ, visc (6) T k −1 v ms k where is the specific discharge in ( ), and Φ the viscous dissipation function. is the average of the reduced permeability tensor ( k ), as we are working with a 2D model: k k xx xz k = . (7) k k zx zz −4 −1 −6 −1 Considering viscosity of water and realistic ranges (10 ms to 10 ms ) of the spe- cific discharge in aquifers, viscous dissipation will have a negligible contribution to the total entropy production. This was also shown by Costa (2006). Thus, in homogeneous porous media, viscous dissipation only has a small effect on the total entropy production in a system. However, if preferred flow paths of a convective system are blocked by zones of low permeability caused by spatially heterogeneous permeability, the importance of viscous dissipation effects can increase. In case of a two-dimensional model like in this work, Eq. (4) can be written as: Niederau et al. Geotherm Energy (2019) 7:35 Page 8 of 27 2 2 ∂T ∂T μ b f ′′ 2 2 dS = + + [v + v ]. (8) x z total T ∂x ∂z kT 0 0 By integrating over the model domain, one obtains the total entropy production S : total ′′ ˙ ˙ S = S dV . total (9) total Non-dimensionalising the total entropy production yields the entropy production Num- ber N (Siavashi et al. 2015), which is used in this study as a representative of the total entropy production in the system. Using a non-dimensionalized version of entropy pro- duction facilitates comparability of our results with the literature: S q total N = with Y = ˙ (10) VY  T , where q is the total specific heat flow (Siavashi et al. 2015). Several thermal rock proper - ties, such as thermal conductivity, are used for calculating the entropy production num- ber. Each grid cell has a bulk thermal conductivity value, which is the geometric mean of thermal conductivities of the rock matrix and pore-filling fluid. The temperature- dependent bulk thermal conductivity value  , used in Eq. (10), is the average of bulk thermal conductivities in the model grid. Temperature-dependent changes of matrix thermal conductivity are small. For water, however, thermal conductivity and viscosity vary strongly with temperature. Water thermal conductivity in the model cells is calcu- lated by the formulation in Zyvoloski et al. (1996), while its viscosity for calculating the entropy production number is calculated based on the IAPWS formulation (Wagner and Kretzschmar 2008). Reservoir model used in this study Different types of models of the Yarragadee Aquifer are presented in this study. Mod - els range from a geometrical geological model, over conceptual stratigraphic and petro- physical models to a numerical model of relevant partial differential equations presented in the methods section. The synthesis of information from all these models is a reservoir model of the Yarragadee Aquifer. In the following, we will briefly describe the geologi - cal model with focus on the main aquifer unit, the stratigraphic model underlying the spatial facies distribution of the main aquifer unit in this model, and the petrophysical model relating porosity and permeability in the aquifer. Finally, we will present the dis- cretized numerical model which combines the features of the previous in a rectilinear hexahedral grid for simulating the heat and mass transport in the Yarragadee Aquifer. Geological model The geological model in this manuscript is a 2D section extracted from a 3D geologi - cal model of the central Perth Basin presented in Niederau et  al. (2017). Its location is indicated in Fig. 4. The Perth Basin is the result of rifting between Australia and Greater India in the late Paleozoic to Mesozoic (Harris 1994) during the breakup of Gondwana. Sedimentary successions filling the Perth Basin can be up to 15 km thick. Thickness of Niederau et al. Geotherm Energy (2019) 7:35 Page 9 of 27 Fig. 4 Map showing well locations and model boundaries in the study area in the central Perth Basin (modified from Niederau et al. (2017)). Location of the 2D cross section is indicated by the dashed line the Yarragadee Aquifer in the central Perth Basin range from 500  m to about 3000  m. More detailed information about the geological and tectonic history of the Perth Basin can be found in Playford et  al. (1976), Backhouse (1984), Harris (1994) and Crostella and Backhouse (2000). Lithological units in the geological model are divided based on their hydraulic properties (see Niederau et al. 2017), i.e., divided into aquifer and aqui- tard units, based on geological and petrophysical data (Crostella and Backhouse 2000; Davidson and Yu 2008; Sheldon et al. 2012; Timms et al. 2015). This is often referred to as “zonation method” (Testerman 1962; Alpay 1972). The main aquifer unit of interest is the jurassic Yarragadee Aquifer, which comprises lithologies of three stratigraphical units, namely Cadda, Yarragadee, and Gage formations. Working groups of the Western Niederau et al. Geotherm Energy (2019) 7:35 Page 10 of 27 Australian Geothermal Center of Excellence (WAGCoE) published several detailed stud- ies covering the petrography of the Yarragadee Aquifer (see Timms et  al. (2015) and references therein), its estimated lateral extension and thickness, and estimated temper- atures in the aquifer (see Sheldon et  al. 2012 and references therein). Our 2D geologi- cal model comprises six geological units, (three aquifer units and three aquitard units), dissected by three major faults. One fault, the Darling Fault, displaces the sedimentary units of the Perth Basin against the Yilgarn Craton, which bounds the Perth Basin to the east. The Darling Fault is the east boundary of the 2D geological model. Two fur - ther faults displace the basement and the Yarragadee Aquifer, visible in the cross sec- tion in Fig.  4. The Yarragadee Aquifer is covered by Post-rift sediments, of which the South-Perth-Shale acts as a thermal blanket for the Yarragadee Aquifer. In Niederau et al. (2017), we showed that reservoir structure has a potentially stronger impact on a convection pattern than spatially heterogeneous permeability. There, we conclude that abrupt changes in reservoir thickness particularly impact a convective pattern. Conse- quently, when studying convection in the Perth Basin, fault displacements should be rep- resented in the model geometry. The two fault displacements considered in this study, are at x ≈ 30, 000 km , and, more pronounced, in the west, at x ≈ 14, 000 km in Fig. 4. However, Niederau et  al. (2017) kept the spatial correlation length of permeability constant; thus, their conclusion relies on a Monte Carlo study with a singular spatial cor- relation length. Next to differences in aquifer thickness due to tectonic activity, spatially heterogene - ous permeability caused by differential sedimentation and compaction was found to affect a convective pattern. In the following, we will present a stratigraphic model of the Yarragadee Aquifer, giving insight in the spatial distribution of lithofacies and cor- responding spatial differences in aquifer permeability. These spatial differences are the heterogeneity, whose spatial continuity is bound to the lithofacies. As their spatial distri- bution and continuity is uncertain, so is that of permeability. Stratigraphical model The Yarragadee Formation was deposited under a predominantly fluvial sedimentary environment, as first inferred by Playford et  al. (1976). Timms et  al. (2015) conducted a detailed study analyzing lithostratigraphy and petrography of multiple core samples from wells in the central Perth Basin (among others the Cockburn 1 well). They identi - fied a variable lithofacies distribution within the Yarragadee Aquifer, with predominant lithofacies corresponding to channel fill deposits and fluvial barforms. Furthermore, they state that formation porosity strongly correlates with lithofacies, which is sup- ported by the previous studies (Delle Piane et al. 2013). Delle Piane et al. (2013) divided core measurements into groups of different flow-zone indicators (FZI). This method is used for identifying classes of rocks which presumably belong to the same hydraulic unit based on porespace geometry and flow properties (Prasad 2003). As lithofacies are cor - related to porosity according to Timms et al. (2015), and FZI deduced from porosity data connect similar hydraulic units (Delle Piane et al. 2013), it may be concluded that litho- facies and hydraulic units correlate. As the depositional environment of the Yarragadee Aquifer varies spatially, types and changes of depositional environment correlate with different lithofacies, and porosity Niederau et al. Geotherm Energy (2019) 7:35 Page 11 of 27 and permeability significantly vary spatially. This may be obvious, considering the flu - vial nature of the Yarragadee Formation, and that hydraulic parameters in a fluviatile sedimentation environment can significantly vary in space. Studies using stratigraphic forward modeling to characterize spatial sediment distribution and distribution of hydraulic properties also show strong spatial variations of porosity and grain size (Cor- bel et  al. 2012a). Vertical variability in porosity can also be determined from wells and geophysical well-log analysis. In Niederau et  al. (2017), we deduced a vertical correla- tion length of around 150 m for porosity, based on porosity measurements on samples from the borehole Cockburn  1. In sedimentary depositional environments, horizontal correlation of lithofacies and connected properties is generally larger than their vertical correlation. Timms et al. (2015) defined facies based on metrics from Gamma ray logs, called GR facies. Those GR facies broadly correspond to lithofacies described in the Yar - ragadee Aquifer. They did not observe clear cross-well correlation of the GR facies and concluded that the corresponding lithofacies are not continuous over distances greater than 5 km. The horizontal correlation length in our study varies between 0.5  and 10 km and is, therefore, significantly larger than the conclusion of Timms et  al. (2015). We extend the horizontal correlation beyond 5 km to assess if larger correlation lengths still affect the occurrence of hydrothermal convection. In the case of the Perth Basin, how - ever, we focus our main analysis on horizontal correlation lengths between 0.5 and 5 km, considering the findings of Timms et al. (2015). Petrophysical model In Niederau et  al. (2017), we deduced porosity from the sonic log of Cockburn  1 after estimating clay content from the Gamma ray log. Using the Wyllie method (Wyllie et al. 1956), porosity was estimated from measured reciprocal interval velocity (“slowness”) Δt: Δt = φΔt + (1 − φ)(V Δt + (1 − V )Δt ) ma (11) f sh sh sh with φ being porosity, V being the shale volume deduced from the GR log, and t , t , sh f sh t the slowness of formation water, shale fraction, and remaining rock matrix, respec- ma tively. Equation (11) was calibrated using available porosity measurements on core samples from Cockburn  1. For relating porosity with permeability, an equation based on fractal theory (Pape et al. 1999) was previously calibrated to the Yarragadee Aquifer, using data from CSIRO (2007) (see Niederau 2014 and Niederau et  al. 2017 for more details): 2.8 8.01 −18 k =[166φ + 19904φ + 388(10φ) ]× 10 , (12) where k is permeability in m . We created a suite of porosity samples (Fig. 5 left) based on mean and standard deviation of the porosity log deduced from Eq. (11). A correlated permeability distribution was generated using Eq.  (12) (Fig.  5 right). We chose to gen- erate a porosity distribution based on mean and standard deviation of the porosity log over using the porosity log directly to have idealized distributions. Interpretation of the results would be significantly more difficult if the underlying distributions of spatial het - erogeneity are not normal, or log-normal distributions, respectively, but asymmetrical or skewed. Niederau et al. Geotherm Energy (2019) 7:35 Page 12 of 27 Fig. 5 Histograms of porosity and permeability samples used for populating the parameter space in the simulation ensembles. Mean and standard deviation of porosity was deduced by log analysis of the borehole Cockburn 1. The permeability distribution was calculated from the porosity distribution by applying Eq. (12) Monte Carlo ensembles Considering the stratigraphic and petrophysical models, we spatially vary porosity and permeability within the Yarragadee Aquifer using Sequential Gaussian Simulation (Deutsch and Journel 1998). This algorithm requires a distribution to draw samples from and a correlation length for defining the spatial variance. The previously described ver - tical correlation length cl of 150  m is kept constant, while we increase the horizontal correlation length cl from 500 to 10,000 m by 500 m steps. For each horizontal correla- tion length, we generate an ensemble comprising 240 members, all with different—but equally likely—porosity and permeability fields, based on the distributions presented in Fig.  5. Ensembles are simulated on the JURECA-HPC (Jülich Supercomputing Cen- tre 2016) cluster and each ensemble comprises up to 240 members to fully utilize the JURECA architecture, i.e., 24 cores per node with 10 nodes allocated for each ensem- ble. Therefore, the particular ensemble size is a result of the architecture of the JURECA cluster, where we ran our simulations on. For comparison, we assess a simulation with homogeneous porosity and permeability (i.e., infinite spatial correlation). Values for the homogeneous model are presented in Table 1. Porosity and permeability of the homoge- neous model resemble mean values of the corresponding distributions. To evaluate the influence of vertical correlation length, we simulate a series of ensembles with a vertical correlation length of 100  m. Again, we increase the horizontal correlation length from 500 to 5000 m by 500 m steps. Numerical model All presented studies on the impact of spatially heterogeneous permeability on the entropy production caused by convection are based on the same model geometry: a two- dimensional cross section taken from the model presented by Niederau et al. (2017). In the following, we describe the essential model parameters in more detail. Niederau et al. Geotherm Energy (2019) 7:35 Page 13 of 27 Fig. 6 Initial temperature (color coded) and boundary conditions for the 2D model of the Perth Basin. Black lines indicate interfaces between different model units Table 1 Geological units implemented in  the  2D model and  corresponding petrophysical properties which were used in the simulation Model unit Porosity Vertical Thermal conductivity Radiogenic −1 −1 permeability (matrix; W m K ) heat production −15 2 −3 (1 0 m ) (μ W m ) Leederville Aquifer 0.3 118.0 3.4 0.6 South Perth Shale 0.1 0.01 1.8 0.82 Yarragadee Aquifer Mean: 0.2 Mean: 11.8 4.3 0.52 Range: 0.14–0.33 Range: 0.3–4100 4.3 0.52 Cattamarra coal measures 0.1 0.018 4.1 0.45 −4 Basement/Yilgarn Craton 0.06 4 4 1.18 10 References 1 2 3 4,5 References: (1) CSIRO (2007); (2) Schilling et al. (2013); (3) HDRPL (2008); (4) Botman (2010); (5) Middleton (2013) Boundary conditions and initial values The lower ( T ) and upper (T ) boundary temperatures of the model domain are kept l u constant (Fig.  6). T equals 95 C , which is shown to be a representative temperature at a depth of 3.5 km in the central Perth Basin (Sheldon et al. 2012; Schilling et al. 2013), and 19 C , the average annual surface temperature around Perth. Lateral boundaries are set as zero Neumann boundary conditions, implying no heat or mass transport across the boundaries. The top boundary for the hydraulic potential was set as a function of topography and sea level respectively, varying around 3500 m, as the onshore part of the model does not have high topography. The initial conditions for temperature are results of a purely conductive steady-state solution, using model parameters stated in Table  1 and displayed in Fig.  6. For the hydraulic head, we applied a constant initial value of 3500 m, representing hydrostatic conditions in the entire model. Spatial discretization The calculated entropy production is affected by the grid discretization of a numerical model. Kleidon and Lorenz (2005) state that the value of entropy production decreases with smaller cell sizes until a threshold value is reached. However, a finer numerical grid, consequently, means more unknowns and longer simulation times, i.e., higher computa- tional cost. To find the optimal discretization with respect to a trade-off between com - putational cost and accuracy, we simulate the 2D model with different discretizations, ranging from 100 to 5  m in cell size. Figure  7 shows the resulting entropy production Niederau et al. Geotherm Energy (2019) 7:35 Page 14 of 27 Fig. 7 Entropy production number N plotted versus cell size for different discretizations of the 2D model shown in Fig. 6. N decreases with smaller cell sizes until it approaches a steady value of around 1.28 number N for different cell sizes. With equidistant sizes of ≤ 25 m , changes in N are ˙ ˙ S S minor compared to the exponential increase in unknowns. Therefore, an equidistant cell size ( dx = dy = dz ) of 25 m is used for our simulations, yielding a total of 196,000 cells. This represents a discretised domain of 35 km in length and 3.5 km in depth (1400 cells in x-direction times 140 cells in z-direction). Transient evolution Each transient simulation in this study starts from the same initial conditions, a purely conductive temperature field. From this conductive state, convection develops over a certain time interval. It is not unusual, that simulations with convection do not reach a steady state, i.e., the numerical solution does not change with time. A reason for this is often oscillation of convection cells (Niederau et  al. 2017), or their superposition with a regional groundwater flow (Schilling et  al. 2013). While a steady-state solution may not be found, simulations with convection can start to show oscillations, e.g., within a temperature interval. Methods exist to assess the timing when such oscillation within a “solution interval” starts. One method is to observe transient convection in a phase space and to see if the solution orbits around the so-called attractors (Lorenz 1963; Vin- cent and Yuen 1988). In Niederau et al. (2016), we showed that the concept of attractors is a viable method to determine the onset of transient oscillation in the Perth Basin. This onset, and the following oscillation, is called “dynamic steady state”. To assess the mini- mal simulation time needed for allowing a convective system in a dynamic steady state to develop, we study the entropy production over time for a small ensemble comprising 12 members with cl equal to 500  m. Generally, the smaller cl is, the longer simula- h h tions need to reach a dynamic steady state. This is why we used cl of 500  m as proxy for assessing the minimum simulation time for the Monte Carlo study. If no significant changes in the entropy production number with time can be observed, a dynamic steady state is reached. Using the entropy production as a measure based on thermodynam- ics, the dynamic steady state can be seen as a thermodynamic steady state. Figure  8 shows the transient evolution of the entropy production number of such an ensemble Niederau et al. Geotherm Energy (2019) 7:35 Page 15 of 27 Fig. 8 Change of entropy production number N with time of a test ensemble comprising 12 members with a horizontal correlation length of 500 m. The average entropy production number stabilizes after around 45,000 years, whereas the majority of ensemble members do not show a strong variation after 60,000 years with a horizontal correlation length of 500 m. The entropy production number increases steeply from the simulation start, suggesting that the whole ensemble is an supercriti- cal state, i.e., above the corresponding critical Rayleigh number for this system. Despite that all realizations show indications of free convection, peak entropy production num- bers differ as well as the transient development until 50,000–60,000 years. This suggests that spatially heterogeneous permeability not only influences whether convection cells develop, but also how long a convective system needs to stabilize. After about 60,000 years, the convective system stabilizes into a (thermo-) dynamic steady state. This behavior for the convective system in the Perth Basin is also described in Schilling et al. (2013) and Niederau et al. (2017). Based on these findings, we limit the simulated time of ensemble members to 60,000 years. Results In this section, we present the simulation results of the aforementioned Monte Carlo ensembles with focus on the influence of changing spatial correlation length, analyzed using the average entropy production number (Eq.  10). Per definition, higher entropy production numbers correspond to more vigorous free convection, i.e., stronger lateral heat transport. Following the hypothesis that increasing the finite horizontal correlation length of permeability coincides with a less likely development of a convective system, we assess the average entropy production number of each ensemble, and relate it to the scenario with homogeneous permeability. Average entropy production with correlation length Figure 9 shows the entropy production number N for ensembles with a horizontal cor- relation length up to 10 km. The average ensemble entropy production number ( N ) decreases with increasing horizontal correlation length ( cl ), suggesting that more ensemble members do not show convection at higher cl . The entropy production of ensembles with heterogeneous permeability is compared to a homogeneous case, where Niederau et al. Geotherm Energy (2019) 7:35 Page 16 of 27 Fig. 9 Variation of the ensemble-average entropy production number N with different horizontal correlation length cl (vertical correlation length cl = 150 m ). Left: including non-convecting ensemble members. Right: h v excluding non-convecting ensemble members aquifer permeability does not vary spatially and is equal to the mean of the permeabil- ity distribution used in the heterogeneous ensembles (Fig.  5 right). Although correla- tion length in a homogeneous case can theoretically be seen as infinite (i.e., cl >> model domain), we plot it at cl equal zero in Fig.  9. The reason is that infinite and finite cor - relation length are different. Both finite correlation lengths cl and cl would have to be h v magnitudes larger than the model extent to mimic a homogeneous case, which is practi- cally the same as not stating a finite correlation length at all. In that case, only the per - meability value, and not also its spatial change, would be a controlling factor whether a convection system forms, similar to a model with homogeneous permeability. However, due to the sequential nature of the SGSim algorithm, the first determined permeability value would be used for the whole domain in that case, and likely not represent the aver- age permeability of the reservoir. In a case with observable heterogeneous permeability, both, the permeability value, and its change in space, are important to consider, making a homogeneous approach not feasible. This means that in Fig.  9 and following visualiza- tions, a finite correlation length cl of zero means that no finite correlation length exists, i.e., a homogeneous case. The homogeneous case, i.e., no finite correlation length, is characterized by a higher entropy production number, associated with an ideal convection system, where cells are not distorted or diverted by lenses with low permeability. The entropy production of a homogeneous case is used as a reference point for emphasizing the importance of considering heterogeneous permeability in an aquifer (here specifically the Yarragadee Aquifer) when modeling hydrothermal convection. While the entropy production of the whole ensemble continuously decreases with increasing cl (Fig.  9 left), it stabilizes at around 2000 m horizontal correlation length if members without convection are ignored when calculating N (Fig.  9 right). No significant influence of cl on the aver- ˙ h age entropy production of an ensemble is observed for correlation lengths larger 4000 m (Fig.  9 right). This suggests that heat transport is not strongly influenced by changes in cl beyond this range. This can further mean that the number of convection cells is rela - tively stable for cl larger 4000 m. We will test this hypothesis later in this manuscript by assessing the relation between number of convection cells and entropy production. Niederau et al. Geotherm Energy (2019) 7:35 Page 17 of 27 Fig. 10 Boxplots of the entropy production number for each ensemble with different horizontal correlation length cl . At the beginning, average entropy production decreases with increasing cl , but stabilizes at larger h h cl Figure 10 shows boxplots for each ensemble, including only members showing convec- tion. While the mean entropy production decreases until cl = 2000 m , ensemble spread increases with increasing cl as indicated by the longer whiskers. As said before, the ensemble average of the entropy production number is constant for larger cl . Howe ver , analysis of all the ensembles shows oscillation of the ensemble spread for larger hori- zontal correlation lengths. Outliers (i.e., above 1.5 times the upper quartile or below 1.5 times the lower quartile) are mostly found at the lower end of the ensembles (diamonds in Fig. 10). This may be linked directly to the top and bottom Dirichlet boundary condi - tions for temperature, meaning that there is a fixed upper limit for heat transport in the model, thus reducing the likelihood of “positive” outliers. Simulation results suggest that the increase in horizontal correlation length does not have a strong influence once a convection system has formed. We do not observe a systematic connection between cl and ensemble variance, which is expressed by the varying whisker lengths for boxplots with increasing cl in Fig.  10. It should be noted, however, that variations in whisker lengths may also be due to the ensemble size, so con- clusions should be made cautiously. The ensemble median, however, is relatively con - stant over this range of cl , indicating no significant shift of the ensemble distribution, but a change in its width, i.e., its standard deviation. Influence of vertical correlation length In a sedimentary environment, lithofacies and associated petrophysical parameters do not only vary laterally, but also vertically. Such vertical variations due to layering have a stronger influence on the resulting vertical heat transport than lateral variations. This can be illustrated by comparison with a simplified circuit model. In case of variations along the transport path, i.e., vertical variations when looking at vertical heat transport, resistors in the electric circuit would be connected in series. That is, adding a resistor with a high resistivity would already change the resulting voltage significantly. For lat - eral variations, the analog circuit would be resistors connected in parallel, and the total Niederau et al. Geotherm Energy (2019) 7:35 Page 18 of 27 Fig. 11 Schematic showing the significant influence of vertical variations along the transport path compared to lateral variations, and analog electrical circuits as examples. Arrows represent the relative proportions of q passing through the material, not its magnitude resistance is lower than the smallest individual resistor in the system. Here, an additional resistor in series with a high resistivity would have a less significant influence. Figure  11 illustrates this analog. As a result, one would use the harmonic mean for averaging vertical variations along the transport path, and the arithmetic mean for averaging the lateral variations. Accordingly, changing the vertical correlation length of permeability should have a more significant influence on the associated convection system than changing the horizontal correlation length. To test this hypothesis, we performed the same Monte Carlo simulations with cl equal to 100  m. Results show that the vertical correlation length ( cl ) of permeability has a significant effect on the likelihood of convective heat transport. Again, increasing horizontal correlation length of permeability correlates with increasing number of realizations showing no convective heat transport for both vertical correlation lengths (Fig.  12). While over 80% of members in ensembles with cl = 150 m show convection, this fraction decreases drastically for ensembles with cl = 100 m down to < 50 % . For both vertical correlation lengths, the number of con- vecting ensemble members mainly decreases between 1500 and 3500 m horizontal correlation length. In ensembles with cl larger than 3500 m, the fraction of members showing convection stabilizes for cl = 100 m , while it continuously slightly decreases for cl = 150 m. Fig. 12 Variation of the fraction of realizations with free convection N with horizontal correlation convection length cl for two different cl . Ensembles with cl = 100 m show a significantly stronger decrease in h v v convecting realizations, compared to ensembles with cl = 150 m v Niederau et al. Geotherm Energy (2019) 7:35 Page 19 of 27 These results emphasize how important knowledge about the (vertical) spatial continuity of permeability in a modeled aquifer is for quantifying the probability of convective heat transport. Furthermore, we observe that cl has a stronger influence on the development of convection cells than cl , which does not influence the heat transfer, expressed by entropy production, over a certain range. This observation has important implications, as the vertical correlation length cl can usually be deduced from well logs. A shorter vertical correlation length means more frequent alterna- tion of different permeabilities, meaning more little permeable zones, which strongly influence the average permeability of the sedimentary column, as it is averaged by the harmonic mean. Number of convection cells and entropy production In the context of this work, one convection cell is understood as the connected pair of zones of upflow and downflow. To test the hypothesis that spatially heterogeneous permeability influences the convection pattern, i.e., the number of convection cells, we count the convection cells in each simulated model. For this purpose, we apply a peak count method to horizontal temperature profiles through the model. Each counted posi - tive peak can be associated with a zone of upflow. Before presenting the results of the cross section of the Perth Basin, we show, as an example, the influence of cl on the number of convective rolls in closed boxes of dimensions 2 km (x-direction) ×1 km (z-direction) in Fig.  13. In these simplified models without complicated reservoir structure, the change in the number of convec- tion cells is solely controlled by the change in horizontal correlation length of perme- ability. Figure 13 shows isotherms for three different cl . At cl = 50 m , three distinct h h convection cells develop. Increasing cl to 250 m and 500 m yields a reduction in the number of convection cells caused by wider zones of low permeability (gray-to-white fields in Fig.  13). Results of this analysis for the Perth Basin ensembles are presented as a joint plot of correlation length and number of convection cells (wave number) in Fig.  14. A joint plot shows the bivariate distribution between two variables, here cl and wave number, as well as their univariate distributions at the margin. The number of convection cells decreases at shorter horizontal correlation lengths, visualized by the darker colors in Fig. 14. Darker blue indicates that a wave number was counted in Fig. 13 Simulation of convection in a closed box with heterogeneous permeability. Horizontal correlation length cl of permeability from left to right: 50 m, 250 m, and 500 m. Permeability is displayed by gray fields, while isotherms are presented as colored lines Niederau et al. Geotherm Energy (2019) 7:35 Page 20 of 27 Fig. 14 Heatmap of wave number (number of convection cells) and horizontal correlation length for the ensembles. Color intensity represents the amount of ensemble members showing the according number of convection cells with lighter colors correlating to more models in a specific bin more ensemble members. At cl = 500 m , the majority of simulations have eight con- vection cells. At cl = 1000 m , the average wave number decreased to seven, and to six at cl = 4000 m . This observation agrees with results from simpler models, where no major changes in the number of convection rolls are observed for larger cl . Simil arly to the entropy production number (Fig. 10), the number of convection cells stabilizes at longer correlation lengths. While 12 peaks are counted in the model with homo- geneous permeability, the number of peaks undulates between 6 and 7 for ensembles with a correlation length greater than 1500 m. That is, at 1500 m, the average number of convection cells in the model is halved compared to a simulation with homogene- ous permeability. The similarity between the distribution of the entropy production number and the wave number suggests that the former can be related to the probable number of convection cells, given a closed model with a fixed size. Changes in the model size, i.e., its aspect ratio, also influence the number of convection cells (Beck 1972; Bejan 1979; Turcotte and Schubert 2014; Bejan 2013). Discussion The presented results are analyzed following the hypothesis that convective systems are less likely to develop when we increase the finite horizontal correlation length in an aquifer, and that this can be measured by the entropy production number. This hypothesis may need further explanation. We assume that all parameters relevant for convection, such as temperature difference across the porous medium, or its thick - ness, are kept constant. A model with homogeneous permeability can be treated as a case with an infinite spatial correlation length. In such a model, only the magni - tude of permeability determines if convection is likely or not. If we consider the same model with heterogeneous permeability, i.e., a finite correlation length, the magni - tude of permeability and its spatial variation would control the onset of free convec- tion. It may appear logical that longer correlation lengths might facilitate the onset of Niederau et al. Geotherm Energy (2019) 7:35 Page 21 of 27 convection, as permeable areas stretch over a longer distance. However, so do areas of low permeability. These work as barriers to the buoyancy-driven convective cur - rent, which is mainly vertical. At shorter finite horizontal correlation lengths, these currents just get diverted around lenses of low permeability, whereas a lateral elonga- tion of those lenses can completely inhibit the onset of convection. And such a lateral elongation corresponds to a longer horizontal correlation length. One controlling factor in a convective system is the temperature difference between top and bottom of the convecting layer. Consequently, the type of boundary conditions applied to a model influences the simulated convection. We apply Dirichlet boundary conditions at the top and bottom boundaries of the model. Thus, the applied boundary conditions supply a constant amount of energy to the system. It needs to be discussed whether these boundary conditions reflect a realistic state of the system. If tempera - tures at depth are known, applying Dirichlet boundary conditions is reasonable. How- ever, as this is usually not the case, thermal Neumann boundary conditions are applied as bottom boundary condition (Smith and Chapman 1983). In the Perth Basin region, several studies use a basal Dirichlet boundary condition for temperature and yield real- istic simulation results (Corbel et  al. 2012b; Schilling et  al. 2013; Niederau et  al. 2017). Their results suggest that, for the bottom boundary, a Dirichlet boundary condition for temperature is suitable for the Perth Basin. Still, it is worthy to discuss the implications which the choice of boundary condition has on the entropy production. Constant force boundary conditions (Dirichlet boundary conditions) yield a maximum of estimated entropy production and constant flux boundary conditions (Neumann boundary condi - tions) define a lower bound. We focused here on maximum conditions to obtain a com - parison to the previous work on cell structures (Börsing et al. 2017). Furthermore, under the scope of this paper, i.e., the possibility to use entropy production to obtain insight in the behavior of convection in heterogeneous systems, we generally expect similar results concerning the geometric interpretation of the entropy production. Additionally to permeability of a reservoir, its structure, i.e., abrupt changes in reser- voir thickness, has a strong influence on the convective system. Such changes are usually caused by fault zones which not only influence reservoir structure, but normally also the spatial permeability distribution, and can thus influence groundwater flow paths (Frey - mark et al. 2019; Inbar et al. 2019). In terms of their permeability relative to that of the surrounding rock, fault zones may be more permeable, similarly permeable, or less per- meable. For instance, fault-related processes, such as clay smear (Vrolijk et al. 2016), can drastically influence permeability of a fault, e.g., by acting as a barrier to fluid flow. The existence of clay smear alone, however, is not a defining indicator on the permeability of the fault, as potential breaching of the clay layer coincides with increased flow across the fault (Kettermann et  al. 2017). Furthermore, experimental data showed that the influ - ence of clay smear on permeability changes with fault activity, i.e., rock deformation (Takahashi 2003). We only consider fault displacement in our simulations which over- simplifies the impact of fault systems. Schilling et al. (2013) observe that zones of upflow occur in fault zones or in their vicinity, depending if the faults are highly permeable or impermeable respectively. While it is logical that upflow is focused within permeable fault zone, it can appear counterintuitive that zones of upflow occur next to fault zones if these are impermeable. However, aquifer thickness of the Yarragadee Aquifer is highest Niederau et al. Geotherm Energy (2019) 7:35 Page 22 of 27 in the hanging walls of normal faults. Thus, formation of convection cells in these areas is facilitated, especially when fluid flow across a fault is inhibited by its low permeability. However, fault zones are not homogeneously permeable over their whole extent, mean- ing that the results presented by Schilling et  al. (2013) reflect extreme conditions. But even if faults are modeled just as displacement, they have an influence on the convective pattern (Niederau et al. 2017). It is not known whether faults in the central Perth Basin are sealing or transmissive. Faults studied in the northern Perth Basin suggest that they may be sealing (Olierook et al. 2014), but this cannot be assumed to be true in the cen- tral Perth Basin. Considering the great uncertainty in fault permeability in general, and also in the Perth Basin, we decided to model only the fault displacement, as our study focuses on the influence of the correlation length of matrix permeability on the convec - tive system, and not on the influence of fault permeability. It would, however, be possible and important, to extend the analysis to simulations which consider fluid flow in and across fault zones, for example by representing faults as 2D elements in a finite-element discretization. The presented diagnostic equations for the entropy production would be equally applicable. Our simulation results suggest that varying the correlation length of permeability in the Yarragadee Aquifer changes the convection pattern in the aquifer. These results agree with the theoretical studies on the influence of the spatial continuity on the onset of convection (Nield et al. 2010; Simmons et al. 2010). For instance, Nield and Simmons (2007) state that weak vertical heterogeneity has presumably little or no effect on the onset of convection. In our simulations, reducing the vertical correlation length by 50 m caused a significant reduction of ensembles showing convection at higher horizon - tal correlation lengths. By reducing the vertical correlation length, we increase the ver- tical heterogeneity. However, the variance of the permeability distribution applied in our model (see Fig.  5 right) is too large to be considered weak heterogeneity after the definition provided by Nield and Simmons (2007) and Simmons et al. (2010). Thus, the influence of vertical heterogeneity is more significant in our study, but this does not con - tradict previous studies. In the case of the Perth Basin, it can further be concluded that a horizontal correla- tion length greater than 2000 m has a small influence on the convective pattern. Such a stabilization in convective systems was also found by Prasad and Simmons (2003) who observed a decrease of overall degree of instability with increasing horizontal correlation length. While the model which they studied is a heterogeneous Elder Problem (Elder 1967), their conclusions can be confirmed to some extent in this study with a geometri - cally more complex model. Furthermore, our findings coincide with those by Prasad and Simmons (2003) that reservoir models using representative averages or flow parameters will often significantly underestimate the importance of their spatial heterogeneities, e.g., their effect on the onset of convection. We use a geostatistical approach for varying horizontal correlation lengths, as we have no real information about the latter in the Yarragadee Aquifer. Cross correlation between wells is a usual approach for assessing a horizontal correlation length. How- ever, the distance between wells in the study area is between 40 and 100 km, thus too long to assess a meaningful horizontal correlation length. The vertical correlation length, however, as well as the value range of porosity and permeability used in our models are Niederau et al. Geotherm Energy (2019) 7:35 Page 23 of 27 representative for the Yarragadee Aquifer, as it is based on core- and log data from wells in the Perth Basin (Niederau et al. 2017). Another approach for modeling the spatial per- meability distribution would be stratigraphic forward modeling. Existing stratigraphic models for the Yarragadee Aquifer by Corbel et al. (2012a) provide a good basis for esti- mating a range of correlation lengths and a different approach for modeling a spatial permeability distribution. In their model, sediment influx in the Perth Basin was from the north and east. Despite yielding a probable permeability distribution, stratigraphic forward modeling in such a sedimentary environment is non-trivial and model results are just one possible realization of the sediment distribution in the subsurface. As such, it is comparable to geostatistical realizations of the same distribution. The key difference between the two approaches is that stratigraphic forward modeling is a process-based approach, while geostatistics is a data-based one. In contrast to analyzing one strati- graphic forward model, a Monte Carlo approach enables us to draw statistically more robust conclusions. Ideally, both approaches should be combined. That is, conditioned geostatistical ensembles of permeability distributions based on stratigraphic forward models could provide a model representative for the Yarragadee Aquifer. A classic approach in studying convecting system is assessing the systems’ Rayleigh number. While the Rayleigh number is very useful for estimating if a system is likely to show convection, it does not provide information about the likely convective pattern. The same Rayleigh number may characterize different convection patterns which have different entropy production numbers (Börsing et al. 2017). Also, assessing a meaning - ful Rayleigh number in a non-idealized setting (e.g., a convecting layer with non-parallel boundaries) is difficult, as the height of the convecting layer, or the temperature differ - ence across it is two input parameters in classical Rayleigh number formulations. Both parameters are not constant in a real-world example. For instance, estimating the Ray- leigh number of a simple model similar to the one shown in Fig.  1a is straightforward, as is the estimation of N . Introducing lenses of low permeability as in Fig.  1b changes the convective system. Such a change in the average permeability is represented by a change in the Rayleigh number and N of the model. If introduction of heterogeneity in permeability is not reflected by a change in the average permeability, a correct Ray - leigh number analysis gets problematic, especially if the harmonic mean of the vertical permeability of the studied aquifer is subcritical, as in Fig. 1c, d. Such a conclusion, that traditional Rayleigh analysis can be misleading, as convection can occur at subcritical Rayleigh numbers (due to subcritical average permeability), was also drawn by Prasad and Simmons (2003), who stated “the traditional Rayleigh stability criterion based upon an average permeability is inadequate for describing solute transport processes in het- erogeneous systems”. A traditional Rayleigh number analysis may only be valid when heterogeneity in permeability is small, i.e., a deduced log-permeability distribution has a small standard deviation. In a system with spatially heterogeneous permeability, a Rayleigh number analysis may lead to incorrect predictions, as low-permeable layers in a system significantly reduce the average vertical permeability (e.g., when using the harmonic mean) and thereby the estimated Rayleigh number. Thus, entropy production can be used as a thermodynam - ically-based, diagnostic metric which may vary in systems of equal Rayleigh number and thus represent different convection patterns. This is one advantage of the entropy Niederau et al. Geotherm Energy (2019) 7:35 Page 24 of 27 production number compared to a Rayleigh number assessment. Concerning the Nus- selt number, entropy production has the advantage that a volumetric estimation of con- vection is possible, due to the formulation by Bejan (2013). Furthermore, the entropy production number can be used for assessing a probable convection pattern in a res- ervoir of more complex geometry. In our simulations, the number of convection cells decreased with increasing horizontal correlation length of permeability, at least at lower correlation lengths. The average entropy production number similarly decreased at lower correlation lengths, suggesting that it correlates with the number of convec- tion cells. It should, however, be noted that entropy production provides the potential to observe whether major changes in a convective system occur. It does not directly help in locating zones of up flow or down flow in a realistic environment. In this sense, we can - not expect the measure of entropy production to indicate regions of upflow in an area. However, it provides an estimation about the expected number, and accuracy, of up flow and down flow regions in the entire simulation domain. Conclusion In this manuscript, we analyze how different correlation lengths in permeability affect a convection pattern by using a Monte Carlo approach. We use the measure of entropy production for quantifying this effect. Simulation results of the MC ensembles show that the number of ensemble members with dominant convective heat transport decreases with increasing correlation length, as expressed by a continuously decreasing entropy production number ( N ). Excluding non-convective ensemble members reveals that N ˙ ˙ S S stabilizes at a horizontal correlation length of around 2000 m. As expected, the num- ber of convection cells also decreases with increasing correlation length, similarly to the entropy production number, also stabilizing at a horizontal correlation length of around 2000 m. Vertical correlation length has a stronger influence on the possibility of convec - tion development in our model. This result has important implications for similar case studies, as deducing the vertical correlation length from well logs is usually feasible in reservoir studies. The measure presents itself as a powerful diagnostic tool when studying hydrothermal convection in real-world reservoir systems, and that in such settings, it can be more use- ful than a classical Rayleigh number assessment. It is not fully clear which conditions in our simulations facilitate the development of a convective system. It is likely that certain combinations of heterogeneous permeability and reservoir structure are an important factor for convection in our model of the Perth Basin. The evaluation of entropy produc - tion shows further that it is important to consider spatially heterogeneous permeability with an accurate correlation length to reliably simulate the probable convection present in the Perth Basin. We consider the use of entropy production as a diagnostic tool and focus here on 2-D analyses, as an extension and application to previous theoretical work (Börsing et  al. 2017). However, one strength of the measure itself is that it can be applied to any more complex setting,  both in the choice of boundary conditions, as well as in the dimen- sionality of the problem. Including faults as own hydraulic units in a finite-element dis - cretization, for instance, the influence of a constant flux boundary condition, or how entropy production is affected by the development of elongated convection rolls, rather Niederau et al. Geotherm Energy (2019) 7:35 Page 25 of 27 than cells, in three-dimensional problems. We already performed the first studies inves - tigating the development entropy production and its interpretability in three-dimen- sional settings (Huang 2017) and see this as an interesting path for future research. Acknowledgements The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium on the supercom- puter JURECA at Forschungszentrum Jülich for the JARA-HPC project 11061. The authors thank the reviewers and editor for their valuable comments which helped improving the manuscript. Authors’ contributions JN sets up the geological model, did the ensemble simulation and majority of analysis. He wrote the majority of the manuscript. JFW was a major contributor in structuring and revising the manuscript. NB significantly contributed to the manuscript, and provided the basis code for the entropy production analysis. All authors read and approved the final manuscript. Funding The conducted research was not part of a funded project. Computing time on the supercomputer JURECA were Granted by JARA-HPC under the project number 11061. Availability of data and materials The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request. Competing interests The authors declare that they have no competing interests. Author details Institute for Applied Geophysics and Geothermal Energy, E.ON Energy Research Center, RWTH Aachen University, 2 3 Aachen, Germany. Jülich-Aachen Research Alliance (JARA-HPC), RWTH Aachen University, Aachen, Germany. Present Address: Geothermal Energy and Geofluids Group, Department of Earth Sciences, ETH Zürich, Zurich, Switzerland. Com- putational Geoscience and Reservoir Engineering, RWTH Aachen University, Aachen, Germany. Institute of Geophysics, Department of Earth Sciences, ETH Zürich, Zurich, Switzerland. Received: 2 August 2019 Accepted: 8 November 2019 References Alpay OA, et al. A practical approach to defining reservoir heterogeneity. J Petrol Technol. 1972;24(07):841–8. Backhouse J. Revised late Jurassic and early Cretaceous stratigraphy in the Perth Basin. Geol Surv Western Aust Rep. 1984;12:1–6. Barbier E. Geothermal energy technology and current status: an overview. Renew Sustain Energy Rev. 2002;6(1):3–65. https ://doi.org/10.1016/S1364 -0321(02)00002 -3. Beck JL. Convection in a box of porous material saturated with fluid. Phys Fluids. 1972;15(8):1377–83. Bejan A. A study of entropy generation in fundamental convective heat transfer. J Heat Transf. 1979;101(4):718–25. Bejan A. Convection heat transfer. Hoboken: Wiley; 2013. Börsing N, Wellmann JF, Niederau J, Regenauer-Lieb K. Entropy production in a box: analysis of instabilities in confined hydrothermal systems. Water Resour Res. 2017;53(9):7716–39. https ://doi.org/10.1002/2017W R0204 27. Botman C. Determining the geothermal properties of the Perth Basin lithology to aid geothermal energy projects. Perth: Curtin University of Technology B.Sc (Geophysics) Honours thesis Unpublished. 2010. Clauser C. Thermal storage and transport properties of rocks, II: thermal conductivity and diffusivity. Dordrecht: Springer; 2011. p. 1431–48. https ://doi.org/10.1007/978-90-481-8702-7_67. Clauser C, Bartels J. Numerical simulation of reactive flow in hot aquifers: SHEMAT and processing SHEMAT. Berlin: Springer; 2003. Corbel S, Griffiths C, Dyt C, Ricard L. Hot sedimentary aquifer characterization using stratigraphic forward modelling, Perth Basin, Australia. In: Proceedings of the thirty-seventh workshop on geothermal reservoir engineering, Stanford University, Stanford, California. 2012a. Corbel S, Schilling O, Horowitz F, Reid L, Sheldon H, Timms N, Wilkes P. Identification and geothermal influence of faults in the Perth metropolitan area, Australia. In: Thirty-seventh workshop on geothermal reservoir engineering, Stanford, CA. 2012b. Costa V. On natural convection in enclosures filled with fluid-saturated porous media including viscous dissipation. Int J Heat Mass Transf. 2006;49(13–14):2215–26. Crostella A, Backhouse J. Geology and petroleum exploration of the central and southern Perth Basin, Western Australia. 57, Geological Survey of Western Australia. 2000. CSIRO Csiro Pressureplot(tm) software and Pressuredb (tm) database. 2007. www.press urepl ot.com. Accessed 08 Apr Davidson WA, Yu X. Perth regional aquifer modelling system (PRAMS) model development: Hydrogeology and ground- water modelling. Western Australia Department of Water. 2008. Delle Piane C, Esteban L, Timms L, Ramesh Israni S. Physical properties of mesozoic sedimentary rocks from the Perth Basin, Western Australia. Aust J Earth Sci. 2013;60(6–7):735–45. Niederau et al. Geotherm Energy (2019) 7:35 Page 26 of 27 Deutsch C, Journel A. GSLIB Geostatistical software library and user’s guide. New York: Oxford University Press; 1998. Diersch HJG, Kolditz O. Variable-density flow and transport in porous media: approaches and challenges. Adv Water Resour. 2002;25(8):899–944. https ://doi.org/10.1016/S0309 -1708(02)00063 -5. Elder J. Transient convection in a porous medium. J Fluid Mech. 1967;27(3):609–23. Freymark J, Bott J, Cacace M, Ziegler M, Scheck-Wenderoth M. Influence of the main border faults on the 3d hydraulic field of the central Upper Rhine Graben. Geofluids. 2019;. https ://doi.org/10.1155/2019/75207 14. Grant M, Donaldson I, Bixley P. Geothermal reservoir engineering. New York: Academic Press; 1982. https ://doi. org/10.1016/B978-0-122-95620 -1.X5001 -2. Harris L. Structural and tectonic synthesis for the Perth Basin, Western Australia. J Petrol Geol. 1994;17(2):129–56. HDRPL Geothermal energy potential in selected areas of Western Australia (Perth Basin). Hot Dry Rocks Pty Ltd. (HDRPL), South Yarra, Australia. 2008. Horowitz F, Regenauer-Lieb K, Wellmann J, Chua H, Wang X, Poulet T. Evidence for hydrothermal convection in the Perth Basin, Australia. In: Gurgenci H, Budd AR, editors. Proceedings of the Sir Mark Oliphant international frontiers of sci- ence and technology Australian geothermal energy conference. 2008;18. Huang P. Investigating different formulations for hydrothermal convection in geothermal systems. Master’s thesis, RWTH Aachen University. 2017. Inbar N, Rosenthal E, Magri F, Alraggad M, Möller P, Flexer A, Guttman J, Siebert C. Faulting patterns in the lower yarmouk gorge potentially influence groundwater flow paths. Hydrol Earth Syst Sci. 2019;23(2):763–71. Irvine DJ, Sheldon HA, Simmons CT, Werner AD, Griffiths CM. Investigating the influence of aquifer heterogeneity on the potential for thermal free convection in the Yarragadee Aquifer, Western Australia. Hydrogeol J. 2014;. https ://doi. org/10.1007/s1004 0-014-1194-1. Jureca Jülich Supercomputing Centre. General-purpose supercomputer at jülich supercomputing centre. J Large-Scale Res Facil. 2016;2:A62. https ://doi.org/10.17815 /jlsrf -2-121. Kettermann M, Urai JL, Vrolijk PJ. Evolution of structure and permeability of normal faults with clay smear: insights from water-saturated sandbox models and numerical simulations. J Geophys Res Solid Earth. 2017;122(3):1697–725. Kleidon A, Lorenz R. Entropy production by earth system processes. In: Kleidon A, Lorenz R, editors. Non-equilibrium thermodynamics and the production of entropy. Heidelberg: Springer; 2005. p. 1–20. Kolditz O, Ratke R, Diersch HJG, Zielke W. Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models. Adv Water Resour. 1998;21(1):27–46. Lipsey L, Pluymaekers M, Goldberg T, van Oversteeg K, Ghazaryan L, Cloetingh S, van Wees JD. Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands. Geothermics. 2016;64:135–51. Lorenz EN. Deterministic nonperiodic flow. J Atmos Sci. 1963;20(2):130–41. Middleton M. Radiogenic heat generation in the Darling Range, Western Australia. Explor Geophys. 2013;44:206–14. Niederau J. Calibration of a fractal model relating porosity to permeability and its use for modeling hydrothermal transport processes in the Perth Basin, Australia. Energy Procedia. 2014;59:293–300. https ://doi.org/10.1016/j.egypr o.2014.10.380. Niederau J, Börsing N, Wellmann F, Clauser C. Entropy production and attractors: Measures to quantify uncertainty and complexity introduced by convection. In: Proceedings European geothermal congress, 19–24 September, Stras- bourg, France. 2016. Niederau J, Ebigbo A, Marquart G, Arnold J, Clauser C. On the impact of spatially heterogenous permeability on free convection in the Perth Basin, Australia. Geothermics. 2017;66:119–33. Nield D, Simmons CT. A discussion on the effect of heterogeneity on the onset of convection in a porous medium. Transp Porous Media. 2007;68(3):413–21. Nield D, Kuznetsov A, Simmons CT. The effect of strong heterogeneity on the onset of convection in a porous medium. Transp Porous Media. 2009;77(2):169. Nield D, Kuznetsov A, Simmons CT. The effect of strong heterogeneity on the onset of convection in a porous medium: 2d/3d localization and spatially correlated random permeability fields. Transp Porous Media. 2010;83(3):465–77. Olierook HK, Timms NE, Hamilton PJ. Mechanisms for permeability modification in the damage zone of a normal fault, northern Perth Basin, Western Australia. Marine Petrol Geol. 2014;50:130–47. https ://doi.org/10.1016/j.marpe tgeo.2013.10.012. Pape H, Clauser C, Iffland J. Permeability prdiction based on fractal pore-space geometry. Geophysics. 1999;64(5):1447–60. Playford PE, Cockbain AE, Low G. Geology of the Perth Basin, Western Australia. Geological Survey of Western Australia, Perth, Australia. 1976. Prasad M. Velocity-permeability relations within hydraulic units. Geophysics. 2003;68(1):108–17. Prasad A, Simmons CT. Unstable density-driven flow in heterogeneous porous media: a stochastic study of the elder “short heater” problem. Water Resour Res. 2003;39(1):SBH 4–1–SBH 4–21. https ://doi.org/10.1029/2002W R0012 90. Pujol M, Ricard LP, Bolton G. 20 years of exploitation of the Yarragadee Aquifer in the Perth Basin of Western Australia for direct-use of geothermal heat. Geothermics. 2015;57:39–55. https ://doi.org/10.1016/j.geoth ermic s.2015.05.004. Rana R, Horne R, Cheng P. Natural convection in a multi-layered geothermal reservoir. J Heat Transf. 1979;101(3):411–6. Rath V, Wolf A, Bücker H. Joint three-dimensional inversion of coupled groundwater flow and heat transfer based on automatic differentiation: sensitivity calculation, verification, and synthetic examples. Geophys J Int. 2006;167(1):453–66. Rayleigh L. On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Lond Edinburgh Dublin Philos Mag J Sci. 1916;32(192):529–46. Rubin Y, Mavko G, Harris J. Mapping permeability in heterogeneous aquifers using hydrologic and seismic data. Water Resour Res. 1992;28(7):1809–16. Schilling O, Sheldon H, Reid L, Corbel S. Hydrothermal models of the perth metropolitan area, Western Australia: implica- tions for geothermal energy. Hydrogeol J. 2013;21(3):605–21. https ://doi.org/10.1007/s1004 0-012-0945-0. Sheldon HA, Florio B, Trefry MG, Reid LB, Ricard LP, Ghori KAR. The potential for convection and implications for geother- mal energy in the Perth Basin, Western Australia. Hydrogeol J. 2012;20(7):1251–68. Niederau et al. Geotherm Energy (2019) 7:35 Page 27 of 27 Siavashi M, Bahrami HRT, Saffari H. Numerical investigation of flow characteristics, heat transfer and entropy generation of nanofluid flow inside an annular pipe partially or completely filled with porous media using two-phase mixture model. Energy. 2015;93:2451–66. Siavashi M, Bordbar V, Rahnama P. Heat transfer and entropy generation study of non-darcy double-diffusive natural convection in inclined porous enclosures with different source configurations. Appl Therm Eng. 2017;110:1462–75. https ://doi.org/10.1016/j.applt herma leng.2016.09.060. Simmons CT, Kuznetsov AV, Nield DA. Eec ff t of strong heterogeneity on the onset of convection in a porous medium: importance of spatial dimensionality and geologic controls. Water Resour Res. 2010;. https ://doi.org/10.1029/2009W R0086 06. Smith L, Chapman DS. On the thermal effects of groundwater flow: 1. Regional scale systems. J Geophys Res Solid Earth (1978–2012). 1983;88(B1):593–608. Takahashi M. Permeability change during experimental fault smearing. J Geophys Res Solid Earth. 2003. https ://doi. org/10.1029/2002J B0019 84. Testerman J, et al. A statistical reservoir-zonation technique. J Petrol Technol. 1962;14(08):889–93. Timms NE, Olierook HK, Wilson ME, Delle Piane C, Hamilton PJ, Cope P, Stütenbecker L. Sedimentary facies analysis, mineralogy and diagenesis of the mesozoic aquifers of the central Perth Basin, Western Australia. Marine Petrol Geol. 2015;60:54–78. Tolman RC, Fine PC. On the irreversible production of entropy. Rev Mod Phys. 1948;20(1):51. Turcotte D, Schubert G. Geodynamics. Cambridge: Cambridge University Press; 2014. Vincent AP, Yuen DA. Thermal attractor in chaotic convection with high-prandtl-number fluids. Phys Rev A. 1988;38:328– 34. https ://doi.org/10.1103/PhysR evA.38.328. Vrolijk PJ, Urai JL, Kettermann M. Clay smear: review of mechanisms and applications. J Struct Geol. 2016;86:95–152. Wagner W, Kretzschmar HJ. Iapws industrial formulation 1997 for the thermodynamic properties of water and steam. In: International steam tables: properties of water and steam based on the industrial formulation IAPWS-IF97 2008. p. 7–150. Wyllie MRJ, Gregory AR, Gardner LW. Elastic wave velocities in heterogeneous and porous media. Geophysics. 1956;21(1):41–70. Zyvoloski G, Robinson B, Dash Z, Trease L. Models and methods summary for the FEHMN application. Tech. rep., Los Alamos National Lab., NM (United States). 1996. https ://www.osti.gov/servl ets/purl/36645 5. Accessed 12 Oct 2016. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Journal

Geothermal EnergySpringer Journals

Published: Nov 21, 2019

References