Open AccessArticle

Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer

Francisco Alhama

¹,

José Antonio Jiménez-Valera

^2,*

and

Iván Alhama

Applied Physics Department, Technical University of Cartagena, C/Doctor Fleming s/n., 30202 Cartagena, Spain

Mining and Civil Engineering Department, Technical University of Cartagena, Paseo Alfonso XIII, n° 52 3020352, 30203 Cartagena, Spain

Author to whom correspondence should be addressed.

Appl. Sci. 2024, 14(2), 922; https://doi.org/10.3390/app14020922

Submission received: 18 November 2023 / Revised: 17 January 2024 / Accepted: 19 January 2024 / Published: 22 January 2024

(This article belongs to the Special Issue Advances and Applications of Numerical Analysis and Intelligent Computing)

Download

Browse Figures

Versions Notes

Abstract

A general and precise protocol that follows the standards of an inverse problem in engineering is proposed to estimate groundwater velocity from experimental lectures of temperature vertical profiles in a 2D aquifer. Several values of error in the temperature measurements are assumed. Since a large quantity of parameters and initial conditions influence the solution of this process, the protocol is very complex and needs to be tested to ensure its reliability. The studied scenario takes into account the input temperature of the water as well as the isothermal conditions at the surface and bottom of the aquifer. The existence of an input region, in which profiles develop to become linear, allows us to eliminate experimental measurements beyond such a region. Once the protocol is developed and tested, it is successfully applied to estimate the regional (lateral) groundwater velocity of the real aquifer and the result compared with estimations coming from the piezometric map.

Keywords:

inverse problem; groundwater flow; temperature profiles; porous media; heat flux

1. Introduction

The steady-state temperature profiles in a depth aquifer, where groundwater penetrates from the left border to a limited region where profiles develop, is strongly dependent on the groundwater velocity and its input temperature, which generally differs from the surface and bottom temperatures of the domain. This problem has been studied decades ago [1,2,3,4] and more recently [5,6,7,8,9].

The horizontal groundwater flow with a constant temperature coming from an inlet border, in an extended aquifer whose surface and bottom are subjected to isothermal (Diritlech) conditions, gives rise to a temperature field within the aquifer that depends on both the depth and the horizontal location. The extension of the region in which the temperature field is developed has been recently studied and characterized dimensionally [10]. These authors provide expressions for the width of the developed region as a function of the depth of the aquifer, groundwater velocity and thermal diffusivity of the porous media, and the latter combines the thermal conductivity of the soil–fluid matrix and the specific heat of the fluid. Based on the developed length of the temperature profile, the vertical dependence of the temperature profiles can be expressed as a function of the relative horizontal location with respect to that length, while the horizontal profile throughout the development region is a function of the relative location with respect to the depth of the aquifer. Beyond this development region, the profiles are constant and do not allow estimations of water flows. This explains the non-dependence of the profiles on the global horizontal extent of the aquifer, a novel result not assumed in these scenarios by other authors who impose the added and unrealistic condition of the existence of a constant horizontal thermal gradient along the entire aquifer, resulting in different profiles throughout the domain [11,12]. After these studies, no other works have been published on the estimation of horizontal fluxes with the hypothesis of the existence of a development region. In addition, all the temperature profiles resulting from the temperature field are very sensitive to the dimensionless group defined from the three temperatures that set the boundary conditions.

An estimation of groundwater velocity in the form of an inverse problem could be carried out using the universal type-curves related to temperature profiles, but this would require a large number of abacuses, because of the many potential values of the dimensionless group related to the three boundary temperatures. The other and more general technique for the velocity estimation is that of setting a classical inverse problem protocol [13], based on the minimization of a functional that compares, in successive steps, the experimental measurements with those coming from the numerical simulation of scenarios that progressively approach the solution [14]. This technique has been successfully applied to inverse problems of similar complexity in different fields [15,16,17,18]. Due to the large number of parameters involved in these problems (depth, boundary temperatures and thermal properties, as well as the experimental thermal profiles), the proposed protocol is complex and requires a sufficient number of iterations. Two functionals are defined, one for each position of the temperature measurements. After numerical iterations, the sum of these functionals converges toward a minimum value corresponding to a scenario that provides the desired estimates.

This document is structured as follows. In Section 3, the physical and mathematical model of the scenario under study is described. The steps and block diagram of the applied inverse problem protocol, which uses a classical technique [13], for flow estimation are presented in Section 4 and the results are verified in Section 5. The verification includes the introduction of temperature deviations in the measurements to simulate the errors of the measuring instruments. Finally, in Section 6, such a protocol is applied to a real aquifer, and the estimates flows are compared with those obtained applying Darcy’s Law.

2. Physical Scenario

The scenario studied is shown schematically in Figure 1. In the saturated 2D rectangular aquifer, with impervious top and bottom surfaces, the groundwater penetrates through the left vertical border with constant velocity (

v_{o}

) which is maintained throughout the aquifer due to the imposed piezometric gradient. Thermal conditions of the first class (Diritlech), with different temperature values, are imposed at the left, upper and bottom boundaries. It is assumed that the aquifer is large enough so that the region of development of the vertical temperature profiles, limited by

l_{T}^{*}

, is less than the length of the aquifer (

L

). Beyond this region, the vertical temperature profiles are straight lines. In order for this scenario to reproduce the behavior of real aquifers, no other conditions or assumptions regarding the natural processes of heat advection–diffusion within the aquifer are added. The imposition of horizontal thermal gradients assumed by other authors [11,12], for example, would radically change the solution of the problem, since the vertical temperature profile would change throughout the entire scenario regardless of its horizontal size.

