Elucidating the relationship between aquifer heterogeneity and thermal dispersivity using stochastic heat transport simulations

Abstract

Homogeneous parametrization of aquifer properties is conventionally adopted due to data limitations and computational constraints, where the determination of effective model parameters and consideration of heterogeneity remains a key challenge. In this study, numerical simulations were developed to investigate the effect of heterogeneity on groundwater heat transport processes (conduction, dispersion and advection), and elucidate the relationship between aquifer heterogeneity and thermal dispersivity (α_L). A simple finite element model was developed to simulate the operation of a geothermal system within a closed heterogeneous aquifer using stochastic permeability realizations. Sensitivity analysis revealed how the α_L value is affected by (1) the scale of heterogeneity (R), (2) the distance between the well pair (L), and (3) the Darcy flux (q_in). The effective permeability of the medium was found to continuously decrease as R increases; however, cold water reaches the outflow side faster through the high-conducting channels; thus, the thermal breakthrough time decreases in heterogeneous media. Calibrated α_L values allowed homogeneous models to accurately fit the breakthrough curves obtained from the heterogeneous simulations. Increase in both R and L corresponded to higher α_L, but R remained the dominant parameter. Modification of the Péclet number revealed that the thermal conduction weakens, yet does not remove, the effect of heterogeneity or α_L. It is concluded that well-calibrated α_L values are necessary for accurate predictions of heat transport processes in groundwater flow systems. This study provides guidance on the estimation of the longitudinal α_L value for use in numerical modelling of practical geothermal problems.

Résumé

La paramétrisation homogène des propriétés aquifère est généralement adoptée en raison des limitations des données et des contraintes de calcul, pour lesquelles la détermination des paramètres effectifs du modèle et la prise en compte de l’hétérogénéité restent une difficulté majeure. Dans cette étude des simulations numériques ont été développées pour investiguer les effets de l’hétérogénéité sur les processus de transport de chaleur dans les eaux souterraines (conduction, dispersion, advection) et élucider la relation entre l’hétérogénéité de l’aquifère et la dispersivité thermique (α_L). Un modèle simple en éléments finis a été mis en œuvre pour simuler l’opération d’un système géothermal au sein d’un aquifère hétérogène fermé en utilisant des réalisations stochastiques de perméabilité. L’analyse de sensibilité révèle comment la valeur α_L est affectée par (1) l’échelle de l’hétérogénéité (R), (2) la distance entre les forages constituant le doublet (L) et (3) le flux de Darcy (q_in). La perméabilité effective du milieu a montré une diminution continue avec l’augmentation de R. Cependant l’eau froide atteint plus rapidement le côté de la sortie en circulant dans des conduits de perméabilité élevée ce qui explique la diminution du temps de percée thermique dans un milieu hétérogène. Les valeurs calibrées α_L permettent aux modèles hétérogènes de reproduire précisément les courbes de restitution obtenues avec les simulations hétérogènes. Une augmentation de R et L correspond à un α_L plus élevé, mais R reste le paramètre le plus important. La modification du nombre de Péclet indique que la conduction thermique faiblit, bien qu’elle ne supprime pas l’effet de l’hétérogénéité ou α_LEn conclusion, des valeurs bien calibrées de α_L sont nécessaires pour des prédictions précises des processus de transport de chaleur dans les systèmes aquifères. Cette étude fournit une orientation pour l’estimation de la valeur α_L longitudinal à utiliser pour la simulation numérique de problèmes géothermiques pratiques.

Resumen

La parametrización homogénea de las propiedades de los acuíferos se adopta convencionalmente debido a las limitaciones de los datos y las restricciones computacionales, donde la determinación de los parámetros efectivos del modelo y la consideración de la heterogeneidad siguen siendo un desafío clave. En este estudio, se desarrollaron simulaciones numéricas para investigar el efecto de la heterogeneidad en los procesos de transporte de calor de las aguas subterráneas (conducción, dispersión y advección), y dilucidar la relación entre la heterogeneidad del acuífero y la dispersividad térmica (α_L). Se desarrolló un sencillo modelo de elementos finitos para simular el funcionamiento de un sistema geotérmico dentro de un acuífero heterogéneo cerrado utilizando realizaciones estocásticas de la permeabilidad. El análisis de sensibilidad reveló cómo el valor de α_L se ve afectado por (1) la escala de heterogeneidad (R), (2) la distancia entre el par de pozos (L), y (3) el flujo Darcy (q_in). Se comprobó que la permeabilidad efectiva del medio disminuye continuamente a medida que aumenta R. Sin embargo, el agua fría alcanza más rápidamente el lado de salida a través de los canales de alta conductividad, por lo que el tiempo de rotura térmica disminuye en los medios heterogéneos. Los valores α_L calibrados permitieron que los modelos homogéneos se ajustaran con precisión a las curvas de penetración obtenidas a partir de las simulaciones heterogéneas. El aumento tanto de R como de L se correspondió con una mayor α_L, pero R siguió siendo el parámetro dominante. La modificación del número de Péclet reveló que la conducción térmica debilita, aunque no elimina, el efecto de la heterogeneidad o de α_L. Se concluye que son necesarios valores de α_L bien calibrados para realizar predicciones precisas de los procesos de transporte de calor en los sistemas de flujo de aguas subterráneas. Este estudio proporciona orientación sobre la estimación del valor longitudinal de α_L para su uso en la modelización numérica de problemas geotérmicos prácticos.

摘要

由于数据限制和计算约束，通常采用含水层属性参数的均质化，其中有效模型参数的确定和异质性的考虑仍然是一个关键挑战。在本研究中，开发了数值模拟来研究异质性对地下水热传输过程（传导、弥散和对流）的影响，并阐明含水层异质性与热弥散性 (α_L ) 之间的关系。开发了一个简单的有限元模型，通过随机渗透率实现，在封闭的异质性含水层中模拟地热系统的运行。敏感性分析揭示了 α_L 值如何受（1）异质性尺度（R），（2）对井之间的距离（L），以及（3）达西流量（q_in）的影响。随着 R 的增加，介质的有效渗透率被发现持续下降。然而，通过高导流通道，冷水更快到达出口侧，因此在异质介质中，热穿透时间减少。校准的 α_L 值使得均质模型能够准确拟合从异质模拟中获得的穿透曲线。R 和 L的增加对应于更高的 α_L，但 R 仍然是主导参数。Péclet 数的修改显示，热传导削弱了，但并未消除异质性或α_L 的影响。结论表明，为了在地下水流动系统中准确预测热传输过程，必须精确校准 α_L 值。本研究为在实际地热问题的数值建模中估算纵向 α_L值提供了指导。

Resumo

A parametrização homogênea das propriedades do aquífero é convencionalmente adotada devido à limitações de dados e restrições computacionais, onde a determinação de parâmetros efetivos do modelo e a consideração da heterogeneidade continuam sendo um desafio fundamental. Neste estudo, simulações numéricas foram desenvolvidas para investigar o efeito da heterogeneidade nos processos de transporte de calor das águas subterrâneas (condução, dispersão e advecção) e elucidar a relação entre heterogeneidade do aquífero e dispersividade térmica (α_L). Um modelo simples de elemento finito foi desenvolvido para simular a operação de um sistema geotérmico dentro de um aquífero heterogêneo fechado usando realizações de permeabilidade estocástica. A análise de sensibilidade revelou como o valor α_L é afetado por (1) a escala de heterogeneidade (R), (2) a distância entre o par de poços (L) e (3) o fluxo de Darcy (q_in). A permeabilidade efetiva do meio diminuiu continuamente à medida que R aumenta. No entanto, a água fria atinge mais rápido o lado de saída através dos canais de alta condução, então o tempo de ruptura térmica diminui em meios heterogêneos. Valores de α_L calibrados permitiram que modelos homogêneos se ajustassem com precisão às curvas de ruptura obtidas das simulações heterogêneas. O aumento em R e L correspondeu a α_L mais alto, mas R permaneceu como o parâmetro dominante. A modificação do número de Péclet revelou que a condução térmica enfraquece, mas não remove, o efeito da heterogeneidade ou α_L. Conclui-se que valores de α_L bem calibrados são necessários para previsões precisas de processos de transporte de calor em sistemas de fluxo de águas subterrâneas. Este estudo fornece orientação sobre a estimativa do valor de α_L longitudinal para uso em modelagem numérica de problemas geotérmicos práticos.

