The Problem of Heat and Mass Transfer in Melting Snow and Ice Cover Taking into Account The External Heat Flow

1. Introduction

Recently, various snow cover models have become increasingly popular. They are used to solve problems related to avalanche movement, the contribution of snow cover to runoff formation in river basins, and the spread of pollutants in melting snow. The thermal balance and regime of the earth’s surface are shaped by a powerful natural factor—solar radiation flux. Accounting for solar radiation flux is important during snowmelt, as it allows for more accurate predictions of the thermal regime of the snow-ice cover and the thermodynamic processes within it. Models typically use surface temperature as a boundary condition, which ignores important physical processes occurring in the uppermost layer of snow. Solar radiation flux significantly influences the thermal balance of the snow and ice surface, where intense processes such as melting, evaporation, and other processes occur in the near-surface layer under the influence of radiation, exerting a decisive influence on the dynamics of the underlying snow cover and the thermal regime. Snow and ice are light-scattering media over a fairly wide range of the solar radiation spectrum, so their optical properties must be taken into account.

One way to account for heat flux in snowmelt modeling is to specify it for the boundary layer on the snow surface. This method for determining heat flux is presented in the work [1], where a condition for heat exchange with the atmosphere was specified at the upper boundary of the soil (or snow cover), while a geothermal heat flux was introduced at the lower boundary. The boundary condition at the snow-soil boundary determined the equality of temperatures and heat fluxes. The work included calculations to assess the influence of various factors on the snow surface temperature and the development of a method for determining the thermal stability of the snow cover and the effective thermal conductivity of snow based on the temperature of the upper soil layer.

Another method is the partial reflection of shortwave radiation from the snow surface and its partial penetration into the snow layer. In many studies, this process is presented in the form of Bouguer’s law [2], which takes into account such important optical properties as the scattering coefficient, the attenuation coefficient, and albedo. The albedo for snow, where the attenuation coefficient reaches 30–60 m⁻¹, reaches limit values close to 100 [2]. Thus, freshly fallen snow is characterized by the highest albedo value A>90[2]. A good description of the reflectivity of snow cover is presented in the work [3], which also proposes a model for calculating the albedo based on air temperature and total solar radiation data, which can be used in the absence of observations. In [4], a one-dimensional multiphase model of snow cover was considered, in which the change in heat content for each phase is described by a heat flux equation accounting for convection, thermal conductivity, solar radiation, and phase transitions.

The heat flux in [5] to describe heat and moisture exchange in snow was defined based on P.P. Kuzmin’s method. The flux accounted for direct and diffuse shortwave radiation, longwave radiation from the atmosphere, longwave radiation from snow, turbulent heat exchange, heat exchange during evaporation, and heat input with liquid precipitation. To model the hydrothermal regime of snow cover at subzero temperatures, a heat source was used due to the absorption of shortwave radiation flux within the snowpack, based on the Bouguer-Lambert law. Problems of heat and mass transfer in melting snow and ice cover without taking into account an external heat source were considered in works [6,7,8,9,10,11].

In the paper [12], a heat and moisture exchange model was developed that calculates all components of the heat and water balance on land, as well as state variables (effective landscape surface temperature, soil temperature, soil moisture content, amount of frozen water in the soil, albedo, etc.). Here, the transfer of shortwave radiation within the snowpack is taken into account in the temperature equation, the radiation intensity varies with depth according to the Bouguer-Lambert law, and an empirical relationship is proposed for determining the attenuation coefficient. For the albedo coefficient, the relationship from [13] was used, taking into account the aging of the upper snow layer.

The process of radiative snowmelt is described in detail in the paper [14]. Experimental studies were conducted to examine the influence of the radiation component on the snowmelt process and changes in snow cover albedo. The paper [15] compares two parameterizations for solar radiation transfer in sea ice. The first assumed exponential radiation decay with constant transmittance and attenuation coefficients. The second employed a two-layer scheme for solar radiation penetration into ice, simulating the near-surface transition layer. Numerical experiments were conducted to simulate the seasonal evolution of the snow-ice cover thickness and the vertical temperature profile within it.

2. Statement of the Problem

Melting snow is considered as a three-phase medium consisting of water ($i=1$), air ($i=2$), and ice ($i=3$). The process is described using the mass conservation equation for each phase, the generalized Darcy law, and the energy conservation equation for melting snow [5,16,17,18]:

$${\displaystyle \frac{\partial {\rho }_i}{\partial t}}+div\left({\rho }_i{\overrightarrow{u}}_i\right)=\sum _{j=1}^3{I}_{ji},i=1,2,3,I_{ji}=-I_{ij},\sum _{i,j=1}^3{I}_{ij}=0;$$

(1)

$$s_i\varphi ({\overrightarrow{u}}_i-{\overrightarrow{u}}_3)=-K_0{\displaystyle \frac{k_{0i}}{{\mu }_i}}(\nabla p_i+{\rho }_i^0\overrightarrow{g}),i=1,2,p_2-p_1=p_c(s_1,\theta ),\sum _{i=1}^2{s}_i=1;$$

(2)

$$\left(\sum _{i=1}^3{\rho }_i^0{c}_i{\alpha }_i\right){\displaystyle \frac{\partial \theta }{\partial t}}+\left(\sum _{i=1}^3{\rho }_i^0{c}_i{\overrightarrow{v}}_i\right)\nabla \theta =div({\lambda }_c\nabla \theta )-\nu I_{31}.$$

(3)

Here t is the time; ${\overrightarrow{u}}_i$ is the velocity of the i-th phase; ${\rho }_i$ - the reduced density related to the true density ${\rho }_i^0$ and the volume concentration ${\alpha }_i$ by the relation ${\rho }_i={\alpha }_i{\rho }_i^0$ (the condition ${\sum }_{i=1}^3{\alpha }_i=1$ is a consequence of the definition of ${\rho }_i$); $I_{ji}$ is the intensity of the transition of mass from the j-th to the i-th component in unit volume per unit time; ${\overrightarrow{v}}_i=\varphi s_i{\overrightarrow{u}}_i$ are the filtration rates of water and air; $\varphi$ is the porosity of snow; $s_1$, $s_2$ – water and air saturation (${\alpha }_1=\varphi s_1$, ${\alpha }_2=\varphi s_2$, ${\alpha }_3=1-\varphi$); $K_0$ – permeability tensor, $k_{0i}$ – relative phase permeabilities ($k_{0i}=k_{0i}\left(s_i\right)\ge 0,k_{0i}{|}_{s_i=0}=0$), ${\mu }_i$ – dynamic viscosity, $p_i$ – phase pressures, $p_c$ – capillary pressure , $\overrightarrow{g}$ is the gravity acceleration vector; $\theta$ – temperature of the environment (${\theta }_i=\theta$, $i=1,2,3$); $c_i=const>0$ – heat capacity of the i -th phase at constant volume; $\nu =const>0$ – specific heat of fusion of ice; ${\lambda }_c$ – thermal conductivity of snow (${\lambda }_c=a_c+b_c{\rho }_c^2$, ${\rho }_c={\sum }_{i=1}^3{\rho }_i^0{\alpha }_i$, $a_c=const>0$, $b_c=const>0$ [5]).

The system (1)–(3) is supplemented by the hypotheses ${\overrightarrow{u}}_3=0$, $I_{12}=0$, $I_{31}=I_{31}\left(\theta \right)$, ${\rho }_i^0={\rho }_i^0\left(\theta \right)$, $i=1,2,3$.

The purpose of the work is to construct an exact solution of the system of equations (1)–(3) under natural boundary conditions.

3. Self-Similar Solution to the Problem

We introduce the final values of temperature ${\theta }^{-}$, ${\theta }_1$ and ${\theta }^{+}$ (ice melting point). Let $0<{\theta }^{-}<{\theta }_1<{\theta }^{+}$. We assume that for all $\theta \in (0,\infty )$ the following relations hold: ${\alpha }_3\left(\theta \right)=0$ when $\theta \ge {\theta }^{+}$, ${\alpha }_3\left(\theta \right)=1-{\varphi }^{-}-{\varphi }_1(\theta -{\theta }_1)$ when ${\theta }_1\le \theta \le {\theta }^{+}$, ${\alpha }_3\left(\theta \right)=1-{\varphi }^{-}$ for $\theta \le {\theta }_1$. Here, ${\varphi }^{-}=\varphi \left({\theta }^{-}\right)\in (0,1)$, ${\varphi }_1=(1-{\varphi }^{-})/({\theta }^{+}-{\theta }_1)$ are the given parameters. In addition, it is assumed that the porous medium is homogeneous ($K_0=const>0$); ${\rho }_i^0=const>0$ (these conditions do not affect the generality of the results), in the coordinate system $xyz$ the vector $\overrightarrow{g}=(0,0,-g)$, the functions included in the system (1) - (3) depend on $z,t$. Eliminating in (1) $I_{31}$ , we obtain the system of equations