Assuming a homogeneous and isotropic 2-D aquifer, the mathematical model is formed via the heat balance governing equation, which involves diffusion, convection and storage processes, plus the boundary thermal conditions:

k (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}}) - ρ_{e, w} c_{e, w} (v_{o} \frac{\partial T}{\partial x}) = ρ_{e} c_{e} \frac{\partial T}{\partial t}

(1)

T_{(x = 0, y, t)} = T_{w}

(2)

T_{(x, y = 0, t)} = T_{s}

(3)

T_{(x, y = H, t)} = T_{b}

(4)

T_{(x, y, t = 0)} = T_{i n i}

(5)

According to [10], the horizontal extent of the development of temperature profiles is

l_{T}^{*} = C_{1} \frac{H^{2} v_{o}}{α_{m}}

(6)

with

C_{1}

as a constant that is the function of the chosen criterion to define this length. So, when the temperature at the right end reaches 99% of its steady-state value at

y = H / 2

C_{1}

= 0.49. As regards the curvature of the profiles within the development region, this depends on the value of the monomial

\frac{T_{w} - T_{s}}{T_{b} - T_{w}}

The mathematical model is numerically simulated using the free software Ngspice [19] via a precise model based on the network simulation method [20], a tool that has demonstrated efficiency and reliability in many other problems of similar or higher complexity [21,22,23,24,25,26]. Ngspice provides the exact solution of the networks due to its powerful algorithms, so that the negligible errors caused by computation are only due to the grid size. In the meshes considered in this article (or the order of 40 × 100), the resulting errors are below 0.1% [19].

3. The Inverse Problem, $v_{o}$ versus $T - y$

This inverse problem consists of determining flow velocities via direct measurements of two sets of N discrete values of temperature–depth at two horizontal positions,

{(T - y)}_{x_{a}, i}

and

{(T - y)}_{x_{b}, i}

, 1 ≤ i ≤ N and

x_{a} > x_{b}

. It is assumed that steady-state temperature conditions have been reached throughout the whole aquifer. The duration of the transient to stabilize the profiles depends essentially on the water velocity, which marks the advection process below the water table of temperature

T_{s}

and which defines, in turn, the time it takes for the water to travel through the aquifer extension (

l_{T}^{*}

) where the temperature–depth profile develops.

The geometric shape of these two discrete temperature profiles, which are affected by the instrument error, must be non-linear to ensure that the positions of the measurements are within the development region of the temperature profiles (

x_{a}

l_{T}^{*}

and

x_{b}

l_{T}^{*}

) and that the proposed inverse problem protocol can, therefore, be applied. The curvature of the experimental temperature profiles ensures that the latter condition is satisfied. The temperature field in the aquifer and, as a consequence, the abovementioned profiles, are determined by the groundwater velocity

v_{o}

. It is assumed that locations

x_{a}

and

x_{b}

with respect to the origin of coordinates (or the water inlet region location) are not known, since such an origin is also unknown. Of course, the value of

l_{T}^{*}

is not known a priori since this value depends on

v_{o}

. If the conditions

x_{a}

l_{T}^{*}

and

x_{b}

l_{T}^{*}

are not satisfied (which means that experimental profiles are linear or quasi-linear), no solution can be obtained. Then, locations must be corrected by choosing instead another (or others) closer to the water inlet region. It is also assumed that the horizontal thermal diffusivity is negligible against advection to make the influence of velocity, as opposed to horizontal diffusion, predominant and that the time necessary to reach the steady state has elapsed. Finally, parameters

α

and

H

are assumed to be known, as well as temperatures

T_{s}

and

T_{b}

T_{i n i}

is irrelevant since we work at steady-state conditions. In short, the input data of the inverse problem are

H

α

T_{s}

T_{b}

{(T - y)}_{x_{a}, i}

{(T - y)}_{x_{b}, i}

and

x_{o}

(the distance

x_{a} - x_{b}

). The input temperature of the water,

T_{w}

, is also an unknown.

The protocol for the solution of the inverse problem synthesizes as follows:

Step 1. A scenario is numerically simulated with the known data. A value of

T_{w, 1}

is proposed that can be initially adjusted from the shape of the measured temperature profiles, determining whether it is an inlet temperature outside the limits of these profiles (upper or lower) or between such limits [10]. A first random velocity

v_{o, I}

(

I

refers to first iteration) is also imposed, greater than the supposedly expected one. The numerical solution provides the temperature field in the aquifer that is named

T {(x, y)}_{s i m, I}

Step 2. The discrete profile

T - y

closest to the water entry boundary

{(T - y)}_{x_{b}}

is compared with those resulting from the simulation in successive horizontal positions, starting with values close to locations where the profile is nearly linear. The following functional is defined:

Ψ_{b, x_{s i m}} = \sum_{i = 1}^{i = N} {[{(T - y)}_{x_{b}, i} - {(T - y)}_{x_{s i m, i}}]}^{2}

(7)

where

{(T - y)}_{x_{s i m, I}}

is the set of N values of the simulation corresponding to the horizontal position. Once obtained,

Ψ_{b, x_{s i m}}

and the new

Ψ_{b, x_{s i m}}

are evaluated from successive locations,

x_{s i m} = x_{s i m} - {Δ x}_{s i m}

, until reaching the location for which