Introduction

Groundwater plays a crucial role in geothermal energy extraction, as an important factor in heat transport processes. Geothermal water extraction from aquifers depends on the rate of replacement; therefore, water injection is necessary to ensure long-term sustainable water production. In geothermal applications (e.g. well doublet systems), the system has to maintain a balance between the water extraction and injection to avoid aquifer depletion; however, the injected water can significantly reduce the temperature of the extracted thermal water through thermal feedback mechanisms (Banks 2009). The latter phenomenon is typically known as a thermal breakthrough (Clyde and Madabhushi 1983), which can be effectively modelled both in homogeneous and heterogeneous media.

The most popular approach is that the groundwater flow system or the geothermal field is characterized by formations, layers, different rock types or hydrostratigraphic units, and that the homogeneous approximation is applied for each layer/domain. This modeling approach has been used to investigate the role of conductive and convective heat transfer in the Rheingraben, Germany (Clauser and Villinger 1990), to study the relationship between the results of a convective groundwater flow model and surface observations beneath Lake Chad, West-Central Africa (Lopez et al. 2016), to characterize the groundwater flow system and heat transport mechanism in a regional-scale semi-arid basin in Central Mexico (Ortega Guerrero 2022), to investigate the origin of deep geothermal anomalies in the German Molasse Basin, Germany (Przybycin et al. 2017), to identify the presence of hot springs in the Kurşunlu geothermal field, Turkey (Akar et al. 2021), to predict hydrothermal processes during the operational lifetime of a deep geothermal reservoir development project in Northeast German Basin, Germany (Blöcher et al. 2010), and to determine the hydrothermal coupling between the confined and unconfined aquifers in the Buda Thermal Karst, Hungary (Havril et al. 2016; Szijártó et al. 2021) etc. In some models, the diversity of the individual categories nearly imitated the large-scale geological heterogeneity—e.g. in the Momotombo geothermal field, Nicaragua (Porras et al. 2007), in the Larderello–Travale geothermal field, Italy (Romagnoli et al. 2010), or in Rotorua Geothermal Field, New Zealand (Ratouis et al. 2016).

Data scarcity is one of the main reasons for developing homogeneous models; however, geostatistical methods (e.g. Goovaerts 1997; Bruckmann and Clauser 2020) can be used to account for geological heterogeneity, even with limited data. Geostatistics-based stochastic modeling provides an established approach to simulate transport processes in heterogeneous media, enabling the quantification of uncertainty in modelled variables (Deutsch and Journel 1997; Pyrcz and Deutsch 2014). For groundwater flow systems, mass (i.e., solute) transport has been studied more intensively (Ferguson 2007) than heat transport. One reason could be that thermal diffusivity exceeds the hydromechanical diffusivity by several orders of magnitude. As such, one expects heat transport to be less sensitive to aquifer heterogeneities than solute transport; however, experimental evidence suggests that heterogeneity can significantly control heat transport in aquifers.

Shook (2001) attempted to predict the thermal breakthrough in heterogeneous media using tracer tests. The results indicated that an injection tracer test can be readily analyzed to determine the velocity (and therefore the breakthrough time) of temperature fronts in heterogeneous, porous, permeable media saturated by a single-phase fluid, provided that thermal conduction and dispersion can be neglected. Notably, the temperature front moves slower than the pore fluid when thermal equilibrium is assumed between the pore fluid and the matrix. This difference is explained by the retardation factor expressing that energy is transported through both the fluid-filled pores and the rock matrix, making the heat plume is retarded relative to fluid particles (Shook 2001). The onset of the thermal breakthrough is defined as the point when the temperature of the produced water begins to decline (Gringarten 1978; Gringarten and Sauty 1975; Clyde and Madabhushi 1983). Irregular flow paths in a heterogeneous medium, as well as faults and fracture zones, can result in premature thermal breakthrough, leading to economic risks and operational challenges. Instances of premature thermal breakthrough have been reported in numerous geothermal reservoirs, including The Geysers in the USA (Beall et al. 1994), Miravalles in Costa Rica (Parini et al. 1996) and Cerro Prieto in Mexico (Ocampo et al. 1998). Fadel et al. (2022) investigated the causes of premature thermal breakthrough in a hydrothermal project within the Upper Jurassic carbonates of the North Alpine Foreland Basin in Germany, concluding that regional reservoir heterogeneity plays a crucial role in early thermal breakthrough.

Stochastic-based model analyses in hydrogeology that investigate the effects of heterogeneity on heat transport are limited. Crooijmans et al. (2016) studied the impact of facies heterogeneity on the performance of well doublets in low-enthalpy geothermal sedimentary reservoirs. Babaei and Nick (2019), Liu et al. (2019) and Di Dato et al. (2022) examined the effect of heterogeneity on the thermal breakthrough time, thereby on the lifetime of well doublets, finding that the thermal breakthrough time decreases due to heterogeneity. Additionally, Di Dato et al. (2022) revealed in shallow geothermal systems, that heterogeneity and thermal dispersivity had a negligible effect on the recirculation ratio (the fraction of water returning to the production well), which quantifies the long-term evolution of thermal feedback. Well patterns were studied by Liu et al. (2020) in heterogeneous media (Dezhou geothermal field, China) to maximize the performance of the geothermal system.

Since heterogeneous transport models embedding spatially variant parametrization can be computationally demanding, several modelers choose to simulate homogeneous models equipped with “effective” or “upscaled” model parameters (e.g. Sanchez-Vila and Fernàndez-Garcia 2016; de Barros and Rubin 2011). Mechanical dispersivity and heterogeneity are known to be intimately related (Dagan 1989; Gelhar et al. 1992; Schulze-Makuch 2005); therefore, for solute transport problems, the most studied effective parameter has been by far the longitudinal dispersivity. At the local scale, dispersion results from velocity variations within pore spaces, while at the macroscale, contrasts in hydraulic conductivities (or permeabilities) increase the plume size, causing a “macrodispersive” effect (e.g. Dagan 1989). Less effort has been placed into finding a link between heterogeneity and thermal dispersivity; hence, the role of thermal dispersivity is often ignored in heat transport modelling because advection and diffusion are considered more dominant processes (e.g. Porras et al. 2007; Tóth et al. 2020; Szijártó et al. 2021, etc.).

Under certain conditions and simplifications, it is understood how to relate the solute transport equation to the heat transport equation (e.g. Shook 2001; Hecht‐Méndez et al. 2010); however, it is questionable whether thermal dispersivity can be equated to hydrodynamic/mechanical dispersivity. Willemsen and Groeneveld (1989) indicated that the use of hydrodynamic dispersivity to characterize the thermal dispersion may lead to an overestimation of the process. Bear (1972), Ingebritsen and Sanford (1998) argued that the effects of thermal dispersion are negligible compared to conduction and set dispersivity to zero. Similar conclusions were presented by Constantz et al. (2003) and Vandenbohede et al. (2008) who found that thermal dispersion in heat transport is relatively unimportant and thermal dispersivity values are significantly smaller than mechanical dispersivity values. Conversely, other studies reached the opposite conclusion. Windqvist and Hyden 1976, Sauty et al. (1982a, b) and de Marsily (1986) found that thermal dispersivity is on the same order of magnitude as mechanical dispersivity. Hidalgo et al. (2009) and Molina-Giraldo et al. (2011) investigated the effect of thermal dispersivity on groundwater heat exchangers (GHEs). Hidalgo et al. (2009) performed steady-state numerical model simulations on temperature plumes and defined a relationship between the transverse thermal dispersivity and the variance of logarithm of the heterogeneous hydraulic conductivity field. Molina-Giraldo et al. (2011) determined a Darcy flux range in their analytical model calculations where the effect of thermal dispersivity on the temperature plume distribution around the GHE is not negligible. That is, the role and the value of thermal dispersion is less clear, and strongly controversial, even today; thus, the question of the relationship between heterogeneity and thermal dispersivity is still open, as pointed out by Anderson (2005), Ferguson and Woodbury (2005), Vandenbohede et al. (2008), Hidalgo et al. (2009) and Molina-Giraldo et al. (2011) Furthermore, the discrepancy between different scales of examination/observation scales alone have not resolved this contradiction.

