1 Introduction

Geothermal energy is one of the most promising alternative power sources. Among the number of indisputable advantages of this type of energy, it can be distinguished that it is renewable; the novel technologies of its extraction do not pollute the environment and it approaches the traditional forms of energy in efficiency. Recent years show a sharp increase in the volume and expansion of spheres for geothermal resources utilization. In a number of countries, geothermal energy technologies become dominant, and the share of geothermal energy in the global energy balance is growing steadily. The geothermal resources of various temperatures are widely used in the power production, district heating and other fields of application.

A lot of monographs, reviews, and papers (Elder 1981; Grant et al. 1982; Tsypkin 1995; Alkhasov et al. 2000; O’Sullivan et al. 2001; Alkhasov 2008; Tsypkin 2009; Alkhasov and Ramazanov 2014) focus on issues related to the extraction and use of geothermal energy. Mathematical simulation of processes associated with the movement of liquids and gases in geothermal systems based on Darcy law and equations of mass and energy balance for two-phase filtration were studied by many authors (Tsypkin 1995; Alkhasov and Ramazanov 2014; Brownell et al. 1977; Verigin and Golubev 1975; Garg and Pritchett 1988; O’Sullivan 1985; Tsypkin and Woods 2004; Richardson and Tsypkin 2004; Ramazanov and Alkhasova 2017; Ramazanov et al. 2017). A brief critical review of these studies is given in the work of Tsypkin (2009), in which the two-phase flow problem with phase transitions in a porous medium is consistently and strictly formulated. In that paper a plane one-dimensional problem is considered which makes it possible to find a relatively simple self-similar solution. However, in the vicinity of the wells, an apparently flat model may not be acceptable, and a radially symmetric model will be more adequate. This assumption has indeed confirmed, and the critical permeability for the flat and radial models, as it turned out, may differ by several orders of magnitude (Alkhasov and Ramazanov 2014).

But we have noted that there is a gap in the mentioned investigations. When simulating the movement of a vapor–water mixture in a reservoir–well system, the authors usually assume that the heat carrier boils in the well, and the influence of the formation is taken into account by the so-called operating characteristic of the reservoir (Alhasov et al. 2000; Chisholm 1983; Kutateladze and Sorokin 1961; Naimanov 1970; Droznin 1971). In the current paper, we turned to poorly studied case of a two-phase water–steam medium flow into a geothermal reservoir–well system when a phase transition occurs in the reservoir itself. In many cases there is a thermodynamic equilibrium of the vapor–water mixture in the geothermal reservoir (Grant et al. 1982). The study of heat and mass transfer in such a mixture is described in Richardson and Tsypkin (2004), O’Sullivan (1980, 1981), Sorey et al. (1980). Tsypkin (2009) provides a diagram of various heat and mass transfer regimes with phase transitions in a geothermal formation within the framework of flat one-dimensional self-similar solutions. A similar diagram for a radially symmetric problem that does not allow a self-similar solution was obtained in Alkhasov and Ramazanov (2014). Ramazanov and Alkhasova (2017) have found the conditions for the existence of heat and mass transfer in the reservoir saturated with vapor–water mixture, in violation of which a region saturated with either pure liquid or pure steam is formed near the well. On the other hand, there are many works where the properties of two-phase flows in pipes and heat exchangers were studied, and in particular in geothermal wells (Chisholm 1983; Tong and Tang 1997; Kakac and Bon 2008; Alhasov et al. 2000).

In contrast to the discussed works, in the current paper we consider the coupled problem of heat and mass transfer in the reservoir around the production well and in the wellbore itself. The properties of the distribution of thermo-mechanical fields and vapor fraction have been studied in these regions up to the wellhead. In this case, the heat transfer between the formation and the well with the surrounding impermeable rock mass was taken into account. The dependence of the results on the specific type of phase permeability for water and vapor when the mixture moves in the reservoir has been investigated. We have considered also the cases of linear and nonlinear dependence of the phase permeability on the water saturation of the reservoir.

In the papers of Churchill (1977), Alkhasov et al. (2009, 2011) the effect of natural convection on heat transfer in a geothermal stratum within the system of a borehole and permeable rock mass was studied. Now, we assume that the conditions are met under which the gravitational effects can be neglected (Tsypkin 2009).

2 Heat and mass transfer in a reservoir

2.1 Formulation of the problem

Suppose that at the initial moment, a geothermal reservoir saturated with a vapor–water mixture has a pressure P0, a saturation temperature T0= F(P0), and water saturation S0. As a result of well operation, its bottom hole pressure drops to P0. This will lead to a change in water saturation in the vicinity of the well; but, according to the calculations, not always in additional vaporization. At the unsteady stage, the water saturation near the well, depending on the parameters of the problem, may increase despite the pressure drop (Ramazanov and Alkhasova 2017). It is required to study the properties of the distribution of thermo-mechanical fields and, above all, the void fraction α and mass dryness fraction x (Chisholm 1983) in the well and the reservoir with various bed parameters.

We write out the set of equations describing the heat and mass transfer in the reservoir filled with the vapor–water mixture being in thermodynamic equilibrium, following Tsypkin (2009):

$$m\frac{\partial }{\partial t}S\rho_{w} + \text{div} \rho_{w} \upsilon_{w} = M,$$
$$m\frac{\partial }{\partial t}\left( {1 - S} \right)\rho_{v} + \text{div} \rho_{v} \upsilon_{v} = - M,$$
$$\frac{\partial }{\partial t}\left( {{\uprho}e} \right)_{m} + \text{div} \left( {\uprho_{w} h_{w}\upupsilon_{w} +\uprho_{v} h_{v}\upupsilon_{v} } \right) = \text{div} \left( {\uplambda_{m} \;\text{grad} \;T} \right),$$
$$\upupsilon_{j} = - \frac{k}{{\mu_{j} }}f_{j} \left( S \right)\;\text{grad} \,P,\quad j = w,v$$
$$\uprho_{w} = \uprho_{w0} \left({1 + \upalpha_{P} \left({P - P_{0}} \right) - \upbeta_{T} \left({T - T_{0}} \right)} \right),$$
(2.1)
$$P = \uprho_{v} RT,\;\ln \frac{P}{{P_{a}}} = A + \frac{B}{T},\;e_{j} = h_{j} - \frac{P}{{\uprho_{j}}}$$
$$dh_{w} = C_{w} dT + \frac{{1-{\upbeta} T}}{{\uprho_{w} }}dP,\;dh_{v} = C_{p} dT,\;de_{s} = C_{s} dT,$$
$$\lambda_{m} = mS\lambda_{w} + m\left( {1 - S} \right)\lambda_{v} + \left( {1 - m} \right)\lambda_{s} ,$$
$$\left( {\rho e} \right)_{m} = mS\uprho_{w} e_{w} + m\left( {1 - S} \right)\uprho_{v} e_{v} + \left( {1 - m} \right)\uprho_{s} e_{s} .$$

