On the Limitations and Implications of Modeling Heat Transport in Porous Aquifers by Assuming Local Thermal Equilibrium

Heat transport in natural porous media, such as aquifers or streambeds, is generally modeled assuming local thermal equilibrium (LTE) between the fluid and solid phases. Yet, the mathematical and hydrogeological conditions and implications of this simplification have not been fully established for natural porous media. To quantify the occurrence and effects of local thermal disequilibrium during heat transport, we systematically compared thermal breakthrough curves from a LTE with those calculated using a local thermal nonequilibrium (LTNE) model, explicitly allowing for different temperatures in the fluid and solid phases. For the LTNE model, we developed a new correlation for the heat transfer coefficient representative of the conditions in natural porous aquifers using six published experimental results. By conducting an extensive parameter study (>50,000 simulations), we show that LTNE effects do not occur for grain sizes smaller than 7 mm or for groundwater flow velocities that are slower than 1.6 m day. The limits of LTE are likely exceeded in gravel aquifers or in the vicinity of pumped bores. For such aquifers, the use of a LTE model can lead to an underestimation of the effective thermal dispersion by a factor of up to 30 or higher, while the advective thermal velocity remains unaffected for most conditions. Based on a regression analysis of the simulation results, we provide a criterion which can be used to determine if LTNE effects are expected for particular conditions.