Ψ_{b, x_{s i m}}

is minimum.

{Δ x}_{s i m}

is a finite value sufficiently small to produce an appreciable change in the functional. Let us name

x_{b, I}

and

Ψ_{b, m i n, I}

the final convergent values of this iteration.

Step 3. The values

{(T - y)}_{x_{a}, i}

are compared with those of the simulated model at

x_{a, I} = x_{b, I} + x_{o}

x_{o} = x_{a} - x_{b}

. The new functional is built whose value is

Ψ_{a, I}

Ψ_{a, j} = \sum_{1}^{N} {[{(T - y)}_{x_{a}, i} - {(T - y)}_{x_{a, j, i}}]}^{2}

(8)

Once the former steps related to the first iteration (

j

= I) are finalized, the partial results

v_{o, I}

x_{b, I}

Ψ_{b, m i n, I}

x_{a, I}

and

Ψ_{a, I}

are retained.

Step 4.

v_{o, I}

is modified by decreasing it by a small value,

{∆ v}_{o, I}

(a small fraction of

v_{o, I}

, for example, is 0.05·

v_{o, I}

). The model is simulated again for the new velocity of

v_{o, I I} = v_{o, I} - {∆ v}_{o, I}

. Steps 2 and 3 are repeated for the new iteration (

j

= II) providing the new values,

v_{o, I I}

x_{b, I I}

Ψ_{b, m i n, I I}, x_{a, I I}

and

Ψ_{a, I I}

Step 5. If

Ψ_{b, m i n, I I} + Ψ_{a, I I}

is less than

Ψ_{b, m i n, I} + Ψ_{a, I}

, go back to Step 4 (which takes one back to the application of steps 2 and 3, with

j = I I I

) to obtain the new values

v_{o, I I I}

x_{b, I I I}

Ψ_{b, m i n, I I I}, x_{a, I I I}

and

Ψ_{a, I I I} .

Continue with new iterations until reaching one,

j

+ 1, for which

Ψ_{b, m i n, j + I} + Ψ_{a, j + I}

is greater than

Ψ_{b, m i n, j} + Ψ_{a, j}

. The estimated velocity is that corresponding to iteration

j

. Experimental profiles

{(T - y)}_{x_{a}, i}

and

{(T - y)}_{x_{b}, i}

are estimated to be at locations

x_{b, j}

and

x_{a, j}

. The characteristic length can then be determined from its expression using Equation (6).

Step 6. If

Ψ_{b, m i n, I I} + Ψ_{a, I I}

is greater than

Ψ_{b, m i n, I} + Ψ_{a, I}

, go back to Step 4, but take a negative increment for the velocity, −

{∆ v}_{o}

Step 7. To estimate

T_{w}

, new successive values of this input temperature,

T_{w, 2}

T_{w, 1}

{Δ T}_{w}

T_{w, 3}

T_{w, 2}

{Δ T}_{w}

…, with

{Δ T}_{w}

, a fraction of

T_{w, 1}

, and of each of these temperatures, steps 1 to 6 are repeated. The final estimates of

v_{o}

x_{b}

and

x_{a}

are those values of these quantities for which the functionals

Ψ_{b}

and

Ψ_{a}

are minimum.

Figure 2 shows a block (flow) diagram of the protocol, which could be optimized for faster convergence through programming routines, selecting non-constant values (dependent on the results of each iteration) for

{∆ v}_{o}

4. Verification of the Protocol

The proposed protocol is verified as follows by a specific test that guarantees its reliability. Firstly, a direct problem is posed (with all known parameters) and numerically simulated. From its solution, we deduce the temperature field of the scenario and, from it, the characteristic length (

l_{T}^{*}

) in which temperature profiles develop. Then, we read temperature values at different depths at two arbitrary positions,

x_{a}

and

x_{b}

, within the profile development region, 0 <

x

l_{T}^{*}

. Let us call

{(T - y)}_{x_{a}, i}

and

{(T - y)}_{x_{b}, i}

, with

x_{a} > x_{b}

i

= 1, 2, … and N, the discrete set of regularly distributed temperatures in depth. To introduce deviations in these readings, and to simulate measurement errors, random errors of maximum value

e

% are applied to each reading. The new sets of temperature values, denoted as

{(T - y)}_{x_{a}, e, i}

and

{(T - y)}_{x_{b}, e, i}

, will be the input data for the inverse

v_{o}

estimation problem.

The parameter values for the direct problem (for which a negligible horizontal diffusivity effect is guaranteed, compared to advection) and the discrete temperature profiles affected by the error, which constitute the input data for the inverse problem, are shown in Table 1. Two random errors, 1 and 2%, are chosen, while the number of measurements regularly distributed along the depth is N = 8.

A refined mesh of 41 (horizontal) × 100 (vertical) is adopted, requiring a computing time of a few seconds for each iteration. As ordered, 30 iterations have been simulated by applying steps 1 to 8 of the protocol. The final estimates and a sample of partial results, with

v_{o}

= 5 × 10⁻⁶ m/s and

{∆ v}_{x, o}

= 0.1 × 10⁻⁶ m/s, are shown in Table 2 and Table 3. The input data are

H

= 1.0 m,

T_{1}

= 0.0 °C,

T_{2}

= 0.2 °C,

T_{3}

= 1.0 °C,

α

= 10⁻⁷ m²/s,

v_{o}

= 5 × 10⁻⁶ m/s,

x_{a}

= 3 m, and

x_{b}

= 9 m.

