Coupled groundwater flow and heat transport simulation for estimating transient aquifer–stream exchange at the lowland River Spree (Germany)

Subsurface flow and heat transport near Freienbrink, NE Germany, was simulated in order to study groundwater–surface water exchange between a floodplains aquifer and a section of the lowland River Spree and an adjacent oxbow. Groundwater exfiltration was the dominant process, and only fast surface water level rises resulted in temporary infiltration into the aquifer. The main groundwater flow paths are identified based on a 3D groundwater flow model. To estimate mass fluxes across the aquifer–surface water interfaces, a 2D flow and heat transport modelling approach along a transect of 12 piezometers was performed. Results of steady‐state and transient water level simulations show an overall high accuracy with a Spearman coefficient ρ = 0.9996 and root mean square error (RMSE) = 0.008 m. Based on small groundwater flow velocities of about 10^?7 to 10^?6 ms^?1, mean groundwater exfiltration rates of 233 l m^?2 d^?1 are calculated. Short periods of surface water infiltration into the aquifer do not exceed 10 days, and the infiltration rates are in the same range. The heat transport was modelled with slightly less accuracy (ρ = 0.8359 and RMSE = 0.34 °C). In contrast to the predominant groundwater exfiltration, surface water temperatures determine the calculated temperatures in the upper aquifer below both surface water bodies down to 10 m during the whole simulation period. These findings emphasize prevailing of heat conduction over advection in the upper aquifer zones, which seems to be typical for lowland streams with sandy aquifer materials and low hydraulic gradients. Moreover, this study shows the potential of coupled numerical flow and heat transport modelling to understand groundwater–surface water exchange processes in detail. Copyright © 2013 John Wiley & Sons, Ltd.     

Relative Recovery of Thermal Energy and Fresh Water in Aquifer Storage and Recovery Systems

This paper explores the relationship between thermal energy and fresh water recoveries from an aquifer storage recovery (ASR) well in a brackish confined aquifer. It reveals the spatial and temporal distributions of temperature and conservative solutes between injected and recovered water. The evaluation is based on a review of processes affecting heat and solute transport in a homogeneous aquifer. In this simplified analysis, it is assumed that the aquifer is sufficiently anisotropic to inhibit density‐affected flow, flow is axisymmetric, and the analysis is limited to a single ASR cycle. Results show that the radial extent of fresh water at the end of injection is greater than that of the temperature change due to the heating or cooling of the geological matrix as well as the interstitial water. While solutes progress only marginally into low permeability aquitards by diffusion, conduction of heat into aquitards above and below is more substantial. Consequently, the heat recovery is less than the solute recovery when the volume of the recovered water is lower than the injection volume. When the full volume of injected water is recovered the temperature mixing ratio divided by the solute mixing ratio for recovered water ranges from 0.95 to 0.6 for ratios of maximum plume radius to aquifer thickness of 0.6 to 4.6. This work is intended to assist conceptual design for dual use of ASR for conjunctive storage of water and thermal energy to maximize the potential benefits.     

Heat transport in the Red Lake Bog,Glacial Lake Agassiz Peatlands

We report the results of an investigation on the processes controlling heat transport in peat under a large bog in the Glacial Lake Agassiz Peatlands. For 2 years, starting in July 1998, we recorded temperature at 12 depth intervals from 0 to 400 cm within a vertical peat profile at the crest of the bog at sub‐daily intervals. We also recorded air temperature 1 m above the peat surface. We calculate a peat thermal conductivity of 0·5 W m^?1 °C^?1 and model vertical heat transport through the peat using the SUTRA model. The model was calibrated to the first year of data, and then evaluated against the second year of collected heat data. The model results suggest that advective pore‐water flow is not necessary to transport heat within the peat profile and most of the heat is transferred by thermal conduction alone in these waterlogged soils. In the spring season, a zero‐curtain effect controls the transport of heat through shallow depths of the peat. Changes in local climate and the resulting changes in thermal transport still may cause non‐linear feedbacks in methane emissions related to the generation of methane deeper within the peat profile as regional temperatures increase. Copyright © 2006 John Wiley & Sons, Ltd.     

Numerical Analysis of Thermal Remediation in 3D Field‐Scale Fractured Geologic Media

