Volume 4 - Year 2021- Pages 157-164
DOI: 10.11159/ijci.2021.019

The Effect of Anisotropic Permeability on the Temperature Profiles Obtained In a River Discharge Scenario to a Deep Aquifer

José Antonio Jiménez-Valera1, Francisco Alhama2

1Technical University of Cartagena, Department of Mining and Civil Engineering
Paseo Alfonso XIII, nº52, Cartagena, Spain 30203
2Technical University of Cartagena, Applied Physics Department

Abstract - The discharge-recharge patterns between rivers and aquifers in soils are strongly influences by the degree of anisotropy in their hydraulic permeability. In consolidated soils, the ratio ‘horizontal/vertical’ permeabilities can reach values up to 10 and even more. As a consequence, as a result of the coupling between fluid flow and heat transport, the last caused by temperature boundary conditions, the temperature profiles turn also dependent of such anisotropy. In this paper, a typical river-aquifer scenario with permeable bottom is studied, searching both the flow patterns for different degree of anisotropy and the effect of the flow field in the temperatures vertical profiles. A harmonic temperature condition is applied to the soil surface while constant and different temperatures are imposed to the bottom of the aquifer and to the river. The coupled effect between flow and heat transport produces the formation of eddies of temperature in the form of thermal waves that displace to the up and down boundaries and to the right one due to the drag of the fluid. The mean temperature profiles tend to be lineal far from the river. The extension of the region where the profiles develop increase with the ratio ‘horizontal/vertical’ permeabilities. The problem is numerically solved by the network simulation method using the free software Ngspice.

Keywords: Recharge-discharge, River-water interactions, Anisotropic permeability, Harmonic temperature profiles, Flow and heat transport

© Copyright 2021 Authors - This is an Open Access article published under the Creative Commons Attribution License terms. Unrestricted use, distribution, and reproduction in any medium are permitted, provided the original work is properly cited.

Date Received: 2021-07-20
Date Accepted: 2021-07-26
Date Published: 2021-08-23

1. Introduction

The interaction between surface and deep water bodies, particularly between rivers and their underlying aquifers (the object of this work), is a topic of great interest in hydrogeology due to its importance for the storage and maintenance of underground water resources. Its study not only covers the mechanical problem of the investigation of 2D patterns of water flow, dependent on the contour conditions of hydraulic potential and soil permeability, but also extends to the temperature patterns determined by the coupling between water flow and heat transport imposed by temperature conditions at the surface and at the bottom of the domain.

Numerous works, most of them dedicated to in situ studies (Duque et al. [1]), have been published in the last decades. There are no analytical solutions other than for 1D cases of constant vertical flow (1D) and harmonic surface temperature conditions (Bredehoeft & Papadopulos [2] and Taniguchi [3]). The study of vertical temperature profiles has recently been used in the form of an inverse problem to determine recharge-discharge rates that are difficult and expensive to obtain by other means (Duque et al. [4]).

In consolidated soils, to a greater or lesser degree, the hydraulic permeability anisotropy can become important, with values of the ‘horizontal/vertical permeability’ ratio greater than 10 (Childs et al. [5] and Chapuis and Gill [6]). As a consequence, the velocity field (decoupled from the thermal problem) and the temperature field (coupled to the previous one) are strongly dependent on this quotient.

In the present work, a symmetric, 2D scenario composed of a river and its underlying aquifer with a thickness less than half the horizontal extension is investigated by means of numerical simulations using the network simulation method (Horno [7] and García-Ros et al. [8]) and the free software Ngspice ([9]).

The boundary conditions of constant piezometry of the river, large permeability at the bottom of the aquifer, and free condition at the right edge produce 2D flow patterns in which the development region of the linear (mean temperature) profile can be clearly distinguished. The extension of this region, clearly influenced by anisotropy in permeability, spreads as the ratio ‘horizontal/vertical’ permeability increases and can reach the right boundary deleting the formation of eddies.

In addition, the combination of the heat transported by dragging and the harmonic variation of the temperature imposed on the up and down surfaces produces the generation of thermal waves or eddies that moves along the domain towards the outlet direction. According to the magnitude of the flow field, which in turn is strongly dependent on the ‘horizontal/vertical’ permeability ratio, the region more influenced by the river temperature extends (as velocity increases) towards the outlet up to reach all the domain. The mean temperature-depth profiles are very changeable and difficult to characterize.

2. Physical, Mathematical and Network Models