Partial results of

Ψ_{b, m i n, j}

Ψ_{a, j}

and

Ψ_{b, m i n, j} + Ψ_{a, j}

are shown for the same pairs of relative positions

x_{a}

and

x_{b}

. There is no direct correlation between the evolution of the values of individual functionals. These values can increase or decrease with successive iterations, but their sum is always monotonously decreasing until it converges to its minimum value. The estimated final value for velocity in both cases verify the efficiency and accuracy of the proposed protocol. These estimates with errors, sufficiently small for this engineering problems, are

v_{o} = 4.9 \times 10^{- 6} m / s e_{v_{o}} = 2 % e = 1 %

v_{o} = 5.1 \times 10^{- 6} m / s e_{v_{o}} = 2 % e = 2 %

Note that the size mesh might influence notably the results for coarse meshes, since the partial locations

x_{a}

and

x_{b}

have to be located according to the chosen thickness of the volume element. This is the reason why the partial locations in Table 2 are the same for the two selected random errors.

5. Application to a Real Aquifer

The above protocol has been applied to the Agua Amarga coastal aquifer (Alicante, SE of Spain) to estimate lateral groundwater velocity. It is a very controlled aquifer due to environmental issues related to the salt marsh ecosystem developed on its surface [27]. This ecosystem interferes with the extraction wells that supply the Alicante I and II desalination plants [28] and the existence of an old salt industry in the coastal region. Agua Amarga is a depressed natural brackish water lagoon separated from the Mediterranean Sea by a sand bar with a length of 2 km and a width of 50–200 m. This Quaternary–Neogene aquifer, connected to the sea, has an average hydraulic gradient between 0.1 and 1.5%. The aquifer has an approximate extension of 1700 hectares, of which 150 hectares are occupied by plots of land where a salt industry was developed, whose activity took place approximately between 1925 and 1975 [29].

Figure 3 shows the aquifer location, its limits, and the position of the measurement wells, P-3 and P-4, separated by a distance of 300 m and far away from the salt marsh and the extraction wells of the desalination plants. The piezometric isolines derived from the aquifer control well network, also shown in Figure 3, allow us to ensure that the groundwater flow between the P-3 and P-4 boreholes is horizontal and directed towards the coastline. A thermal diffusivity of 3.6 × 10⁻⁶ m²/s has been estimated [30]. Figure 4 shows the geologic cross-section through the P-3 and P-4 well line to the sea. The aquifer, formed by an approximately 20 m thick layer of silty sands and conglomerates from the Pliocene–Pleistocene, is located under a thin layer of clayey mud with gypsum and salt (Quaternary). Beneath the aquifer, a limestone marl formation (Upper Messinian) acts as an impermeable barrier.

The temperature profiles of these wells, corresponding to March 2019, are shown in Figure 5. Six temperatures have been taken at positions separated by 1 m, between depths

y

= 3.0 to

y

= 8.0 m measured from the water table, thus avoiding the influence on temperature profiles due to changes in ambient temperature during the day. According to Figure 5, temperatures in the upper and lower boundary conditions at the domain are 18.2 and 20.2 °C, respectively. The simulated domain is 900 m, an extension large enough to ensure that, at the right edge, the profiles are fully developed and not dependent on the horizontal position. Partial results of the application of the protocol for successive iterations in velocity are as follows: functionals,

Ψ_{b, m i n, j}

Ψ_{a, j}

and its sum

Ψ_{b, m i n, j}

Ψ_{a, j}

. Estimated positions of the wells, in relation to the water inlet source,

x_{a, j}

and

x_{b, j}

, are shown in Table 4. The convergence of these results provides the following estimates for velocity and well positions:

v_{o}

= 8.0 × 10⁻⁵ m/s,

x_{a}

= 487.5 m (P-3) and

x_{b}

= 187.5 m (P-4). The protocol for estimating the temperature of the inlet groundwater flow, the partial results of which are shown in Table 5, gives a value of

T_{w}

= 20.0 (°C). This value corrects the positions of wells P-3 and P-4 to

x_{a}

= 424.5 m and

x_{b}

= 124.5 m, respectively. These results do not signify the existence of a regional groundwater source at 124.5 m from P-4. The regional flow of this aquifer comes mainly from the Campo de Elche region, an agricultural area in the province of Alicante (SE of Spain) located about 15 km from the piezometers. What rather confirms the estimation is that the water temperature of the regional flow is 20.0 °C when reaching the horizontal location of 124.5 m upstream P-4. It is from this location that the temperature profiles that have been used as data to estimate the groundwater velocity start to develop.

The above results satisfactorily approximate the empirical estimates deduced using Darcy’s law and hydraulic potential isoline profiles [29]. To determine the hydraulic conductivity, field tests were carried out in the P-3 and P-4 piezometers. The results of these tests are as follows:

P - 4 Gilg - Gavard test K = 5.1 \times 10^{- 3} m / s

P - 3 Lefranc test K = 2.5 \times 10^{- 3} m / s

Undoubtedly, this is a heterogeneous aquifer with hydraulic conductivity values closer to more permeable soils (clean sands) in the region under study. With the piezometric data,

h_{P - 4}

= −2.15 m and

h_{P - 3}

= −4.00 m, and the distance between wells,

∆ x

= 300 m, the resulting hydraulic gradient is

\frac{h_{P - 4} - h_{P - 3}}{∆ x} =

0.0062. Furthermore, the range of Darcy velocities is

v_{D a r c y}