Thermal methods are promising for remediating fractured geologic media contaminated with volatile organic compounds, and the success of this process depends on the coupled heat transfer, multiphase flow, and thermodynamics. This study analyzed field‐scale removal of trichloroethylene (TCE) and heat transfer behavior in boiling fractured geologic media using the multiple interacting continua method. This method can resolve local gradients in the matrix and is less computationally demanding than alternative methods like discrete fracture‐matrix models. A 2D axisymmetric model was used to simulate a single element of symmetry in a repeated pattern of extraction wells inside a large heated zone and evaluate effects of parameter sensitivity on contaminant recovery. The results showed that the removal of TCE increased with matrix permeability, and the removal rate was more sensitive to matrix permeability than any other parameter. Increasing fracture density promoted TCE removal, especially when the matrix permeability was low (e.g., <10^?17 m²). A 3D model was used to simulate an entire treatment zone and the surrounding groundwater in fractured material, with the interaction between them being considered. Boiling was initiated in the center of the upper part of the heated region and expanded toward the boundaries. This boiling process resulted in a large increase in the TCE removal rate and spread of TCE to the vadose zone and the peripheries of the heated zone. The incorporation of extraction wells helped control the contaminant from migrating to far regions. After 22 d, more than 99.3% of TCE mass was recovered in the simulation.     

Estimation of vertical water and heat fluxes in the semi-confined aquifers in tokyo metropolitan area,japan

Groundwater circulation is known to be one of the agents responsible for the redistribution of geothermal energy by acting as a source or sink in the course of its movement through porous media. Heat transport in groundwater systems is considered to be a coupled process and the theory based on this was used to analyse temperature profiles of 30 thermally stable observation wells in a deep, semi-confined aquifer system in the Tokyo Metropolitan area. Vertical water fluxes in the semi-confined aquifers and the associated upward heat fluxes were estimated from a heat flux equation that describes convection and conduction processes of heat transport in one dimension. The vertical downward water fluxes in Shitamachi lowland, Musashino and Tachikawa terraces were 0.69.26.91 × 10^?9, 1.46-70.92 × 10^?9 and 2.61.2204 × 10^?9 m/s, respectively. A vertical upward water flux of 1.80-33.60 × 10^?9 m/s was estimated in Shitamachi lowland. The water flux generally decreased with increasing depth for observation wells which intercepted more than one semi-confining layer. The estimated upward heat fluxes for Shitamachi lowland, Musashino and Tachikawa terraces were 0.32-1.12, 0.49-1.21 and 1.00-11.62 W/m², respectively. The heat flux was highest in Tachikawa terrace where a major fault, the Tachikawa fault, is located. Generally, the estimated heat flux was higher in the semi-confining layers than in the aquifers. Areas with heat sources and sinks as well as groundwater flow patterns in the semi-confined aquifers were revealed by heat flux and temperature distributions in the study area.     

Modeling axisymmetric flow and transport

Unmodified versions of common computer programs such as MODFLOW, MT3DMS, and SEAWAT that use Cartesian geometry can accurately simulate axially symmetric ground water flow and solute transport. Axisymmetric flow and transport are simulated by adjusting several input parameters to account for the increase in flow area with radial distance from the injection or extraction well. Logarithmic weighting of interblock transmissivity, a standard option in MODFLOW, can be used for axisymmetric models to represent the linear change in hydraulic conductance within a single finite-difference cell. Results from three test problems (ground water extraction, an aquifer push-pull test, and upconing of saline water into an extraction well) show good agreement with analytical solutions or with results from other numerical models designed specifically to simulate the axisymmetric geometry. Axisymmetric models are not commonly used but can offer an efficient alternative to full three-dimensional models, provided the assumption of axial symmetry can be justified. For the upconing problem, the axisymmetric model was more than 1000 times faster than an equivalent three-dimensional model. Computational gains with the axisymmetric models may be useful for quickly determining appropriate levels of grid resolution for three-dimensional models and for estimating aquifer parameters from field tests.     

Groundwater temperature evolution in the subsurface urban heat island of Cologne,Germany