In principle, this set of equations can be written as a closed set of two evolutionary equations: the equation of the first order for water saturation and the second order on pressure. Therefore, the initial and boundary conditions can be written as follows:

$$t = 0:\;P = P_{0} ,\;S = S_{0} ,$$
$$r = r_{c} :\;P = P^{0} ,$$
(2.2)
$$r = L:\,P = P_{0} ,\,S = S_{0} .$$

The temperature distribution, including at the boundaries, is determined by the pressure distribution from the Clapeyron–Clausius equation. The distribution of water saturation, in particular its value in the well, is found by solving the problem.

The first two equations in the set (1) are the equations of continuity for water and vapor, respectively. The next is the heat transfer equation, then the generalized Darcy equations for the liquid and vapor phases, the equations of state of the phases, and the Clapeyron–Clausius equation, which determines the phase equilibrium.

We introduce the function Q, proportional to the mass flow rate of the mixture

$$\rho_{w} v_{w} + \rho_{v} v_{v} = - \frac{Q}{r},$$

and rewrite the set (2.1) in a form

$$\frac{\partial P}{\partial r} = \frac{{\upmu_{w}\upmu_{v} }}{{k\left( {\uprho_{w} f_{w}\upmu_{v} +\uprho_{v} f_{v}\upmu_{w} } \right)}}\frac{Q}{r},$$
$$m\frac{\partial}{\partial t}\left[{S\uprho_{w} + \left({1 - S} \right)\rho_{v}} \right] = \frac{1}{r}\frac{\partial Q}{\partial r},$$
(2.3)
$$\begin{aligned} & \frac{\partial }{\partial t}\left[ {\left( {1 - m} \right)\uprho_{s} C_{s} T + m\left( {S\uprho_{w} + \left( {1 - S} \right)\uprho_{v} } \right)\bar{h}_{w} + m\left( {1 - S} \right)\uprho_{v} q} \right] - m\frac{\partial P}{\partial t} + \\ & \quad + {\text{div}}\left[ { - \frac{{Q\bar{h}_{w} }}{r} +\uprho_{v} v_{v} q} \right] = {\text{div}}\left( {\uplambda_{m} \,{\text{grad}}\,T} \right)\,,\,\,\,\,\, \\ \end{aligned}$$
$$\ln \frac{P}{{P_{a} }} = A + \frac{B}{T},$$
$$\uprho_{v} v_{v} = - \frac{{\uprho_{v}\upmu_{w}\, f_{v} }}{{\uprho_{w}\, f_{w}\upmu_{v} +\uprho_{v}\, f_{v}\upmu_{w} }}\frac{Q}{r},\;\uprho_{w} v_{w} = - \frac{{\uprho_{w}\upmu_{v}\, f_{w} }}{{\uprho_{w}\, f_{w}\upmu_{v} +\uprho_{v}\, f_{v}\upmu_{w} }}\frac{Q}{r},$$
$$\uprho_{w} =\uprho_{w0} \left( {1 + \alpha \left( {P - P_{0} } \right) -\upbeta\left( {T - T_{0} } \right)} \right),$$
$$P =\uprho_{v} RT,\;\bar{h}_{w} (P,T) = h_{w} (P,T) - h_{w} (P^{0} ,T^{0} ),$$
$$\uplambda_{m} = mS\uplambda_{w} + m\left( {1 - S} \right)\uplambda_{v} + \left( {1 - m} \right)\uplambda_{s} .$$

Here q is the specific heat of the phase transition.

Note that the third equation in the set (2.3) is obtained from the energy equilibrium equation of the set (2.1) and a mass phase conservation equation, and the equality

$$\bar{h}_{w} \equiv h_{w} - h_{w} (P^{0} ,T^{0} ),\;dh_{w} = d\bar{h}_{w} ,\;h_{v} = h_{w} + q.$$

Taking into account the dependence of pressure on temperature described by the Clapeyron–Clausius equation, we can write

$$\bar{h}_{w} (P,T) = h_{w} (P,T) - h_{w} (P^{0} ,T^{0} ) = C_{w} (T)\left( {T - T^{0} } \right),$$

where through \(C_{w} (T)\) the effective heat capacity is denoted.

The first four equations of the set (2.3) form a closed system of equations regarding four unknowns Q, S, T, P [we assume that the remaining variables are expressed in terms of these four ones with the help of the remaining equalities in the set (2.3)].

2.2 Exact stationary solution

Let us find the stationary solution of the set (2.3). In this case, it will be written in the form

$$P^{\prime } = \frac{{\upmu_{w}\upmu_{v} }}{{k\left( {\uprho_{w}\, f_{w}\upmu_{v} +\uprho_{v}\, f_{v}\upmu_{w} } \right)}}\frac{Q}{r},\;Q = {\text{const}},$$
$$\text{div} \left[ { - \frac{{Q\bar{h}_{w} }}{r} +\uprho_{v} v_{v} q} \right] = \text{div} \left( {\uplambda_{m} \text{grad} T} \right),$$
$$\ln \frac{P}{{P_{a}}} = A + \frac{B}{T},$$
$$\uprho_{v} v_{v} = - \frac{{\uprho_{v} \mu_{w}\, f_{v} }}{{\uprho_{w}\, f_{w}\upmu_{v} +\uprho_{v}\, f_{v}\upmu_{w} }}\frac{Q}{r},\;\uprho_{w} v_{w} = - \frac{{\uprho_{w}\upmu_{v}\, f_{w} }}{{\uprho_{w}\, f_{w}\upmu_{v} +\uprho_{v}\, f_{v}\upmu_{w} }}\frac{Q}{r},$$
$$\rho_{w} = \rho_{w0} \left( {1 + \alpha \left( {P - P_{0} } \right) - \beta \left( {T - T_{0} } \right)} \right),$$
$$P =\uprho_{v} RT,\;\bar{h}_{w} (P,T) = C_{w} \left( {T - T^{0} } \right),$$
$$\uplambda_{m} = mS\lambda_{w} + m\left( {1 - S} \right)\uplambda_{v} + \left( {1 - m} \right)\uplambda_{s} .$$

The boundary conditions are

$$r = r_{c} :\;P = P^{0} ,$$
$$r = L:\;P = P_{0} ,\;S = S_{0} .$$

Integrating the heat transfer equation, we have

$$- \frac{Q}{r}\left( {\bar{h}_{w} + \frac{{\rho_{v} \mu_{w}\, f_{v} q}}{{\rho_{w}\, f_{w} \mu_{v} + \rho_{v} f_{v} \mu_{w} }}} \right) - \lambda_{m} T^{\prime } = \frac{\text{const}}{r}.$$