The physical scheme of the problem (symmetric) is depicted in Figure 1. Temperature and flow boundary conditions are shown in this figure. A river with constant hydraulic (height) head discharges on the sub-soil aquifer whose bottom and vertical right boundary are permeable (outflow water) while the left boundary is impermeable (no flow). Vertical and horizontal components of the velocity are very different from one point to another depending on the anisotropy in the hydraulic conductivity.

A seasonal (annual) harmonic temperature is imposed at the upper region (soil surface) while the temperature of the river-water is constant. Vertical left boundary is adiabatic to reflect the symmetry of the scenario, the bottom maintains a constant value temperature and finally, the right vertical boundary is free to flow (heat is dragged by the water according to its temperature at the boundary).

Figure 1. Physical scheme of the problem. Hydraulic and heat boundary conditions.

As for the mathematical model, rectangular coordinates with the origin at the upper-left corner of the aquifer are chosen (see Figure 1). The list of equations is as follows



Eq. 1, Darcy constitutive law, is decoupled from the heat transfer. In this equation, kx and ky (m/s) are the horizontal and vertical hydraulic conductivities respectively, while h (J/N) represents the piezometric head.


Eq. 2, is the balance of the different heat flows (convection, conduction and storage) that locally counteract one another (heat conservation flow). Fourier law is used for the conduction flow. The term ceρe (calg−1°C−1) is the specific heat of the water-soil matrix, km(cals-1m-1°C-1) is the thermal conductivity of the water-soil matrix, cwρw (calg−1°C−1) is the specific heat of water, vx and vy are the horizontal and vertical components of the water velocity, and T (°C)  and t (s) the temperature and time variables.




As regards the flow boundary conditions, Eq. 3a is a Dirichlet condition for the head variable along the river. In this equation, lr (m) is the width of the river and hr (J/N) its piezometric head. Eq. 3b is the no-flow condition at the ground surface and at the symmetric vertical wall, and Eq. 3c is again a first class condition that emulates the constant head ho (J/N) at the outflow regions (bottom and right limits) of the aquifer. Temperature boundary conditions are described in Eqs. 4a, 4b and 4c.





A first class condition is imposed at the river-water, Eq. 4a, where Tr (°C) is the constant river temperature and at the ground surface. In Eq. 4b, surface soil temperature is dependent according to annual harmonic variation. A pure sinusoidal dependence has been chosen for simplifying. Tmed (°C) is the average soil surface temperature while ω (rad/s) is calculated as 2π/τ, with τ the period of the oscillation. Homogeneous second class or Neumann condition is imposed at the symmetric (adiabatic) line (no heat flow), Eq. 4c, while the right vertical boundary is assumed free, which means that flow dragged the heat according to the local temperature at the boundary. Temperature at the bottom (Eq. 4d) To is constant (°C).


Finally, Eq. (5) represent the initial temperature, Tini (°C).

This scenario has been adopted as typical without prejudice to the fact that real scenarios may be slightly or substantially different. Thus, it is common to find aquifers with an impermeable bottom and shallow aquifers, for which the boundary conditions are different. The initial temperature does not affect the final solution of the thermal field (steady state). The problem described by the mathematical model represents a unique case in which the most important deviations of the temperature field occur in a region of the aquifer close to the river bank. Beyond this region, whose extension would be interesting to determine, there are vertical temperature profiles where the recharge effect is not appreciated. These profiles represent the maximum or minimum value of the temperature that is reached at each point. From the mathematical point of view, the application of any other boundary condition does not present special difficulties from the point of view of numerical calculation.

The problem parameter values are listed in Table 1. Anisotropy in the hydraulic permeability is imposed by increasing its horizontal value 5 and 10 times in relation with the vertical value which retains a same value.

The design of the numerical model follows the rules of the network simulation method, a tool largely employed in the literature to solve many types of coupled, non-lineal similar or even more difficult scenarios (such as Bénard, Yusa and many others) that the one here studied. Really, the numerical model is a simple modification of the software ‘Fahet’ (Alhama et al. [10]), developed by the research group of the UPCT and applied for the solution of some benchmark problems (Alhama et al. [11]).

As we consider the network model as an auxiliary tool, we refer to the cited references to deepen in its design. A mesh grid of 100 X 20 cells has been used. The time step varies according to the smoothness of the solution and its values are ruled by the software according to the tolerance chosen by the user. Other parameters of the simulation are pointed out later.

Table 1. Parameter values of the problem

3. Numerical Solution

3. 1. Patterns of equipotential and flow lines

Preserving the geometric scale of the domain, Figures 2 and 3 show the patterns of iso-potentials and streamfunction, respectively.