Long‐term heating of shallow urban aquifers is observed worldwide. Our measurements in the city of Cologne, Germany revealed that the groundwater temperatures found in the city centre are more than 5 K higher than the undisturbed background. To explore the role of groundwater flow for the development of subsurface urban heat islands, a numerical flow and heat transport model is set up, which describes the hydraulic conditions of Cologne and simulates the transient evolution of thermal anomalies in the urban ground. A main focus is on the influence of horizontal groundwater flow, groundwater recharge and trends in local ground warming. To examine heat transport in groundwater, a scenario consisting of a local hot spot with a length of 1 km of long‐term ground heating was set up in the centre of the city. Groundwater temperature‐depth profiles at upstream, central and downstream locations of this hot spot are inspected. The simulation results indicate that the main thermal transport mechanisms are long‐term vertical conductive heat input, horizontal advection and transverse dispersion. Groundwater recharge rates in the city are low (<100 mm a^?1) and thus do not significantly contribute to heat transport into the urban aquifer. With groundwater flow, local vertical temperature profiles become very complex and are hard to interpret, if local flow conditions and heat sources are not thoroughly known. Copyright © 2014 John Wiley & Sons, Ltd.     

Determining thermal diffusivity using near‐surface periodic temperature variations and its implications for tracing groundwater movement at the eastern margin of the Tibetan Plateau

The thermal diffusivity is the key parameter that controls near‐surface temperature where periodic temperature variation is progressively attenuated and delayed with depth. This article presents the results of apparent thermal diffusivity using temperatures recorded by a bedrock temperature measurement network in the fault zones of western Sichuan. High sensitivity temperature sensors (10^?4 K) were installed at a maximum depth reaching 30 m. The apparent thermal diffusivities were deduced from both amplitude damping and phase shifting of annual temperature variations between two different depths. Under pure conduction, the thermal diffusivity determined through the phase method (α_Φ) should be equivalent to that determined through the amplitude method (α_A), whereas effects of the upward (downward) water flow are evidently reflected in the amplitude decay to make α_Φ larger (lesser) than α_A. The discrepancy between α_Φ and α_A can thus be a tracer of water movement or convective heat transfer. The calculated α_Φ of the measurement stations varies from 1.22 × 10^?6 to 3.00 × 10^?6 m²/s, and the estimated α_A ranges from 0.93 × 10^?6 to 2.41 × 10^?6 m²/s. Two regimes of heat transfer underground were suggested from the results. Conductive heat transport prevails over the nonconductive processes at five stations, which is characterized by α_Φ coincident with α_A for the same depth pair. On the contrary, the values of α_Φ differ from α_A at six stations in the intersection area of the Y‐shaped fault system, implying that convective heat transfer also plays a comparably important role. This finding is consistent with the hot springs distribution of the area. The results also indicate that water moves upward with an average Darcy velocity of approximately ?1 × 10^?7 m/s in this region. Our research provides new evidence for the hydrothermal activity in the fault zones at the eastern margin of the Tibetan Plateau.     

The dynamical and thermal structure of deep mantle plumes

If the interpretation of the D″ layer at the base of the mantle as a thermal boundary layer, with a temperature increment in the order of 800 K, is correct, then the formation of deep-mantle plumes to vent material from it appears inevitable. We demonstrate quantitatively that the strong temperature dependence of viscosity guides the upward flow into long-lived chimneys that are ～ 20 km in diameter near the base of the mantle and decrease in width with progressive upward softening and partial melting of plume material. The speed of flow up the axis of the plume is correspondingly fast; 1.6 m y^?1 at the base and 4.8 m y^?1 at 670 km depth. Thermal diffusive spreading of a heated plume is compensated by a slow horizontal convergence of mantle material toward the chimney in response to the lower pressure there. This convergence, which contributes only a small increment to the flux of material up the plume, also serves to throttle the flow in the chimney. The global plume mass flux necessary to transport 1.6 × 10¹² W of core heat upward through the mantle is 1.8 × 10⁶ kg s^?1. At its base, plume material is probably still significantly below its solidus or eutectic temperature, but substantial partial melting is very likely as it rises. We speculate that a small fraction of this fluid component eventually emerges at the surface in “hot spots”, with the fate of the remainder being unknown. The behaviour and properties of D″ and of plumes are closely coupled. Not only are plumes a necessary consequence of a thermal boundary layer, but their existence is impossible without that layer.     