Generally, local thermal equilibrium (LTE) is assumed when modeling heat transport. LTE is defined as an instant exchange of thermal energy between the solid and the fluid phases of a porous medium, an assumption which neglects potential temperature differences (Quintard et al., 1997). The LTE assumption allows for merging of the two separate energy equations describing temperature in the fluid and solid phases into one, simplifying the modeling procedure by volumetric averaging (Whitaker, 1991). Previous work about the validity of the LTE assumption in porous media focused on engineering applications such as packed bed reactors (Al-Nimr & Abu-Hijleh, 2002;Al-Sumaily et al., 2013;Amiri & Vafai, 1998;Khashan & Al-Nimr, 2005). While several criteria exist to examine the appropriateness of the LTE assumption for such applications (Amiri & Vafai, 1998;Hamidi et al., 2019;Kim & Jang, 2002;Minkowycz et al., 1999;Zanoni et al., 2017;Zhang et al., 2009), a thorough investigation of the conditions under which LTE is valid for flow in porous aquifers is currently lacking. In fact, adopting the criteria developed for engineering purposes is not straightforward, as the experimental conditions and reference parameters are not representative for those typical of flow in porous aquifers. For example, the required "characteristic length" in technical applications is commonly defined as channel length or width and cannot easily be determined for heat transport in subsurface environments.
Some of the few studies focusing on geotechnical applications showed that the LTE assumption can fail in geothermal systems hosted in fractured rocks Shaik et al., 2011), in partly saturated systems, for example, during infiltration of rain or melt water in frozen soil (Heinze & Blöcher, 2019), in fractured rocks as heterogeneous porous media with large permeability contrasts (Hamidi et al., 2019), and under streambed conditions with very low Reynolds numbers . In these cases, a local thermal nonequilibrium (LTNE) model may be applied. LTNE models apply two separate energy equations to describe the temperature in the fluid and solid phases (Sözen & Vafai, 1990) that are coupled through a heat transfer term which is dependent on the heat transfer coefficient and the heat transfer area. The heat transfer coefficient is a fundamental physical parameter whose value is determined experimentally (Kaviany, 1995). Again, available experimental work stems from the engineering field, and the derived parameter values and correlations are valid only under certain conditions (Singhal et al., 2017a(Singhal et al., , 2017bSun et al., 2015;Tavassoli et al., 2013Tavassoli et al., , 2015Zhu et al., 2019). For instance, these correlations are mostly derived from air or gas as the fluid, using high Reynolds numbers and wider porosity ranges that are commonly found in porous aquifers (Gunn, 1978;Tavassoli et al., 2015;Zhu et al., 2019). In addition, most of these models assume perfect spheres with a single size (Sun et al., 2015) that cannot be assumed for a heterogeneous and anisotropic aquifer. Furthermore, it is difficult to identify and quantify the errors induced by using an LTE model with the currently available LTNE measures. Standard criteria such as the temperature difference between fluid and solid phases (Abdedou & Bouhadef, 2015;Al-Sumaily et al., 2013;Khashan & Al-Nimr, 2005;Khashan et al., 2006) or the difference between the local fluid temperatures of both models (Hamidi et al., 2019) have a technical meaning, and the transition to a natural environment is difficult. The implications for the heat transport modeling in the field are mostly unknown.
In this study we examine the conditions of forced convective heat transport in natural porous aquifers under which the assumption of LTE applies and determine when it is expected to fail. To achieve this, we systematically compare thermal breakthrough curves (BTCs) calculated using LTE and LTNE models for a range of parameter values that are typical of realistic aquifer conditions. We further develop a new Nusselt number (Nu) correlation to determine the heat transfer coefficient representing the conditions expected in porous aquifers using all available data sets. Finally, we present discrete limits for the conditions under which the commonly used LTE model is applicable when heat transport in porous aquifers is quantified.

Method
To investigate the occurrence and impact of LTNE effects on heat transport in porous media, a parameter study with one-dimensional (1-D) numerical LTNE and LTE models is conducted. The notations are given in Appendix A. In a first step, a suitable Nu correlation for porous aquifers to determine the heat transfer coefficient between fluid and solid is derived based on appropriate literature data. Then a numerical model with a step input is used to create thermal BTCs using the LTE and LTNE models. In the last step, the resulting thermal BTCs are evaluated with different methods to quantify LTNE effects. The dominant LTNE parameters are identified through a global sensitivity analysis. More than 50,000 modeling runs are conducted to cover the whole range of possible parameter combinations representative of natural porous aquifers. To evaluate a criterion which allows to easily assess whether the occurrence of LTNE effects is likely, a regression analysis with the sensitive parameters as explanatory variables is conducted.

Heat Transfer Coefficient
The heat transfer coefficient is a vital parameter for LTNE models, as it describes the heat transfer rate between the fluid and solid phases. As seen above, h sf depends on Nu, the particle size d p , and the thermal conductivity of the fluid λ f 10.1029/2020WR027772

Water Resources Research
GOSSLER ET AL. (Wakao et al., 1979). While λ f and d p are usually known in ideal packed beds, Nu must be estimated using an appropriate correlation which commonly depends on the particle Reynolds number, Re, and the Prandtl number, Pr (see supporting information Table S1 for examples). Generally, such a correlation is empirically derived from laboratory experiments (Achenbach, 1995;Collier et al., 2004;Naghash et al., 2016;Nie et al., 2011;Shent et al., 1981;Wakao & Kaguei, 1982;Zanoni et al., 2017) or, more recently, by particle resolved direct numerical modeling (Chen & Müller, 2019;Singhal et al., 2017b;Sun et al., 2015;Tavassoli et al., 2015;Zhu et al., 2019). As natural porous aquifers are typified by Re < 50 and a porosity (fluid volume fraction) n < 0.5, the correlations established by previous studies are covering much broader ranges. Table S1 shows an overview of these correlations and demonstrates that no single model appropriately describes the conditions expected in a porous aquifer. Therefore, a new regression using only the experimental data representative for conditions in porous aquifers (Re < 50 and n < 0.5) is carried out. We found only one study (Kunii & Smith, 1961) inspecting Nu values for very low Re (Re < 0.1), and their study was disregarded because the mathematical model was criticized by several authors (Gunn & De Souza, 1974;Littman et al., 1968;Wakao et al., 1979). Because the raw data of all of the published experiments were determined with a Pr of 0.7-1 (gas), we corrected the Nu values based on the correlations of Wakao et al. (1979) and Zhu et al. (2019) (Table S1) using the Pr value of 9 (water, 10°C) for each Nu of the respective Re value z: Nu orig, Re = z gives the originally measured Nu value at a Re of z. To correct for the higher Pr of water, the amount of increase of the Nu due to the higher Pr is added. This amount (Nu Pr = 9, Re = z − Nu Pr = 0.7, Re = z ) is calculated with the Nu correlation of Wakao et al. (1979) and Zhu et al. (2019) (section S1.1 and Table S1).
Theoretical considerations for heat transfer between a sphere and a fluid flowing around it show that Nu must have a minimum value of 2 (Ranz & Marshall, 1952;Shent et al., 1981). As argued, for example, by Nelson and Galloway (1975), lower Nu numbers may occur in a packed bed of spheres. In our work, due to the low range of Pr expected in groundwater conditions (Pr 1/3~2 .08 at 10°C and Pr 1/3~1 .94 at 20°C), the common dependency of Nu on Pr 1/3 is not resolved. Therefore, following the common dependency of Nu on Re with a lower limit (see, e.g., Table S1), a correlation of is fitted to all available data using a and b as free parameters.
In addition to the parameter study of the new correlation, a parameter study is conducted for five (Achenbach, 1995;Singhal et al., 2017b;Wakao et al., 1979;Zanoni et al., 2017;Zhu et al., 2019) other Nu correlations, to evaluate the influence of the Nu correlation choice.

LTE and LTNE Models
The heat transport of the fluid and solid phases can be described as a continuous single average temperature field (Equation 4) (Nield & Bejan, 2017) assuming that the temperature difference between the fluid and solid phases within the representative elementary volume is negligible (LTE assumption). The equations are as follows: The first part on the right side of Equation 4 is the effective thermal dispersion coefficient D l,eff (Equation 5) consisting of the thermal diffusivity and thermal mechanical dispersion coefficient D l (Equation 6) (Rau et al., 2012a). The second part describes the advective heat transport (Equation 7).
When the LTE assumption does not fit, then heat transport can be described with a dispersion particle-based two equation model (LTNE assumption; Equations 8 and 9) (Amiri & Vafai, 1994). The two temperatures are interlinked by a coupling term (Wakao et al., 1979) based on the heat transfer coefficient (Equation 1) and the specific surface area, here assumed for a bed of spherical particles a sf = (6(1 − n))/d p (Dullien, 1979). The equations are as follows: In the simulations conducted as part of our work, the LTE (Equation 4) and LTNE (Equations 8 and 9) models with the following initial (Equation 10a) and boundary conditions (Equation 10b) were used:

Numerical Solution and Parameter Study
The LTE and LTNE models explained above were solved using the MATLAB function "pdepe" which utilizes a finite-element, piecewise nonlinear Galerkin/Petrov-Galerkin method with second-order accuracy in space (Skeel & Berzins, 1990). To avoid boundary effects, the model domain was setup with 15 m length. The spatial discretization was 0.0025 m for distances less than 4.5 m (1.5 · maximum distance of interest = 3 m) and increased along the remaining distance. The fine numerical discretization and the chosen input parameter range lead to simulations in which the nodal distance was smaller than the used grain size. While the representative elementary volume is usually larger than the mean grain size in the volume averaging approach (Rau et al., 2014), we do not consider this problematic, as we consider homogenous parameter conditions within the whole model domain. This computationally demanding configuration was chosen as conservative setup to assure mesh independency, but a sensitivity analysis revealed that the results are still acceptable at spatial discretization up to 0.01 m and when reducing the model domain length to 5 m. The time discretization is automatically adapted. A Linux Cluster (at the Leibniz Supercomputing Centre of the Bavarian Academy of Science and Humanities) was utilized to run the simulations in parallel. The numerical solution was validated with analytical solutions (van Genuchten & Alves, 1982;Wakao & Kaguei, 1982) and illustrates good agreement for the LTE and LTNE models (see section S2 for further information). Table 1 shows the range of values used in our investigation. The parameters are independent of each other, for example, the porosity does not constrain the seepage velocity and so forth. Dependent fluid properties, such as the thermal conductivity λ f , the density ρ f , and the specific heat capacity c f , were calculated from the fluid temperature T f (Furbo, 2015). Nu was computed through the new correlation given by Equation 3 (Figure 2b). The resulting Re ranges between 0.015 and 70 with 95% of the simulations within the range of 1 to 50.
To investigate the entire parameter space, first, a uniformly distributed random sample for each independent parameter (Table 1) was created. To decrease the computational cost of the sensitivity analysis, the random parameter sets were extended using Saltelli's (2002) scheme. Then, the dependent fluid parameters were calculated. The number of parameter sets was varied from 48 to 32,000 resulting in 50,608 different parameter  (Banks, 2012) Parameter Range Thermal conductivity solid (λ s ) Vol. heat capacity solid (ρ s c s ) 1.5 * 10 6 -3.5 * 10 6 J m

10.1029/2020WR027772
Water Resources Research sets for each Nu correlation. Each parameter set was used in a model run where a unit step boundary condition was applied to create thermal BTCs for the LTE and LTNE models. The thermal BTCs were analyzed for the distances of 0.5, 1, 2, and 3 m.

Quantification of LTNE
Three methods were used to quantify the magnitude of the effect of LTNE. Method 1 was newly developed in this study and compared to two existing methods (Methods 2 and 3): Method 1: Comparison of the difference of effective thermal dispersion and advective thermal velocity between the LTE and LTNE models via the following approach: 1. Simulation of fluid temperature BTC with the LTNE model; 2. Estimation of the effective thermal dispersion coefficient D l, eff, fit and the advective thermal velocity v t, fit from the LTNE calculation (previous step) by fitting the analytical LTE model (van Genuchten & Alves, 1982) (Equation S1) using a standard nonlinear least square fitting procedure; 3. Comparison of the estimated D l, eff, fit and v t, fit from the LTE model with the used D l, eff (Equation 5) and v t, LTE (Equation 7) in Step 1.
Method 2: Maximum and average temperature difference between the normalized temperature obtained using the LTE model and the fluid temperature obtained using the LTNE model (Hamidi et al., 2019) for each evaluated distance x (N = number of timesteps). and Method 3: Maximum and average normalized temperature difference between fluid and solid (Abdedou & Bouhadef, 2015;Al-Sumaily et al., 2013;Khashan et al., 2006;Khashan & Al-Nimr, 2005) obtained using the LTNE model for each distance x: and

Global Parameter Sensitivity Analysis
To identify the dominant LTNE parameters a global parameter sensitivity was calculated using the function "sobolSalt" as part of the R package "sensitivity" (Iooss et al., 2019). This implements the Monte Carlo estimation of the variance-based sensitivity indices (Sobol indices) of the first and total order for each independent parameter.

New Nusselt Correlation
The available Nu values found in literature with conditions expected in natural porous aquifers (Re < 50, n < 0.5) are shown in Figure 1. The original Nu values determined with gas as fluid (Pr = 0.7) were adapted for water as fluid (Pr = 9) using the Nu correlations and Equation 2 (see section S1 for further information).

Water Resources Research
GOSSLER ET AL.

result in very low
Nu numbers for the investigated Re range (Figure 2a). A reason for the discrepancy in Nu values between some publications (Naghash et al., 2016;Nelson & Galloway, 1975;Nie et al., 2011;Zanoni et al., 2017) (first group) and other studies (Achenbach, 1995;Singhal et al., 2017b;Sun et al., 2015;Wakao et al., 1979;Zhu et al., 2019) (second group) is not obvious, and therefore, the unrealistic data for porous aquifers from the first group were not included in our regression. The new regression based on the available data and Equation 3 is shown in Figure 2. This correlation with best fit values a = 3.1 and b = 0.57 is used in all of the presented analysis.

Comparison of Thermal BTCs Obtained From the LTE and LTNE Models
Significant differences between the modeled thermal BTCs calculated using the LTE and LTNE models can be found within the assumed parameter ranges representative of the conditions in porous aquifers. As an example, Figure 3 highlights the thermal BTCs obtained from the LTE and LTNE models for particle sizes (<0.5 cm) with slow flow velocities (<2 m day −1 ) at the lower end of the investigated parameter range (Figures 3a-3c) and for large particles sizes (>7.5 cm) with high flow velocities (>20 m day −1 ) (Figures 3d-3f). Furthermore, the three different methods used to evaluate the degree of LTNE are  illustrated. For low flow velocities and small particle sizes, the LTE and LTNE models result in similar (visually the same) thermal BTCs regarding effective dispersion and advective thermal velocity (Method 1).
The computed normalized temperature differences between the LTE and LTNE fluid temperature (Method 2) and LTNE fluid and solid temperature (Method 3) are very low (<0.005). In contrast, in the simulations for large particle sizes and high flow velocities, the thermal BTCs from the two models differ significantly. The solid and fluid temperature of the LTNE model is noticeably more dispersed (Method 1), and significant temperature differences between the LTE temperature and LTNE fluid temperature (Method 2) as well as solid and fluid temperature (Method 3) can be observed.

Influence of the Nusselt Correlation
To investigate the influence of the Nu correlation choice, the parameter study was extended to the Nu correlations suggested by Zanoni et al. (2017) as well as four others (Achenbach, 1995;Singhal et al., 2017b;Wakao et al., 1979;Zhu et al., 2019 Zanoni et al. (2017) shows no significant retardation compared to the fluid front. By contrast, the solid phase temperature increases much slower due to the very low heat exchange between the fluid and solid phases as caused by the low Nu values from this correlation (see Table S2 for values of Nu). Both fluid and solid BTCs obtained from the LTNE model for the Nu correlation of Zanoni et al. (2017) differ significantly from the temperature calculated using the LTE model. In fact, fitting the analytical LTE model to the fluid temperature BTC obtained from the LTNE model fails to achieve a satisfactory result (see Figure S6). All other Nu correlations lead to similar BTCs. This shows that Nu correlations leading to very low Nu numbers are not suitable for natural conditions when simply adjusting the Pr to a value appropriate for water.

Water Resources Research
In conditions with large particle sizes and high seepage velocities (Figures 4c and 4d), all BTCs of the different Nu correlations differ from the LTE model showing higher dispersion while slightly differing from each other. Therefore, the general outcome that LTNE conditions lead to higher thermal dispersion is nearly independent of the choice of the Nu correlation (see also Figure S7).

Analysis of the Parameter Sensitivity
The parameter sensitivity analysis reveals that when evaluating all three methods, the most sensitive parameters are the seepage velocity, the particle size, and the porosity. The least sensitive parameters are found to be the fluid temperature, the volumetric heat capacity of the solid, and the thermal dispersivity. The full details of the sensitivity analysis are given in section S4. In many engineering applications, the ratio of fluid to solid phase thermal conductivities is reported as an important LTE criterion (Dehghan et al., 2014;Lee & Vafai, 1998;Minkowycz et al., 1999). The ratio between the thermal conductivity of water (Huber et al., 2012) and common aquifer materials (Côté & Konrad, 2005) is usually in the range~0.1-0.5. The differences between the thermal conductivities in engineering applications with air as fluid and metals as a solid phase span a much broader range. We believe that our analysis did not find this to cause sensitivity because of the limited range of thermal conductivities expected in aquifer settings, as was the focus in our study.

The Influence of LTNE Effects on the Advective Thermal Velocity
The modeled against the fitted thermal velocity is shown in Figure 5a, highlighting an excellent fit (R 2 > 0.99) with generally very small differences (for 95% of the simulations, the deviation is smaller

10.1029/2020WR027772
Water Resources Research than 5%). A small number of simulations at low Peclet numbers and large particle sizes deviate from the expected thermal advective velocity (Figure 5b). Aside from these conditions, the influence of LTNE effects on the advective velocity is generally very low. The fitted flow velocities exceed the modeled velocities slightly (<5%) at very high velocities. The results also illustrate that this is not flow distance dependent, as the correlations are similar for all investigated distances ( Figure S10). These observations are in agreement with a recent study (Gossler et al., 2019) which experimentally investigated possible LTNE effects. No significant influence on the advective thermal velocity could be observed within the analyzed range of seepage velocities expected in gravel aquifers (5-50 m day −1 ). In a numerical study, Roshan et al. (2014) investigated the influence of LTNE effects on velocity estimates from the damping and phase shifting of the diel temperature signal with depth in a streambed for flux conditions leading to 0.001 < Re < 0.01. They come to an opposite conclusion, stating that LTNE effects are limited to very slow flow velocities but can lead to velocity deviations up to a factor of 150. A possible explanation of these findings is their choice of the Nu correlation which leads to very low Nu numbers at low flow velocities but increases significantly at higher velocities. Furthermore, the mathematical model used to derive the Nu values (Kunii & Smith, 1961) on which the correlation of Roshan et al. (2014) is based was criticized by several authors (Gunn & De Souza, 1974;Littman et al., 1968;Wakao et al., 1979).

The Influence of LTNE Effects on Thermal Dispersion
In contrast to the advective thermal velocity, the effective thermal dispersion can be significantly influenced by LTNE effects. We evaluate that an increase by a factor of over 30 (LTNE method1; Dl ) is possible within the range of parameters investigated in this study ( Figure 6). Increasing flow velocities and particle sizes also lead to increasing LTNE effects and, consequently, to a higher effective thermal dispersion. Similar results are obtained for all investigated distances ( Figure S11). To elucidate the conditions under which this increase in thermal dispersion can be expected, the LTNE effects on thermal dispersion have been categorized in Table 2. As the thermal dispersion is generally an uncertain parameter in modeling, we considered an increase up to 50% as within the usual uncertainty range. An increase above a factor of 10 is considered as highly influenced by LTNE effects. Following this categorization and the grain Figure 5. (a) The differences between the fitted advective thermal velocity and the modeled advective thermal velocity are very small for most simulations. For 95% of the simulations, the deviation is smaller than 5%. The fitted velocity deviates at very high velocities by a small amount (<5%) from the modeled advective thermal velocity. This shows that the influence of LTNE effects on the advective thermal velocity is generally very small. (b) For a small part of the simulations at low Peclet numbers in combination with large particle sizes, the fitted velocity is smaller than the modeled thermal velocity. This deviation at low Peclet numbers, representing conduction dominated situations at small flow distances, is contributed to the differences in the boundary conditions immanent in the chosen LTE and LTNE models.

10.1029/2020WR027772
Water Resources Research size-based aquifer types according to Table 2, we reveal that LTNE effects can be expected in gravel aquifers with high Darcy fluxes (Figure 7). Furthermore, we identify the following threshold values considering all combinations of parameters for which no significant LTNE effects (LTNE method1; Dl < 1.5) could be observed: If either the particle size is <7 mm or the seepage velocity is <1.6 m day −1 , no significant LTNE effects (LTNE method1; Dl < 1.5) should be expected. Two previous laboratory studies (Bandai et al., 2017;Gossler et al., 2019) observed increasing thermal dispersion values with increasing particle sizes and suspected LTNE effects as the cause for this. Our results confirm the assumption that LTNE effects can significantly enhance the effective thermal dispersion coefficient.

Comparison of the Methods to Quantify LTNE Effects
The comparison of the results of the three methods for all simulations shows a significant relationship between them ( Figure 8). Yet, for some parameter settings (low Pe and large d p ), the differences between the LTE and fluid LTNE temperatures (LTNE method2 ) are large, while the differences between solid and fluid temperature (LTNE method3 ) are close to zero. These findings are in accordance with the results of some other studies (Hamidi et al., 2019;Rees et al., 2008) which showed that LTE and LTNE models can differ even if the temperature differences between the fluid and solid phase are nearly equal. Interestingly, while LTNE method1,Dl does not capture these LTNE effects, LTNE method1,v shows a strong correlation with LTNE method2, max (Figure 8b). These deviations are limited to conditions of low Pe numbers (Figures 5b and 8b).
A possible explanation for these discrepancies at low Pe is the differences in the boundary conditions between the LTE and LTNE models. When trying to model a scenario in which a hot or cold fluid is injected into a porous medium, the boundary condition necessary for the LTE model and the assumption inherent to the LTE model induce that under conduction-dominated conditions (low Pe), the thermal conductivity of the solid phase also contributes significantly to heat transport at the boundary (see Equations 4 and 5). The LTNE model allows to set more appropriate boundary conditions, allowing to only affect the fluid phase at the boundary.

Criterion to Estimate LTNE Conditions
To estimate if LTNE effects are likely to occur in porous aquifers, we use the four most sensitive parameters, particle size d p , porosity n, seepage velocity v a , and thermal conductivity of the solid phase λ s , as explanatory variables to derive a new equation in a two-step regression procedure. First, to find the best model formula, an exhaustive screening approach with the R package "glmulti" (Calcagno & de Mazancourt, 2010) was applied on the normalized (ordered quantile normalization; Peterson & Cavanaugh, 2019) LTNE method1; Dl value and the scaled explanatory variables (min-max normalization) (for further details, refer to

Water Resources Research
section S6). In other words, the abovementioned explanatory variables are used to predict the value of LTNE method1; Dl . The value of LTNE method1; Dl , predicted by the regression, is called LTNE cat . All possible models including pairwise interactions were considered up to a maximum number of seven terms. The selection was based on Bayesian information criterion which penalizes models with higher number of parameters to prevent overfitting (Schwarz, 1978). The model formula in Equation 17 with six terms shows the best results concerning limited complexity and adequate accuracy (adjusted R 2 = 0.89). In order to enable a practical application, the regression coefficients were determined by applying this model to the original data (not normalized and scaled) (adjusted R 2 = 0.64) resulting in Figure 7. Categorized LTNE effects based on Darcy flux and particle size. Increased thermal dispersion due to LTNE effects is mainly expected for conditions with high flow velocities and large grain sizes like gravel aquifers.
The predicted value LTNE cat (Equation 17) was compared with the associated LTNE category (Table 2) for each simulation in a cumulative distribution plot (Figure 9). The boundaries of the categories (Table 3 and Figure 9) are based on the 5% and 95% quantiles of this cumulative distribution. This can be viewed as a probability measure, for example, if LTNE cat is <1.4, it is likely that no LTNE effects occur. Ninety-five percent of the simulations which show low LTNE effects are above this value. No simulation with medium or high LTNE effects results in an LTNE cat < 1.4. If LTNE cat is in the low LTNE-medium LTNE range (3.3-4.6), both categories are equally possible, as simulation with parameter sets resulting in this range can be within both categories (Figure 9). The accuracy of this approach is reasonable, as LTNE cat gives only a qualitative assessment of the LTNE conditions.

Limitations of the Used Approach
Our study focuses on a 1-D setup with homogenous parameter settings in each simulation to investigate LTNE effects. Even though most of the parameter settings expected in porous aquifers are covered, the influence of macroscopic heterogeneity on LTNE effects remains unclear. A recent publication investigated LTE and LTNE approaches for fractured porous media (Hamidi et al., 2019). They showed that high permeability and porosity contrasts can lead to significant LTNE effects resulting in a difference of up to 7% in local fluid temperatures. This indicates that macroscopic heterogeneity could additionally increase LTNE effects.
Furthermore, we used the arithmetic mean particle size as a referential particle size in our simulations. The sensitivity analysis showed that the particle size is a highly influential parameter for LTNE effects. Similarly, Heinze and Blöcher (2019) revealed that the particle size is a crucial parameter for LTNE effects during infiltration processes. For transferring the results to aquifers with unsorted grain size distributions, further research is necessary to reveal if, for instance, a representative grain size would be suitable. Candidates are statistical grain size values such as the geometric or harmonic means as common for calculation of thermal or hydraulic conductivity.   (Table 3).

Water Resources Research
Another assumption of our study is that we consider most parameters independent of each other (see section 2.3). Certainly, some parameters like particle size and porosity can be correlated (e.g., Urumović & Urumović, 2014), but these correlations depend on many influences like grain size distribution, shape of the grains, and compaction. Therefore, we have not elaborated further on which conditions are more or less frequent in natural aquifers and rather provide the full range of parameter combinations.

Implications for Heat Transport Modeling and Interpretation
The increased thermal dispersion caused by LTNE under certain conditions can have significant consequences for the modeling and management of geothermal groundwater use, such as groundwater-based heat pump systems, which cycle groundwater by operating extraction and injection well doublets. With a steady increase of such systems leading to densely used aquifers (Pophillat et al., 2020), there is a growing need for reliable prediction of the induced thermal plumes (Böttcher et al., 2019;Epting et al., 2017;Ferguson, 2009;Hähnlein et al., 2013;Herbert et al., 2013). To delineate a thermal range or a thermally affected zone in an aquifer, the value for thermal dispersion is of crucial importance (Molina-Giraldo et al., 2011;Park et al., 2015Park et al., , 2018Pophillat et al., 2020). Consequently, the processes that determine its magnitude must be well understood. For example, we expect that LTNE effects can occur even in aquifers with low natural seepage water velocities as the flow velocity increases significantly toward injection and extraction wells (Park et al., 2015). Furthermore, the highly transient operation typical for groundwater-based heat pump systems (Muela Maya et al., 2018) can enhance the effects of LTNE (Minkowycz et al., 1999). Therefore, an increase in effective thermal dispersion due to LTNE effects should be carefully considered under the conditions depicted in Figure 7. Thermal nonequilibrium potentially also affects the storage of thermal energy in aquifers composed of large particle sizes. The increased spreading of the thermal front in the fluid phase due to an incomplete storage in the solid phase leads to an initially lower energy density in the area of interest than is expected when a LTE model is used. Finally, when using heat as a tracer, the aim is often to quantify fluxes (Irvine et al., , 2017Rosenberry et al., 2015) or use the advective part for an inversion of the hydraulic conductivity field (Somogyvári & Bayer, 2017;Somogyvári et al., 2016). In dynamic interfaces with surface-water bodies, such as when aquifers interact with rivers, groundwater can reach velocities up to 10 m day −1 and higher (e.g., Angermann et al., 2012;Cremeans et al., 2018;Rau et al., 2014). These systems are possibly influenced by LTNE effects. As the groundwater velocities for interactions with lakes are usually well below 1 m day −1 (Rosenberry et al., 2015), the LTE assumption is mostly applicable for these conditions. Only at low Pe numbers in combination with large particle sizes an LTE model can overestimate the advective thermal velocity. Even in conditions with strong LTNE effects, for some of these cases, the LTE model is appropriate (Figure 7) because the LTNE effects do not significantly influence the advective thermal velocity at advection-dominated conditions.

Conclusion
The limitations of the LTE assumption in heat transport modeling for natural porous aquifers have been investigated by an extensive 1-D parameter study comparing thermal BTCs of LTE and LTNE models. A new correlation to determine the Nusselt number based on available literature data tailored to porous aquifers was derived. As all available data on Nu are determined with air as fluid, more laboratory studies measuring the heat transfer with water as a fluid are required.
While LTNE effects do not occur for grain sizes smaller than 7 mm or for flow velocities that are slower than 1.6 m day −1 , LTNE effects can occur within the conditions expected in porous aquifers like gravel aquifers with high flow velocities and large grain sizes. The effective thermal dispersion can be increased by a factor of over 30 caused by LTNE effects. The advective thermal velocity is not significantly influenced even in conditions with strong LTNE effects. Only at low Pe conditions in combination with large particle sizes, the LTE model can lead to an overestimation of the advective thermal velocity. A criterion to predict if LTNE conditions will occur is developed based on the most influential properties, the porous media particle size, seepage velocity, porosity, and thermal conductivity of the solid phase. Our results can be used as a guide toward more accurate modeling of heat transport in natural porous media.