The primary aim of this study was twofold. First, the study aimed to elucidate, through detailed numerical model calculations, how the spatial variability of aquifer permeability affects coupled groundwater flow and heat transport processes, particularly in relation to thermal breakthrough. A two-dimensional (2D) closed aquifer–aquiclude system was built up to simulate heat transport in a synthetic well doublet, reproducing the re-injection of cooler water into the aquifer along the upstream boundary and the extraction of warmer water from the downstream boundary. Sequential Gaussian simulations (SGS) were used to represent random spatial permeability fields for model parametrization. The effects of different permeability heterogeneities, model length and Darcy flux were systematically assessed using control parameters, such as outlet temperature, thermal breakthrough time (defined as the time at which 10% of cooling occurs) and effective permeability to demonstrate the impact of heterogeneity on heat transport and the flow system. The second objective was to correlate the scale of heterogeneity with effective thermal dispersivity. To achieve this, homogeneous model calculations were carried out with varying longitudinal dispersivities to match the results from heterogeneous permeability models, allowing for the determination of effective (upscaled) α_L values. The numerical models were solved both with and without thermal conduction to isolate the effects and better understand interactions in heat transport processes.

Model description

Reference area

For the construction of the synthetic model, which studies a deep geothermal system, model parameters were derived from a reference area in the Békés Basin, located in the southeastern part of the Pannonian Basin, Hungary, where Quaternary and Neogene sediments were deposited in a thickness of more than 2 km (Horváth et al. 2015). The Pannonian Basin is known for its excellent geothermal potential, with heat fluxes and geothermal gradients significantly exceeding the continental average (Lenkey et al. 2021). In the Békési Basin, the surface heat flux and geothermal gradient reach values of 100 mW/m² and 50 °C/km, respectively. The Újfalu Formation is a major sedimentary geothermal aquifer in the reference area, which is often the main target of thermal water exploration in Hungary due to its high permeability (~10^–12 m²), temperature (~120 °C) and large thickness (several hundred meters) (Buday et al. 2015; Horváth et al. 2015). The Újfalu Formation is a typical heterogeneous sedimentary aquifer with alternating semiconsolidated sediments, consisting of sandstone with limited interbeds of poorly permeable muds and other fine-grained beds. This property makes it suitable for numerical modelling of geothermal reservoirs using the stochastic method. Typical well spacing for geothermal systems in the region ranges from 500 to 1000 m, with cold water injected through screened sections ~100 m long. It should be noted that the reference geothermal area serves only to frame the problem, as the main objective of the study is to uncover a more general relationship between heterogeneity and thermal dispersivity through synthetic simulations.

Permeability fields

The effect of heterogeneity on heat transport was analyzed within a stochastic modeling framework. Lognormally distributed permeability fields k (m²) were generated using unconditional Sequential Gaussian Simulations (SGS) with the SGeMS code (Remy et al. 2009), which is a widely used algorithm for the stochastic characterization of geological sites (Deutsch and Journel 1997; Delbari et al 2009; Manchuk and Deutsch 2012; Nguyen et al. 2016). SGS randomly visits a cell on the grid and determines its value by kriging, then stores the value and uses it for calculating the value of the next randomly visited cell. The process continues like this until the entire grid is generated. The result of the process is countless realizations with the same probability depending on the order of visited cells. Conditioning data were not used, so the distributions had standard normal distribution (zero mean and unit standard deviation). In the next step, these were transformed to a mean of μ = –12, where μ = log₁₀(k), but the standard deviation remained the same.

Each distribution represents an independent and identically distributed (IID) permeability field with a given mean (10^–12 m²), variance (one order of magnitude) and isotropic spatial correlation length denoted by R. The different realizations were statistically uniform, while only the correlation length (heterogeneity scale) was varied. Figure 1 illustrates the lognormal distribution and the variogram of permeability data for one realization where R = 10 m. The heterogeneity scale (R) was defined as the correlation length of the variogram.

Fig. 1

a Histogram and b variogram of the lognormally distributed permeability field of one realization with R = 10 m. Blue is the horizontal direction, red is the vertical direction

Specifically, four isotropic correlation lengths were tested: R = 5, 10, 20 and 50 m, adjusted according to the numerical resolution and model size. Preliminary test simulations were conducted to determine the number of realizations required to obtain a statistically stable solution at a given correlation length. It was established, that the variability of the mean and the standard deviation of the numerical solutions stabilized above 10 realizations, so 10 realizations were generated for each correlation length. Figure 2 presents one permeability realization for different correlation lengths, while Fig. 3 shows four different permeability realizations at R = 20 m.

Fig. 2

Heterogeneous permeability fields using the same mean (–12) and variance (1) with different correlation lengths: a R = 5 m, b R = 10 m, c R = 20 m and d R = 50 m. The logarithm of the permeability is shown

Fig. 3

a-d Heterogeneous permeability realizations using the same mean (–12) and variance (1) at R = 20 m. The logarithm of the permeability is shown

Physical model

The solution of heat transport in porous media requires the coupled handling of the continuity equation (Eq. 1), Darcy’s law (Eq. 2) and the heat transport equation (Eq. 3) governing the mass conservation, the momentum conservation without thermal buoyancy and the energy conservation by advection, thermal diffusion and thermal dispersion (Nield and Bejan 2017),

$$\nabla \mathbf{q}=0,$$

(1)

$$\mathbf{q}=\boldsymbol{ }-\frac{k}{\upeta } \nabla p,$$

(2)

$$\left(\phi\rho_\text{w}c_\text{w}+\left(1-\phi\right)\rho_\text{m}c_\text{m}\right)\frac{\partial T}{\partial t}=-\rho_\text{w}c_\text{w}\mathbf q\nabla T+\nabla\left\{\left[\phi\lambda_\text{w}+\left(1-\phi\right)\lambda_\text{m}+\boldsymbol\Lambda^{\mathrm{disp}}\right]\nabla T\right\},$$

(3)

where q, p and T denote the unknown Darcy flux, pressure and temperature, respectively, while Φ, k, η, t, Λ^disp denote the porosity, permeability, dynamic viscosity of water, time and the thermal dispersion tensor, respectively. The density, specific heat and thermal conductivity of the water and the matrix are represented by ρ, c and λ with subscripts w or m, respectively. The thermal dispersion tensor, which is used only in homogeneous models, is defined by

$$\boldsymbol\Lambda_{\boldsymbol i\boldsymbol j}^{\mathrm{disp}}=\rho_\text{w}c_\text{w}\left[\alpha_\text{T}V_\text{q}\delta_{ij}+\left(\alpha_\text{L}-\alpha_\text{T}\right)\frac{q_iq_j}{V_\text{q}}\right],$$

(4)