How Does Subsurface Characterization Affect Simulations of Hyporheic Exchange?

We investigated the role of increasingly well‐constrained geologic structures in the subsurface (i.e., subsurface architecture) in predicting streambed flux and hyporheic residence time distribution (RTD) for a headwater stream. Five subsurface realizations with increasingly resolved lithological boundaries were simulated in which model geometries were based on increasing information about flow and transport using soil and geologic maps, surface observations, probing to depth to refusal, seismic refraction, electrical resistivity (ER) imaging of subsurface architecture, and time‐lapse ER imaging during a solute tracer study. Particle tracking was used to generate RTDs for each model run. We demonstrate how improved characterization of complex lithological boundaries and calibration of porosity and hydraulic conductivity affect model prediction of hyporheic flow and transport. Models using hydraulic conductivity calibrated using transient ER data yield estimates of streambed flux that are three orders of magnitude larger than uncalibrated models using estimated values for hydraulic conductivity based on values published for nearby hillslopes (10^?4 vs. 10^?7 m²/s, respectively). Median residence times for uncalibrated and calibrated models are 10³ and 10⁰ h, respectively. Increasingly well‐resolved subsurface architectures yield wider hyporheic RTDs, indicative of more complex hyporheic flowpath networks and potentially important to biogeochemical cycling. The use of ER imaging to monitor solute tracers informs subsurface structure not apparent from other techniques, and helps to define transport properties of the subsurface (i.e., hydraulic conductivity). Results of this study demonstrate the value of geophysical measurements to more realistically simulate flow and transport along hyporheic flowpaths.     

Evaluating MT3DMS for Heat Transport Simulation of Closed Geothermal Systems

Owing to the mathematical similarities between heat and mass transport, the multi-species transport model MT3DMS should be able to simulate heat transport if the effects of buoyancy and changes in viscosity are small. Although in several studies solute models have been successfully applied to simulate heat transport, these studies failed to provide any rigorous test of this approach. In the current study, we carefully evaluate simulations of a single borehole ground source heat pump (GSHP) system in three scenarios: a pure conduction situation, an intermediate case, and a convection-dominated case. Two evaluation approaches are employed: first, MT3DMS heat transport results are compared with analytical solutions. Second, simulations by MT3DMS, which is finite difference, are compared with those by the finite element code FEFLOW and the finite difference code SEAWAT. Both FEFLOW and SEAWAT are designed to simulate heat flow. For each comparison, the computed results are examined based on residual errors. MT3DMS and the analytical solutions compare satisfactorily. MT3DMS and SEAWAT results show very good agreement for all cases. MT3DMS and FEFLOW two-dimensional (2D) and three-dimensional (3D) results show good to very good agreement, except that in 3D there is somewhat deteriorated agreement close to the heat source where the difference in numerical methods is thought to influence the solution. The results suggest that MT3DMS can be successfully applied to simulate GSHP systems, and likely other systems with similar temperature ranges and gradients in saturated porous media.     

Thermal structure and energy of the Hakone volcano,Japan

The thermal energy balance and the temperature profile of the Hakone volcano are considered quantitatively. Across the Hakone volcano and its surroundings the heat flow values vary from 10^–1 to 10³ mW/m², due to thermal conduction and mass flow involving volcanic steam and hot spring discharge. An area with extremely low heat flow is observed in the western side of the caldera showing the presence of percolating meteoric water. The hydrothermal activity is intense in the eastern half of the caldera.The total heat discharge from the high temperature zone (discharge area) of the Hakone volcano amounts to 11.0×10⁷ W. The magmatic steam energy discharge is 95.0×10⁶ W. The thermal energy by redistribution of the terrestrial heat flow by the lateral deep ground water flow is calculated to be 9.00×10⁶ W. For the model having the vertical vent in the volcano's central part up to 1 km depth below the ground surface from a magma reservoir the computed temperature distribution is consistent with the observed values. The depth of the magma reservoir is 7 km below the ground surface and the diameter is 5 km.     

Estimation of heat discharge rates using infrared measurements by a helicopter-borne thermocamera over the geothermal areas of Unzen volcano, Japan