From here

$$\frac{N(S,T)Q}{{\uplambda_{m} r}}T + T^{\prime } = \frac{{N_{0} Q}}{{\uplambda_{m} r}}T_{0} + \frac{{M_{0} Q}}{{\uplambda_{m} r}},$$
$$P^{\prime } = \frac{{\mu_{w} \mu_{v} }}{{k\left( {\rho_{w} f_{w} \mu_{v} + \rho_{v} f_{v} \mu_{w} } \right)}}\frac{Q}{r},\quad Q = {\text{const}},$$
$$N(S,T) = C_{w} + \frac{{\uprho_{v}\upmu_{w} f_{v} }}{{\uprho_{w} f_{w}\upmu_{v} +\uprho_{v} f_{v}\upmu_{w} }}\frac{q}{T},\;N_{0} = N(S_{0} ,T_{0} ),$$
$$\begin{aligned} &M(S,T) = \lambda_{m} \frac{dT}{dP}\frac{{\mu_{w} \mu_{v} }}{{k\left( {\rho_{w} f_{w} \mu_{v} + \rho_{v} f_{v} \mu_{w} } \right)}},\\& M_{0} = M(S_{0} ,T_{0} ).\end{aligned}$$

Transforming the first equation with the help of the second one, we obtain the following first integral connecting the required quantities:

$$N(S,T)T + M(S,T) = N_{0} T_{0} + M_{0} ,$$
(2.4)
$$N(S,T) = C_{w} + \frac{{\uprho_{v}\upmu_{w} f_{v} }}{{\uprho_{w} f_{w}\upmu_{v} +\uprho_{v} f_{v}\upmu_{w} }}\frac{q}{T},\;N_{0} = N(S_{0} ,T_{0} ),$$
$$M(S,T) = \lambda_{m} \frac{dT}{dP}\frac{{\upmu_{w}\upmu_{v} }}{{k\left( {\uprho_{w} f_{w}\upmu_{v} +\uprho_{v} f_{v}\upmu_{w} } \right)}},\;M_{0} = M(S_{0} ,T_{0} ).$$

Next, we consider the simplest form of the functions for relative permeability of the phases

$$f_{w} (S) = S,\;f_{v} (S) = 1 - S.$$
(2.5)

Using (2.5) and the equation of phase equilibrium, we express from (2.4) water saturation through temperature. Thus, we obtain the following set of equations:

$$\frac{N(S,T)Q}{{\lambda_{m} r}}T + T^{\prime } = \frac{{N_{0} Q}}{{\lambda_{m} r}}T_{0} + \frac{{M_{0} Q}}{{\lambda_{m} r}},$$
$$S(T) = \frac{{\left( {N_{0} T_{0} + M_{0} } \right)\rho_{v} \mu_{w} - \rho_{v} \mu_{w} C_{p} T - \frac{dT}{dP}\frac{{\mu_{w} \mu_{v} \left[ {\left( {1 - m} \right)\lambda_{s} + m\lambda_{v} } \right]}}{k}}}{{\left( {\rho_{w} \mu_{v} C_{w} - \rho_{v} \mu_{w} C_{p} } \right)T + \frac{dT}{dP}\frac{{\mu_{w} \mu_{v} m\left( {\lambda_{w} - \lambda_{v} } \right)}}{k} - \left( {N_{0} T_{0} + M_{0} } \right)\left( {\rho_{w} \mu_{v} - \rho_{v} \mu_{w} } \right)}},$$
$$\begin{aligned} & N(S,T) = C_{w} + \frac{{\rho_{v} \mu_{w} \left( {1 - S} \right)}}{{\rho_{w} S\mu_{v} + \rho_{v} \left( {1 - S} \right)\mu_{w} }}\frac{q}{T},\\&N_{0} = N(S_{0} ,T_{0} ), \end{aligned}$$
$$\begin{aligned} & M(S,T) = \lambda_{m} \frac{dT}{dP}\frac{{\mu_{w} \mu_{v} }}{{k\left( {\rho_{w} S\mu_{v} + \rho_{v} \left( {1 - S} \right)\mu_{w} } \right)}},\\&M_{0} = M(S_{0} ,T_{0} ),\end{aligned}$$
$$\ln \frac{P}{{P_{a}}} = A + \frac{B}{T},\;\rho_{v} = P/\left({RT} \right)a,\;\rho_{w} = \rho_{w0} \left({1 + \alpha \left({P - P_{0}} \right) - \beta \left({T - T_{0}} \right)} \right).$$

This system is integrated in quadratures. Considering the boundary conditions and assuming the radius of the well per unit of length, we get the following exact solution:

$$\int\limits_{{T^{0} }}^{T} {\frac{dT}{{\frac{{N(S(T),T)T\uplambda_{m0} }}{{N_{0}\uplambda_{m} (T)}} - \frac{{\uplambda_{m0} }}{{\uplambda_{m} (T)}}T_{0} - \frac{{M_{0}\uplambda_{m0} }}{{N_{0}\uplambda_{m} (T)}}}} = - \left( {b - 1} \right)\ln r,} \quad b = \frac{{N_{0} Q}}{{\uplambda_{m0} }} + 1,$$
$$S(T) = \frac{{\left( {N_{0} T_{0} + M_{0} } \right)\rho_{v} \mu_{w} - \left( {C_{w} T + q} \right)\rho_{v} \mu_{w} - \frac{dT}{dP}\frac{{\mu_{w} \mu_{v} \left[ {\left( {1 - m} \right)\lambda_{s} + m\lambda_{v} } \right]}}{k}}}{{C_{w} T\left( {\rho_{w} \mu_{v} - \rho_{v} \mu_{w} } \right) - \rho_{v} \mu_{w} q + \frac{dT}{dP}\frac{{\mu_{w} \mu_{v} m\left( {\lambda_{w} - \lambda_{v} } \right)}}{k} - \left( {N_{0} T_{0} + M_{0} } \right)\left( {\rho_{w} \mu_{v} - \rho_{v} \mu_{w} } \right)}},$$
$$\ln \frac{P}{{P_{a} }} = A + \frac{B}{T},\;\uprho_{v} = P/\left( {RT} \right),$$
(2.6)
$$\uprho_{w} =\uprho_{w0} \left( {1 +\upalpha\left( {P - P_{0} } \right) - \beta \left( {T - T_{0} } \right)} \right),$$
$$\begin{aligned} & N(S,T) = C_{w} + \frac{{\uprho_{v}\upmu_{w} \left( {1 - S} \right)}}{{\uprho_{w} S\upmu_{v} +\uprho_{v} \left( {1 - S} \right)\upmu_{w} }}\frac{q}{T},\\&N_{0} = N(S_{0} ,T_{0} ),\end{aligned}$$
$$\begin{aligned} &M(S,T) =\uplambda_{m} \frac{dT}{dP}\frac{{\upmu_{w}\upmu_{v} }}{{k\left( {\uprho_{w} S\upmu_{v} +\uprho_{v} \left( {1 - S} \right)\upmu_{w} } \right)}},\\&M_{0} = M(S_{0} ,T_{0} ),\end{aligned}$$
$$\uplambda_{m} = mS\uplambda_{w} + m\left( {1 - S} \right)\uplambda_{v} + \left( {1 - m} \right)\uplambda_{s} .$$