where α_L, α_T V_q, δ denote the longitudinal and transversal thermal dispersivity, the magnitude of the Darcy flux and the Kronecker delta, respectively. The parameters for both the heterogeneous and the homogeneous models are summarized in Table 1. The temperature dependence of the water density and viscosity was neglected during the simulation to concentrate on the effect of heterogeneity. Consequently, the equation system of Eqs. (1)–(3) is only partially coupled. While the Darcy flux influences the time-dependent temperature field, the temperature does not influence the steady-state Darcy flux, due to the constant water density and the viscosity. In homogeneous simulations, the transversal dispersivity was fixed to the longitudinal one, α_T = α_L/10.

Table 1 Parameter values for the base model and the studied parameter range (italic). Thermal dispersivity is present only in homogeneous aquifer

A 2D numerical model was developed based on the reference area described in section ‘Reference area’. The heterogeneous aquifer–aquiclude system was characterized by a length of L = 500 m and a thickness of d = 300 m. The aquifer with a thickness of 100 m was bounded above and below by two aquicludes (Fig. 4a). Heterogeneous permeability distributions were assumed for the aquifer (as described in detail in ‘Permeability fields’), while the permeability of the two aquicludes was set to zero. The heat transport process included both advection and conduction in the aquifer, but only conduction in the aquicludes. Temperature boundary conditions were fixed along the top and the bottom boundaries as T_top = 115.5 °C and T_bot = 124.5 °C, respectively. The initial condition for the temperature was calculated assuming only conduction. Cold water at T_in = 50 °C was injected through the left side of the aquifer, with an average horizontal Darcy flux of q_in = 10^–7 m/s, that leaves through the right side of the aquifer. Boundary and initial conditions for the temperature and the flow are presented in Fig. 4b.

Fig. 4

Base model a illustrating the permeability field and b presenting the boundary conditions for the Darcy flux (blue) and the temperature (red) as well as the initial thermal condition

The coupled flow and heat transport problem (Eqs. 1–3) was solved using COMSOL Multiphysics v5.3, a finite element numerical software package (Zimmermann 2006). The model domain was discretized by triangle elements with a maximum size in the aquifer of 1 m to ensure the appropriate resolution for the small-scale heterogeneities (R = 5 m). In the aquiclude, the maximum element size increased gradually from 1 to 10 m away from the aquifer. Thus, the mesh consisted of 155,782 finite elements for the base model, which varied between 63,060–558,376 as the model length was increased from 200 to 2000 m. The pressure and temperature fields within the elements were approximated by quadratic and linear polynomials, respectively. Simulations continued until a stationary solution was reached, typically requiring 500 years, with a range from 50 to 5000 years, depending on model length and Darcy flux. Initial time steps were increased exponentially from 0 to 2 years to stabilize transition effects, after which the maximum time step was fixed, generally at 0.5 years, ranging from 0.05 to 5 years (depending on model length and the Darcy flux). Time-dependent simulations were run on an Intel Workstation with eight cores, requiring ~3.5–12 h of CPU time and 1.5–4 GB memory per model. In total, 348 time-dependent coupled flow and heat transport simulations were conducted.

Model parametrization

The influence of several key parameters characterizing heat transport in heterogeneous porous media was investigated with and without the consideration of thermal conduction. These parameters include:

• 1. The scale of heterogeneity (the range/correlation length of the permeability field), R

• 2. The model length (distance between the injection and production well), L

• 3. The inlet Darcy flux (volume flux of water through the aquifer), q_in

First, the heterogeneity is present at all spatial scales in the groundwater flow systems, making it essential to understand its impact on the heat transport processes (Ferguson 2007; Di Dato et al. 2022). The smaller the scale of heterogeneity, the closer the numerical result is to the homogeneous solution; therefore, the lowest value was chosen as R = 5 m. The upper limit was R = 50 m, which is comparable to the thickness of the aquifer (100 m). However, it is important to note that R = 50 m may exceed the model’s representative elementary volume (REV) size, as defined by Bear (1972) and Freeze and Cherry (1979). Second, it is known that mechanical longitudinal dispersivity increases with the observation scale (Gelhar et al. 1992; Schulze-Makuch 2005). To reveal the relationship between the heterogeneity scale and the observation scale, model lengths representing typical distances between injection and production wells (L = 200–2000 m) were used (Babaei and Nick 2019). Third, significant variability in well yield exists; thus, the inlet Darcy flux was varied by two orders of magnitude (q_in = 10^–8–10^–6 m/s). Finally, in order to demonstrate the importance of thermal conduction, it was turned on and off during the modeling. Neglecting thermal conduction made model heat transport behavior more analogous to mass transport, allowing solutions to be rescaled by using the retardation factor—see the electronic supplementary material (ESM). Additionally, the role of thermal conduction (related to advection) can change with different yields, model lengths and heterogeneity scales.

The influence of these model parameters on groundwater flow and heat transport was quantified by using control parameters, including:

• Outlet temperature (T_out)

• Time of 10% cooling (t_10%)

• Stationary temperature (T_stat)

• Effective permeability (k_eff)

The outflow temperature time series averaged on the right side of the aquifer (T_out) were calculated for each model realization. The time of cooling by 10% on the outlet side (i.e., the 10% temperature drop between the initial and stationary temperature) is a well-defined and characteristic point of the outlet breakthrough curves (BTCs), which were used as a reference for the model interpretation. Additionally, the time of 10% cooling can be important information for the operation of the geothermal system. The value of the stationary temperature effectively characterizes the long-term operation of the well-doublet and it involves the ratio of the advection and conduction (Péclet number) in the model calculation. The effective permeability quantifies the effect of heterogeneity on the flow system. It was determined from the linearized form of the Darcy’s law (see Eq. 2), which presents how much effective permeability k_eff would cause Δp pressure difference between the vertical boundaries of the model if the medium were homogeneous. The average Δp pressure difference, which is independent of time, was calculated during the simulations, and effective permeability was estimated by

$${k}_{\text{eff}}=\boldsymbol{ }\frac{{q}_{\text{in} }\eta L}{\Delta p}$$

(5)

To elucidate the relationship between the heterogeneity and the thermal dispersivity, firstly the time of 10% cooling in heterogeneous media was determined. (Heterogeneous models do not comprise thermal dispersivity, the dispersion is represented purely by the heterogeneity of the permeability.) Secondly, homogeneous model calculations with different thermal dispersivity values (in a range of α_L = 0–50 m) were performed, and the time of 10% cooling was computed based on the T_out curves. Thirdly, the time of 10% cooling obtained from heterogeneous and homogeneous model results were compared to reveal the relationship between heterogeneity and thermal dispersivity.

Results

Representative simulations

Figure 5 shows the visual impact of the permeability heterogeneity on flow and temperature at small-scale (R = 5 m, left) and large-scale (R = 20 m, right) heterogeneities (correlation lengths). Irregular flow paths form within the medium as the groundwater flows through the network of high-permeability zones (Fig. 5c,d). At larger scales, the number of distinct “flow channels” decreases but the channels become wider. The magnitude of the Darcy flux is very sensitive to the local permeability values, resulting in a significant difference between the flow channels and the low-permeability zones up to four orders of magnitude. It is important to note that the fluid density is constant; thus, the flow is stationary. Qualitatively, as the scale of heterogeneity decreases, the flow pattern becomes more similar to the homogeneous solution.

Fig. 5

Permeability fields with heterogeneity scales of a R = 5 m and b R = 20 m. Logarithm of the stationary Darcy flux magnitude at c R = 5 m and d R = 20 m. Temperature snapshots after 50 years e–f in the absence and g–h in the presence of thermal conduction at R = 5 m and R = 20 m, respectively. Isotherm of T = 85 °C is denoted by the green contour