$${\displaystyle \frac{\partial }{\partial t}}(\varphi s_1{\rho }_1^0+{\rho }_3^0(1-\varphi ))+{\displaystyle \frac{\partial }{\partial z}}\left({\rho }_1^0{v}_1\right)=0;$$

(4)

$${\displaystyle \frac{\partial }{\partial t}}\left(\varphi s_2{\rho }_2^0\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\rho }_2^0{v}_2\right)=0;$$

(5)

$$v_i=-K_0{\displaystyle \frac{k_{0i}}{{\mu }_i}}({\displaystyle \frac{\partial p_i}{\partial z}}-{\rho }_i^{0}g),i=1,2,p_2-p_1=p_c(s_1,\theta ),s_1+s_2=1;$$

(6)

$$\left(\sum _{i=1}^3{\rho }_i^0{c}_i{\alpha }_i\right){\displaystyle \frac{\partial \theta }{\partial t}}+\left(\sum _{i=1}^2{\rho }_i^0{c}_i{v}_i\right){\displaystyle \frac{\partial \theta }{\partial z}}={\displaystyle \frac{\partial }{\partial z}}\left({\lambda }_c{\displaystyle \frac{\partial \theta }{\partial z}}\right)-\nu {\rho }_3^0{\displaystyle \frac{\partial \varphi }{\partial t}}.$$

(7)

For the system (4) - (7) consider the following problem: snow occupies the area of $(-\infty ,ct)$, $t>0$. At $z=-\infty$, there is no water ($s_1=0$, $v_1=0$), the air is stationary ($v_2=0$) and the temperature is set to $\theta ={\theta }^{-}$ (below the melting point ice); at $z=ct$ the velocities of water and air are known ($v_1=v_1^{+}$, $v_2=v_2^{+}$), the air pressure is ($p_2=p^{+}$) and the boundary condition is set ${\lambda }_c{\displaystyle \frac{d\theta }{dz}}=Q-\alpha (\theta -{\theta }^0)$ ($\alpha$ is the heat transfer coefficient, ${\theta }^0$ is the atmospheric air temperature, Q is the set heat flow). Assuming that all the sought functions depend only on the variable $\xi =z-ct$ (c is an unknown constant) from (4) - (7) we obtain

$$-c{\displaystyle \frac{d}{d\xi }}(\varphi s_1{\rho }_1^0+{\rho }_3^0(1-\varphi ))+{\displaystyle \frac{d}{d\xi }}\left({\rho }_1^0{v}_1\right)=0;$$

(8)

$$-c{\displaystyle \frac{d}{d\xi }}\left(\varphi s_2{\rho }_2^0\right)+{\displaystyle \frac{d}{d\xi }}\left({\rho }_2^0{v}_2\right)=0;$$

(9)

$$v_i=-K_0{\displaystyle \frac{k_{0i}}{{\mu }_i}}({\displaystyle \frac{d{p}_i}{d\xi }}-{\rho }_i^{0}g),p_2-p_1=p_c(s_1,\theta ),s_1+s_2=1;$$

(10)

$$\left(\sum _{i=1}^3{\rho }_i^0{c}_i(v_i-c{\alpha }_i)\right){\displaystyle \frac{d\theta }{d\xi }}-c\nu {\rho }_3^0{\displaystyle \frac{d\varphi }{d\xi }}={\displaystyle \frac{d}{d\xi }}\left({\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}\right);$$

(11)

$$s_1{|}_{\xi \to -\infty }=0,\theta {|}_{\xi \to -\infty }={\theta }^{-},{\displaystyle \frac{d\theta }{d\xi }}{|}_{\xi \to -\infty }=0,v_i{|}_{\xi \to -\infty }=0;$$

(12)

$$p_2\left(0\right)=p^{+},({\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}+\alpha (\theta -{\theta }^0)-Q){|}_{\xi =0}=0,v_i\left(0\right)=v_i^{+},i=1,2.$$

(13)

The sought functions are $s_1\left(\xi \right)$, $v_i\left(\xi \right)$, $p_i\left(\xi \right)$, and the constant c. The solution to problem (8) - (13) is constructed as follows. We find the constant c and obtain representations for the filtration rates and temperature by integrating equations (8), (9) and (10). Using these representations and (10), we arrive at an equation for the saturation of $s_1\left(\xi \right)$. Investigation of the solvability of the problem for $s_1\left(\xi \right)$ completes the construction of the solution to the problem (8)-(13). In the following, we use the notation $s_1\left(\xi \right)\equiv s\left(\xi \right)$.

Definition 1. A weak solution to the problem (8)–(13) is the function $\theta \left(\xi \right)$, $s\left(\xi \right)$, $v_i\left(\xi \right)$, $p_i\left(\xi \right)$ and a fixed parameter c if:

• 1) the function $\theta \left(\xi \right)$ has a continuous derivative and satisfies the equation

  $${\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}=c{\rho }_3^0(c_1-c_3)M\left(\theta \right)+(A_1{c}_1+A_2{c}_2)(\theta -{\theta }^{-})-\nu {\rho }_3^{0}c(\varphi -{\varphi }^{-})\equiv f_1\left(\theta \right),$$

  $${\displaystyle M\left(\theta \right)\equiv \underset{{\theta }^{-}}{\overset{\theta }{\int }}(1-\varphi \left(\zeta \right))d\zeta =\left\{\begin{array}{c}{M}_1=(1-{\varphi }^{-})(\frac{1}{2}{\theta }^{+}-{\theta }^{-}+\frac{1}{2}{\theta }_1),\theta \ge {\theta }^{+},\hfill \\ (1-{\varphi }^{-})(\theta -{\theta }^{-})-\frac{{\varphi }_1}{2}{(\theta -{\theta }_1)}^2,{\theta }_1\le \theta \le {\theta }^{+},\hfill \\ (1-{\varphi }^{-})(\theta -{\theta }^{-}),\theta \le {\theta }_1.\hfill \end{array}\right.}$$

  and the conditions $({\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}+\alpha (\theta -{\theta }^0)-Q){|}_{\xi =0}=0$, $\theta {|}_{\xi \to -\infty }={\theta }^{-}$, ${\displaystyle \frac{d\theta }{d\xi }}{|}_{\xi \to -\infty }=0$;
• 2) the function $s\left(\xi \right)$ has a continuous derivative with weight $a_0\left(s\right)$ ($a_0\left(0\right)=a_0\left(1\right)=0$), which satisfies the equation

  $$a_0\left(s\right){\displaystyle \frac{ds}{d\xi }}={\phi }_1{\phi }_2\gamma {\displaystyle \frac{p_0^{\prime }}{p_0}}{\displaystyle \frac{f_1}{{\lambda }_c}}+{\displaystyle \frac{1}{p_0}}\overline{g}{\phi }_1{\phi }_2+{\displaystyle \frac{1}{p_0}}\left|c\right|\varphi As-{\displaystyle \frac{1}{p_0}}\left|c\right|(\varphi -{\varphi }^{-})B\equiv f_2(s,\theta ),$$

  and the conditions $s\left(0\right)=s^{+}$, $s{|}_{\xi \to -\infty }=0$;
• 3) the functions $v_i\left(\xi \right)$ satisfy the equalities

  $$v_1=c\varphi s+c{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}({\varphi }^{-}-\varphi ),v_2=c\varphi (1-s)-c{\varphi }^{-},$$

  and the conditions $v_i\left(0\right)=v_i^{+}$, $v_i{|}_{\xi \to -\infty }=0$;