The mass flow rate of the mixture 2πQ can be found from the condition

$$b = 1 - \frac{1}{\ln L}\int\limits_{{T^{0} }}^{{T_{0} }} {\frac{dT}{{\left[ {\frac{{N(S(T),T)T\uplambda_{m0} }}{{N_{0}\uplambda_{m} (T)}} - \frac{{\uplambda_{m0} }}{{\uplambda_{m} (T)}}T_{0} - \frac{{M_{0}\uplambda_{m0} }}{{N_{0}\uplambda_{m} (T)}}} \right]}}} ,$$
(2.7)
$$Q = \frac{{\uplambda_{m0} }}{{N_{0} }}\left( {b - 1} \right).$$

Here L is the radius of the reservoir area around the production well saturated with vapor–water mixture, on the boundary of which the values of P0, T0, S0are set.

Let us express from (2.4) the dependence of the stratal permeability on the water saturation at the inflow to the well:

$$k(S^{0} ) = \frac{{M(S^{0} ,T^{0} ) - M_{0} }}{{N_{0} T_{0} - N(S^{0} ,T^{0} )T^{0} }}.$$
(2.8)

SuggestingS0 = 1 in this formula, we find (if the right side is positive) one of the critical values of permeability, when the vapor–water mixture at the well condenses completely forming pure water. Similarly, assuming that S0 = 0 in (2.8), in the case of the positivity of the right-hand side, we find the second critical permeability, when the vapor–water mixture at the well completely turns into vapor. If any of these equations has no solution, i.e., the right side of (2.8) is negative, the corresponding phase in its pure form cannot be obtained at any permeability for given values of the remaining parameters.

2.3 Quasi-stationary solution

We will proceed from the non-stationary set of Eq. (2.3). The estimates show that after a time much longer than \(t_{1} = r_{c}^{2} \left( {\rho C} \right)_{m} /\lambda_{m}\), but rather short compared with the characteristic time of the problem under consideration (\(t < < L^{2} \left( {\rho C} \right)_{m} /\lambda_{m}\)), in the vicinity of the well, an almost stationary solution is established in the region that expands with time (Alkhasov 2008).

Denote the stationary solution (2.6), and (2.7) by

$$T_{c} (r,L),\;P_{c} (r,L),\;S_{c} (r,L),\;Q_{c} (L).$$

Here it is emphasized that the solution was obtained in the interval 1 ≤ r ≤ L.

Let RT(t) is the law of the disturbance front motion of the initial temperature distribution in the reservoir. Present a quasi-stationary solution in the form