On the other hand, the temperature field is time-dependent as the cold water propagates from the inlet side (left) toward the outlet. Figure 5e–h illustrates the temperature snapshots after 50 years at different scales, both in the absence (Fig. 5e,f) and in the presence (Fig. 5g,h) of thermal conduction. Neglecting thermal conduction, the thermal front is very complex, it follows the margin of the permeable network. It points out that the only heat transport mechanism is the advection. For R = 20 m, 50 years is sufficient for cold water to reach the outlet, whereas for R = 5 m, the thermal front is less channeled, and cold water does not reach the outlet boundary, suggesting that the thermal breakthrough time is reduced at larger heterogeneity scales (Molnár and Galsa 2022). This behavior is analogous to the early arrival times observed for solute particles migrating in heterogeneous porous media (e.g. Pedretti et al 2013; Bianchi and Pedretti 2017). Allowing thermal conduction (Fig. 5g,h), the heat plume front is less complex, as thermal diffusion smooths the sharp margins. It seems that the role of heterogeneity is suppressed by conduction. Thus, conduction controls the impact of heat plume channeling in preferential flow zones, although the difference between the temperature snapshots of R = 5 and 20 m is observable. Finally, it was noticed that cooling appears in the bedding aquicludes only in the presence of thermal conduction, which is the only mechanism that can transfer heat from the aquicludes.

Parametric analysis

This section presents the effects of three model parameters, including the scale of heterogeneity (R), the model length (L) and the inlet Darcy flux (q_in) on the control parameters. For each set of model parameters in the heterogeneous permeability field (see ‘Permeability fields’), 10 realizations were generated and analysed; however, only the averaged values with standard deviations and the averaged time series of the numerical simulations are shown here. The interpretation of numerical results will first be provided for cases without thermal conduction as a reference state, followed by cases that include thermal conduction.

Figure 6a shows the averaged outlet temperature curves as a function of the time for the homogeneous model and for the models with different scales of heterogeneity, both without thermal conduction (solid line) and with thermal conduction (dotted line). As the scale of heterogeneity increases, the cold water arrives earlier at the right side of the model, indicating the cooling starts earlier, so the thermal breakthrough time decreases. This is consistent with qualitative observations shown in Fig. 5e,f. The time needed for cooling by 10% along the outlet side affirms that increasing the scale of heterogeneity decreases the breakthrough time (Fig. 6b). While the time of 10% cooling is 73 years at the range of R = 5 m (by ~23 years less than in the homogeneous model), it is only 55 years at R = 20 m. However, the standard deviation increases with the scale of heterogeneity, because t_10% is sensitive to the connectivity of larger zones with high permeability. Although the cold water arrives earlier at the outlet side of the model, complete cooling of the aquifer requires more time due to the presence of low permeability zones. Consequently, the slope of the temperature curves of high R decreases—that is the interval of the temperature change becomes wider. Actually, this phenomenon refers to the presence of thermal dispersion occurring in a heterogeneous medium. In the absence of thermal conduction, the stationary outlet temperature is 50 °C.

Fig. 6

Average outlet temperature time series at a different scales of heterogeneity, c model lengths and e inlet Darcy fluxes. The time of 10% cooling with its standard deviations plotted against b the scale of heterogeneity, d the model length and f the inlet Darcy flux. d and f Shows the nondimensional t_10% times for better comparison. Solid and dotted lines denote the models without and with conduction, respectively. Standard deviation is calculated from the solution of the 10 permeability realizations

In the presence of thermal conduction (dotted line), the effect of permeability heterogeneity becomes smaller. The difference among the time series is less considerable and the stationary outlet temperature stabilizes at a higher value (~75 °C). This indicates that thermal conduction weakens the influence of dispersion in heterogeneous media. This effect is quantitatively confirmed in Fig. 6b where the decrease in the time of 10% cooling is significantly smaller than in the case without conduction. The decrease in breakthrough time (t_10%) is ~40 years between the homogeneous model and the R = 20 m model without conduction. However, this value reduces to 13 years in the presence of conduction. It is noted that in the homogeneous model without conduction, the outlet temperature drops instantaneously (solid black curve in Fig. 6a). Slight deviation from the perfectly vertical direction is the consequence of numerical diffusion which is strongly exceeded by the physical thermal diffusion. It is noted that the cooling appears the latest in the homogeneous model without conduction (~97 years); therefore, the curves cross each other in Fig. 6b.

It is evident that the cooling of the outlet water starts later as the model length increases, and the stationary temperature is stabilized at 50 °C in the absence of thermal conduction (Fig. 6c). However, allowing conduction, the resulting steady-state temperature increases substantially and the difference between the curves (with and without conduction) becomes larger with increasing model length. The longer the model, the conduction has more time to heat up the fluid, resulting in a slighter cooling due to conduction as L increases.

In order to compare the results of the 10% cooling obtained from numerical models with different lengths, t_10% was nondimensionalized by using the advective time scale (L/q_in) and plotted against model length. Figure 6d proves again that the cool water arrives earlier at the outlet side without conduction at all model lengths. Nevertheless, the nondimensional time for 10% cooling seems larger in the longer model domain, which suggests that the role of heterogeneity/dispersion is reduced by time and/or distance, because the thermal front is less sharp, the temperature gradient is lower (Eq. 3). The standard deviation decreases with increasing model length due to the fact that the larger model domain results in statistically more robust numerical solution.

At higher inlet Darcy flux (q_in = 10^–6 m/s), the breakthrough time decreases both with and without thermal conduction (Fig. 6e). The stationary outlet temperature is only marginally influenced by conduction. Slower groundwater flow enhances the role of conduction relative to advection (i.e., reduces the Péclet number) resulting in a much higher outlet temperature. The nondimensional t_10% curve (Fig. 6f) remains constant when neglecting thermal conduction, demonstrating that only heat advection influences the outlet temperature; however, in the presence of thermal conduction, the picture is more complex. At a low Darcy flux (q_in = 10^–8 m/s), the role of horizontal conduction becomes relatively more important reducing the nondimensional t_10%, which approximates zero. At high Darcy flux (q_in = 10^–6 m/s), the role of horizontal conduction is insignificant, as the breakthrough time approaches that of the model without conduction (blue). The standard deviations are almost constant, except for simulations where the role of conduction is large compared to advection (q_in ≤ 2 × 10⁻⁸ m/s with thermal conduction). Simulations with a low Péclet number result in smaller standard deviations.

The stationary outlet temperature is independent of the scale of heterogeneity both in the absence (T_stat = 50 °C, not shown) and in the presence of thermal conduction (Fig. 7a). Although T_stat does not vary with R, the increasing standard deviation indicates the higher variability among different realizations. This is consistent with the analysis presented in Fig. 5, increasing R renders the random spatial fields (i.e., the k fields) less well-mixed and more spatially ordered, favoring preferential flow and heat transport channeling.

Fig. 7

Stationary outlet temperature plotted against a the scale of heterogeneity, b the model length and c the inlet Darcy flux. Blue and red colors denote the model results without and with conduction, black dashed line represents the purely conductive solution (b and c)

Figure 7b indicates that larger distances between the wells increase the stationary outlet temperature if the thermal conduction is taken into account, as T_stat tends toward 120 °C (the average temperature at this depth without injection). On the contrary, enhancing the Darcy flux (injection rate) reduces the stationary outlet temperature (Fig. 7c) by decreasing the role of thermal conduction to advection (increases Péclet number). In the absence of conduction, the stationary temperature is constant, T_stat = 50 °C, which corresponds to the injected water temperature.

Finally, the effect of the three model parameters on the value of the effective permeability was examined. Since the temperature dependence of the water density was disregarded, the effective permeability computed from the flow field is the same for the models of including or ignoring thermal conduction. Although the geometric mean of permeability is 10^–12 m² = 1 Darcy for all permeability realizations, the effective permeability continuously decreases as the scale of heterogeneity increases (Fig. 8a). Approximately 30% decrease was noticed in the simulations; thus, the effective permeability is lower than the average permeability of the heterogeneous media and the difference grows with the scale of the heterogeneity. The reason for this behavior can be explained as follows. Recall that prescribed flux boundary conditions are set up at the model inlet and outlet, and no flux conditions are applied on the remaining aquifer sides. Then, the pressure difference in the moving fluid increases between the two sides of an obstacle, represented here by low-permeability zones. Synthetic model examples in the Appendix represent how the effective permeability calculated from the pressure difference depends on the radius of a circular obstacle. The model length does not significantly alter the effective permeability at R = 10 m (Fig. 8b), as a slight increase with L remains within the standard deviation. It is worth noting that the standard deviation is higher in models with larger correlation lengths or shorter model lengths. The inlet Darcy flux does not affect the effective permeability, since the flow pattern is independent of q_in (Fig. 8c).