Unzen volcano is situated on Shimabara Peninsula, western Kyushu, Japan. On the flank of the volcano, the Obama, Unzen and Shimabara hot springs are aligned in a direction from the southwest to the northeast across the peak. At Obama and Shimabara, heat is transferred mainly by water flow. But at Unzen heat is transferred by the discharge of natural steam and by conduction as well as water flow. In order to estimate the heat discharge by mechanisms other than water flow, infrared measurements by a helicopter-borne thermocamera were conducted over the Unzen hot spring area. The heat discharge was calculated from the thermal image by a method based on heat balance of the ground surface resulting in a value of 1.9 × 10⁶ cal/s (7.9 MW). The heat discharged by all mechanisms including that by water flow is estimated to be 5.0 × 10⁶ cal/s (21 MW). Similar preliminary estimates have been made for heat discharge at the Obama and Shimabara hot springs giving values of about 1.2 × 10⁷ cal/s (50 MW) and 1.0 × 10⁵ cal/s (0.4 MW), respectively. These values indicate that the heat discharge decreases with distance from the southwest to the northeast direction across the volcano. The total heat discharge from three hot spring areas on Unzen volcano is about 1.7 × 10⁷ cal/s (71 MW).The heat balance method appears useful for quantitative analysis of regional trends but its accuracy may not be always sufficient for detailed surveys. Several methods of determining heat flow, including the heat balance method, were compared at a test field in the Unzen hot spring area. The values obtained by the heat balance method coincide roughly with the other results but more detailed analysis is necessary to improve the accuracy of current methods of geothermal measurements.     

Efficient flow and transport simulations in reconstructed 3D pore geometries

Upscaling pore-scale processes into macroscopic quantities such as hydrodynamic dispersion is still not a straightforward matter for porous media with complex pore space geometries. Recently it has become possible to obtain very realistic 3D geometries for the pore system of real rocks using either numerical reconstruction or micro-CT measurements. In this work, we present a finite element–finite volume simulation method for modeling single-phase fluid flow and solute transport in experimentally obtained 3D pore geometries. Algebraic multigrid techniques and parallelization allow us to solve the Stokes and advection–diffusion equations on large meshes with several millions of elements. We apply this method in a proof-of-concept study of a digitized Fontainebleau sandstone sample. We use the calculated velocity to simulate pore-scale solute transport and diffusion. From this, we are able to calculate the a priori emergent macroscopic hydrodynamic dispersion coefficient of the porous medium for a given molecular diffusion D_m of the solute species. By performing this calculation at a range of flow rates, we can correctly predict all of the observed flow regimes from diffusion dominated to convection dominated.     

Impact of Groundwater Flow and Energy Load on Multiple Borehole Heat Exchangers

The effect of array configuration, that is, number, layout, and spacing, on the performance of multiple borehole heat exchangers (BHEs) is generally known under the assumption of fully conductive transport. The effect of groundwater flow on BHE performance is also well established, but most commonly for single BHEs. In multiple‐BHE systems the effect of groundwater advection can be more complicated due to the induced thermal interference between the boreholes. To ascertain the influence of groundwater flow and borehole arrangement, this study investigates single‐ and multi‐BHE systems of various configurations. Moreover, the influence of energy load balance is also examined. The results from corresponding cases with and without groundwater flow as well as balanced and unbalanced energy loads are cross‐compared. The groundwater flux value, 10^?7 m/s, is chosen based on the findings of previous studies on groundwater flow interaction with BHEs and thermal response tests. It is observed that multi‐BHE systems with balanced loads are less sensitive to array configuration attributes and groundwater flow, in the long‐term. Conversely, multi‐BHE systems with unbalanced loads are influenced by borehole array configuration as well as groundwater flow; these effects become more pronounced with time, unlike when the load is balanced. Groundwater flow has more influence on stabilizing loop temperatures, compared to array characteristics. Although borehole thermal energy storage (BTES) systems have a balanced energy load function, preliminary investigation on their efficiency shows a negative impact by groundwater which is due to their dependency on high temperature gradients between the boreholes and surroundings.     

Uncertainty in 1D Heat‐Flow Analysis to Estimate Groundwater Discharge to a Stream