For the isotropic case (Figures 2a and 3a), it can be seen that, as expected, most of the flow is directed towards the bottom of the aquifer and in the area under the river. In fact, 90% of the total flow from the river is located in the region 0 <x <lcarac However, it is important to note that neither the vertical nor the horizontal velocity of the fluid (which are relatively similar) in the left central region of the aquifer can be neglected, these give rise to strong couplings between the water and heat flows in this region. Beyond, to the right of the aquifer, the same phenomenon occurs, although with already negligible velocities, which induces standard vertical profiles practically independent of flow. On the other hand, it is noted that the discharge velocity at the bottom of the river changes with position, which confirms as the most suitable choice of the boundary condition that of constant hydraulic potential instead of that of constant velocity.

The effect of anisotropy on permeability can be seen in Figures 2b and 3b (for Kx = 5ky) and 2c and 3c (for kx = 10ky). As expected, the equipotential lines expand as the permeability ratio increases, although the line with the lowest potential remains within the domain due to the null potential condition at the right edge. The flow lines are more illustrative since they show that much of the water flow exits the right edge of the domain for increasing horizontal permeabilities. In fact, for kx = 10ky, more water comes out of the right border of the domain than at the bottom. The vertical and horizontal components of velocity can be seen through the slope of these lines

(a) kx = ky

(b) kx = 5ky

(c) kx = 10ky
Figure 2. Lines of constant hydraulic potential.

(a) kx = ky

(b) kx = 5ky

(c) kx = 10ky
Figure 3. Stream lines pattern.

From a negligible horizontal velocity almost throughout the domain for the case kx = ky, it passes to the emergence of wide regions where the horizontal velocity is appreciable for substantially higher horizontal conductivities. However, a central domain seems to be maintained in which horizontal and vertical velocities are comparable.

It can be affirmed that, in domains that are horizontally extensive enough, the isotropic case gives rise to a well-defined extension region in which the thermal profiles (of average temperature) will develop until a linear profile is reached. Beyond this region, velocities can be considered negligible and such thermal profiles invariant and linear.

The anisotropic case extends this characteristic length and can reach the right vertical border, which eliminates the linear thermal profile region for the mean temperature.

3.2 Patterns of isotemperature-time

To illustrate the complex variation of the time-temperature patterns in which the harmonic displacement of the heat islands or eddies is represented, Figure 4 shows these patterns for eight different times within the same annual period, once the initial transient has been overcome. In all cases these figure shows the emergence, development and fading of the islands of hot and cold temperatures which displace throughout the domain. For kx = ky (Figure 4), there is negligible horizontal water flow and the islands move vertically towards the surface and bottom of the aquifer while for kx >> ky, the appreciable horizontal flow force the islands to move also towards the right boundary (Figures 5 and 6).

Figure 4. Isotemperature pattern corresponding to nine successive time intervals over a period, for the isotropic case

Figure 5. Isotemperature pattern corresponding to nine successive time intervals over a period, for kx = 5ky

Figure 6. Isotemperature pattern corresponding to nine successive time intervals over a period, for kx = 10ky

3.3 Temperature Profiles

Figures 7 represents the maximum and minimum values of the temperature along different vertical columns within the aquifer. Three locations, all far from the river have been chosen. In the case of kx = ky, Figure 7a, the profiles are nearly symmetric, and the small difference between the location x = 44.75 m and x = 49.75 m is justified in that these locations are beyond the characteristic region in which the profiles develop.

For kx = 5ky (Figure 7b) the profiles have lost their symmetry because they are within the region of their development. The same occurs for kx = 10ky. The temperature Tmin is close to that of the river, an effect more accused for kx = 10ky, Figure 7c. As regards Tmax, the profile is more difficult to explain because the locations x= 34,75 m and 44,75 m fall within the profile development region where the coupling is strong.

These explanations are coherent with the profiles of mean temperature, Figure 8. For kx = ky, the profile is quite lineal just at the boundary and curves smoothly below such locations. The same occurs for kx = 5ky and kx = 10ky, but with profiles that are clearly curved. Due to the temperature boundary condition at z=H, the mean temperature converges in all cases to the bottom temperature forcing the profile to be bent near the bottom. Really, boundary equation (4d) can be considered as an idealized condition since it would be applied below the aquifer bottom in most cases.

Finally, Figure 9 shows the horizontal temperature profiles in two typical positions, y = 5.5 m (dashed line) and y = 15.5 m (continuous line). Again, this region shows the existence of a horizontal developing region, small than the length of the domain for kx = ky (where the temperature tends to its mean value whatever be the depth coordinate) and larger for kx > 5ky or 10ky, (where the mean temperatures tend to a different value according the depth because the developing region is larger than the length of the domain).