[1.76 × 10⁻⁷, 2.81 × 10⁻⁶] m/s.

Since this is a coastal aquifer, it is necessary to investigate whether density-driven flows caused by concentration changes due to salinity changes produce fluid velocities comparable to the estimated velocity. The expression for these velocities is

v_{d d} = \frac{κ g (∆ ρ_{e, w})}{ε μ}

, with

∆ ρ_{e, w}

indicating the saltwater density changes, κ the intrinsic permeability, μ the dynamic viscosity and ε the porosity of the aquifer (

ε

= 0.6) [31]. For calculations, Figure 6 shows the salinity-depth profiles of wells P-3 and P-4, measured in March 2019. Although this is a complex coastal aquifer, the salt concentrations in the aquifer caused by marine intrusion are negligible compared to the remaining concentrations caused by former salt mining, as shown in these curves.

Regarding vertical density gradients, in the worst case, it only occurs in a vertical region near the surface in P-4. These density changes induce fluxes of the value

v_{d d} = \frac{κ g (∆ ρ_{e, w})}{ε μ}

, with the values

κ =

4.6 × 10⁻¹¹ m² [32], g = 9.8 m/s²,

∆ ρ_{e, w}

= 6.40 kg/m³, and

μ

= 10⁻³ kg/(m·s).

v_{d d}

= 4.81 × 10⁻⁶ m/s, which is 6.0% of the estimated velocity. These calculations allow the effects of density-driven flows in the problem to be neglected overall.

6. Final Comments and Conclusions

An estimation of the horizontal groundwater flow from discrete temperature profile measurements at two aquifer locations, within the characteristic region in which the profiles develop to become independent of the horizontal position, has been possible by means of a protocol in the form of an inverse problem. The application of such a technique allows us to obtain, at the same time, an estimation of the water temperature at the aquifer entrance as well as the limits of the development region of the thermal profile.

Through successive iterations in which the flow velocity and inlet temperature are conveniently modified, a suitable number of numerical simulations of the direct problem are performed with known physical and thermal parameters of the aquifer, until the sum of the functionals that compare the experimental measurements with those obtained from the direct problem is minimized. The successive values of flow velocity and inlet temperature imposed in the direct problem reduce the sum of these functionals to a minimum value that allows us to obtain the final estimates. The initial temperature chosen to start the process can be approximated in view of the profiles measured in the field, whose shape allows us to establish the range in which this temperature is found. As for the initial velocity, any choice is valid.

The verification of the proposed protocol is performed following a classical technique that starts from a direct scenario with all known data and to which random errors are applied to the temperature profile measurements. Two errors are assumed, 1 and 2%. Once the direct problem is solved, an iterative method is applied to a large number of simulations of new scenarios until we retrieve the convergence of the solutions and provide successfully the underground flow estimates for both errors.

Once the protocol has been verified, it is applied to a real 2D scenario in a region around two wells 300 m apart. This is a complex aquifer with high levels of salt concentration due to the activities of the now extinguished salt industry. It is ensured that the flow in the region of the wells is horizontal by the data of the piezometric network of the aquifer control, shown in Figure 3. It is also verified that the density-driven flows, both horizontal and vertical, are negligible compared to the velocity estimates deduced via the application of the protocol, which, in turn, are consistent with the calculations extracted from the hydraulic conductivity of the ground.

Author Contributions

Conceptualization, J.A.J.-V. and F.A.; methodology, J.A.J.-V., I.A. and F.A.; software, J.A.J.-V.; validation, J.A.J.-V. and I.A; formal analysis, J.A.J.-V. and F.A; investigation, J.A.J.-V., I.A and F.A.; re-sources, J.A.J.-V. and I.A; data curation, J.A.J.-V. and I.A; writing—original draft preparation, J.A.J.-V., I.A. and F.A.; writing—review and editing, J.A.J.-V., I.A. and F.A.; visualization, J.A.J.-V., I.A. and F.A.; supervision, I.A. and F.A.; funding acquisition, J.A.J.-V. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Fundación Séneca, Agencia de Ciencia y Tecnología, Región de Mur-cia, grant number 21271/FPI/19.

Data Availability Statement

All data used in the study appear in the manuscript.

Acknowledgments

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)”. Financial support for monitoring field surveys and access to the information in this research were provided by “Mancomunidad de los Canales del Taibilla”.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

