Evolution of a Geothermal System with Fluid Re-Injection

: In this paper, we study the stationary conﬁgurations of a vapor-dominated geothermal system. We propose a general mathematical model, which turns out to be a free boundary problem, and we focus on the stationary solutions. The analysis, developed considering the data of the geothermal ﬁeld in Larderello, Tuscany (Italy), shows that steady-stable states may exist both when the re-injection takes place in depth and when the natural recharge is distributed both in the area occupied by vapor and in the one occupied by liquid.


Introduction
In volcanic areas, geothermal energy represents an important resource. For instance, the Tuscan basins of Larderello and Monte Amiata provide almost 30% of the energetic need of the whole region (see, for example, [1,2]). Several studies have been devoted to geothermal reservoirs with the aim of providing mathematical models-some of them very sophisticated-for extraction wells [3][4][5][6][7][8][9][10].
Here we refer to a specific situation-a vapor-dominated system [11,12]-in which the reservoir is sealed at the top by a nearly impermeable rock layer, under which water is in the vapor state up to a depth at which the pressure and temperature are such that it is present as liquid. Wells are drilled through the impermeable layer to reach the high-pressure vapor zone. The vapor then leaves the reservoir, and extraction is sustained by the vaporization of deep water (sometimes called plutonic or juvenile water). This is the case for the Larderello geothermal system (one of the largest in the world), which is characterized by a reservoir pressure that is much lower than hydrostatic pressure. This is a consequence of its natural evolution from an initial liquid-dominated to the current vapor-dominated system [3,10].
The genesis of vapor-dominated reservoirs has intrigued researchers for many years. In 1973, Truesdell and White [6] reviewed all the principal conceptual models proposed until then for the reservoir and possible mechanisms for superheating the steam. These models aimed to give a qualitative description of the physical processes leading to the development of a vapor-dominated reservoir. However, all efforts notwithstanding, many open questions remain.
A reservoir is always formed by a complex of fractured permeable rocks confined between impermeable layers and saturated by water, in a liquid and/or vapor phase. As far as the Larderello geothermal field is concerned, the heat source is a magmatic intrusion located at a relatively small depth (5-10 km) producing temperatures of 500 ÷ 700°C in the rocks surrounding the reservoir. Both the lower and upper boundaries of the reservoir depend on the local lithology. The so-called K horizon (a seismic reflecting surface), strictly related to a temperature of about 450°C, is considered to represent the bottom of the geothermal reservoir [10,13,14]. The top of the productive reservoir is set in correspondence to the first fractured level in the central areas of the field, where a great amount of well data are available. In the peripheral areas, where well data are lacking, the top is set in correspondence to the 250°C isotherm [13][14][15][16]. Usually, the shallow carbonate-anhydrite reservoir shows a temperature of 220 ÷ 250°C and pressure of about 15 ÷ 20 bar at 1000 m depth, while the deep metamorphic reservoir has a temperature of 300 ÷ 350°C and pressure around 70 bar at 3000 m depth [13].
Vapor-dominated geothermal fields are still little understood, and their production histories are so short that premonitory indications of exhaustion are not yet generally recognized. In the 1970s, Sestini [30] concluded that the production of steam from the Larderello fields at that time was approaching a steady state and could have been maintained for many years. However, Truesdell and White [6] came to quite different conclusions, saying that the Larderello field was approaching exhaustion. This appears to have been sustained 30 years later, as the time series of production shows, in the papers by Romagnoli et al. [10] and Barelli [9]. In this context, the study of water levels reported in 1978 by Celati et al. [31] is of great importance.
In the last 20 years, due to the increase of the power demand from renewable sources, many studies have investigated the Larderello-Travale geothermal field using fully three-dimensional models [10,14,[32][33][34]. Clearly, the reliability of the forthcoming conclusions of these approaches depends heavily upon the reliability of the boundary and initial conditions selected. In general, these studies are validated against the current rates of production with the aim of predicting future scenarios.
The present contribution focuses on the modeling of geothermal systems in which extraction, natural and artificial recharge mechanisms may act simultaneously. Indeed, the extracted geothermal fluid, after its utilization in a power plant, can be injected back into the reservoir through specific injection wells. For instance, in Larderello, re-injection activities started in the 1970s to contrast the production decline [2,13]. In this way, the natural recharge of the reservoir (by meteoric or juvenile water) is artificially integrated, drastically reducing the impact on the environment of power-plant operations [35,36].
The aim of this study is to investigate the existence and stability of possible equilibrium configurations that can guarantee sustainable production. We develop a model that is able to give some basic information. Of course, approximations have been introduced in order to develop a model which is mathematically treatable. For example, the fluid in the reservoir is assumed to be pure water; this is a very schematic description, as, in practice, many chemicals are carried by the water, which makes the real thermodynamics much more complex. Next, we consider a region with a cluster of wells distributed so that the vapor motion underneath can be reasonably considered to be one-dimensional. This may appear to be an excessive approximation, but actually it is largely justified by the fact that problems such as this are characterized by a great level of data uncertainty. Our model is a one-dimensional parabolic free boundary problem, where both extraction and recharge are treated as source terms in the continuity equation. In this framework, we develop a model describing the behavior of the geothermal field, the location of the water vaporization front and its stability with respect to perturbations. More precisely, we examine the steady (or quasi-steady) state attained by recharging the reservoir with water, looking for those recharging rates compatible with a stable state. In particular, we discuss the effects of both extraction and recharge (both natural and artificial). We determine the steady solution by solving a non-linear equation which, when peculiar conditions are met, has more than a single solution. A linear stability analysis of the stationary liquid/vapor interface (if it exists) is then used to select the physically acceptable solution. In particular, we show that selecting the parameters according to the field data allows us to obtain a stable interface that is qualitatively in agreement with the borehole extraction/injection data. In this respect, our model, although extremely simple, allows us to obtain some insights regarding the state of the reservoir and to predict the location of the phase transition surface. On one hand, the simplicity of the model has the advantage that one can deal with an exact solution; on the other hand, the price to pay is that a classical validation remains questionable,as a sharp liquid/vapor interface is very unlike to exist in a real field. However, it should also be taken into account that that even highly sophisticated three-dimensional models depend significantly upon the boundary and initial conditions imposed. In the majority of cases (and for Larderello in particular), these conditions suffer from a great degree of uncertainty that greatly influences the reliability of the model predictions.
The paper continues as follows. Section 2 is devoted to the general 1D model. The analysis of the stationary solutions is illustrated in Sections 3 and 4. In Section 5, we relate the results to some historical field data for a better interpretation and justification of the model. The conclusions are drawn in the last section.

The Mathematical Model
We consider a geothermal fluid composed of water in a gaseous and liquid state. We denote by x the vertical axis pointing upward. The bottom and the top of the geothermal reservoir correspond to x = L i (L i < 0) and to x = L s (L s < 0), respectively. The ground surface corresponds to x = 0, so that the domain of interest is [L i , L s ]. We assume that the width of the capillary fringe is negligible with respect to the vertical dimensions of the reservoir, and thus that there is a clear separation between the liquid and the vapor phase through a sharp interface, x = s(t), that may vary in time. So, the lower part of the domain L i ≤ x ≤ s (t), is saturated with liquid, while in the upper part of the domain, s (t) ≤ x ≤ L s , there is only vapor. At x = L s the vapor pressure is known, and we denote it as P s . The latter is a measured datum and essentially depends on the exploitation of the basin.
In Table 1, we list the values of the constants involved in the model. These parameters are compatible with those available for some areas of the geothermal field of Larderello, Tuscany [10,14]. Table 1. List of parameters referring to the reservoir of Larderello [10,14]. All dimensionless constants used in the present work are evaluated by using these values. A full list of additional simbols is reported in Table 3.

Symbol Description
Value Units According to the experimental data [37][38][39], the temperature within the reservoir is linear: where γ > 0, constant in time, is the known thermal gradient (see Table 1). We also set L = L s − L i (= 2500 m) .
(2) Remark 1. The given value of γ in Table 1 implies that the temperature at depth L i = −3500 m is 650 K.
Assuming a surface temperature of 300 K, we have a geothermal gradient in the caprock of γ cap = 0.25 K/m. This may be surprising, since this value is about three times the value commonly indicated in textbooks as the typical gradient in a geothermal field (generally 0.08 K/m) and actually requires that the average thermal conductivity of the caprock be rather low-say 0.32 W/(mK). A simple explanation is that the mean gradient of the field is indeed based upon the estimated depth of the K horizon (in the present case, around −5000 m) and the isotherm which identifies the K horizon itself (about 750 K), leading, correctly, to γ mean = 0.08 K/m. On the other hand, the question is much more complex. In a very recently published article [39], the great variability of the thermal conductivity K T of rocks at Larderello and their strong dependence on temperature and, to a much lesser extent, on pressure were clearly emphasized. Reported values vary from less that 2 W/(mK) to more then 4.5 W/(mK) depending on the core sample depth and exploration well. Furthermore, the evaluated geothermal gradient has a great variability with depth, reaching values, for some layers, close to γ cap = 0.2 K/m-not far from the value proposed in Table 1.
The domain s(t) ≤ x ≤ L s is saturated by vapor whose continuity equation is where φ (x) is medium porosity (see Table 1), v is the velocity of the vapor, ρ is its density and d v is the vapor source/sink density rate (with dimensions kg m −3 s −1 ). In general, d v is due to both the recharge vapor (generated from the vaporization of meteoric recharge water), which mixes with the vapor inside the reservoir, and the extraction/re-injection due to external devices [40]. In general, the rate d v depends on the local pressure, or rather on the difference between the external pressure and the local pressure. More specifically, we split d v into two contributions: a recharge generated from the vaporization of recharge water (which may depend on the local pressure) and a point-wise recharge due to N injection/extraction wells whose extraction/injection rates depend on the difference between the local pressure and the re-injection pressure: where the subscripts "met" and "inj" mean, respectively, "of meteoric origin" and "due to injection", P denotes the vapor pressure, δ denotes the Dirac delta, x i is the location of the ith injection head, d is the order of magnitude of the density rate, where J (i) v (dimensionless) gives its dependence on the pressure difference, and P (i) inj is the injection pressure. The vapor is assumed to behave as a perfect gas; i.e., where r is the molar vapor constant (see Table 1). Using the generalized Darcy's law [41], where µ is the vapor viscosity, which depends on both temperature and pressure, and κ (x) is the permeability of the medium (see again Table 1), ρ is as in Equation (5) and g is the gravity acceleration; thus, we rewrite (3) as where T is given by (1). In the liquid region L i ≤ x ≤ s (t), the continuity equation is written as ∂φρ ∂t where ρ is the density of the liquid, v is the liquid velocity and d is the liquid water source density rate (with a dimension of s −1 ), due to both re-injection and natural recharge fluid. The solution to Equation (8) is where, for simplicity, we have neglected v (L i , t), and where the pressure P is the hydrostatic pressure. The continuity of the mass flux across s (t) implies that ρv − ρ v = −ṡ (ρ − ρ) [42]. As ρ ρ , we can approximate this relation withṡ Thus, recalling Equations (6) and (9), we obtaiṅ where which (with a dimension of m s −1 ) is the total recharge rate in the liquid domain. In particular, we assume that we can identify two distinct effects: (i) Recharge due to the meteoric waters, denoted by d ,met , which, as d v,met , may depend on the local pressure; (ii) Recharge due to clusters of re-injection wells which, for simplicity, we assume to be concentrated in M locations.
Thus, setting where, as in (4), the injection rates depend on the difference between the re-injection pressure P inj and the local pressure. Indeed, d ,inj , i = 1, 2, ..., M, gives the injection rate order of magnitude, inj , which is dimensionless, gives the dependence on the pressure difference. We thus have where A second free boundary condition is provided by Clapeyron's law, which gives the relationship between P and T at the liquid/vapor interface. In the range of temperature relevant to our problem, this law takes the following empirical form: where T is in Kelvin, P in Pascal and a and b have values of a = 9.617 × 10 2 Pa, b = 17.35.
Consequently, the second condition on the free boundary is and (11) rewrites asṡ We remark that (16) requires the following condition be fulfilled: in order to obtain a physically consistent solution.
To complete the information, we need to give the pressure at the reservoir top, P (L s , t) = P s (see Table 1) and the initial conditions (18) and (16) Thus, the mathematical problem to be solved is the following: find a time τ > 0, a function , so that, for T (x) given by Equation (1), the following system is fulfilled: , t ∈ (0, τ) , According to the results in [43,44], the problem has one unique solution provided that

Steady-State Solutions
We consider a steady state (meaning that neither d v nor d depend on time) and normalize space, temperature and pressure in the following way: where P c is a characteristic pressure, which still needs to be selected. Recalling Equation (1), we consider T to be an independent variable; namely where We remark that, because of the change of the variable in Equation (20), T = 1 corresponds to the top of the basin while 1 + Γ to the bottom. Next, we introduce as the temperature at the interface, so that Then, we set We also set ρ T = ρ ,o ρ T , µ P, T = µ o µ P, T , κ T = κ o κ T , with µ, ρ and κ being dimensionless; φ T = φ o φ T , with φ o giving the porosity order of magnitude, so that and with where L (i) has been defined in Equation (14) and Setting P c = P s , the steady dimensionless version of problem (19) where η = g rγ ≈ 0.533, f (σ) is given by Equation (24) and and Here and in the following, we omit the " " to keep the notation as clear as possible.
The steady (dimensionless) version of (19), 4 is where with In Section 4, we see that the set of steady solutions (Equation (32)) depends upon a set of n suitable parameters ζ 1 , . . . ζ n . Thus, Equation (32) writes, more precisely, where ζ 1 and . . . ζ n are dimensionless. Our approach consists first of verifying that F = 0 can be made explicit, under some limitations, with respect to a specific parameter-say ζ 1 -that is Then, we look for the existence of stable, steady solutions to the system (Equations (29) and (32)) by varying the other parameters. Thus, for a given steady state σ st , the study of its stability to small perturbations is reduced to investigate the following:

Numerical Simulations
In this section, we analyze some specific cases which represent some of the most significant physical situations. For the sake of simplicity, we take ρ T = µ P, T = κ T = φ T = 1, so that h, given by Equation (25), simplifies to h (σ) = f (σ) /σ. Next, we assume that the injection rate is independent of the difference between the injection pressure and local pressure, and therefore There are many cases, and it is impossible to analyze them all. We therefore focus on these issues: To see the explicit form of Equation (36), it is more convenient to proceed by subsequent steps. Suppose, for a moment, that β v = 0. This will help us to understand (and also visualize graphically) the effect of β v = 0. Proposition 1. For β v = 0, Equation (33) can be written in the form and has a maximum in (σ o , ω 2,o ) ≈ 1.066, 1.29 × 10 −6 .
Proof. The solution to (29) is where In this case, so that Ω (σ) = Ω 1 (σ), with Next, assuming that the liquid recharge occurs only in a part of the domain, D given by (27) becomes for some (given) T 1 < 1 + Γ and F, as defined by Equation (33), can be rewritten as with and Ω is given by Equation (42). Clearly, Equation (44) is meaningful only for 1 < σ < T 1 < 1 + Γ, showing that T 1 remains in the liquid domain. Function Ω 1 monotonically increases and is negative in (1, 1 + Γ). Therefore, we have −∞ < Ω 1 (σ) ≤ Ω 1 (1 + Γ) < 0, meaning that we can express F| β v =0 = 0 in the form of Equation (38). With the values in Table 1, the function F o has a maximum value in (σ o , ω 2,o ) ≈ 1.066, 1.29 × 10 −6 . (38) is a limit case in which the reservoir receives mass only in the liquid layer. This does not seem to be the case for Larderello, where the vapor region is also fed by steam of meteoric origin, as reported by many authors [40,45]. However, the artificial re-injection of geothermal wastewater is practiced in many geothermal fields as a means for reservoir pressure support and cannot be ignored [46][47][48]. We consider these more significant situations in the following three subsections.

Case 1. Recharge Due to Liquid Groundwater Recharge and a Single Vapor Extraction/Injection Well
We now prove that, also in this case, the equation F = 0 can be solved explicitly as in Equation (38) provided that the injection rate and the well depth fulfill suitable conditions.
In this case, D is still given by Equation (43) Integrating Equation (29), 1 once, we obtain where C 1 is an integration constant, θ is the Heaviside function and T o is the injection well depth. Setting z = P 2 /2, Equation (47) becomes whose general solution is where θ(T − T o ) is the Heaviside function. By imposing boundary conditions (29), 2, and (29), 3, we get with A (σ) given by Equation (40), which reduces to Equation (39) whenβ v = 0. In particular, so that ∂P ∂T σ tends to (41) whenβ v → 0 and when T o → 1 (as physically expected). To prove that the equation F = 0 can be solved explicitly, notice first that, because of (49), Ω has a different form, namely where Ω 1 is given by Equation (42) and To use the explicit form in Equation (36), we need Ω(σ, T o ) = 0. It is convenient now to fix an upper limit to the variability of T o -say T o . Then, it can be seen that As a consequence, condition Ω = 0 is fulfilled provided the variability ofβ v is bounded below, namelyβ Remark 3. The critical value (52) tends to −∞ for T o going to 1. Thus, we always consider only values of T o close to but distinct from 1. We notice also that we have to take care of this lower bound only when we consider the case of extraction; that is,β v < 0. For example, for T o = 1.0005, we are allowed to consider onlŷ In the present case, steady states are described by and Ω is given by Equation (50). Recalling Equation (37), we need to study the sign of Figure 1 provides a synthetic picture of the varying influence ofβ v in determining stable/unstable states. We have (linearly) stable steady solutions only for 0 < ω 2 ≤ ω 2,cr and σ ≥ σ cr . Figure 1 shows that by decreasingβ v from positive values (injection) to negative ones (extraction), ω 2,cr increases and σ cr decreases (i.e., the physical free boundary shifts towards the top of the reservoir). This means that the liquid region increases due to the augmented meteoric liquid recharge, as physically expected. Figure 2 shows that solutions represented in Figure 1 satisfy the compatibility condition (18).

Case 2. Groundwater Recharge Occurring in the Whole Domain and a Single Vapor Extraction/Injection Well
The procedure to obtain the explicit form of F = 0 is similar to that in Section 4.1, but function Ω is more complex. Calculations are presented in such a way to emphasize the differences with respect to Section 4.1.

Remark 4.
Recalling what we emphasized in Remark 3, in order to obtain physically meaningful solutions to Equation (38), we need Ω to remain negative when we changeβ v and B v . For fixed T o ∈ (1, T o ) and    We have stable steady solutions only for 0 < ω 2 ≤ ω 2,cr and σ ≥ σ cr . Figure 4 shows that by increasingβ v and keeping B v > 0 fixed, ω 2,cr decreases while σ cr increases. On the contrary, increasing B v and keepingβ v < 0 fixed causes ω 2,cr to decrease and σ cr to decrease. This behavior also agrees with physical expectations. Indeed, increasing the vapor extraction with a fixed vapor recharge (i.e., from the blue curve to the red one) implies that σ cr decreases, and the liquid natural recharge needs to be adjusted at higher rates to maintain stability. On the other hand, a fixed extraction rate and a decreasing natural vapor recharge (i.e., from the red curve to the black one) imply that, to maintain stability, the liquid natural recharge must increase and the critical interface must increase. Figure 5 shows that the compatibility condition (18) is satisfied.

Case 3. Recharge Due to a Single Injection Well Localized in the Liquid Domain
In this case, β v = 0 and d ,met = 0, so that Equation (26) can be rewritten as which is the injection well located in the liquid domain. The solution to (29) is given by Equation (39), and so F as defined by Equation (33) becomes with Ω 1 (σ) given by Equation (42). Thus, setting , the function F acquires this form: where ω 1 has been defined in Equation (45), 1 without any constraints on ω 3 since Ω 1 is always strictly negative. Figure 6 shows the stable steady states in the present case and Figure 7 shows that the compatibility condition (18) is fulfilled.

Sustainability of Equilibria and Re-Injection
The main environmental and industrial issue concerning a geothermal field is the long-term sustainability of its power production capability. Since the middle of the last century, this question has been deeply discussed in the specialized literature [40,45,47,49,50], and the primary role of re-injection has been clearly emphasized. It is interesting to express the equilibrium conditions investigated in this study in dimensional variables and to compare them with some historical field data referring to the Larderello reservoir. Due to the severe approximations considered, such a comparison does not have any specific intention but to meet, at least, the order of magnitude of the most significant and available data.
An estimate of the volume V field of the Larderello reservoir [10,11,14] is about 8 × 10 11 m 3 . Assuming the porosity of Table 1, the volume V fluid of fluid in the reservoir is on the order of 1.6 × 10 10 m 3 . The meteoric recharge from the shallow aquifer surrounding the reservoir is not constant in time; it depends on a variety of factors and can only be roughly estimated [40]. The rate of steam extraction Q extr has changed considerably during the last decades; for example, it was about 300 kg/s in the 1950s and declined to about 150 kg/s in the 1980s [11]. Precise historical series are available from ENEL, the company holding the rights of exploitation of the Larderello geothermal fields [2]. In the 1990s, the natural recharge Q nat from the shallow aquifers in the area was estimated to be about one quarter of the extraction rate Q extr [40]. Due to re-injection, it is likely that the ratio of extraction/natural recharge has not changed very much [10,13]. To solve these ürpnöe,s, we assume Q nat = 0.25 Q extr at any time. Then, we define Consider, now, "case 2" as discussed in Section 4.2. According to the model, the total discharge per unit volume (that is, in kg/(m 3 s)) is given by where d ,met , d v,met and d v,inj are determined by the dimensionless parameters ω 2 ,β v , and B v through (45), 2, (46), 2, and (53), 2, i.e., where, by using data from Table 1, For fixed values ofβ v < 0 and B v > 0, any stable steady state is represented by a point (σ, ω 2 ), with 0 < ω 2 < ω crit and σ > σ crit , on the curve ω 2 = F(σ;β v , B v ). Therefore the real location of the interface depends upon the effective value of ω 2 . At any time during the evolution of the reservoir, we can relate q model with Q field since and exploiting (57) and (58), we have where C 2 is given by Equation (59). To evaluate the left-hand side of Equation (61), we extracted some data from the historical series reported by Pallotta [2]. Notice that Q field V fluid ≈ −10 −8 kg/(m 3 s), so that Q field / (C 2 V fluid ) = O(10 2 ) (see Table 2). Using B v andβ v as free parameters, Equation (61) and Table 2 provide a simple criterion to relate the Larderello field data with the equilibrium interface. For example, selecting Q field / (C 2 V fluid ) = −69 (the 2015 datum),β v = −7.5 × 10 4 , B v = 1.6 × 10 4 , Equation (61) is fulfilled with σ ≈ 1.088, (i.e., s ≈ −2220 m), which is stable because, in case 2, with this choice of parameters, all interfaces σ > 1.018 are stable (see Figure 4). If we keep B v = 1.6 × 10 4 and decreaseβ v up to −1.2 × 10 5 , the new stable position of the interface is around −2560 m. According to Figures 22 and 23 in [10], at a depth of about −2000 m, the fluid pressure in 2010 was 10 bar and the temperature was around 550 K. The real fluid present in the reservoir is not pure water but a complex mixture of water, carbon dioxide and several other chemical components with a great variability of molar fractions. Even for a simple binary water-carbon dioxide system, the estimated saturation pressure would be greatly above the values reported in [10]. It must be emphasized, however, that for a real geothermal fluid, precise indications about the vapor saturation pressure are still lacking, and only some models are supplied [51][52][53]. This suggests, in any event, that for the Larderello field, the vapor/liquid interface must be deeper than −2000 m a.s.l. [30]. All this considered, our model, besides its indisputable simplicity and considerable approximations, locates the interface below the above limit provided that the parameters are selected consistently with the field data. This is certainly a very encouraging result. Nomenclature and other symbol definitions can be found in Table 3.

Conclusions
The main conclusion of our investigation is that stable stationary states of the liquid/vapor interface are possible either with liquid or vapor recharge combined with vapor extraction. The main parameters ω 2 (natural liquid recharge), β v < 0 (vapor extraction) and B v > 0 (vapor meteoric recharge) need to be combined conveniently to maintain a stable interface. Increasing the vapor extraction moves the interface towards the bottom of the reservoir, which at the limit is shown to be increasingly vapor-dominated. Interestingly, stable interfaces can generally be obtained only if ω 2 is less than a critical value; that is, only if the natural liquid recharge is sufficiently large.
The results obtained considering uniform porosity and permeability provide a good qualitative description of the possible steady configurations of a vapor-dominated basin such as Larderello.
An important characteristic of the model is that stationary stable states are compatible with the presence of natural and artificial recharging mechanisms. Such a conclusion agrees with the results of other geological studies [18], despite the drastic simplifications introduced in our model: 1D geometry, an absence of chemicals dissolved in water and uniform geological parameters (which, in a real geothermal field, are far from being homogeneous and may change significantly due to local temperature and pressure). The differential Equation (19), 4, has been proven to be quite satisfactory for modeling the vapor drawdown from a reservoir and for recharging from both aquifers and injection wells. We also think that our analysis is likely to produce physically meaningful results also for equilibrium configurations of water-dominated fields. In conclusion, we believe that, despite the many simplifications made in the model, this study can be an useful step towards improving our understanding of the dynamics of geothermal fields.