(a) kx = ky

(b) kx = 5ky

(c) kx = 10ky
Figure 7. Maximum and minimum temperature depth-profiles

Figure 8. Average temperature depth-profiles

Figure 9. Horizontal temperature profiles

4. Final Comments and Conclusions

The interaction between rivers and underlying aquifers in order to determine the recharge-discharge between them, knowing the piezometric conditions of the contour, to determine the velocity field is already a complex problem that requires numerical solutions.

Flow patterns are basically determined by the permeable or non-permeable nature of the aquifer bottom and the anisotropy in hydraulic permeability that can be pronounced in consolidated soils. The flow field is coupled to the heat transport problem, resulting in temperature patterns that can still be higher under the condition of harmonic temperature on the ground surface.

For cases of aquifers that are more or less extensive horizontally and with permeable bottom, the numerical simulations carried out show the emergence of a horizontal extension or characteristic length in which a linear profile of average temperatures that prevails beyond this length progressively develops. The increasing anisotropy of the soil, defined by the 'horizontal / vertical permeability' ratio, which can reach values of several tens, increases the size of this region of development of the linear profile and gives rise to more or less regular profiles of mean temperature dependent on the anisotropy within this region.

However, the maximum and minimum temperature profiles, of more complex interpretation, apparently lack regularity thanks to the strong coupling between flow and heat transport. These profiles curve and converge to the idealized boundary condition imposed by the temperature at the bottom. Thus, the Tmed profiles are the only ones capable of being used in the inverse problem of the determination of velocities and flows from the experimental knowledge of the profiles

The strong coupling of flow and heat transport, on the other hand, translates into the emergence of heat islands or eddies that evolve from the central region of the aquifer towards the horizontal boundaries and flow outflow due to diffusive and drag effects. Finally, the horizontal temperature profile highlights the existence of a profile development length that, eventually, can go beyond the right border for large anisotropies.

5. Acknowledgements

We would like to thank “Fundación Séneca” for the scholarship awarded to José Antonio Jiménez Valera, which will allow us to continue this investigation. “21271/FPI/19. Fundación Séneca. Región de Murcia (Spain)”


[1] C. Duque, M. L. Calvache and P. Engesgaard, "Investigating river-aquifer relations using water temperature in an anthropized environment (Motril-Salobreña aquifer)," Journal of Hydrology, vol. 381, no. 1-2, pp. 121-133, 2010. View Article

[2] J. D. Bredehoeft and I. S. Papadopulos, "Rates of vertical groundwater movement estimated from the Earth's thermal profile," Water Resources Research, vol. 1, no. 2, pp. 325-328, 1965. View Article

[3] M. Taniguchi, "Evaluation of vertical groundwater fluxes and thermal properties of aquifers based on transient temperature-depth profiles," Water Resources Research, vol. 29, no. 7, pp. 2021-2026, 1993. View Article

[4] C. Duque, S. Müller, E. Sebok, K. Haider and P. Engesgaard, "Estimating groundwater discharge to surface waters using heat as a tracer in low flux environments: The role of thermal conductivity," Hydrological Processes, vol. 30, no. 3, pp. 383-395, 2016

[5] E. C. Childs, N. Collis-George and J. W. Holmes, "Permeability Measurements in the Field as an Assessment of Anisotropy and Structure Development," Journal of Soil Science, vol. 8, no. 1, pp. 27-41, 1957. View Article

[6] R. P. Chapuis and D. E. Gill, "Hydraulic anisotropy of homogeneous soils and rocks: influence of the densification process," Bulletin of the International Association of Engineering Geology-Bulletin de l'Association Internationale de Géologie de l'Ingénieur, vol. 39, no. 1, pp. 75-86, 1989. View Article

[7] J. Horno, "Network simulation method," Research Signpost, 2002.

[8] G. García-Ros, I. Alhama, J. L. Morales, "Numerical simulation of nonlinear consolidation problems by models based on the network method," Applied Mathematical Modelling, vol. 69, pp. 604-620, 2019. View Article

[9] Ngspice [Software]. Open Source Mixed Mode, Mixed Level Circuit Simulator (Based on Berkeley's Spice3f5). [Online].Available: View Article

[10] I. Alhama, A. Soto, M. Cánovas and F. Alhama, "FAHET: Flow and Heat Transport simulator", © 2004 2009 UPCT, 2011.

[11] I. Alhama, M. Cánovas and F. Alhama, "Simulation of fluid flow and heat transport coupled processes using FAHET software", Journal of porous media, vol. 18, no. 5, pp. 537-546, 2015. View Article