$C_{1}$	constant (dimensionless)
$c_{e}$	volumetric heat capacity of the soil–fluid matrix (cal·m⁻³·°C⁻¹)
$c_{e, w}$	volumetric heat capacity of the water, (cal·m⁻³·°C⁻¹)
$e$	relative random error (%)
$g$	acceleration of gravity (m/s²)
$H$	height (total depth) of the domain (m)
$k$	thermal conductivity of the soil–fluid matrix (cal·s⁻¹·m⁻¹ °C⁻¹)
$L$	length of the aquifer (m)
$l_{T}^{*}$	characteristic length in which temperature profiles develop (m)
$N$	number of experimental measurements for the inverse problem
P-3, P-4	piezometers in Agua Amarga coastal aquifer
$t$	time (s)
$T$	temperature (°C)
$T_{s}$	temperature at the soil surface (°C)
$T_{b}$	temperature at the bottom of the aquifer (°C)
$T_{w}$	temperature of the water at the inlet border (°C)
$T_{i n i}$	initial soil temperature (°C)
$(T - y)$	vertical temperature–depth profile (°C)
$v_{o}$	horizontal groundwater flow velocity (m/s)
$v_{D a r c y}$	Darcy velocity (m/s)
$x, y$	spatial coordinates (m)
$x_{a}$ , $x_{b}$	horizontal locations within the interval [0, $l_{T}^{*}$ ] (m)
$x_{o}$	$x_{a}$ − $x_{b}$ (m)
$x_{s i m}$	horizontal location in simulation (m)
$α$	thermal diffusivity of the soil–fluid matrix (m²/s), $α = k / ρ_{e} c_{e}$
$α_{m}$	$α_{m} = k / ρ_{e, w} c_{e, w}$ (m²/s)
$∆ T_{w}$	temperature increment at the inlet boundary (°C)
Δ $v_{o}$	velocity increment used in inverse problem (m/s)
Δ $x$	increment of the horizontal location (m)
ε	porosity (dimensionless)
$κ$	intrinsic permeability (m/s)
$μ$	fluid dynamic viscosity (g/(m·s))
$ρ_{e}$	wet bulk density of the soil–fluid matrix (g/m³)
$ρ_{e, w}$	fluid density of the water (g/m³)
Ψ	mathematical functional

Subscripts

$a, b$	refers to x_a and x_b, respectively
$d d$	refers to density-driven flow
$e$	refers to random error (%)
$h, v$	refers to horizontal and vertical components of v
$i$	$i$ = 1, 2, … N index of a particular temperature of the profile
$j$	$j$ = I, II, III, … index of the iteration for calculate the functional
$m i n$	refers to minimum value of the functional
$s i m$	refers to simulated results
$x, y$	related to spatial directions $x$ and $y$ , respectively
$x_{a}$ , $x_{b}$	refers to horizontal locations (m)
$x_{s i m}$	location simulated (m)
$1, 2, 3$	refers to the successive values of inlet temperature of the water

References