• 4) the functions $p_i\left(\xi \right)$ satisfy the equalities

  $$p_2\left(\xi \right)=p\left(\xi \right)+p_0\left(\theta \right)b\left(s\right),p_1\left(\xi \right)=p_2\left(\xi \right)-p_c(s\left(\xi \right),\theta \left(\xi \right)),$$

  $${\displaystyle p\left(\xi \right)=p^{+}-p_0\left({\theta }^{+}\right)b\left(s^{+}\right)-\underset{\xi }{\overset{0}{\int }}{f}_3(s\left(\zeta \right),\theta \left(\zeta \right))d\zeta ,}$$

  and the condition $p_2\left(0\right)=p^{+}$.

Determination of filtration rates

From (8) and (9) it follows that

$${\rho }_1^0{v}_1-c(\varphi s_1{\rho }_1^0+(1-\varphi ){\rho }_3^0)=A_1=const,$$

(14)

$${\rho }_2^0{v}_2-c\varphi s_2{\rho }_2^0=A_2=const.$$

(15)

From (14) and (12) we have $A_1=-c{\rho }_3^0(1-{\varphi }^{-})$, ${\varphi }^{-}\equiv \varphi \left({\theta }^{-}\right)$. From (15) and (12) we find $A_2=-c{\rho }_2^0{\varphi }^{-}$.

From (14), (15) we obtain representations for the filtration velocities:

$$v_1=c\varphi s+c{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}\left(\varphi \right),v_2=c\varphi (1-s)-c{\varphi }^{-}.$$

(16)

Representation for Temperature

Using (16) and conditions (12), we obtain

$${\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}=c{\rho }_3^0(c_1-c_3)M\left(\theta \right)+(A_1{c}_1+A_2{c}_2)(\theta -{\theta }^{-})-\nu {\rho }_3^{0}c(\varphi -{\varphi }^{-})\equiv f_1\left(\theta \right),$$

(17)

where