Fig. 8

Effective permeability plotted against a the scale of heterogeneity, b the model length and c the inlet Darcy flux

Effective thermal dispersivity

Model calculations with homogeneous permeability fields were performed using different longitudinal thermal dispersivity values to compare the effect of heterogeneity and thermal dispersivity. Based on homogeneous model results, calibration curves were set up between the longitudinal thermal dispersivity and the time for 10% cooling.

Figure 9 illustrates the dispersivity values as a function of the time of the 10% cooling nondimensionalized by the advection time scale. Figure 9a presents the homogeneous model results at different model lengths in the absence (solid line) and in the presence (dotted line) of thermal conduction. Without thermal conduction and dispersion (α_L = 0), the curves tend toward the same point of 0.629, the retardation factor (Molina-Giraldo et al. 2011; Rau et al. 2012). This means that the 10% cooling occurs twice as late in a model which is twice as long; however, increasing the dispersivity separates the curves and shortens t_10% times, being an additional horizontal heat transport mechanism. In longer models, this time shortening is not so effective (ca. Fig. 6d). Including thermal conduction, the starting points of the curves (dotted line) separate because the conduction influences the heat transfer in a different way. The conduction delays the cooling in most model calculations, especially in longer model domains; however, in a short model domain (L = 200 m) with low dispersivity, the conduction facilitates the horizontal heat transport, resulting in lower t_10%.

Fig. 9

Calibration curves between the longitudinal thermal dispersivity and the nondimensional time of the 10% cooling at different a model lengths and b inlet Darcy fluxes without (solid line) and with (dotted line) thermal conduction. Quadratic polynomials were fitted on numerical results from homogeneous model calculations

In contrast with model length, the longitudinal thermal dispersivity is independent of the inlet Darcy-flux when there is no conduction (Fig. 9b). The curves are perfectly identical, which means that the inlet Darcy flux, being higher by one order of magnitude, results in cooling faster by one order of magnitude at each dispersivity value. In the case of slow flow (q_in = 10^–8 m/s), thermal conduction has a significant role in horizontal heat transport, decreasing by t_10%. In the case of fast flow (q_in = 10^–6 m/s), the role of thermal conduction is reduced, and the dispersivity curve tends to the pure advection curves (without conduction).

Knowing the time of 10% cooling in heterogeneous media (Fig. 6b,d,f), the calibration curves can be used to establish the relationship between the scale of the heterogeneity and the longitudinal thermal dispersivity. Figure 10 shows the calculated longitudinal dispersivity values for different heterogeneity scales, model lengths and inlet Darcy fluxes. It was found that the increase in both the heterogeneity scale and the model length corresponds to higher thermal dispersivity (Fig. 10a,b). However, the slope of the curves seems to decrease with R and L. At higher values of R, the variability among different permeability realizations is enhanced resulting in a larger standard deviation. It should be noted that R = 50 m is likely to be too large relative to the aquifer thickness of the numerical model (100 m).

Fig. 10

Calculated longitudinal thermal dispersivity at different values of a scale of heterogeneity, b model length and c inlet Darcy flux. Blue solid and red dotted lines denote results without and with conduction, respectively. Standard deviation is calculated from the 10 permeability realizations

Figure 10c illustrates the influence of the inlet Darcy flux on the calculated thermal dispersivity values. By neglecting thermal conduction, the water and the heat travel along the same path independently of the inlet Darcy flux; thus, α_L is insensitive to q_in. However, thermal conduction decouples the fluid flow and heat transport. The impact of conduction becomes dominant at low Darcy flux resulting in a relevant decrease in dispersivity by smoothing the thermal front. The effect of conduction becomes minor at high Darcy flux approximating the purely advection-driven system (without conduction) and resulting in higher dispersivity values.

Based on simulation results, it is unequivocal that the presence of thermal conduction reduces the role of thermal dispersion—for example, the heterogeneity scale of R = 10 m corresponds to α_L = 29.66 m neglecting thermal conduction, but it is decreased to α_L = 7.55 m including conduction (L = 500 m, q_in = 10^–7 m/s). The quantitative relation between R and α_L is summarized in Table 2. Figure 11 illustrates the average outlet temperature curve in the heterogeneous base model (R = 10 m, L = 500 m, q_in = 10^–7 m/s) and that in the homogeneous model using the calculated longitudinal thermal dispersivity values in the absence and in the presence of thermal conduction. This emphasizes that the effect of heterogeneous permeability on the outlet temperature can be approximated very well by applying suitable longitudinal thermal dispersivity. The relative deviation between the outlet temperature curves from the heterogeneous model and the homogeneous model with fitted thermal dispersivity is quantified by ε which is based on the area of the difference between heterogeneous and homogeneous BTCs:

$$\varepsilon =\frac{{\int }_{0}^{{t}_{\text{stat}}}\left|{T}_{\text{out}}^{\text{het}}-{T}_{\text{out}}^{\text{hom}}\right|\text{d}t}{{\int }_{0}^{{t}_{\text{stat}}}\left|{T}_{\text{out}}^{\text{het}}-{T}_{\text{stat}}\right|\text{d}t}100$$

(6)

where t_stat denotes the time, which is needed to reach the stationary outlet temperature, and het and hom indexes indicate the heterogeneous and homogeneous models, respectively. The calculated relative deviations in Table 2 indicate that BTCs from homogeneous models using properly-tuned dispersivity values approximate the BTCs from heterogeneous model simulations well.

Table 2 Relationship between the scales of heterogeneity and longitudinal thermal dispersivity neglecting and including thermal conduction in the base model (L = 500 m and q_in = 10^–7 m/s). ε shows the relative deviation between the heterogeneous and homogeneous BTCs. SD denotes standard deviation

Fig. 11

Comparison between the averaged outlet temperature curve from heterogeneous permeability realizations and the outlet temperature curve from homogeneous permeability model with the suitable longitudinal thermal dispersivity value a in the absence and b in the presence of thermal conduction. Parameters correspond to the base model: R = 10 m, L = 500 m, q_in = 10^–7 m/s. The temperature difference between the heterogeneous and homogeneous BTCs, ΔT is plotted in red

Discussion

The interaction of groundwater flow and heat transport in heterogeneous porous aquifers was systematically investigated using synthetic 2D numerical models, with a focus on thermal dispersion. It was established that injected cold water reaches the outflow side of the model faster in the heterogeneous medium (compared to the homogeneous one); thus, the permeability heterogeneity decreases the thermal breakthrough time, which agrees with the results of Liu et al. (2019) and Di Dato et al. (2022). A larger scale of heterogeneity resulted in a continuous decrease of the thermal breakthrough time until R became comparable to the model size—at R = 50 m, the size (i.e. diameter) of the heterogeneity is close to equal to the aquifer thickness. This effect can be explained by the phenomenon of thermal dispersion; additionally, a larger heterogeneity scale increased the standard variation of t_10%. Otherwise, the heterogeneity does not affect the stationary temperature, which shows the long-term effect, as Di Dato et al. (2022) pointed out in the case of a shallow geothermal well doublet. Moreover, Sbai and Larabi (2023) analyzed the relationship between doublet lifetime and thermal breakthrough using homogeneous approximation at the Dogger geothermal exploitation site in the Paris basin. They found that the doublet lifetime was twice as long as the breakthrough time. The results of this study indicate that this ratio might be larger in heterogeneous media.