Suzuki, S. Percolation measurements based on heat flow through soil with special reference to paddy fields. J. Geophys. Res. 1960, 65, 2883–2885. [Google Scholar] [CrossRef]
Stallman, R.W. Computation of ground-water velocity from temperature data. USGS Water Supply Pap. 1963, 1544, 36–46. [Google Scholar]
Bredehoeft, J.D.; Papadopulos, I.S. Rates of vertical groundwater movement estimated from the Earth’s thermal profile. Water Resour. Res. 1965, 1, 325–328. [Google Scholar] [CrossRef]
Lapham, W.W. Use of Temperature Profiles Beneath Streams to Determine Rates of Vertical Ground-Water Flow and Vertical Hydraulic Conductivity; Paper 2337; U.S. Geological Survey Water-Supply: Reston, VA, USA, 1989.
Holzbecher, E. Inversion of temperature time series from near-surface porous sediments. J. Geophys. Eng. 2005, 2, 343–348. [Google Scholar] [CrossRef]
Constantz, J. Heat as a tracer to determine streambed water exchanges. Water Resour. Res. 2008, 44, W00D10. [Google Scholar] [CrossRef]
Szymkiewicz, A.; Tisler, W.; Burzyński, K. Examples of numerical simulations of two-dimensional unsaturated flow with VS2DI code using different interblock conductivity averaging schemes. Geologos 2015, 21, 161–167. [Google Scholar] [CrossRef]
Duque, C.; Müller, S.; Sebok, E.; Haider, K.; Engesgaard, P. Estimating groundwater discharge to surface waters using heat as a tracer in low flux environments: The role of thermal conductivity. Hydrol. Process. 2016, 30, 383–395. [Google Scholar] [CrossRef]
Bendjoudi, H.; Cheviron, B.; Guérin, R.; Tabbagh, A. Determination of upward/downward groundwater fluxes using transient variations of soil profile temperature: Test of the method with Voyons (Aube, France) experimental data. Hydrol. Process. Int. J. 2005, 19, 3735–3745. [Google Scholar] [CrossRef]
Jiménez-Valera, J.A.; Alhama, F. Dimensionless Characterization to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in Aquifers. Mathematics 2022, 10, 2717. [Google Scholar] [CrossRef]
Lu, N.; Ge, S. Effect of horizontal heat and fluid flow on the vertical temperature distribution in a semiconfining layer. Water Resour. Res. 1996, 32, 1449–1453. [Google Scholar] [CrossRef]
Reiter, M. Using precision temperature logs to estimate horizontal and vertical groundwater flow components. Water Resour. Res. 2001, 37, 663–674. [Google Scholar] [CrossRef]
Beck, J.V.; Blackwell, B.; Clair, C.R.S. Inverse Heat Conduction: Ill-Posed Problems; Wiley-Interscience: New York, NY, USA, 1985. [Google Scholar]
Chantasiriwan, S. Inverse heat conduction problem of determining time-dependent heat transfer coefficient. Int. J. Heat Mass Transf. 1999, 42, 4275–4285. [Google Scholar] [CrossRef]
Zueco, J.; Bég, O.A. Numerical Determination of the Temperature Dependent Thermophysical Properties in Solid Materials: Experimental Instrumentation. In Proceedings of the International Heat Transfer Conference, Washington, DC, USA, 8–13 August 2010; pp. 285–291. [Google Scholar]
Sánchez-Pérez, J.F.; García-Ros, G.; Castro, E. Simultaneous determination of the position, release time and mass release rate of an unknown gas emission source in short-term emissions by inverse problem. Chem. Eng. J. 2022, 445, 136782. [Google Scholar] [CrossRef]
Wang, S.; Zhang, L.; Sun, X.; Jia, H. Solution to two-dimensional steady inverse heat transfer problems with interior heat source based on the conjugate gradient method. Math. Probl. Eng. 2017, 2017, 2861342. [Google Scholar] [CrossRef]
Lebedev, L.P.; Vorovich, I.I.; Gladwell, G.M.L. Functional Analysis: Applications in Mechanics and Inverse Problems; Springer Science & Business Media: Dordrecht, The Netherlands, 2012; Volume 41. [Google Scholar]
Ngspice. Open Source Mixed Mode, Mixed Level Circuit Simulator (Based on Berkeley’s Spice3f5). 2016. Available online: http://ngspice.sourceforge.net/ (accessed on 10 November 2023).
Horno, J. Network Simulation Method; Research Signpost: Trivandrum, India, 2002. [Google Scholar]
Bég, O.A.; Takhar, H.S.; Zueco, J.; Sajid, A.; Bhargava, R. Transient Couette flow in a rotating non-Darcian porous medium parallel plate configuration: Network simulation method solutions. Acta Mech. 2008, 200, 129–144. [Google Scholar] [CrossRef]
Caravaca, M.; Sanchez-Andrada, P.; Soto, A.; Alajarin, M. The network simulation method: A useful tool for locating the kinetic–thermodynamic switching point in complex kinetic schemes. Phys. Chem. Chem. Phys. 2014, 16, 25409–25420. [Google Scholar] [CrossRef]
Fernández-Gracía, M.; Sánchez-Pérez, J.F.; del Cerro, F.; Conesa, M. Mathematical Model to Calculate Heat Transfer in Cylindrical Vessels with Temperature-Dependent Materials. Axioms 2023, 12, 335. [Google Scholar] [CrossRef]
Garratón, M.C.; del Carmen García-Onsurbe, M.; Soto-Meca, A. A new Network Simulation Method for the characterization of delay differential equations. Ain Shams Eng. J. 2023, 14, 102066. [Google Scholar] [CrossRef]
Sánchez-Pérez, J.F.; Marín-García, F.; Castro, E.; García-Ros, G.; Conesa, M.; Solano-Ramírez, J. Methodology for Solving Engineering Problems of Burgers–Huxley Coupled with Symmetric Boundary Conditions by Means of the Network Simulation Method. Symmetry 2023, 15, 1740. [Google Scholar] [CrossRef]
Solano, J.; Balibrea, F.; Moreno, J.A.; Marín, F. Dry Friction Analysis in Doped Surface by Network Simulation Method. Mathematics 2023, 11, 1341. [Google Scholar] [CrossRef]
Alhama, I.; García-Ros, G.; González-Alcaraz, M.N.; Álvarez-Rogel, J. Long-term artificial seawater irrigation as a sustainable environmental management strategy for abandoned solar salt works: The case study of Agua Amarga salt marsh (SE Spain). Catena 2022, 217, 106429. [Google Scholar] [CrossRef]
González-Alcaraz, M.N.; Aránega, B.; Tercero, M.C.; Conesa, H.M.; Álvarez-Rogel, J. Irrigation with seawater as a strategy for the environmental management of abandoned solar saltworks: A case-study in SE Spain based on soil–vegetation relationships. Ecol. Eng. 2014, 71, 677–689. [Google Scholar] [CrossRef]
Alhama Manteca, I. Simulation and consequences of successive anthropogenic activity in the Agua Amarga coastal aquifer (southeast Spain). Hydrol. Sci. J. 2013, 58, 1072–1087. [Google Scholar] [CrossRef]
Jiménez-Valera, J.A.; Alhama, I.; Duque, C. Advanced Type Curves for Vertical Groundwater Flow Estimation from Temperature Profiles; Applications to Real Scenarios; Mining and Civil Engineering Department, Technical University of Cartagena: Cartagena, Spain, 2023. [Google Scholar]
Chapuis, R.P. Estimating the in situ porosity of sandy soils sampled in boreholes. Eng. Geol. 2012, 141, 57–64. [Google Scholar] [CrossRef]
Bense, V.F.; Kooi, H. Temporal and spatial variations of shallow subsurface temperature as a record of lateral variations in groundwater flow. J. Geophys. Res. 2004, 9, B04103. [Google Scholar] [CrossRef]

Figure 1. Physical scheme and boundary conditions.

Figure 2. Block diagram of the inverse problem protocol.

Figure 3. Location and limits of Agua Amarga coastal aquifer with wells P-3 and P-4 (Google Earth Pro).

Figure 4. Cross-section of the aquifer on the line defined by wells P-3 and P-4.

Figure 5. Real temperature–depth profiles measured in wells P-3 and P-4 of the Agua Amarga aquifer in March 2019.

Figure 6. Vertical profiles of electrical conductivity recorded in March 2019 in boreholes P-3 and P-4.

Table 1. Parameters of the scenario and sets of input temperature profiles at locations

x_{a}

and

x_{b}

, for e = 1 and 2%.

Table 1. Parameters of the scenario and sets of input temperature profiles at locations

x_{a}

and

x_{b}

, for e = 1 and 2%.