$${\displaystyle M\left(\theta \right)\equiv \underset{{\theta }^{-}}{\overset{\theta }{\int }}(1-\varphi \left(\zeta \right))d\zeta =\left\{\begin{array}{c}{M}_1=(1-{\varphi }^{-})(\frac{1}{2}{\theta }^{+}-{\theta }^{-}+\frac{1}{2}{\theta }_1),\theta \ge {\theta }^{+},\hfill \\ (1-{\varphi }^{-})(\theta -{\theta }^{-})-\frac{{\varphi }_1}{2}{(\theta -{\theta }_1)}^2,{\theta }_1\le \theta \le {\theta }^{+},\hfill \\ (1-{\varphi }^{-})(\theta -{\theta }^{-}),\theta \le {\theta }_1.\hfill \end{array}\right.}$$

Considering (14), (15) under the conditions (13), for unknown parameters c, $s^{+}=s_1\left(0\right)$, ($s_1\equiv s$, $s_2\equiv 1-s$), $\theta \left(0\right)$ we obtain the following nonlinear system of equations

$$v_1^{+}=c(\varphi \left(\theta \left(0\right)\right)s^{+}+{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}({\varphi }^{-}-\varphi \left(\theta \left(0\right)\right))),$$

(18)

$$v_2^{+}=c(\varphi \left(\theta \left(0\right)\right)(1-s^{+})-{\varphi }^{-}),$$

(19)

$$({\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}+\alpha (\theta -{\theta }^0)-Q){|}_{\xi =0}=0.$$

(20)

From the equation for temperature: ${\lambda }_c{\displaystyle \frac{d\theta }{d\xi }}\equiv f_1\left(\theta \right),$ then the system can be rewritten as:

$$\left\{\begin{array}{c}{v}_1^{+}=c[\varphi \left(\theta \left(0\right)\right)s\left(0\right)+{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}({\varphi }^{-}-\varphi \left(\theta \left(0\right)\right))],\hfill \\ v_2^{+}=c[\varphi \left(\theta \left(0\right)\right)(1-s\left(0\right))-{\varphi }^{-}],\hfill \\ (f_1\left(\theta \right)+\alpha (\theta -{\theta }^0)-Q){|}_{\xi =0}=0,\hfill \end{array}\right.$$

here $f_1\left(\theta \right)=c\{{\rho }_3^0(c_1-c_3)M\left(\theta \right)-(\theta -{\theta }^{-})({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0(\varphi -{\varphi }^{-})\}.$

The solution to the problem (12), (13), (17) can be represented as

$${\displaystyle I\left(\theta \right)\equiv \underset{\theta }{\overset{{\theta }^{+}}{\int }}\frac{d\zeta }{f_1\left(\zeta \right)}=\underset{\xi }{\overset{0}{\int }}\frac{d\zeta }{{\lambda }_c\left(\zeta \right)}\equiv \psi \left({\lambda }_c\left(\xi \right)\right),\theta \left(\xi \right)=I^{-1}\left(\psi \left({\lambda }_c\left(\xi \right)\right)\right).}$$

(21)

For $\theta \in [{\theta }_1,{\theta }^{+}]$ from (21) we obtain

$${\displaystyle I\left(\theta \right)=\frac{1}{b_1}\underset{\theta }{\overset{{\theta }^{+}}{\int }}\frac{d\zeta }{{(\zeta -{\theta }_1+\alpha )}^2-{\beta }^2}=\psi \left({\lambda }_c\left(\xi \right)\right),\theta \in [{\theta }_1,{\theta }^{+}].}$$

(22)

Here $d_1=c{\rho }_3^0(c_1-c_3)(1-{\varphi }^{-}+A_1{c}_1+A_2{c}_2)-\nu c{\rho }_3^0{\varphi }_1$, $b_1=-c{\varphi }_1{\rho }_3^0(c_1-c_3)/2$, $\alpha =d_1/\left(2b_1\right)$, ${\beta }^2=(d_1^2-4a_1{b}_1)/\left(4b_1^2\right)$.

For $\theta \in [{\theta }^{-},{\theta }_1]$ from (21) we obtain

$${\displaystyle \theta \left(\xi \right)={\theta }^{-}+({\theta }_1-{\theta }^{-})exp(-b_2\underset{\xi }{\overset{{\xi }_1}{\int }}\frac{d\zeta }{{\lambda }_c\left(\zeta \right)}).}$$

(23)

Here $b_2=c{\rho }_3^0(c_1-c_3)(1-{\varphi }^{-}+A_1{c}_1+A_2{c}_2).$The solvability of this problem is proved in detail in [6]. Thus, for a given function ${\lambda }_c(s,\theta )$, the representation (21) and its special cases (22), (23) determine the temperature for all $\xi \in (-\infty ,0)$.

Algorithm for solving the system (18-20)

The system solution for unknown parameters $\theta \left(0\right)$, $c<0$, $s^{+}$ is divided into three stages:

1. Let $\varphi \left(\theta \left(0\right)\right)(1-s\left(0\right))-{\varphi }^{-}=0$, then $v_2^{+}=0$, $v_1^{+}<0$ and $s\left(0\right)\ne 1$.

$$s\left(0\right)=1-{\displaystyle \frac{{\varphi }^{-}}{\varphi \left(\theta \right(0\left)\right)}},c={\displaystyle \frac{v_1^{+}}{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(1-{\rho }_3^0/{\rho }_1^0)}}<0,0\le s\left(0\right)<1.$$

To find $\theta \left(0\right)$, we consider the options and use the flow condition on the boundary.

If $\theta \ge {\theta }^{+}$, then $M\left({\theta }^{+}\right)=M_1=const$, $\varphi \left({\theta }^{+}\right)=1$ - the solution is nonphysical.

If $\theta \le {\theta }_1$ , then $M\left({\theta }^{-}\right)=0$ and $c=0$ - the solution is nonphysical.

For the case ${\theta }_1\le \theta \left(0\right)\le {\theta }^{+}$, we substitute c into $f_1$ and get:

$F_1\equiv v_1^{+}\left({\rho }_3^0(c_1-c_3)\right((1-{\varphi }^{-})(\theta \left(0\right)-{\theta }^{-})-{\displaystyle \frac{{\varphi }_1}{2}}{(\theta \left(0\right)-{\theta }_1)}^2)-(\theta \left(0\right)-{\theta }^{-})({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0{\varphi }_1(\theta \left(0\right)-{\theta }_1))-{\varphi }_1(\theta \left(0\right)+{\theta }_1)(1-{\rho }_3^0/{\rho }_1^0)(\alpha (\theta \left(0\right)-{\theta }^0)-Q)=0.$

This is the quadratic equation for $\theta \left(0\right)$:

$A_1{\theta }^2+B_1\theta +C_1=0,$ where

$A_1=\alpha {\varphi }_1(1-{\rho }_3^0/{\rho }_1^0)-v_1^{+}{\rho }_3^0{\displaystyle \frac{{\varphi }_1}{2}}(c_1-c_3)$,

$B_1=v_1^{+}\left({\rho }_3^0(c_1-c_3)\right(1-{\varphi }^{-}+{\varphi }_1{\theta }_1)-({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0{\varphi }_1)-{\varphi }_1(1-{\rho }_3^0/{\rho }_1^0)(\alpha ({\theta }^0+{\theta }_1)+Q)$,

$C_1=v_1^{+}\left({\rho }_3^0(c_1-c_3)\right({\displaystyle \frac{{\varphi }_1}{2}}{\theta }_1^2-{\theta }^{-}(1-{\varphi }^{-}))+{\theta }^{-}({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})+\nu {\rho }_3^0{\varphi }_1{\theta }_1)+{\varphi }_1{\theta }_1(1-{\rho }_3^0/{\rho }_1^0)(\alpha {\theta }^0+Q)$. Here ${\rho }_1^0,{\rho }_2^0,{\rho }_3^0,c_1,c_2,c_3,{\theta }_1,{\theta }^{-},\nu ,{\varphi }_1,{\varphi }^{-}$ are given parameters, and the parameters $v_1^{+},\alpha ,{\theta }^0,Q$ must be subject to conditions so that the following relations are satisfied:

• If ${\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1<0$ , then there are no real roots.
• If ${\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1=0$ and $B_1/A_1>0$, then there are no solutions.
• If ${\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1=0$ and $B_1/A_1<0$, then the solution is:
  $\theta \left(0\right)=-B_1/A_1,$ here $\theta \left(0\right)$ must satisfy the condition ${\theta }_1\le \theta \left(0\right)\le {\theta }^{+}$.
• If ${\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1>0$, then there are two roots:
  ${\theta }_1\left(0\right)=-B_1/A_1+{({\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1)}^{1/2}$, which

  –

  for $B_1/A_1<0$ must satisfy the condition ${\theta }_1\le {\theta }_1\left(0\right)\le {\theta }^{+}$;

  –

  for $B_1/A_1>0$ must satisfy the inequalities

  ${({\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1)}^{1/2}>-B_1/A_1$ and ${\theta }_1\le {\theta }_1\left(0\right)\le {\theta }^{+}$.

  ${\theta }_2\left(0\right)=-B_1/A_1-{({\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1)}^{1/2}$, which

  –

  for $B_1/A_1>0$ is not a solution;

  –

  for $B_1/A_1<0$ must satisfy the inequalities

  ${({\left({\displaystyle \frac{B_1}{2A_1}}\right)}^2-C_1)}^{1/2}<-B_1/A_1$ and ${\theta }_1\le {\theta }_2\left(0\right)\le {\theta }^{+}$.

2. Let $\varphi \left(\theta \left(0\right)\right)s\left(0\right)+{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}({\varphi }^{-}-\varphi \left(\theta \left(0\right)\right))=0$, then $v_1^{+}=0$, $v_2^{+}<0$.

$$s\left(0\right)={\displaystyle \frac{{\rho }_3^0/{\rho }_1^0(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})}{\varphi \left(\theta \right(0\left)\right)}},c={\displaystyle \frac{v_2^{+}}{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(1-{\rho }_3^0/{\rho }_1^0)}}<0,0\le s\left(0\right)\le 1.$$

To find $\theta \left(0\right)$, we consider the options and use the flow condition on the boundary.

If $\theta \ge {\theta }^{+}$, then $M\left({\theta }^{+}\right)=M_1=const$, $\varphi \left({\theta }^{+}\right)=1$ - the solution is nonphysical.

If $\theta \le {\theta }_1$ , then $M\left({\theta }^{-}\right)=0$ and $c=0$ - the solution is nonphysical.

For the case ${\theta }_1\le \theta \left(0\right)\le {\theta }^{+}$, we substitute c into $f_1$ and get:

$F_2\equiv v_2^{+}\left({\rho }_3^0(c_1-c_3)\right((1-{\varphi }^{-})(\theta \left(0\right)-{\theta }^{-})-{\displaystyle \frac{{\varphi }_1}{2}}{(\theta \left(0\right)-{\theta }_1)}^2)-(\theta \left(0\right)-{\theta }^{-})({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0{\varphi }_1(\theta \left(0\right)-{\theta }_1))-{\varphi }_1(\theta \left(0\right)+{\theta }_1)(1-{\rho }_3^0/{\rho }_1^0)(\alpha (\theta \left(0\right)-{\theta }^0)-Q)=0.$

This is the quadratic equation for $\theta \left(0\right)$:

$A_2{\theta }^2+B_2\theta +C_2=0,$ where

$A_2=\alpha {\varphi }_1(1-{\rho }_3^0/{\rho }_1^0)-v_2^{+}{\rho }_3^0{\displaystyle \frac{{\varphi }_1}{2}}(c_1-c_3)$,

$B_2=v_2^{+}\left({\rho }_3^0(c_1-c_3)\right(1-{\varphi }^{-}+{\varphi }_1{\theta }_1)-({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0{\varphi }_1)-{\varphi }_1(1-{\rho }_3^0/{\rho }_1^0)(\alpha ({\theta }^0+{\theta }_1)+Q)$,

$C_2=v_2^{+}\left({\rho }_3^0(c_1-c_3)\right({\displaystyle \frac{{\varphi }_1}{2}}{\theta }_1^2-{\theta }^{-}(1-{\varphi }^{-}))+{\theta }^{-}({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})+\nu {\rho }_3^0{\varphi }_1{\theta }_1)+{\varphi }_1{\theta }_1(1-{\rho }_3^0/{\rho }_1^0)(\alpha {\theta }^0+Q)$.

Here ${\rho }_1^0,{\rho }_2^0,{\rho }_3^0,c_1,c_2,c_3,{\theta }_1,{\theta }^{-},\nu ,{\varphi }_1,{\varphi }^{-}$ are given parameters, and the parameters $v_2^{+},\alpha ,{\theta }^0,Q$ must be subject to conditions so that the following relations are satisfied:

• If ${\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2<0$ , then there are no real roots.
• If ${\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2=0$ and $B_2/A_2>0$, then there are no solutions.
• If ${\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2=0$ and $B_2/A_2<0$, then the solution is:
  $\theta \left(0\right)=-B_2/A_2,$ here $\theta \left(0\right)$ must satisfy the condition ${\theta }_1\le \theta \left(0\right)\le {\theta }^{+}$.
• If ${\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2>0$, then there are two roots:
  ${\theta }_1\left(0\right)=-B_2/A_2+{({\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2)}^{1/2}$, which

  –

  for $B_2/A_2<0$ must satisfy the condition ${\theta }_1\le {\theta }_1\left(0\right)\le {\theta }^{+}$;

  –

  for $B_2/A_2>0$ must satisfy the inequalities

  ${({\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2)}^{1/2}>-B_2/A_2$ and ${\theta }_1\le {\theta }_1\left(0\right)\le {\theta }^{+}$.

  ${\theta }_2\left(0\right)=-B_2/A_2-{({\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2)}^{1/2}$, which

  –

  for $B_2/A_2>0$ is not a solution;

  –

  for $B_2/A_2<0$ must satisfy the inequalities

  ${({\left({\displaystyle \frac{B_2}{2A_2}}\right)}^2-C_2)}^{1/2}<-B_2/A_2$ and ${\theta }_1\le {\theta }_2\left(0\right)\le {\theta }^{+}$.

3. Let $\varphi \left(\theta \left(0\right)\right)(1-s\left(0\right))-{\varphi }^{-}\ne 0$ and $\varphi \left(\theta \left(0\right)\right)s\left(0\right)+{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}({\varphi }^{-}-\varphi \left(\theta \left(0\right)\right))\ne 0$, then $v_1^{+}<0$, $v_2^{+}<0$.

$$s\left(0\right)={\displaystyle \frac{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(v_1^{+}+\frac{{\rho }_3^0}{{\rho }_1^0}{v}_2^{+})}{\varphi \left(\theta \left(0\right)\right)(v_1^{+}+v_2^{+})}},c={\displaystyle \frac{v_1^{+}+v_2^{+}}{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(1-{\rho }_3^0/{\rho }_1^0)}}<0,0\le s\left(0\right)\le 1.$$

To find $\theta \left(0\right)$, we consider the options and use the flow condition on the boundary.

If $\theta \ge {\theta }^{+}$, then $M\left({\theta }^{+}\right)=M_1=const$, $\varphi \left({\theta }^{+}\right)=1$ - the solution is nonphysical.

If $\theta \le {\theta }_1$ , then $M\left({\theta }^{-}\right)=0$ and $c=0$ - the solution is nonphysical.

For the case ${\theta }_1\le \theta \left(0\right)\le {\theta }^{+}$, we substitute c into $f_1$ and get:

$F_3\equiv (v_1^{+}+v_2^{+})\left({\rho }_3^0(c_1-c_3)\right((1-{\varphi }^{-})(\theta \left(0\right)-{\theta }^{-})-{\displaystyle \frac{{\varphi }_1}{2}}{(\theta \left(0\right)-{\theta }_1)}^2)-(\theta \left(0\right)-{\theta }^{-})({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0{\varphi }_1(\theta \left(0\right)-{\theta }_1))-{\varphi }_1(\theta \left(0\right)+{\theta }_1)(1-{\rho }_3^0/{\rho }_1^0)(\alpha (\theta \left(0\right)-{\theta }^0)-Q)=0.$

This is the quadratic equation for $\theta \left(0\right)$:

$A_3{\theta }^2+B_3\theta +C_3=0,$ where

$A_3=\alpha {\varphi }_1(1-{\rho }_3^0/{\rho }_1^0)-{\rho }_3^0{\displaystyle \frac{{\varphi }_1}{2}}(c_1-c_3)(v_1^{+}+v_2^{+})$,

$B_3=(v_1^{+}+v_2^{+})\left({\rho }_3^0(c_1-c_3)\right(1-{\varphi }^{-}+{\varphi }_1{\theta }_1)-({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})-\nu {\rho }_3^0{\varphi }_1)-{\varphi }_1(1-{\rho }_3^0/{\rho }_1^0)(\alpha ({\theta }^0+{\theta }_1)+Q)$,

$C_3=(v_1^{+}+v_2^{+})\left({\rho }_3^0(c_1-c_3)\right({\displaystyle \frac{{\varphi }_1}{2}}{\theta }_1^2-{\theta }^{-}(1-{\varphi }^{-}))+{\theta }^{-}({\rho }_3^0{c}_1(1-{\varphi }^{-})+{\rho }_2^0{c}_2{\varphi }^{-})+\nu {\rho }_3^0{\varphi }_1{\theta }_1)+{\varphi }_1{\theta }_1(1-{\rho }_3^0/{\rho }_1^0)(\alpha {\theta }^0+Q)$.

Here ${\rho }_1^0,{\rho }_2^0,{\rho }_3^0,c_1,c_2,c_3,{\theta }_1,{\theta }^{-},\nu ,{\varphi }_1,{\varphi }^{-}$ are given parameters, and the parameters $v_1^{+},v_2^{+},\alpha ,{\theta }^0,Q$ must be subject to conditions so that the following relations are satisfied:

• If ${\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3<0$ , then there are no real roots.
• If ${\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3=0$ and $B_3/A_3>0$, then there are no solutions.
• If ${\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3=0$ and $B_3/A_3<0$, then there is a solution:
  $\theta \left(0\right)=-B_3/A_3,$ where $\theta \left(0\right)$ must satisfy the condition ${\theta }_1\le \theta \left(0\right)\le {\theta }^{+}$.
• If ${\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3>0$, then there are two roots:
  ${\theta }_1\left(0\right)=-B_3/A_3+{({\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3)}^{1/2}$, which

  –

  for $B_3/A_3<0$ must satisfy the condition ${\theta }_1\le {\theta }_1\left(0\right)\le {\theta }^{+}$;

  –

  for $B_3/A_3>0$ must satisfy the inequalities

  ${({\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3)}^{1/2}>-B_3/A_3$ and ${\theta }_1\le {\theta }_1\left(0\right)\le {\theta }^{+}$.

  ${\theta }_2\left(0\right)=-B_3/A_3-{({\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3)}^{1/2}$, which

  –

  for $B_3/A_3>0$ is not a solution;

  –

  for $B_3/A_3<0$ must satisfy the inequalities

  ${({\left({\displaystyle \frac{B_3}{2A_3}}\right)}^2-C_3)}^{1/2}<-B_3/A_3$ and ${\theta }_1\le {\theta }_2\left(0\right)\le {\theta }^{+}$.

Determination of saturation and pressure

From (10) and (16) we have

$$v_1=c\varphi s-c{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}(\varphi -{\varphi }^{-})=-K_0{\displaystyle \frac{k_{01}}{{\mu }_1}}({\displaystyle \frac{d{p}_1}{d\xi }}-{\rho }_1^{0}g),$$

$$v_2=c\varphi (1-s)-c{\varphi }^{-}=-K_0{\displaystyle \frac{k_{02}}{{\mu }_2}}({\displaystyle \frac{d{p}_2}{d\xi }}-{\rho }_2^{0}g).$$

Eliminating $p_1$ and $p_2$ from these relations using the second equation in (10), we obtain

$$-k_{01}{k}_{02}{K}_0{\displaystyle \frac{d{p}_c}{d\xi }}=-{\mu }_1{v}_1{k}_{02}+{\mu }_2{v}_2{k}_{01}+K_0{k}_{01}{k}_{02}g({\rho }_1^0-{\rho }_2^0).$$

Following [6], we set: $p_c(s,\theta )=p_0\left(\theta \right)\gamma \left(s\right)$, ${\displaystyle \frac{d\gamma }{ds}}\equiv {\gamma }^{\prime }<0$. For $s\in (-\infty ,+\infty )$, the relative phase permeabilities $k_{0i}$ are defined as follows: $k_{0i}=0$ for $s_i\le 1$, $k_{0i}={\overline{k}}_{0i}(s_i,\theta )s_i^{n_i}$ for $0\le s_i\le 1$, $k_{0i}={\overline{k}}_{0i}(1,\theta )$ for $s_i\ge 1$. Here we have the constants $n_i>1$.

The function $a_0\left(s\right)=-{\gamma }^{\prime }{k}_{01}{k}_{02}>0$ for $s\in (0,1)$ and $a_0\left(s\right)=0$ for $s\le 0$ and $s\ge 1$ is assumed. Using the equation for temperature (17), we represent the equation for saturation in the form

$$\begin{array}{c}{a}_0\left(s\right){\displaystyle \frac{ds}{d\xi }}={\phi }_1{\phi }_2\gamma {\displaystyle \frac{p_0^{\prime }}{p_0}}{\displaystyle \frac{f_1}{{\lambda }_c}}+{\displaystyle \frac{1}{p_0}}\overline{g}{\phi }_1{\phi }_2+{\displaystyle \frac{1}{p_0}}\left|c\right|\varphi As-\\ -{\displaystyle \frac{1}{p_0}}\left|c\right|(\varphi -{\varphi }^{-})B\equiv f_2(s,\theta ),\end{array}$$

(24)

where ${\phi }_i=0$ for $s_i\le 0$, ${\phi }_i=s_i^{n_i}$ for $0\le s_i\le 1$, ${\phi }_i=1$ for $s_i\ge 1$,

$$a_0=-{\phi }_1{\phi }_2{\gamma }^{\prime },A={\overline{\mu }}_1{\phi }_2+{\overline{\mu }}_2{\phi }_1,B={\overline{\mu }}_1{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}}{\phi }_2+{\overline{\mu }}_2{\phi }_1,\overline{g}=g({\rho }_1^0-{\rho }_2^0),$$

$$p_0^{\prime }={\displaystyle \frac{d{p}_0}{d\theta }},{\overline{\mu }}_i={\displaystyle \frac{{\mu }_i}{K_0{\overline{k}}_{0i}}},i=1,2.$$

Equation (24) is considered for $\xi <0$ and the condition $s\left(0\right)=s^{+}$, that is, for $s\left(\xi \right)$ the Cauchy problem is considered, and the condition $s{|}_{\xi =-\infty }=0$ must be justified.

To find the pressures $p_1$ and $p_2$, we consider the equality following from Darcy’s laws (the first equations in (10))

$$\sum _{i=1}^2({\overline{\mu }}_i{v}_i-g{\rho }_i^0{\phi }_i)=-({\phi }_1+{\phi }_2)({\displaystyle \frac{d{p}_2}{d\xi }}-{\displaystyle \frac{{\phi }_1}{{\phi }_1+{\phi }_2}}{\displaystyle \frac{d{p}_c}{d\xi }}).$$

Let’s put [6]

$${\displaystyle p\left(\xi \right)\equiv p_2\left(\xi \right)-p_0\left(\theta \right)b\left(s\right),b\left(s\right)=\underset{0}{\overset{s}{\int }}\frac{{\phi }_1\left(\zeta \right){\gamma }^{\prime }\left(\zeta \right)}{{\phi }_1\left(\zeta \right)+{\phi }_2\left(\zeta \right)}d\zeta .}$$

Then

$$\begin{array}{c}{\displaystyle \frac{dp}{d\xi }}=-{\displaystyle \frac{1}{({\phi }_1\left(s\right)+{\phi }_2\left(s\right))}}[\sum _{i=1}^2({\overline{\mu }}_i{v}_i\left(\xi \right)-g{\rho }_i^0{\phi }_i)-p_0^{\prime }\left(\theta \right){\displaystyle \frac{f_1}{{\lambda }_c}}({\phi }_1\left(s\right)\gamma \left(s\right)+\\ +b\left(s\right)({\phi }_1\left(s\right)+{\phi }_2\left(s\right))\left)\right]\equiv f_3(s,\theta ),\end{array}$$

(25)

$$p\left(0\right)=p^{+}-p_0\left({\theta }^{+}\right)b\left(s^{+}\right).$$

From (25) it follows

$${\displaystyle p\left(\xi \right)=p^{+}-p_0\left({\theta }^{+}\right)b\left(s^{+}\right)-\underset{\xi }{\overset{0}{\int }}{f}_3(s\left(\zeta \right),\theta \left(\zeta \right))d\zeta ,\left[5mm\right]}$$

$$p_2\left(\xi \right)=p\left(\xi \right)+p_0\left(\theta \right)b\left(s\right),p_1\left(\xi \right)=p_2\left(\xi \right)-p_c(s\left(\xi \right),\theta \left(\xi \right)).$$

(26)

If the functions $s\left(\xi \right)$ and $\theta \left(\xi \right)$ are found, then the filtration rates $v_i\left(\xi \right)$ and pressures $p_i\left(\xi \right)$ are determined by formulas (16), (26).

Let us assume that the functions ${\overline{\mu }}_i(s,\theta )$, ${\overline{k}}_{0i}(s,\theta )$, $p_0\left(\theta \right)$, $\gamma \left(s\right)$, $a_0\left(s\right)$ satisfies the following conditions (${\overline{\Omega }}^{*}=[0,1]\times [{\theta }^{-},{\theta }^{+}]$):

$$a)0<{\nu }_0^{-1}\le ({\mu }_i(s,\theta ),{\overline{k}}_{0i}(s,\theta ),p_0\left(\theta \right),|{\displaystyle \frac{d\gamma \left(s\right)}{ds}}\left|\right)\le {\nu }_0,{\displaystyle \frac{a_0\left(s\right)}{s}}{|}_{s=0}=0,$$

$$b\left)\right(\left|\right|{\displaystyle \frac{d\gamma }{ds}}{\left|\right|}_{C[0,1]},\left|\right|{\displaystyle \frac{d{p}_0}{d\theta }}{\left|\right|}_{C[{\theta }^{-},{\theta }^{+}]},\left|\right|{\mu }_i(s,\theta ),{\overline{k}}_{0i}{(s,\theta )\left|\right|}_{C\left({\overline{\Omega }}^{*}\right)})\le {\nu }_0.$$

Theorem 1. 

If conditions “a”, “b” and condition $s^{+}\in (0,1]$ are satisfied, there exists at least one weak solution of the problem (8)–(13). This solution (in addition to Definition 1) has the following properties:

$0\le s\left(\xi \right)\le 1$, ${\theta }^{-}\le \theta \left(\xi \right)\le {\theta }^{+}$ and $c={\displaystyle \frac{v_1^{+}+v_2^{+}}{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(1-{\rho }_3^0/{\rho }_1^0)}}<0$ (found from solving the system (18)–(20)).

The function $\theta \left(\xi \right)$ is monotone in $\xi$ and ${\theta }_1=\theta \left({\xi }_1\right)$, ${\xi }_1<0$. There exists a point ${\xi }_{*}\in (-\infty ,{\xi }_1]$ such that $s\left(\xi \right)=0$ for all $\xi \le {\xi }_{*}$.

Consider the problem

$$a_0\left(s\right){\displaystyle \frac{ds}{d\xi }}=f_2(s,\theta ),{\displaystyle \frac{d\theta }{d\xi }}={\displaystyle \frac{f_1\left(\theta \right)}{{\lambda }_c(s,\theta )}},\xi <0,s\left(0\right)=s^{+},\theta \left(0\right)={\theta }^{+}.$$

(27)

The solvability of this problem is proved in detail in [6].

Let $\epsilon \in (0,1)$, $a_{\epsilon }\left(s\right)\equiv a_0\left(s\right)+\epsilon >0$. For $\xi <0$, instead of (27), we consider the problem

$$a_{\epsilon }\left(s^{\epsilon }\right){\displaystyle \frac{d{s}^{\epsilon }}{d\xi }}=f_2(s^{\epsilon },{\theta }^{\epsilon }),{\displaystyle \frac{d{\theta }^{\epsilon }}{d\xi }}={\displaystyle \frac{f_1\left({\theta }^{\epsilon }\right)}{{\lambda }_c(s^{\epsilon },{\theta }^{\epsilon })}},s^{\epsilon }\left(0\right)=s^{+},{\theta }^{\epsilon }\left(0\right)={\theta }^{+}.$$

(28)

The local solvability of the problem (28) for each $\epsilon >0$ follows from the known results [19], p. 21.

Lemma 1. If $s^{\epsilon }\left(\xi \right)$ is a solution problem (28) and $s^{+}\in [0,1]$, then $0\le s^{\epsilon }\left(\xi \right)\le 1$.

Proof of Lemma 1. 

From (28) we get

$${\displaystyle \frac{d{s}^{\epsilon }}{d\xi }}-R{s}^{\epsilon }=-Q,$$

Where

$$R={\displaystyle \frac{1}{p_0{a}_{\epsilon }}}\left(\right|c|\varphi A+\overline{g}{\displaystyle \frac{{\phi }_1\left(s^{\epsilon }\right)}{s^{\epsilon }}}{\phi }_2+{\displaystyle \frac{{\phi }_1\left(s^{\epsilon }\right)}{s^{\epsilon }}}{\phi }_2\gamma {\displaystyle \frac{f_1}{{\lambda }_c}}{\left(p_0^{\prime }\right)}^{+}),$$

$$Q={\displaystyle \frac{1}{p_0{a}_{\epsilon }}}\left(\right|c|(\varphi -{\varphi }^{-})B+{\phi }_1{\phi }_2\gamma {\displaystyle \frac{f_1}{{\lambda }_c}}{\left(p_0^{\prime }\right)}^{-}),$$

$${\left(p_0^{\prime }\right)}^{+}=max(0,p_0^{\prime }),{\left(p_0^{\prime }\right)}^{-}=-min(0,p_0^{\prime }),p_0^{\prime }={\left(p_0^{\prime }\right)}^{+}-{\left(p_0^{\prime }\right)}^{-}.$$

Since $A>0$, $B>0$ and $R>0$, $Q\ge 0$. Therefore

$${\displaystyle s^{\epsilon }\left(\xi \right)=(s^{+}+\underset{\xi }{\overset{0}{\int }}Q\left(x\right)e^{{\displaystyle -\underset{0}{\overset{x}{\int }}R\left(\zeta \right)d\zeta }}dx)e^{{\displaystyle \underset{0}{\overset{\xi }{\int }}R\left(\zeta \right)d\zeta }}\ge 0.}$$

For the function $(1-s^{\epsilon }\left(\xi \right))$ from (28) we find

$${\displaystyle \frac{d(1-s^{\epsilon })}{d\xi }}-R_1(1-s^{\epsilon })=-Q_1,$$

where

$$R_1={\displaystyle \frac{1}{p_0{a}_{\epsilon }}}\left(\right|c|\varphi A+{\displaystyle \frac{{\phi }_2\left(s^{\epsilon }\right)}{1-s^{\epsilon }}}{\phi }_1\gamma {\displaystyle \frac{f_1}{{\lambda }_c}}{\left(p_0^{\prime }\right)}^{-})>0,$$

$$Q_1={\displaystyle \frac{1}{p_0{a}_{\epsilon }}}\left(\right|c|(\varphi {\overline{\mu }}_1{\phi }_2(1-{\displaystyle \frac{{\rho }_3^0}{{\rho }_1^0}})+{\varphi }^{-}B)+\overline{g}{\phi }_1{\phi }_2+{\phi }_1{\phi }_2\gamma {\displaystyle \frac{f_1}{{\lambda }_c}}{\left(p_0^{\prime }\right)}^{+})>0.$$

Then

$${\displaystyle 1-s^{\epsilon }\left(\xi \right)=(1-s^{+}+\underset{\xi }{\overset{0}{\int }}{Q}_1\left(x\right)e^{{\displaystyle -\underset{0}{\overset{x}{\int }}{R}_1\left(\zeta \right)d\zeta }}dx)e^{{\displaystyle \underset{0}{\overset{\xi }{\int }}{R}_1\left(\zeta \right)d\zeta }}\ge 0,}$$

which completes the proof. □

The function $s\left(\xi \right)$ is continuous on the interval $[{\xi }_1,0]$ and, therefore, there exists a value $s^1\equiv s\left({\xi }_1\right)\in [0,1]$. Therefore, we can consider the problem

$$a_0\left(s\right){\displaystyle \frac{ds}{d\xi }}=f_2(s,\theta ),{\displaystyle \frac{d\theta }{d\xi }}={\displaystyle \frac{f_1\left(\theta \right)}{{\lambda }_c(s,\theta )}},\xi <{\xi }_1,s\left({\xi }_1\right)=s^1,\theta \left({\xi }_1\right)={\theta }_1,$$

(29)

where

$$f_2={\phi }_1{\phi }_2\gamma {\displaystyle \frac{p_0^{\prime }}{p_0}}{\displaystyle \frac{f_1}{{\lambda }_c}}+{\displaystyle \frac{1}{p_0}}\overline{g}{\phi }_1{\phi }_2+{\displaystyle \frac{1}{p_0}}\left|c\right|{\varphi }^{-}As.$$

Lemma 2. Let $s\left(\xi \right)$ be a solution to problem (29) and let the condition ${\displaystyle \frac{{\phi }_1\left(s\right)}{s}}\gamma \left(s\right)\le {\nu }_0$, $s\in [0,1]$ be satisfied. Then there exists a point ${\xi }_{*}\le {\xi }_1$ such that $s\left(\xi \right)\equiv 0$ for all $\xi \le {\xi }_{*}$. If $s^1=0$ and additionally ${\displaystyle \frac{{\phi }_1\left(s\right)}{s}}\gamma \left(s\right){|}_{s=0}=0$, then ${\xi }_{*}={\xi }_1$.

Proof of Lemma 2. 

From (29) it follows

$${\displaystyle \frac{du}{d\xi }+D_1(s,\theta )=D_2(s,\theta ),u=\underset{0}{\overset{s}{\int }}\frac{a_0\left(\zeta \right)}{\zeta }d\zeta .}$$

(30)

Here

$$D_1(s,\theta )={\displaystyle \frac{{\phi }_1\left(s\right)}{s}}{\phi }_2\gamma {\left(p_0^{\prime }\right)}^{-}{\displaystyle \frac{f_1}{p_0{\lambda }_c}},f_1=b_2(\theta -{\theta }^{-})\ge 0,$$

$$D_2(s,\theta )={\displaystyle \frac{1}{p_0}}\left(\right|c|{\varphi }^{-}A+\overline{g}{\displaystyle \frac{{\phi }_1}{s}}{\phi }_2+{\displaystyle \frac{{\phi }_1\left(s\right)}{s}}{\phi }_2\gamma {\left(p_0^{\prime }\right)}^{+}{\displaystyle \frac{f_1}{{\lambda }_c}}).$$

We have

$$D_2(s,\theta )\ge \left|c\right|{\varphi }^{-}A\underset{s,\theta }{min}{\displaystyle \frac{{\overline{\mu }}_1}{p_0}}min(1,\alpha ,{\pi }_{\mu })\equiv D_2^0>0,$$

where $\alpha =\underset{s,\theta }{min}{\displaystyle \frac{{\overline{\mu }}_1}{{\overline{\mu }}_2}}>0$, $\theta \in [{\theta }^{-},{\theta }^{+}]$, $s\in [0,1]$. Using the representation (26) for $\theta \left(\xi \right)\in [{\theta }^{-},{\theta }_1]$, we obtain

$${\displaystyle |D_1|\le D_1^0{e}^{-\frac{1}{{\lambda }_c^{+}}({\xi }_1-\xi )},\underset{\xi }{\overset{{\xi }_1}{\int }}\left|D_1\left(\zeta \right)\right|d\zeta \le \frac{1}{{\lambda }_c^{+}}{D}_1^0,}$$

where

$$D_1^0={\displaystyle \frac{1}{{\lambda }_c^{-}}}{b}_2({\theta }_1-{\theta }^{-})\underset{s,\theta }{max}\left({\displaystyle \frac{1}{p_0}}{\displaystyle \frac{{\phi }_1\left(s\right)}{s}}{\phi }_2\gamma {\left(p_0^{\prime }\right)}^{-}\right)\le {\displaystyle \frac{1}{{\lambda }_c^{-}}}{b}_2({\theta }_1-{\theta }^{-}){\nu }_0^3,$$

$${\lambda }_c^{-}=a_c+b_c{\left({\rho }_2^0{\varphi }^{-}\right)}^2,{\lambda }_c^{+}=a_c+b_c{({\rho }_1^0{\varphi }^{-}+{\rho }_3^0(1-{\varphi }^{-}))}^2.$$

Integrating equation (30) over $\xi$ from an arbitrary value $\xi$ to ${\xi }_1$ , we obtain

$$u\left(s\left({\xi }_1\right)\right)+{\displaystyle \frac{1}{{\lambda }_c^{+}}}{D}_1^0\ge D_2^0({\xi }_1-\xi )+u\left(s\left(\xi \right)\right).$$

(31)

Here $u\left(s\left({\xi }_1\right)\right)=0$ with $s\left({\xi }_1\right)=0$ and ${\displaystyle u\left(s\left({\xi }_1\right)\right)\le \underset{0}{\overset{1}{\int }}\frac{a_0\left(\zeta \right)}{\zeta }d\zeta }$ for $s\left({\xi }_1\right)>0$. The last integral converges by the assumptions of Lemma 2, so $u\left(s\right(\xi \left)\right)\le 0$ for all $\xi \le {\xi }_{*}$, where ${\xi }_{*}$ satisfies the condition

$${\displaystyle D_2^0{\xi }_{*}=D_2^0{\xi }_1-\frac{1}{{\lambda }_c^{+}}{D}_1^0-\underset{0}{\overset{1}{\int }}\frac{a_0\left(\zeta \right)}{\zeta }d\zeta .}$$

Then it follows from the definition of $u\left(s\right)$ that $u\left(s\right(\xi \left)\right)\equiv 0$ for $\xi \le {\xi }_{*}$. Let $s\left({\xi }_1\right)=0$, ${\displaystyle \frac{{\phi }_1\left(s\right)}{s}}\gamma \left(s\right){|}_{s=0}=0$ (in this case $s\left(\xi \right)=0$ satisfies the first equation in (29)). If $s\left(\xi \right)$ is a solution of (29), then, by Lemma 1, the functions $u\left(\xi \right),s\left(\xi \right)$ are continuous in $\xi$. Let us consider a small neighborhood of the point ${\xi }_1$, assuming that at the point $\xi ={\xi }_1-\delta$, $\delta >0$ the inequality $s\left(\xi \right)>0$ holds. For $\xi \in [{\xi }_1-\delta ,{\xi }_1]$ from equation (30) we obtain

$${\displaystyle \frac{du}{d\xi }}\ge D_2^0-\underset{s,\theta }{min}\left({\left(p_0^{\prime }\right)}^{-}{\displaystyle \frac{f_1}{p_0{\lambda }_c}}\right){\displaystyle \frac{{\phi }_1\left(s\right)}{s}}{\phi }_2\gamma \ge {\displaystyle \frac{1}{2}}{D}_2^0$$

by choosing $\delta$ accordingly. Then $0=u\left(s\left({\xi }_1\right)\right)\ge {\displaystyle \frac{1}{2}}{D}_2^0\delta +u\left(s\left(\xi \right)\right)$, i.e. $u\left(s\right(\xi \left)\right)<0$ and, therefore, $s\left(\xi \right)=0$. Repeating this process, at the k -th step we obtain $s\left({\xi }_k\right)=0$, ${\xi }_k={\xi }_1-k\delta$, $k>1$. When k reaches a value for which the inequality $D_1^0\le {\lambda }_c{D}_2^{0}k\delta$ is satisfied, using (31), we obtain $s\left({\xi }_k\right)=0$, $\xi \in (-\infty ,{\xi }_1]$. The lemma is proved. □

4. Numerical Investigation

In the one-dimensional self-similar case, the system (1-3) is written in the form (8-13). Equations (7) for $\theta$ and (24) for s, which will be calculated numerically, will take the form:

$${\displaystyle \frac{d\theta }{d\xi }}=(c{\rho }_3^0(c_1-c_3)M\left(\theta \right)+(A_1{c}_1+A_2{c}_2)(\theta -{\theta }^{-})-\nu {\rho }_3^{0}c(\varphi -{\varphi }^{-}))/{\lambda }_c\equiv F_1\left(\theta \right),$$

(32)

$$\begin{array}{c}{\displaystyle \frac{ds}{d\xi }}=({\phi }_1{\phi }_2\gamma {\displaystyle \frac{p_0^{\prime }}{p_0}}{\displaystyle \frac{f_1}{{\lambda }_c}}+{\displaystyle \frac{1}{p_0}}\overline{g}{\phi }_1{\phi }_2+{\displaystyle \frac{1}{p_0}}\left|c\right|\varphi As-{\displaystyle \frac{1}{p_0}}\left|c\right|(\varphi -{\varphi }^{-})B)/a_0\left(s\right)\equiv F_2(s,\theta ).\end{array}$$

(33)

Since the coefficient $a_0\left(s\right)$ can turn to zero, then when dividing by it, $\epsilon$ is added and $max(a_0\left(s\right),\epsilon )$ is taken.

Algorithm for numerical solution of self-similar problem

Let us introduce a grid $\xi =ih$, $i=0,...,N$, h — step along the spatial coordinate.

We approximate the Newton’s boundary condition:

$$({\lambda }_c{\displaystyle \frac{{\theta }_1-{\theta }_0}{h}}+\alpha ({\theta }_0-{\theta }^0)-Q)=0.$$

(34)

Equation (32) is approximated by an implicit Runge-Kutta scheme of fourth order accuracy [20]

$${\theta }_{i+1}={\theta }_i+(h/6)(k_1+2k_2+2k_3+k_4),$$

(35)

where $k_1=h{F}_1\left({\theta }_i\right)$, $k_2=h{F}_1({\theta }_i+0.5k_1)$, $k_3=h{F}_1({\theta }_i+0.5k_2)$,

$k_4=h{F}_1({\theta }_i+k_3)$.

From here

$${\theta }_1={\theta }_0+(h/6)(k_1+2k_2+2k_3+k_4).$$

(36)

To find $\theta \left(0\right)$ the nonlinear equations (32-33) were solved by successive approximations with an inner cycle over k:

1) assuming $k=0$ for the first approximation, ${\theta }_0$ was taken equal to ${\theta }^{+}$, then ${\theta }_1$ was found from (36);

2) in the next iteration $k=k+1$, ${\theta }_1$ was substituted into (34), where ${\theta }_0$ was found from the boundary condition;

3) steps 1-2 are repeated until the conditions $\left|\right|{\theta }_0^{k+1}-{\theta }_0^k\left|\right|\le \epsilon$ and $\left|\right|{\theta }_1^{k+1}-{\theta }_1^k\left|\right|\le \epsilon$, where $\epsilon ={10}^{-6}$, are satisfied.

To plot the graphs of s and $\theta$ from the first two equations of the system (18-20), we express

$$s\left(0\right)={\displaystyle \frac{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(v_1^{+}+\frac{{\rho }_3^0}{{\rho }_1^0}{v}_2^{+})}{\varphi \left(\theta \left(0\right)\right)(v_1^{+}+v_2^{+})}},c={\displaystyle \frac{v_1^{+}+v_2^{+}}{(\varphi \left(\theta \left(0\right)\right)-{\varphi }^{-})(1-{\rho }_3^0/{\rho }_1^0)}},$$

that is, case 3 of the solution to the system is considered (18-20) for $v_1^{+}<0$, $v_2^{+}<0$.

The velocity c is substituted in the equation for the temperature (32), which is considered to be the Runge-Kutta method of the 4th order of accuracy.

The numerical algorithm for the temperature equation was verified by computational experiments. To simplify the calculations, $F_1=\theta$ was taken as a trial temperature function. An exact solution for this equation was analytically found. Runge’s rule was applied: three calculations were performed on grids with steps $h_1=h=0.1$, $h_2=h/2=0.05$, $h_3=h/4=0.025$.

The equation for water saturation (33) is also approximated by an implicit Runge-Kutta scheme of fourth order accuracy

$$s_{i+1}=s_i+(h/6)(k_1+2k_2+2k_3+k_4),$$

where $k_1=h{F}_2(s_i,{\theta }_i)$, $k_2=h{F}_2(s_i+0.5k_1,{\theta }_i)$, $k_3=h{F}_2(s_i+0.5k_2,{\theta }_i)$,

$k_4=h{F}_2(s_i+k_3,{\theta }_i)$.

$s_0$ is calculated using the formula $s_0={\displaystyle \frac{(\varphi \left({\theta }_0\right)-{\varphi }^{-})(v_1^{+}+\frac{{\rho }_3^0}{{\rho }_1^0}{v}_2^{+})}{\varphi \left({\theta }_0\right)(v_1^{+}+v_2^{+})}}$, taking into account the found ${\theta }_0$.

Figure 1 shows the results of the numerical modeling of temperature and water saturation with depth for the following parameters: $h=0.01$, $Q=7000$, $\alpha =0.2$, ${\lambda }_c=40$, ${\varphi }^{-}=0.62$, $v_1^{+}=v_2^{+}=-0.000001$, ${\theta }^{-}=268.15$, ${\theta }_1=270.15$, ${\theta }^{+}=273.15$, ${\theta }^0=278.15$. The graph of the porosity function is shown in Figure 2.

Figure 3 shows the effect of a given heat flux Q on the temperature at the boundary. It is evident that with increasing Q, the temperature of the upper snow layer also increases. Water saturation remains unchanged.

The results of modeling temperature variations with depth as a function of thermal conductivity ${\lambda }_c$ are presented in Figure 4. The graph shows that thermal conductivity directly influences snow heating: the higher the value of ${\lambda }_c$, the more intense the snow heating. Minor changes in water saturation are observed.

The specified porosity ${\varphi }^{-}$ influences water saturation at the boundary (Figure 5), as is evident from the formula for $s\left(0\right)$. The calculations were performed with the specified parameters ${\lambda }_c=35$, $Q=7000$, and three different values of ${\varphi }^{-}$. The temperature remains unchanged.

5. Conclusions

An exact self-similar solution to the problem of heat and mass transfer in snow and ice cover is constructed, taking into account variable porosity and phase transitions. The solution takes into account the influence of solar radiation flux at the boundary. For water saturation, the maximum principle and the lemma on the finite propagation velocity of disturbances are proven. A theorem on the existence of a weak solution to the problem is formulated. A numerical study of the problem of heat and mass transfer in snow is conducted, taking into account the external heat flux. The results of numerical modeling of temperature and water saturation are presented. The influence of a given heat flux at the boundary and the thermal conductivity coefficient on the temperature distribution over depth is demonstrated. The program is written in C++.

The work was carried out in accordance with the State Assignment of the Russian Ministry of Science and Higher Education entitled ‘Modern models of hydrodynamics for environmental management, industrial systems and polar mechanics’ (Govt. contract code: FZMW-2024-0003).