A relationship was identified between the scale of heterogeneity and the value of longitudinal thermal dispersivity, based on the time required for 10% cooling along the outlet side. The results showed that both the scale of heterogeneity and the model length (observation scale) influence thermal dispersivity obtained from homogeneous model calculations (Fig. 10a,b). The latter process is well-known, at least in terms of mechanical dispersivity (e.g. Gelhar et al. 1992; Schulze-Makuch 2005), but it is much more controversial in terms of thermal dispersivity (Vandenbohede et al. 2008). Here, the results show that the primary factor of the presence of thermal dispersion is the scale of heterogeneity either ignoring or including thermal conduction. The secondary factor, which might be even dominant in some geological situations, is the observation scale. It is worth noting that the thermal dispersivity is independent of the inlet Darcy flux (injection rate) in the absence of thermal conduction (Fig. 10c).

Thermal conduction was found to weaken the effects of both the heterogeneity scale and the observation scale by suppressing the difference in heat transport between the low- and high-permeability zones. The thermal cooling in the pumping well appears later; thus, lower thermal dispersivity values are needed. Many authors state, because the heat propagates not only in the pore space but also in the matrix, that heterogeneity has a smaller effect on the heat transport than on the solute transport. Consequently, the role of the heterogeneity and also the thermal dispersion in the heat transport is negligible (e.g. Bear 1972; Vandenbohede et al. 2008). However, the results of this study indicate that the heterogeneity/dispersivity can have an important influence on heat transport, which is similar to the results of Yuan et al. (1991) and Hidalgo et al. (2009), especially in high Péclet number flow systems, and its neglect must be justified.

Previous studies (Windqvist and Hyden 1976; Sauty et al. 1982a, b; de Marsily 1986) have stated that the value of thermal dispersivity and mechanical dispersivity are similar. When neglecting thermal conduction, the results are similar to the solute transport. Windqvist and Hyden (1976), cited in Sauty et al. (1982b) pointed out that the ratios of the mechanical and thermal dispersivity values range from 0.8 and 1.8, which are all close to unity. Based on the numerical investigation in this study, the ratio of the dispersivity values (without and with thermal conduction) is higher, ranging from 2.2 to 4.1—for example, this ratio approximates 4 at R = 5 and 10 m, and 3 at R = 20 m in the base model (L = 500 m and q_in = 10^–7 m/s). In a word, the thermal dispersivity is lower than the mechanical dispersivity, as pointed by Constantz et al. (2003) and Vandenbohede et al. (2008), but the ratio is not infinite, and it depends on the model parameters.

The longitudinal thermal dispersivity α_L ranges from 10 to 50 m in the absence of thermal conduction and 3–20 m in the presence of thermal conduction. These values are mostly higher than the results of Sommer et al. (2013), where α_L is below 10 m. However, Sommer et al. (2013) investigated the relationship between the longitudinal thermal dispersivity and heterogeneity scale in a shallow ATES (aquifer thermal energy storage) system with a low temperature difference (~4 °C) and small well distances (52–156 m). In contrast, this paper focuses on a deep geothermal system with a significantly higher temperature difference (~70 °C) and larger well distances (200–2000 m), which suggests that α_L can be significantly larger in deep geothermal systems. In addition, this study uses isotropic heterogeneity scales that do not account for the effect of anisotropy. It is important to note that α_L was calculated from different physical parameters: this paper uses t_10% to calculate α_L, whereas Sommer et al. (2013) used thermal recovery. Despite these methodological differences, in the case of a small model length (L = 200 m), in the presence of thermal conduction, both studies report similar longitudinal thermal dispersivity values (α_L ≈ 5 m).

Nondimensional numbers are conventionally used for making comparisons between numerical outcomes and different model parameters. The Péclet number, which is the ratio of thermal advection and conduction (after the original definition), was modified to take into account the role of thermal dispersion. In the modified Péclet number, the thermal advection competes with the thermal conduction and dispersion, analogously to solute transport modelling problems (Simmons and Narayan 1997; Wen et al. 2018; Galsa et al. 2022),

$${\text{Pe}}_{\text{m}}=\frac{{\rho }_{\text{w}}{c}_{\text{w}}{q}_{\text{in}}\frac{{\Delta T}_{\text{x}}}{L}}{\left[\Phi {\lambda }_{\text{w}}+\left(1-\Phi \right){\lambda }_{\text{m}}\right]\left(\frac{{\Delta T}_{\text{x}}}{{L}^{2}}+\frac{{\Delta T}_{\text{y}}}{{d}^{2}}\right)+{\rho }_{\text{w}}{c}_{\text{w}}{q}_{\text{in}}R\frac{{\Delta T}_{\text{x}}}{{L}^{2}}}$$

(7)

where ΔT_y = T_bot–T_top and ΔT_x = T_stat–T_in denotes the stationary vertical and horizontal temperature difference across the model domain. Equation (7) indicates that the role of vertical advection and vertical dispersion is secondary and can be neglected; however, the effect of vertical conduction is notable as it results in the warming-up of the injected water. Furthermore, the dispersivity in heterogeneous models is scaled by R. In simulations where thermal conduction is neglected (λ_m = λ_w = 0), Eq. (7) simplifies to:

$${\text{Pe}}_{\text{m}}=\frac{L}{R}$$

(8)

Figure 12a summarizes the numerical model results studying the scale of heterogeneity, the model length and the inlet Darcy flux in the absence and presence of thermal conduction as functions of the calculated longitudinal thermal dispersivity and the modified Péclet number. As the scale of heterogeneity increases in the absence of thermal conduction (solid line), the modified Péclet number decreases (see Eq. 7), while the thermal dispersivity increases, indicating the relationship between R and α_L (Fig. 10a). In a longer model domain, both the Péclet number and the dispersivity increase, since the observation scale enhances the role of thermal dispersion (Fig. 10b). Without conduction, both Pe_m (Eq. 8) and α_L (Fig. 10c) are independent of the inlet Darcy flux, because the fluid propagates along the same path. In the presence of thermal conduction (dotted line), the trends are similar in the models with varying R and L, but conduction reduces the values of Pe_m and α_L by weakening the relative effect of dispersion and advection. At low Darcy flux, both Pe_m and α_L tend toward zero (advection and dispersion become negligible), while at high Darcy flux, Pe_m and α_L tend toward the base model without conduction (conduction becomes negligible). Figure 12b shows the same solutions, but the longitudinal thermal dispersivity is scaled by the scale of heterogeneity. The graph suggests that a quasi-logarithmic relation appears between α_L/R and Pe_m without conduction (Pe_m axis is logarithmic), where the ratio of α_L/R ratio ranges from 1 to 4, depending on R and L. In the presence of thermal conduction, the α_L/R ratio is suppressed below 1, with the exception of q_in = 10^–6 m/s, where the role of thermal conduction is marginal (Fig. 7c).

Fig. 12

Stationary numerical solutions from heterogeneous model calculations. a The longitudinal thermal dispersivity and b the longitudinal thermal dispersivity scaled by the heterogeneity scale plotted against the modified Péclet number (Eqs. 7 and 8). The results from the parameter test of the scale of heterogeneity (blue), the model length (green) and the inlet Darcy flux (red) are represented in the absence (solid line) and presence (dotted line) of thermal conduction. Base model is characterized by R = 10 m, L = 500 m and q_in = 10^–7 m/s without thermal conduction