$$T(r,R_{T}) = \left\{{\begin{array}{*{20}l} {T_{c} (r,R_{T}),} \hfill & {1 \le r \le R_{T} (t),} \hfill \\ {T_{0},} \hfill & {r > R_{T} (t),} \hfill \\ \end{array}} \right.$$
$$P(r,R_{T} ) = \left\{ {\begin{array}{*{20}l} {P_{c} (r,R_{T} ),} \hfill & {1 \le r \le R_{T} (t),} \hfill \\ {P_{0} ,} \hfill & {r > R_{T} (t),} \hfill \\ \end{array} } \right.$$
(2.9)
$$\begin{aligned} & S(r,R_{T} ) = \left\{ {\begin{array}{*{20}l} {S_{c} (r,R_{T} ),} \hfill & {1 \le r \le R_{T} (t),} \hfill \\ {S_{0} ,} \hfill & {r > R_{T} (t),} \hfill \\ \end{array} } \right.\\&Q = Q_{c} (R_{T} ),\;b = b_{c} (R_{T} ). \end{aligned}$$

In order to find the front motion law RT(t), we multiply the heat transfer equation in the system (2.3) by \(r\) and integrate it over a segment \([1,L]\). For the right side, we take into account that the derivative of the temperature in the unperturbed region is equal to zero, and in the left, according to (2.9), we replace the upper limit \(L\) by the front radius \(R_{T} (t)\). As a result, we get

$$\dot{R}_{T} \int\limits_{1}^{{R_{T} }} {\frac{\partial }{\partial T}\left[ {\left( {1 - m} \right)\uprho_{s} C_{s} T + m\left( {S\uprho_{w} + \left( {1 - S} \right)\uprho_{v} } \right)\bar{h}_{w} + m\left( {1 - S} \right)\uprho_{v} q - mP} \right]} \frac{\partial T}{{\partial R_{T} }}rdr - \left[ {Q\bar{h}_{w} -\uprho_{v} v_{v} rq} \right]\left| \begin{aligned} r = R_{T} \hfill \\ r = 1 \hfill \\ \end{aligned} \right. = - \left[ {\uplambda_{m} rT^{\prime } } \right]_{r = 1} ,$$
(2.10)
$$\uprho_{v} v_{v} = - \frac{{\uprho_{v}\upmu_{w} f_{v} }}{{\uprho_{w} f_{w}\upmu_{v} +\uprho_{v} f_{v}\upmu_{w} }}\frac{Q}{r}.$$

Thus, Eqs. (2.9), (2.10) solve the problem of finding a quasi-stationary solution of the system (2.3).

In conclusion, we note that, using solution (2.9), Eq. (2.10) for the perturbation front can be written in a more convenient for calculations form

$$R_{T} \dot{R}_{T} = 2\gamma (R_{T} ),$$
(2.11)
$$\gamma (R_{T} ) = \frac{1}{2}\left( { - \left[ {Q\bar{h}_{w} -\uprho_{v} v_{v} rq} \right]\left| \begin{aligned} r = R_{T} \hfill \\ r = 1 \hfill \\ \end{aligned} \right. +\uplambda_{m} T^{\prime } (1)} \right)/\left( {F_{1} (T_{0} ) - \int\limits_{{T^{0} }}^{{T_{0} }} {\frac{{F_{1} (T)F_{2} (T)r^{2} \left( {2\ln r + 1} \right)}}{{R_{T}^{2} \ln R_{T} \left( {b - 1} \right)}}} dT} \right),$$
$$T^{\prime } (1) = \frac{b - 1}{{F_{2} (T^{0} )}},\;\uprho_{v} v_{v} = - \frac{{\uprho_{v}\upmu_{w} f_{v} }}{{\uprho_{w} f_{w}\upmu_{v} +\uprho_{v} f_{v}\upmu_{w} }}\frac{Q}{r},$$
$$F_{1} (T) = \left( {1 - m} \right)\uprho_{s} C_{s} \left( {T - T^{0} } \right) + m\left( {S\uprho_{w} + \left( {1 - S} \right)\uprho_{v} } \right)\bar{h}_{w} + m\left( {1 - S} \right)\uprho_{v} q - m\left( {P - P^{0} } \right),$$
$$\begin{aligned} &F_{2} (T) = - \frac{1}{{\frac{{N(S(T),T)T\uplambda_{m0} }}{{N_{0}\uplambda_{m} (T)}} - \frac{{\uplambda_{m0} }}{{\uplambda_{m} (T)}}T_{0} - \frac{{M_{0}\uplambda_{m0} }}{{N_{0}\uplambda_{m} (T)}}}},\\&\bar{h}_{w} = C_{w} (T - T^{0} ).\end{aligned}$$

2.4 Discussion of the results

A high-temperature geothermal formation can often be saturated with a vapor–water mixture in a thermodynamic equilibrium (Grant et al. 1982). In (Tsypkin 2009; Alkhasov and Ramazanov 2014), the linearized models for extracting a vapor–water mixture from a reservoir were considered, respectively, in the flat one-dimensional and radially symmetric cases. These studies have shown that the linearization of the equations leads to a strict limitation of the domain of parameters, where the corresponding model is applicable. Beyond this parameter area, the linear models result in physical contradictions (in negative values water saturations near the well or values water saturations greater than unity ones) and can significantly distort the real values of the sought-for quantities (Alkhasov and Ramazanov 2014). The nonlinear model significantly expands the range of parameters where the considered heat and mass transfer regime is realized. However, its applicability is limited by the positivity condition of the right side of the Eq. (2.8). Otherwise, the whole area under consideration cannot be in a two-phase state, and it is necessary to take into account the presence of single-phase parts in it.

The formulas (2.6) and (2.7) give the exact stationary solution of the problem in the range of distances from the well 1 ≤ r ≤ L (L is the radius of the reservoir region around the well, on the boundary of which unperturbed initial formation conditions are specified). As calculations show, over a sufficiently short time compared with the characteristic one in a certain expanding with time area around the well a quasi-stationary field distribution is established (Alkhasov and Ramazanov 2014). Outside this region, the thermo-mechanical fields have almost unperturbed values. Comparison of the quasi-stationary solution obtained in this way with the exact solution for the unsteady heat conduction equation around the well showed very good agreement of the results (Alkhasov and Ramazanov 2014). In (Alkhasov et al. 2009), a comparison of a quasi-stationary solution is given obtained by numerical methods with the numerical solution of non-stationary equations for the problem of convective heat exchange of a vertical well with an aquifer. And in this case, the quasi-stationary method showed a combination of simplicity with high efficiency and good accuracy. Using this approach, an analytical quasi-stationary solution has been obtained in the current paper.

Distribution of temperature, pressure, and water saturation related to the unperturbed values, depending on the ratio of the distance from the well to its radius is shown in Fig. 1 (with indication of parameters). Attention is drawn to the essentially nonlinear nature of changes in water saturation near the well.

Fig. 1
figure 1

Distribution of reduced thermo-mechanical fields around a production well at P0 =15 MPa, P0 =0.45 MPa, S0 =0.6, k =0.6·10−16 m2: 1P/P0,2S/S0, 3Fi= T/T0

Figure 2 illustrates the different nature of changes in water saturation for different values of formation permeability. At low its values, a sharp increase in the water saturation of the heat carrier occurs near the well. As can be seen from Fig. 2 (line 4), if the permeability is sufficiently small, the water saturation becomes more than one, which, naturally, has no physical meaning. This means that the considered regime of heat and mass transfer is disturbed and it is necessary to expand the model by introducing a region with clean water around the well, coupled with the area of vapor–water mixture. On the other hand, at high reservoir permeability, a sharp increase in the vapor fraction of the heat carrier near the well occurs, i.e., water saturation is reducing. Within the framework of this model, if the initial water saturation of the reservoir is sufficiently small, the water saturation near the well for a sufficiently high permeability or low pressure in the well becomes negative. Thus, with these parameters, the considered regime of heat–mass transfer is also violated and the model must be expanded by introducing an additional region of clean vapor around the well, which is interfaced with the region of the vapor–water mixture.

Fig. 2
figure 2

Water saturation distribution around a production well at P0 = 15 MPa, P0 = 0.45 MPa, S0 = 0.6: 1 k = 10−15м2, 2 10−16, 3 3·10−17, 4 2·10−17

This is confirmed by Fig. 3 representing the diagram of various modes of heat and mass transfer. If the pressure in the well is greater than the critical one (with other parameters fixed), forming the pure vapor near the well is impossible under any permeability.

Fig. 3
figure 3

Critical values of permeability depending on the ratio of pressure in the well to reservoir pressure, k0 = 10−19 m2, P0 =15 MPa:S0 =0.2 (solid lines), S0 =0.25 (dashed lines)

Figure 4 shows the dependences of water saturation near the well on the initial water saturation of the formation at different values of permeability. The presented results confirm the possibility of violation of the considered heat and mass transfer regime and the formation of areas saturated with clean water or pure vapor near the well.

Fig. 4
figure 4

Dependence of water saturation near the well on the value of water saturation far from it, i.e., from the initial water saturation of the reservoir, P0 =15 MPa, P0 =0.45 MPa: 1 k =10−14 m2, 2 10−15, 3 10−16, 4 2·10−17

The obtained quasi-stationary solution depends on time through the effective radius RT (t), which determines the law of motion of the temperature disturbance front, and, accordingly, pressure and water saturation. The dependence RT(t) is obtained from the non-stationary heat transfer equation in the set (2.3). It could seem that in the set (2.3) another non-stationary equation (the equation of continuity for the vapor–water mixture) remains unused, but it is notso. The explanation to (2.3) says that the heat transfer equation was previously transformed using the continuity equation for the mixture and it was taken into account when finding RT(t).

Figure 5 illustrates the evolution of the temperature distribution around the well as the disturbance front radius RT(t) moves, obtained in the framework of the quasi-stationary solution.

Fig. 5
figure 5

Quasi-stationary temperature distribution for different values of the temperature disturbance front radius P0 =15 MPa, P0 =0.45 MPa, k =10−16 m2: 1RT= 3 m, 2 5 m, 3 10 m, 4 100 m

The dependence RT(t) is determined by the numerical solution of the nonlinear ordinary differential Eq. (2.11). However, as the calculations showed, the function \(\gamma (R_{T} )\) in (2.11) depends on its argument rather weakly. Assuming \(\gamma (R_{T} ) \approx \gamma (\bar{R}_{T} ) = {\text{const}}\), where \(\bar{R}_{T}\) is a certain mean value, and integrating Eq. (2.11), we obtain the root law for RT(t)

$$R_{T} = R_{T} (0) + 2\sqrt {\gamma (\bar{R}_{T} )t} .$$
(2.12)

The pattern of the disturbance front RT(t) movement obtained by numerical integration of Eq. (2.11) with two different values of the initial water saturation is illustrated in Fig. 6. The dashed lines show approximate dependences found by formula (2.12). It is clear that the root law (2.12) sufficiently accurately determines the law of motion of the disturbance front.

Fig. 6
figure 6

Motion law for the front of the perturbation of the fields initial distribution T0, P0, S0 in reservoir (RT (m), t (day)), P0 =15 MPa, P0 =0.45 MPa, k = 10−16м2:1S0 =0.6, 2 0.2

3 Heat and mass transfer in the system of well and geothermal reservoir

The calculation for the well was carried out based on the equations from Chisholm (1983).

$$- \frac{dP}{dz} = F_{fr} + \frac{g}{{\text{v}_{mix} }} + G^{2} \frac{{d\varvec{v}_{e} }}{dz},$$
$$\frac{d}{dz}\left[ {xh_{v} + \left( {1 - x} \right)h_{w} } \right] + \frac{d}{dz}\left[ {x\frac{{u_{\varvec{v}}^{2} }}{2} + \left( {1 - x} \right)\frac{{u_{w}^{2} }}{2}} \right] + g = - \frac{{\lambda_{s} }}{{Gr_{c}^{2} }}{\text{Nu}}\left( {T - T_{s0} } \right),$$
(3.1)
$$P =\uprho_{v} RT,\;\ln \frac{P}{{P_{a} }} = A + \frac{B}{T},$$
$$F_{fr} = \frac{{\lambda G^{2} \varvec{v}_{w} }}{{4r_{c} }}{\varphi }_{fr.w.0}^{2} .$$

Here, \({\varphi }_{fr.w.0}^{2}\) is the semi-empirical function that takes into account the two-phase nature of the medium and depends on x and the physical properties of the phases (Chisholm 1983); \(\lambda\) is the coefficient of friction, which was determined by the Churchill formula (Churchill 1977; Chisholm 1983); ve is the effective specific volume defined as

$$MF = G^{2} \varvec{v}_{e} .$$
(3.2)

From (3.2) you can get the formula

$$\varvec{v}_{e} = \frac{1}{{\uprho_{e} }} = \left[ {x\varvec{v}_{v} + K(1 - x)\varvec{v}_{w} } \right]\left[ {x + \frac{1 - x}{K}} \right],\quad K = \frac{{u_{v} }}{{u_{w} }}.$$
(3.3)

Here K is the velocity coefficient (the ratio of the true phase velocities).

The solutions of the problem for the reservoir, Eqs. (2.12.2), and for the well, Eqs. (3.13.3), were “sewn” at the bottom hole. Namely, it was assumed that when the heat carrier passes from the formation to the well, pressure and temperature, as well as mass and energy flax, remain the same. In this case, x and \(\alpha\) can have a discontinuity. Note that \(\alpha = 1 - S\) in the reservoir, but in the well it is expressed in terms of x, K and specific volumes of the phases.

When solving the Eqs. (3.13.3), the simplest semi-empirical formula was used (Chisholm 1983):

$$K = \left( {\frac{{\varvec{v}_{\hom } }}{{\varvec{v}_{w} }}} \right)^{{\frac{1}{2}}} ,\;\varvec{v}_{\hom } = x\varvec{v}_{v} + (1 - x)\varvec{v}_{w} .$$
(3.4)

In the stationary and quasi-stationary cases when heat exchange with surrounding rock mass is not taken into account, the problem (1–2) in the reservoir admits the exact solution obtained in (Ramazanov and Alkhasova 2017). We would not repeat this solution here, but it was used as a test one. In the absence of heat transfer, an analytical solution of the problem (3–6) in the well is permitted also to be in quadratures. Indeed, the equation of energy balance in the set (3) can be solved analytically with respect to x. At the same time, a quadratic equation is obtained for a homogeneous theory, and a cubic equation for a non-homogeneous one. Both solutions give similar values. Expressing x explicitly as a function of pressure (temperature is also expressed in terms of pressure using the Clapeyron–Clausius equation) and substituting this expression into the first equation of the set (3), we obtain an ordinary first-order differential equation for pressure that is integrated in quadratures.

3.1 The problem with taking into account the heat exchange

In the case of heat exchange considering, the problem, both in the reservoir and in the well, was solved numerically by the difference method under the assumption that the heat and mass transfer regime is quasi-stationary. The calculations in the reservoir were performed for the functions of relative permeability of the phases of the simplest type

$$f_{w} (S) = S,\;f_{v} (S) = 1 - S.$$

For thermo-physical coefficients, the following values were used:

\(A = 12.512;\;B = -\,4611.73;\;\uprho_{w} = 800;\;\uprho_{s} = 2000;\;\uplambda_{w} = 0.58;\;\lambda_{s} = 2;\;\lambda_{v} = 0.02;\;\upmu_{v} = \,1.72 \cdot 10^{ - 5} ;\) \(\upmu_{v} = 1.72 \cdot 10^{ - 5} ;\) \(C_{w} = 4.2 \cdot 10^{3} ;\) \(C_{s} = 0.9 \cdot 10^{3} ;\) \(C_{p} = 2.3 \cdot 10^{3} ;\) \(q = 1.85 \cdot 10^{6} ;\) \(R = 461;\) \(H = 2000.\)

Figure 7a gives the characteristic distributions of void fraction \(\alpha\) and mass vapor fraction \(x\) in the well and formation. Through \(\alpha_{0}\) we denote the initial uniform distribution of \(\alpha\) in the reservoir. On the abscissa, the dimensionless length coordinate \(\xi\) is plotted, defined as

Fig. 7
figure 7

a Distribution of void fraction \(\alpha\) and mass dryness fraction \(x\) in the well \(\left( {0 < \xi < 10} \right)\) and reservoir \(\left( {\xi > 10} \right)\) at \(\alpha_{0} = 0.1\) and \(k = 0.3 \cdot 10^{ - 16} ,\,\,10^{ - 16} ,\,\,9 \cdot 10^{ - 16} \,\,(1 - 3)\). b Distribution of void fraction \(\alpha\) and mass dryness fraction \(x\) in the well \(\left( {0 < \xi < 10} \right)\) and formation \(\left( {\xi > 10} \right)\) at \(\alpha_{0} = 0.5\) and \(k = 0.3 \cdot 10^{ - 16} ,\,\,10^{ - 16} ,\,\,9 \cdot 10^{ - 16} \,\,(1 - 3)\)

$$\upxi = \left\{ {\begin{array}{*{20}l} {0 <\upxi < 10} \hfill & {in \, the \, well,} \hfill \\ {\upxi > 10} \hfill & {in \, the \, reservoir} \hfill \\ \end{array} .} \right.$$

The coordinate \(\xi\) grows from the mouth to the bottom and grows further in the reservoir as it moves away from the well.

Three options corresponding to three different values of permeability of the reservoir have been considered. These curves show that near the well in the reservoir, the void fraction has a sharp maximum, which is associated with a significant non-linearity of the problem. With increasing permeability this maximum increases; it becomes steeper and approaches the well. At the same time, the vapor fraction at the entrance to the well also increases with its maximum. This effect of the presence of a maximum of vapor fraction (minimum of water fraction) in a production well was previously shown in (Richardson and Tsypkin 2004; Ramazanov and Alkhasova 2017), respectively, for plane self-similar and radially symmetric problems. In areas of high void fraction, i.e., in the vicinity of the indicated maximum, the possibility of salt precipitation should be expected, worsening the reservoir permeability in this zone (Richardson and Tsypkin 2004; Ramazanov and Alkhasova 2017). As Fig. 7 shows, mass vapor fraction \(x\) from a qualitative point of view behaves similarly \(\upalpha\).

The pressure drop at the bottom hole when operating the well leads to a drop in temperature. This, in turn, with intensive extraction of the heat carrier, leads to the thermal boundary layer forming. As a result, along with the convective component, a substantial part of the heat passes from the reservoir to the well by heat conduction. This results in additional vapor generation in the area close to the bottom of the well and accordingly to the fact that the mass vapor fraction increases abruptly. Then \(x\) (and \(\upalpha\)) almost linearly grows in the wellbore from the bottom to the mouth, due to decrease in the saturation pressure (an increase in the specific gas volume). Figure 7b represents the same curves as Fig. 7a for the case when the initial vapor fraction in the reservoir is uniform and equal \(\upalpha_{0} = 0.5\), while in Fig. 7a \(\upalpha_{0} = 0.1\). In this case, the magnitude of the jump \(\Delta x\) with permeability rises sharply (Fig. 8).

Fig. 8
figure 8

Distribution of \(\alpha\) and \(x\) in the well \(\left( {0 < \xi < 10} \right)\) and reservoir \(\left( {\xi > 10} \right)\) at \(\alpha_{0} = 0.1\). Solid line taking into account heat exchange with surrounding rocks, dashed line without taking into account heat transfer

Figure 9 demonstrates the qualitative nature of heat exchange influence of the reservoir–well system with surrounding rock mass on the distribution of \(\upalpha\) and \(x\). Our estimates show that with the characteristic mass flow rates of heat carrier, the effect of heat exchange with rock mass is relatively small. In this figure, for clarity, the Nusselt numbers were noticeably overestimated compared to the real ones.

Fig. 9
figure 9

Distribution of reduced temperature and pressure accounting (dashed) and excluding (line) heat exchange with surrounding rock in the well \(\left( {0 < \xi < 10} \right)\) and reservoir \(\left( {\xi > 10} \right)\)

In Fig. 9 we can see that heat exchange leads to an increase in the pressure (temperature) gradient in the reservoir in the vicinity of the well. In this case, according to the calculations, the mass flow rate of the mixture decreases, and for vapor itself it rises.

Figure 10 illustrates the nature of the distribution of void fraction in the well and in its surroundings for an initially homogeneous reservoir with different initial vapor saturation. In this case, the value of permeability is chosen for clarity; for larger values, the picture does not change qualitatively. We can see that a sufficiently big difference in the initial distribution results in a relatively small variation of this value in the well, especially at its mouth.

Fig. 10
figure 10

Distribution of void fraction \(\alpha\) in the well \(\left( {0 < \xi < 10} \right)\) and reservoir \(\left( {\xi > 10} \right)\) at \(k = 0.3 \cdot 10^{ - 16}\) and \(\alpha_{0} = 0.1,\,\,0.5,\,\,0.9\)\(\,(1 - 3)\)

Figure 11 shows similar curves for mass dryness fraction \(x\).

Fig. 11
figure 11

Distribution of the mass vapor fraction \(x\) in the well \(\left( {0 < \xi < 10} \right)\) and reservoir \(\left( {\xi > 10} \right)\) at \(k = 0.3 \cdot 10^{ - 16}\) and \(\alpha_{0} = 0.1,\,\,0.5,\,\,0.9\)\(\,(1 - 3)\)

Figure 12a, b show the dependence of void fraction and mass dryness fraction at the wellhead, the bottom hole, and in the reservoir at the well inlet, respectively, on the reduced pressure at the bottom hole (1–3) for two different values of permeability. First of all, from these Figures it follows that the parameters \(\upalpha\) and \(x\) have a gap at the entry of the heat carrier from the reservoir into the well. Namely, these parameters abruptly increase. Further, when moving along the wellbore, these values continue to grow. Figure 12b presents that at relatively low bottom pressure, the gap at the entrance to the well \(\Delta\upalpha\) is small, but \(\Delta x\) is significant. As this pressure increases, \(\Delta\upalpha\) increases and \(\Delta x\) decreases. These effects are determined by the ratio of convective and conductive heat transfer from the reservoir to the well.

Fig. 12
figure 12

a Dependence of void fraction and mass dryness fraction at the wellhead, the bottom of the well, and in the formation at the entrance to the well from the reduced pressure at the bottom of the well, respectively (1–3) at \(k = 10^{ - 16}\). b Dependence of void fraction and mass dryness fraction at the wellhead, the bottom of the well, and in the formation at the entrance to the well from the reduced pressure at the wellhead, respectively (1–3) at \(k = 10^{ - 15}\)

As for the dependence of these parameters on the reduced pressure at the bottom of the well, it is not linear and depending on the parameters (the permeability, in particular) can have a various character.

From the energy point of view, the values of the parameters \(\upalpha\) and \(x\) at the wellhead are of the greatest interest. The curves for these parameters are shown in Fig. 13. At least in the considered range of permeability and pressure and at the wellhead, they change relatively little; although in the wellbore, and especially in the reservoir, these parameters change significantly. However, the total vapor consumption itself at the mouth is proportional to the mass flow rate of the heat carrier and, as Fig. 14 shows, will strongly (approximately directly proportional) depend on the permeability.

Fig. 13
figure 13

Dependence of the volumetric steam content and mass flow rate of steam at the wellhead on the reduced pressure at the bottom of the well for different permeability of the reservoir \(k = 0.3 \cdot 10^{ - 16} ;\,\,\,\,\,10^{ - 16} ;\,\,\,\,\,9 \cdot \,10^{ - 15}\) (1–3)

Fig. 14
figure 14

Dependence of the reduced total vapor discharge in the well on the reduced pressure at the wellhead at different permeability values \(k = 0.3 \cdot 10^{ - 16} ;\,\,\,\,\,10^{ - 16} ;\,\,\,\,\,9 \cdot \,10^{ - 15}\) (1–3)

3.2 Calculations for linear and nonlinear functions of the relative phase permeability in the reservoir

The problem, both in the reservoir and in the well, was solved by a numerical difference method under the assumption that the heat and mass transfer regime is stable. Heat exchange with the impermeable rock mass surrounding the system is not taken into account; the contribution of heat transfer is discussed in Sect. 3.1 (Alkhasov et al. 2017). In the previous section, reservoir calculations were performed for the linear function of the relative permeability of the phases. This section provides a comparative analysis for two types of phase permeability:

  1. 1.

    Linear phase permeability

    $$f_{w} (S) = S,\;f_{v} (S) = 1 - S.$$
  2. 2.

    Non-linear phase permeability

$$f_{w} (S) = \left\{ {\begin{array}{*{20}l} {0,} \hfill & {0 \le S \le S_{wk} } \hfill \\ {\left[ {\left( {S - S_{wk} } \right)/\left( {1 - S_{wk} } \right)} \right]^{7/2} ,} \hfill & {S_{wk} \le S \le 1} \hfill \\ \end{array} } \right.$$
$$f_{v} (S) = \left\{ \begin{aligned} \left[ {\left( {S_{vk} - S} \right)/S_{vk} } \right]^{7/2} ,\,\,\,\,\,\,0 \le S \le S_{vk} \hfill \\ 0,\,\,\,\,\,\,\,S_{vk} \le S \le 1 \hfill \\ \end{aligned} \right.$$

For thermo-physical coefficients, the same values were used as in the Sect. 3.1.

4 Discussion of the results

Two variants, namely linear and nonlinear, of phase permeability have been investigated for two values of permeability in the reservoir \(k = 10^{ - 14}\) (Fig. 15a) and \(k = 10^{ - 16}\) (Fig. 15b). The comparison of the curves shows, in general, the agreement of the qualitative distribution of the parameters for both types of phase permeability. From a quantitative point of view, this agreement can only be called satisfactory.

Fig. 15
figure 15

Distribution of \(\alpha\) and \(x\) in the well \(\left( {0 < \xi < 10} \right)\) and the reservoir \(\left( {\xi > 10} \right)\) with \(\alpha_{0} = 0.2\)a at \(k = 10^{ - 14}\); b at \(k = 10^{ - 16}\); \(\alpha_{1} ,\,x_{1}\) is related to linear phase permeability; \(\alpha_{2} ,\,x_{2}\) is related to nonlinear phase permeability

It follows from the Figures that \(\upalpha\) can have a maximum in the neighborhood of the well. It is associated with a substantial nonlinearity of the problem even for linear phase permeability.

Figure 16 shows that heat transfer leads to an increase in the pressure gradient (and temperature one) in the rock around the well.

Fig. 16
figure 16

Distributions of the reduced temperature \(T/T_{0}\) and the reduced pressure \(P/P_{0}\) in the well \((0 < \xi > 10)\) and in the formation \((\xi > 10)\) at \(\upalpha_{0} = 0.2\): a at \(k = 10^{ - 14}\); b at \(k = 10^{ - 16}\); 1 and 3 refers to linear phase permeability; 2 and 4 refers to nonlinear phase permeability

As the computations show, when increasing the initial void fraction \(\upalpha_{0}\) in the reservoir, the pattern is until some time qualitatively similar to the described one. However, starting from a certain critical value, pure vapor is formed around the well for given values of the other parameters. These issues are discussed in detail in (Alkhasov et al. 2017). Figure 17 illustrates that at sufficiently high \(\upalpha_{0}\) when the pressure and the associated temperature fall below the certain critical value, the pure vapor will form at the well’s mouth. The indicated critical values of \(\alpha\), as shows Fig. 17, are different for linear and nonlinear phase permeability, but the overall picture is qualitatively identical.

Fig. 17
figure 17

Distribution of void fraction in the reservoir in a stationary filtration regime: 1\(\alpha_{0} = 0.8\) (linear phase permeability); 2\(\alpha_{0} = 0.6\) (nonlinear phase permeability)

5 Conclusions

  • The coupled problem has been considered of heat and mass transfer in the vicinity of the production well and in the wellbore itself when extracting the vapor and water mixture. The reservoir is believed to be saturated with vapor–water mixture initially. The problem is solved by the numerical–difference method.

  • Near the well in the reservoir, the distribution of void fraction \(\upalpha\) in a certain range of reservoir permeability has a maximum, which can result in precipitation of salts and deterioration of permeability.

  • In a general case, when the heat carrier passes from the reservoir to the well, the distribution of \(\upalpha\) and mass dryness fraction \(x\) undergo a gap and increase abruptly; this is associated with the forming a thermal boundary layer adjacent to the well. In the wellbore, these values grow from the bottom to the mouth approximately linearly.

  • The dependence of \(\upalpha\) and \(x\) at the wellhead on the parameters of the problem was studied. The qualitative character of the effect of reservoir and well heat exchange with impermeable rock on the structure of thermo-mechanical fields has been shown.

  • The pressure and temperature distributions are found to be rather weakly dependent on the specific type of phase permeability of water and vapor in the reservoir. However, the distribution of vapor fraction is in a good agreement only qualitatively. The quantitative agreement can be called only satisfactory. It means that, at least in rough preliminary calculations, simple linear dependences of phase permeability can be used.

  • The paper first considers the radially symmetric flow of a vapor–water mixture in a geothermal formation taking into account phase transitions. The corresponding mathematical model admits the exact solution of stationary and quasi-stationary nonlinear problems. The distribution of thermo-mechanical fields around the production well is presented.

  • There are critical permeability determining the conditions for the existence of the specified way of heat and mass transfer. We constructed a corresponding diagram.

  • The law of motion is found of the disturbance front for the initial distributions of temperature, pressure, and water saturation in the reservoir. The nature of these fields evolution has been shown.