Temperature measurements have been used by a variety of researchers to gain insight into groundwater discharge patterns. However, much of this research has reduced the problem to heat and fluid flow in one dimension for ease of analysis. This approach is seemingly at odds with the goal of determining spatial variability in specific discharge, which implies that the temperature field will vary in more than one dimension. However, it is unclear how important the resulting discrepancies are in the context of determining groundwater discharge to surface water bodies. In this study, the importance of these variations is examined by testing two popular one‐dimensional analytical solutions with stochastic models of heat and fluid flow in a two‐dimensional porous medium. For cases with low degrees of heterogeneity in hydraulic conductivity, acceptable results are possible for specific discharges between 10^?7 and 10^?5 m/s. However, conduction into areas with specific discharges less than 10^?7 m/s from adjacent areas can lead to significant errors. In some of these cases, the one‐dimensional solutions produced estimates of specific discharge of nearly 10^?6 m/s. This phenomenon is more likely in situations with greater degrees of heterogeneity.     

The thermal boundary-layer interpretation of D″ and its role as a plume source

The “anomalous” layer in the lowermost mantle, identified as D″ in the notation of K.E. Bullen, appears in the PREM Earth model as a 150 km-thick zone in which the gradient of incompressibility with pressure, $\text{d}\text{K}\text{d}\text{P}$, is almost 1.6, instead of 3.2 as in the overlying mantle. Since PREM shows no accompanying change in density or density gradient, we identify D″ as a thermal boundary layer and not as a chemically distinct zone. The anomaly in $\text{d}\text{K}\text{d}\text{P}$ is related to the temperature gradient by the temperature dependence of K_s, for which we present a thermodynamic identity in terms of accessible quantities. This gives the numerical result (?K_s/?T)_P=?1.6×10⁷ Pa K^?1 for D″ material. The corresponding temperature increment over the D″ range is 840 K. Such a layer cannot be a static feature, but must be maintained by a downward motion of the lower mantle toward the core-mantle boundary with a strong horizontal flow near the base of D″. Assuming a core heat flux of 1.6 × 10¹² W, the downward speed is 0.07 mm y^?1 and the temperature profile in D″, scaled to match PREM data, is approximately exponential with a scale height of 73 km. The inferred thermal conductivity is 1.2 W m^?1 K^?1. Using these values we develop a new analytical model of D″ which is dynamically and thermally consistent. In this model, the lower-mantle material is heated and softened as it moves down into D″ where the strong temperature dependence of viscosity concentrates the horizontal flow in a layer ～ 12 km thick and similarly ensures its removal via narrow plumes.     

Modelling Feedback of Chemical Reactions on Flow Fields in Hydrothermal Systems

Coupled reactive transport models of hydrothermal systems provide new insights and deeper understanding of the processes occurring due to fluid flow, heat transfer, solute transport, and chemical reactions. Basic concepts of species transport (diffusion, dispersion, and advection) and chemical precipitation and dissolution reactions are discussed, and five end-member types of reactive transport environments are introduced. One of these reactive transport environments, named ‘reactions within thermal gradients’, is used to demonstrate how free thermal convection can lead to redeposition of minerals and, due to the feedback of reaction on the flow field, a change of the convection pattern. The direct consequence of changing the flow field is a significant variation of the temperature distribution within the modelled area. With the example it is shown how reactive transport simulation can be applied for the detailed study of fossil and recent hydrothermal systems.     

Curie-point Depth from Spectral Analysis of Magnetic Data in Erciyes Stratovolcano (Central TURKEY)

Curie-point depth and heat flow values of the Erciyes region are determined to identify the thermal regime of the Central Anatolia by applying the spectral analysis method to the magnetic anomaly data. To compute the spectrum of the data, the magnetic anomaly of the region is transformed into 2-D Fourier domain to attain the average Curie depth. This method is useful in determining the top boundary of magnetic anomaly sources and reveals the Curie depth as 13.7 km in the study area. The obtained results imply a high thermal gradient (42.3°C km^?1) and corresponding heat flow values (88.8 mWm^?2) in the research area. Using the temperature value measured at borehole drilled by the General Directorate of Mineral Research and Exploration of Turkey (MTA), the values for the thermal gradient and heat flow value were computed as 50.7°C km^?1, 106.5 mWm^?2. From the heat flow value, the Curie-point depth was determined as 11.4 km in this region. It is concluded from the obtained values that the region has very high geothermal potential caused by partial melting of the lower crust.