It is well-known, particularly in solute transport modeling, that the solution is sensitive only regarding the ratio of L/R, and not their individual values (Dagan 1989; Bianchi and Pedretti 2017). In the simulations of this study, ignoring thermal conduction, L/R = 50 works for both model M1 (L = 500 m, R = 10 m, q_in = 10^–7 m/s) and M2 (L = 1000 m, R = 20 m, q_in = 10^–7 m/s). The nondimensional t_10% times for M1 and M2 are t_10% = 0.397±0.026 and 0.402±0.027, respectively, supporting the general assumption. In contrast, when thermal conduction is included in the heat transport simulation, the nondimensional t_10% = 0.51±0.026 and 0.561±0.021 for models M1 and M2, respectively. It is logical, because the role of conduction increases with the model length, but it is independent of the scale of heterogeneity. Moreover, the thermal dispersivity cannot be the same for fixed values of L/R, since α_L increases with both the scale of heterogeneity and the model length (Fig. 10a,b)—for example, the thermal dispersivity is α_L = 29.66±8.08 m for M1, while α_L = 53.89±13.62 m for M2, ignoring thermal conduction. This clearly draws attention to the fact that in the presence of heat conduction, the possibility of rescaling between the mass and heat transport problems is limited, in particular due to the different boundary conditions (see ESM).

To determine the relationship between permeability heterogeneity and thermal dispersivity, only one characteristic point, the time needed for 10% cooling, was used instead of developing an inversion method for all the BTCs. However, quantitative analysis showed that the relative deviation between BTCs is small, with ε ranging from 0.48 to 11.44% with a mean of 4.11% for all simulations. It is noted that better model fitting can be achieved when transversal thermal dispersivity is treated as a free parameter rather than fixed as α_T = α_L/10.

Finally, it is worth noting that the effective permeability calculated from the injection rate and the pressure difference is less than the geometric mean of permeability in the heterogeneous aquifer; hence, the average value does not necessarily characterize the permeability of the medium correctly; it has purely mechanical reason, as pointed out in the Appendix. Additionally, the deviation depends on both the scale of heterogeneity and the model length. Here, in this model set-up, the effect of R was stronger. It would be perspective to redefine the average permeability of heterogeneous porous medium which is independent of the scale of heterogeneity (R).

The numerical models developed in this study, similarly to each model, apply simplifications of the real hydrogeological environment. Two-dimensional approximation entails that the effect of Darcy flux on the dispersivity is underestimated in the surroundings of the well/inlet/outlet, while it is overestimated in the middle zones of the aquifer. During the synthetic simulation, the heterogeneity of permeability was focused on as the primary factor resulting in dispersion; however, other parameters, such as thermal conductivity, specific heat, and porosity, were assumed to be homogeneous. Nevertheless, the variance of these parameters is much less, so these could affect the value of thermal dispersivity to a lower degree. Permeability anisotropy (e.g. Willems et al. 2017; Poulet et al. 2023) was neglected in this simple model to concentrate on the effects of the scale of heterogeneity and observation as well as the injected volume rate as primary factors. Finally, the temperature- (and pressure) dependence of the water density and viscosity was also ignored. Based on preliminary numerical calculations, these effects are marginal at typical Darcy flux values; however, the appearance of free thermal convection in heterogeneous porous media should be an active research topic in the near future. Of course, thermal conduction is always present in a real hydrogeological environment. Turning it off in the model primarily serves a didactic purpose in order to understand and demonstrate basic physical phenomena and to emphasize the effect of thermal conduction.

The results draw attention to another practical issue. In geothermal exploration, when the heterogeneity of the reservoir can be mapped, e.g. by using seismic attributes (Mahrez et al. 2023), and permeability data can be derived from well logging, the scale of heterogeneity can be estimated. By using geostatistical methods, the value of thermal dispersivity characterizing a geothermal reservoir and operation system can be calculated. The applied SGS can be a useful tool for characterizing sedimentary systems with alternating sand and clay lenses such as the Újfalu Formation in Hungary. The results of this paper can provide guidance towards a preliminary estimation of the value of longitudinal thermal dispersivity for numerical modelers, which is relevant, since properly scaled thermal dispersivity can be used to more accurately predict the thermal breakthrough times without the need for a large number of complex heterogeneous model calculations.

The number of heat transport modelling studies is increasing for both shallow and deep geothermal systems. Several studies have focused on shallow ATES systems, assuming homogeneous (Sommer et al. 2015; Bloemendal and Olsthoorn 2018; Stemmle et al. 2023) and heterogeneous (Sommer et al. 2013; Bridger and Allen 2014) reservoirs. It is important to note that the behaviour and modelling of deep geothermal systems differ in many aspects from shallow ones (Babaei and Nick 2019; Liu et al. 2020). Nevertheless, understanding the role of heterogeneity is important for the sustainable operation of geothermal systems and has become a significant focus of current research (Zhao et al. 2021; Wang et al. 2022, 2023, etc.).

Summary and conclusions

Based on the synthetic simulations of a closed heterogeneous aquifer (imitating an injection-pumping well doublet), the effects of the permeability heterogeneity scale (R), model length (L) and the inlet Darcy flux (q_in) were investigated regarding the control parameters, such as the outlet temperature (T_out), the time for 10% cooling (t_10%), the stationary outlet temperature (T_stat), the effective permeability (k_eff) and the longitudinal thermal dispersivity (α_L). In each model parameter set-up, ten heterogeneous permeability realizations were created to characterize the stochastic medium. Numerical modeling of the time-dependent flow and heat transport in heterogeneous porous media was accomplished by both neglecting and including the effect of thermal conduction. The main conclusions of this study include:

• 1. As the scale of heterogeneity increases, injected cold water is preferentially channeled through high permeability zones, reducing the thermal breakthrough time (t_10%) until the scale of heterogeneity becomes comparable to the model size. However, the complete cooling of the pumped water requires more time. This phenomenon can be described by thermal dispersion.

• 2. A numerical relationship between the parameters of the heterogeneous model (scale of heterogeneity, model length and inlet Darcy flux) and the thermal dispersivity of the homogeneous model was established through the thermal breakthrough time. Thus, the effect of heterogeneity is well imitated with an appropriately scaled thermal dispersivity. The scale of heterogeneity was linked to the value of longitudinal thermal dispersivity. Although only one characteristic point of the BTC (t_10%) was used to determine the relationship, the results demonstrated a good fit, with a mean relative deviation of ε ≈ 4.1%.

• 3. The results affirm that the value of longitudinal thermal dispersivity depends on both the scale of heterogeneity (R) and the model length (L) (i.e. observation scale), but the dominant factor is the scale of heterogeneity.

• 4. These simulations highlight that thermal conduction moderates the role of dispersion/heterogeneity, which is an important difference between heat transport and mass transport processes. Still, the thermal dispersion cannot be ignored in accurate heat transport calculations of porous media, especially in deep geothermal exploration.

• 5. The effective permeability (k_eff) of the heterogeneous media is less than the geometric mean of the permeability in all cases, and it depends on the scale of heterogeneity.

Active research is being conducted on the use of temperature-dependent water density to explore the role of thermal buoyancy in groundwater flow and heat transport processes, particularly its impact on thermal dispersivity. Another goal is to investigate the effects of heterogeneity at the basin scale. A key question remains regarding how the basin scale influences the calculated longitudinal dispersivity values.

Appendix

The effective permeability (Eq. 5) decreases as the scale of heterogeneity increases. To illustrate this behaviour, a circular obstacle with a permeability of zero was positioned into a homogeneous medium with k = 10^–12 m² (Fig. 13). The radius of the obstacle, r, was increased, and the effective permeability was calculated from the inlet Darcy flux, q_in = 10^–7 m/s and the pressure difference between the two model sides. It was found that the pressure behind the barrier decreases (Fig. 14a); thus, the pressure difference increases and, hence, the effective permeability decreases (Fig. 14b). The pressure decrease/difference is greater at the larger radius of the barrier. Zones of low permeability can be interpreted as barriers of the flow in a heterogeneous medium.

Fig. 13

Pressure around a circular obstacle with a radius of a 4 m, b 10 m, c 20 m and d 40 m. Streamlines are illustrated by black lines

Fig. 14

a The pressure along the x-axis through the center of the model for different obstacle radii. b The calculated effective permeability as a function of the obstacle radius