y (m)		$T_{x_{a}}$ (°C)			$T_{x_{b}}$ (°C)
y (m)	Real	e = 1%	e = 2%	Real	e = 1%	e = 2%
0.11	0.797	0.798	0.790	0.868	0.876	0.863
0.23	0.595	0.597	0.583	0.725	0.729	0.719
0.35	0.431	0.434	0.425	0.588	0.583	0.590
0.48	0.309	0.308	0.306	0.459	0.458	0.454
0.60	0.219	0.221	0.222	0.340	0.342	0.346
0.72	0.148	0.149	0.151	0.230	0.232	0.228
0.84	0.083	0.083	0.083	0.128	0.128	0.128
0.96	0.019	0.019	0.019	0.029	0.029	0.029

Table 2. Partial and final results for e = 1%.

$v_{o, i}$ (m/s)	$Ψ_{b, m i n, j}$	$Ψ_{a, j}$	$Ψ_{b, m i n, j}$ $+ Ψ_{a, j}$	$x_{a, j}$ (m)	$x_{b, j}$ (m)
2.0 × 10⁻⁶	0.0244	0.1059	0.1303	1.20	7.20
4.0 × 10⁻⁶	0.0066	0.0327	0.0939	2.40	8.40
4.5 × 10⁻⁶	0.0148	0.0156	0.0304	2.60	8.60
4.9 × 10⁻⁶	0.0047	0.0110	0.0157	3.00	9.00
5.1 × 10⁻⁶	0.0091	0.0127	0.0217	3.00	9.00

Table 3. Partial and final results for e = 3%.

$v_{o, i}$ (m/s)	$Ψ_{b, m i n, j}$	$Ψ_{a, j}$	$Ψ_{b, m i n, j}$ $+ Ψ_{a, j}$	$x_{a, j}$ (m)	$x_{b, j}$ (m)
2.0 × 10⁻⁶	0.0084	0.1106	0.1190	1.20	7.20
4.0 × 10⁻⁶	0.0129	0.0374	0.0502	2.40	8.40
4.5 × 10⁻⁶	0.0088	0.0196	0.0284	2.60	8.60
4.9 × 10⁻⁶	0.0209	0.0134	0.0344	3.00	9.00
5.1 × 10⁻⁶	0.0120	0.0115	0.0235	3.00	9.00

Table 4. Partial results for

T_{w, 1}

= 20.20 °C.

Table 4. Partial results for

T_{w, 1}

= 20.20 °C.

$v_{o, i}$ (m/s)	$Ψ_{b, m i n, j}$	$Ψ_{a, j}$	$Ψ_{b, m i n, j}$ $+ Ψ_{a, j}$	$x_{a, j}$ (m)	$x_{b, j}$ (m)
2.0 × 10⁻⁵	0.0382	0.4397	0.4779	346.50	46.50
3.0 × 10⁻⁵	0.0421	0.3363	0.3785	373.50	73.50
4.0 × 10⁻⁵	0.0385	0.1410	0.1795	391.50	91.50
5.0 × 10⁻⁵	0.0385	0.0835	0.1220	415.50	115.5
6.0 × 10⁻⁵	0.0384	0.0505	0.0889	436.50	136.50
7.0 × 10⁻⁵	0.0384	0.0406	0.0790	472.50	172.50
8.0 × 10⁻⁵	0.0312	0.0431	0.0743	487.50	187.50
8.5 × 10⁻⁵	0.0385	0.0504	0.0889	493.50	193.50
9.0 × 10⁻⁵	0.0385	0.0577	0.0962	505.50	205.50

Table 5. Partial and final results for

v_{o}

= 8.0 × 10⁻⁵ m/s.

Table 5. Partial and final results for

v_{o}

= 8.0 × 10⁻⁵ m/s.

$T_{w, k}$ (°C)	$Ψ_{b, m i n, j}$	$Ψ_{a, j}$	$Ψ_{b, m i n, j}$ $+ Ψ_{a, j}$	$x_{a, j}$ (m)	$x_{b, j}$ (m)
20.20	0.0312	0.0431	0.0743	487.50	187.50
20.10	0.0312	0.0427	0.0739	454.50	154.50
20.00	0.0311	0.0406	0.0717	424.50	124.50
19.90	0.0532	0.0384	0.0916	400.50	100.50

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Alhama, F.; Jiménez-Valera, J.A.; Alhama, I. Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer. Appl. Sci. 2024, 14, 922. https://doi.org/10.3390/app14020922

AMA Style

Alhama F, Jiménez-Valera JA, Alhama I. Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer. Applied Sciences. 2024; 14(2):922. https://doi.org/10.3390/app14020922

Chicago/Turabian Style

Alhama, Francisco, José Antonio Jiménez-Valera, and Iván Alhama. 2024. "Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer" Applied Sciences 14, no. 2: 922. https://doi.org/10.3390/app14020922

APA Style

Alhama, F., Jiménez-Valera, J. A., & Alhama, I. (2024). Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer. Applied Sciences, 14(2), 922. https://doi.org/10.3390/app14020922

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer

Abstract

1. Introduction

2. Physical Scenario

3. The Inverse Problem, $v_{o}$ versus $T - y$

4. Verification of the Protocol

5. Application to a Real Aquifer

6. Final Comments and Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Nomenclature

Subscripts

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Article Menu

Inverse Problem Protocol to Estimate Horizontal Groundwater Velocity from Temperature–Depth Profiles in a 2D Aquifer

Abstract

1. Introduction

2. Physical Scenario

3. The Inverse Problem, v o versus T − y

4. Verification of the Protocol

5. Application to a Real Aquifer

6. Final Comments and Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Nomenclature

Subscripts

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

3. The Inverse Problem, $v_{o}$ versus $T - y$