Studying the model of air and water filtration in a melting or freezing snowpack

The great interest to mathematical modelling of thermomechanical processes in and around snow is motivated by a large demand for an adequate description of snow behavior in order to calculate and forecast spring flood hydrographs and water quality in receiving reservoirs, to make assessment of the risk of avalanches in the mountains, etc. To date, there exists a number of works devoted to modelling the filtration of water and air in snow taking into account phase transitions, which use observational data and empirical dependencies, and which are based on various approaches in continuum mechanics and thermodynamics. A fairly extensive review on this topic can be found in [1]. The present article is devoted to a study of the mathematical model of air and water filtration in a snowpack in the presence of “ice-water” phase transitions. The snowpack is modeled as a three-phase continuum consisting of water, air and porous ice skeleton. In the article, the full model is posed in Section 1 and is called Model A. In Section 2, the reduction of Model A is given in the spatially one-dimensional case under some additional physical assumptions and the initial-boundary value problem is formulated for it, which is called Problem B. Section 3 is devoted to a study of this problem. Because of nonlinearity and high inter-connection between equations, the rigorous mathematical results on existence and uniqueness of classical or generalized solutions to Problem B are unavailable, at least, so far. Therefore, in this article we concentrate on building a reasonable approximation. At first, we construct the Rothe scheme for the problem and prove its local in time well-posedness, see Theorem 1. In fact, the classical solution of the Rothe scheme exists and is unique until the moment when the volumetric saturation of water either exceeds unity or drops below zero at some physical position, which is unacceptable from the physical viewpoint, see Remark 2 in Section 3. Secondly, we introduce the globally (in time) well-posed regularized version of the Rothe scheme and give a formal justification of the Rothe method for this scheme, see Theorems 2 and 3. In order to cimplete this introduction, it is worth to notice that the present study is rather close to a set of recent works devoted to the mathematical modelling of processes in snowpack, see [2–4], and can be regarded to as their extension.

1.    Mathematical Model of Snow: Multi-Dimensional Setting

At temperatures close to the freezing point of water, snow can be described as a three-phase medium consisting of water, air (water vapor) and ice. In this case, ice is a solid porous skeleton, and a mixture of water (in liquid form) and air is a two-component continuous medium filtering through the pores. We introduce a mathematical description of the balance of mass, momentum and heat based on the fundamental conservation laws in continuum mechanics, following the presentation from the monographs [5, 6].

We denote the phases of snow such that index i = 1 corresponds to water (in liquid form), index i = 2 corresponds to air (water vapor), and index i = 3 corresponds to ice.

The mass balance equation for each of the phases has the form

∂ρ i                    ³

"df + div x (p i U i ) = ^ I ji ,    i = 1, 2, 3,

j =1

where p i = p i ( x ,t ) is the reduced density of the i th phase, u i = u i ( x ,t ) is the velocity of the i th phase and I j _i is the intensity of mass transition from the j th phase to the i th phase. The reduced density ρ _i is related to the genuine density ρ _i ⁰ and volumetric concentration 3

a i by the formulas p i = a i p 0 , a i >  0 , i = 1, 2, 3 ; ^^a i = 1 . The intensities of phase i =1

transitions Iji meet the requirements Iji = — Iij. Further, we consider that filtration of water and air obeys to Darcy’s law and that the porous ice skeleton is immovable, i.e., vi = —Ko(ф)kOi(Si)(VxPi + p0g), i = 1,2,   U3 = 0,   Vt> 0,           (2)

µi where vi = vi(x, t) is the velocity of filtration, si = si(x, t) is saturation, i.e., the part of porous space occupied by the i th phase, pi = pi (x, t) is the hydraulic pressure, k0i = k0i(si) > 0 is the phase permeability such that k0i(0) = 0, ^i = consti > 0 is the dynamical viscosity of the i th phase; K0 = K0 (ф) is a nonnegative symmetric porous skeleton permeability tensor such that K0(0) = 0, ф = ф(х, t) is the snow porosity, i.e., the volumetric part of pores in the specific volume of snow, and g = —ge1 is the gravitational acceleration (g = const > 0, e1 = (1, 0, 0)t). The velocity of filtration vi is related to the medium velocity ui by the formula vi = ai ui and the volumetric concentration ai is related to the saturation si and the porosity ф by the formula ai = фsi, i = 1, 2. Introducing the notation s := s1 and taking into account that ai = 1, we write out these relations

i =1

in the form v1 = фsu1, v2 = ф(1 — s)u2, a1 = фs, and a2 = ф(1 — s). Here note that s2 = 1 — s, a3 = 1 — ф, 0 < s < 1, and 0 < ф < 1. The difference between the pressures of the filtering phases is determined by Laplace’s law p2 — p1 = pc(s), where the capillary pressure pc = pc(s) is the given function that has the properties pc(s) > 0, pc(0) = to, pc(1) = 0, and pC(s) < 0 [5, Ch. 5, Sec. 1.1].

Equations (2) describe the balance of momentum in the snowpack. Finally, assuming that the temperatures in all three phases coincide, i.e., 6 i ( x ,t) = 6 ( x ,t ) (i = 1, 2, 3 ), we write out the equation of balance of heat in the snowpack as follows:

(^ P i c i ) dt ⁺ (^ P i c i u i ) i =1                 i =1

V x O - diV x ( A c (s^) V x d ) = V V'" 6I ji . i =1 j =1

Here c _i = const _i > 0 is the specific heat capacity of the i th phase at constant volume and A _c (s, ф) is the heat conductivity of snow such that A _c (s, ф) >  A - = const > 0 , V s, ф G R, which is given according to experimental data.

The result of the ab ove considerations is the formulation of the following model of heat and mass transfer in the snowpack.

Model A. Equations of balance of mass (1), momentum (2) and heat (3) together with the set of the ab ove stated physical hypotheses and restrictions on coefficients constitute the basic multidimensional model of dynamics of the snowpack .

Remark 1. Note that the formulations based on modifications of Model A were considered by various authors. For example, in [2], a numerical study of a one-dimensional (in x) formulation of an initial-boundary value problem for a system of the form of Model A was carried out, in which, on the right-hand side of equation (3), the sum -νI13 takes place 33

instead of the sum

EEciBIji, where v = const > 0 is the given specific latent heat of i=1 j=1

the “ice-water” phase transition and the intensity of the phase transition I13 = I13(ф, B, s) is a given function of a very special form. Looking ahead, we note that in the setting considered in the present article, the value I13 is a sought function rather than a given one. The same model (with minor changes) as in [2] was considered earlier in [7], where the unique solvability was proved for the one-dimensional self-similar setting.

In the next section, based on the available phenomenological studies, we introduce the main assumptions regarding the coefficients and nonlinearities of Model A, make a reduction of this model in the spatially one-dimensional case, and formulate the initialboundary value problem for the reduced 1D model.

2.    Reduction of Model A. The 1D Model of Snowpack

Further, we restrict our considerations to the one-dimensional case such that Ox 1 -axis becomes the only spatial axis; we denote x := x 1 . Accordingly, all vector quantities in (1) - (3) become scalar ones, in particular, g = — g e 1 = — g . Matrix K ₀ in (2) 1 also becomes scalar, K ₀ = K 0 , and the operators V _x and div x become the derivative d/dx.

We accept some additional phenomenological hypotheses. Namely, taking into account natural observations (see, for example, [6, page 105]), we neglect sublimation, i.e., we impose the conditions I21 = I12 = 0 and I32 = I23 = 0. We postulate the model law of dependence of porosity on temperature following the article [7]: we assume that ф = ф(х, 9) is a given continuous piece-wise differentiable function such that

ф ( х, 9) = ф for 9 < 9 ,   ф' _е (x, 9 ) >  0 for 9 E [ 9 , 9 + ],   ф(х, 9 ) = ф + for 9 >  9 + , (4)

where Ф , Ф' = const E (0,1] , 9 - ,9 + = const , 0 <9 - <  9 + and 9 + is the temperature of ice melting. We impose the physically reasonable requirements ρ ⁰ ₂ < ρ ⁰ ₃ < ρ ⁰ ₁ on the genuine densities of phases and c ₁ > c ₃ on specific heat capacities. Denote by I the intensity of the “ice-water” phase transition I ₃₁ . We assume that p 0 ( i = 1, 2, 3 ) are constants.

Now, following the ideas from [5, Ch. V, § 1, Sec. 1], using the above hypotheses on I ij , φ , ρ i and c i , we reduce Model A to a system of five differential equations for the five sought functions s , θ , v , I , and p , where v is the total filtration velocity determined by the formula v := v 1 + v _o = sфu ₁ + (1 — s)фu _o , and p is the reduced pressure, which is defined below after formula (5e).

We introduce the new physical characteristics of a snowpack by the formulas

₍ м _v МЛ     ^k 01 ^(s)k 02 ^{(1 — s)}       ‘/x , ( \ _ ^k02 ^{(1 — s)}

(5a)

(,ф)'     0 (ф) М01« + MoW — s) pc():    Mok(s) ’ f M):= — Ko (ф)[ (k01^ p0 +       s p0) L, k(s):= '     + s, (5b)

µ 1 ¹ µ 2       ²                      µ 1           µ 2

K (»,ф) := k(s)K o (Ф), F кФ) := —    ⁽i_, ^ф f кФ) — KW) [ ^M1 s’ p 0 pl ,

(5c)

M o k(s)                       L    M o          J

1 n⁽s) := У s

^k 01 ⁽ (<) ~ M i k(<)    ^’

Q(s, Ф) := C 1 Sфp O + c o (1 — s)ФP ⁰ + С з (1 — Ф)Р ⁰ ’

(5d)

V^(v’ s’ ^Ф ’ ^Z ) ^:= c l P ⁰ [^{(1 — b(s))v — a(s}’ ^{ф)z —} F^(s’ ^ф)]+c o P ⁰ H^s’ ^ф)z+ ^{b(s)v + F (s}’ ^ф)] • ^(5e)

With the help of (5d), the reduced pressure is defined by the formula p := p _o + n(s) .

Obviously, the introduced functions a , b , f , k , K , F , π , Q , and V are uniquely defined by the given coefficients and free terms of Model A. Therefore, they are the given functions of their arguments rather than the sought ones. In turn, since p 2 is an unknown function in Model A, p = p(x, t) is an unknown function as well.

Similarly to [5, Ch. 5, Secs. 1.1–1.2], by rather lengthy but straightforward transformations, from Model A we derive the reduced one-dimensional model of a snowpack in terms of the unknown (sought) functions s , θ , v , I , and p and the given characteristics a , b , f , k , K , F , Q , and V , see equations (6a) – (6e) below.

Finally, we suppose that Q = (0,l) is an open interval in R _x , where the coordinate x = 0 corresponds to the bottom surface of a snowpack bordering on a frozen base (ground, ice, rooftop, etc.), and x = l is the upper surface bordering on an open air.

Now we are in a position to formulate the initial-boundary value problem for the reduced dynamical one-dimensional model of a snowpack .

Problem B. In the space-time domain QT = Q x (0,T), where T = const > 0 is a given time moment, find the water saturation in pores s = s(x, t), the reduced pressure p = p(x/t), the intensity of the “ice-water” phase transition I = I(x,t), the temperature of snow 6 = 6(x,t), and the velocity of filtration v = v(x,t), satisfying the reduced mass balance equations

^d[⁽¹ - s^)ф(x^,6)]

∂t

+     r(^s’ Ф(х, 6 )) It + b ( s ) v + F (s, Ф(x, 6))] = 0, (x, t ) G Q t , (6a)

∂x           ∂x

K^(s,^ф(x,⁶⁾⁾dx - f (s^cx,#)) = (i - pi) ^дф ( ^{х 6 )} ,    ^{( x,t ) g} Q t ,            ^(6b)

дф(x, 6)    I dt = p ’

(x, t ) G Q _t ,

(6c)

the reduced Darcy law v = -K(s,ф(x,6))dx + f (s,ф(x,6)),   (x,t) G Qt,               (6d)

the heat balance equation

^{Q ( s} , ^{ф(x, 6 ))}   + V fv, s, ^ф(x,6), Is) I⁶ -     fMs,^{ф(x, 6))} I⁶) =

∂t                ∂x ∂x ∂x             ∂x

-

⁽c i

-

С з ) 61,

(x,t) G Q T , (6e)

and the boundary and initial conditions s|x=0 = s0 (0, t),

P | x =0 = P 0 ^{( t )} ,

^{6 |} x =0 = ⁶ 0 ^(0,t), t ^{G (0,}T],

(6f)

^s | x = l      ^s 0 ^{( 1,} t),

∂p dx x=i

= 0, 6 | x = i = 6 0 (1,t), t G (0,T],

(6g)

s | t =0 = s o (x, 0),

6 | t =0 = 6 0 (x, 0),   x G [0,1].

(6h)

Here s0 = s0(x, t) and 60 = 60(x, t) are given functions on LIT and p0 = p0(t) is a given function on [0, T] satisfying the boundedness and regularity conditions s0 G [s-,s+], 60 G [6 ,6+] for (x, t) G LT,   s0, 60 G C 2+a (LT), p0 G C1+a[0,T ],   (6i)

where a G (0,1) , 6 — and 6 + are the same as in (4), s - , s ₊ G (0,1) are given constants, and L T = ({ 0 <  x <  1 } x { t = 0 }) U ({ x = 0, x = 1 } x (0, T] ) is the L -shaped part of d Q _t .

Note that, having the five sought functions s , p , I , θ and v found, we can determine separately the velocity of filtration of air v ₂ and the velocity of filtration of water v ₁ by the formulas

∂s

V 2 = a(s, ф(x, 6))— + b(s)v + F(s, ф(x, 6)),   v i = v - V 2 ,

∂x the air pressure p2 and the hydraulic pressure p1 in pores by the formulas

P 2 = p ^{- n(s)}, P i = P 2 ^- P c ^(s),                              ⁽⁸⁾

and the reduced densities of water (ρ1), air (ρ2) and ice (ρ3) by the formulas pi = sфp1, P2 =(1 - s)фp2, рз =(1 - ф)р3.                    (9)

3.    Rothe Scheme for Problem B 3.1.    Construction of Rothe Scheme

Problem B is an initial-boundary value problem for a system of essentially coupled nonlinear differential equations having no definite type, which makes it very difficult to study. In this regard, we study this system approximately, namely, we use the Rothe method, also known as the line method. The essence of the method is that the temporal interval [0, T] is divided into a set of subintervals [0,т] , [ t, 2 т ],..., [(M — 1) т, T] , where т is a small value ( M = T/т E N), derivative дФ/dt (for some Ф = Ф(t) ) is formally replaced by a finite-difference relation of the standard form A _T Ф(t) = (Ф^ф + т ) — Ф(t )) /т , t E [кт, (к + 1) t ] , к = 0,...- M , and on each of the intervals [кт, (к + 1)т] the coefficients in the equations are ‘frozen’ in a suitable way, that is, instead of the original coefficients, we take approximate ones that do not depend on t . The resulting system of equations with boundary conditions is called the Rothe scheme . It is approximated with respect to the original problem and depends on τ as on a parameter. Our closest aim is to build a Rothe scheme for Problem B.

We start by finding the initial data for the pressure function p ⁰ (x) = p(x, 0) . We substitute I from equation (6c) into (6e), then combine the resulting relation with equation (6b), namely, we express the derivative ∂θ/∂t in (6b) using (6c) and (6e), and take into account representation (5e). After quite lengthy, but simple technical transformations, we derive the equation

—^(s^»))u+pi(s-e-dd)dx=P2и^)- ^eqt.  <ю> where we denote

-'i-5 1""

^— ^) ( c i P i ^{(1 — b(s)) + c} 2 P 0 ^b(s) ) ^{K(s- ф(х- 9))ф} б ^{(x- 9)}I⁹ ρ ₁                                                ∂x

^{Q(s- ф(x- 9))} + ^(c i ^{— с} з ⁾⁹Р з ^ф е ^{(x- 9 )}

p 2 (s » -It dr) ^:=

∂x ∂x

^{df (s -ф(x} ,⁹⁾ ∂x

+

( 1 — 4 ) ф ; (x -9)

ρ ₁

^{Q ( s}, ^ф(x , б)) + ^(c i ^{— с} з ^{) 9} р З ^ф' е ^{(x- 9 )}

। (c 1 P ^{0 (1 — b(s))} + ^c 2 P 2 ^b(s) ) ^{f(s -ф(x -9))}

+ ⁽c 2 p 0 - c i p l ) ( a⁽s. ^ф(x. ⁹⁾⁾Is + F⁽s. ^ф(x. ⁹⁾⁾ )] I⁹ - -^ ( A c ^(s . ^ф(x. ⁹⁾⁾I⁹ ) L (12) ∂x              ∂x ∂x           ∂x

Inserting s = s ₀ (x- 0) and 9 = 9 ₀ (x- 0) into (10) and supplementing the resulting equation with boundary conditions (6f) 2 and (6g) 2 , we arrive at the formulation of the boundary value problem for the linear ordinary differential equation for the sought function p ⁰ (x) . Due to the classical theory of ordinary differential equations [8, Ch. XI, Sec. 2], this problem has a unique smooth solution. Indeed, due to (5d) ₂ , (4) and inequalities P 0 < Р з < P i and c _i > c ₃ , we have that Q(s- ф(x- 9)) + (c _i — с ₃ )9р 0 ф д (x- 9) >  Q(s- ф(x- 9)) >  min { c ₂ p 0 - c ₃ p 0 } > 0 for s-9 E [0-1] . Hence, according to (6i), both coefficients P _i | _t =₀ and

P ₂ | t =o are well defined and coefficient K | _t =0 is bounded from above and below by some positive constants independent of x . Inserting p ⁰ (x) , s ₀ (x, 0) and 9 ₀ ( x, 0) into (6d), we find v ⁰ (x) = v(x, 0) . Further in the Rothe scheme we consider that functions p ⁰ (x) and v ⁰ (x) are given alongside s ₀ (x, 0) and 9 ₀ ( x, 0) . Now we are in a position to construct the Rothe scheme itself.

We fix a small value т = T/M , M E N, insert I from (6c) into (6e) and fulfill the partial discretization of (6a), (6b), (6d), and (6e) following Rothe method. Therefore, we construct the following scheme for finding an approximate solution to Problem B.

On the n th step ( n = 0,..., M — 1 ), on the temporal segment { пт < t <  ( n + 1)т } , we find the approximate temperature 9 ⁿ⁺¹ = 9 ⁿ⁺¹ (x) as a solution to the problem

(Q (s ” ,ф( х ,9 " )) + (c i

. fln+1 - fl

— с о ) 9 ” р 0 Ф 6 ( x , 9 ⁿ ))----------

+ V lv " ,s " ,ф( х ,9 " )^1     '' — -d (a. (s " ,ф( x ,9 " ))' ' I

= 0, x E (0,l),

(13a)

dx dx dx                dx

9 ” ⁺¹L =0 = 9 0 (0, (n + 1)r), 9 " ⁺¹ | x = i = 9 0 (l, (n + 1)r),              (13b)

where s ⁿ = s ⁿ (x) , 9 ⁿ = 9 ⁿ (x) and v ⁿ = v ⁿ (x) are either the initial data (for n = 0 ), or the solutions of the Rothe scheme found on the preceding temporal segment (for n >  1 ).

Next, on { пт < t <  (n + 1)т } we find the approximate saturation s ⁿ⁺¹ = s ⁿ⁺¹ (x) as a solution to the problem

ф(х,9”)sn+L^n = AL,",ф(х,9”))ds_ + b(sn)„" + F(sn,ф(х,в"))1 — τ dx               dx n+1 n

— (1 — s ⁿ )ф 6 ( x ,9 ⁿ )-----------, x E (0,l), (14a)

τ s”+1|x=0= s0(0, (n + 1)т),   s”+1|x=l = s0(l,(n + 1)т),                (14b)

where 9n+1 is given as a solution to problem (13a) - (13b), and sn, 9” and v” are the given functions found on the (n — 1) th step. After this, on {пт < t < (n + 1)т} we find the approximate reduced pressure pn+1 = pn+1 (x) as a solution to the problem n+1

-K (sn,ф(x,9n)) dp— — f (sn,ф(x,9n)) = dxdx

0              n+1

1 — pj ф'6(x, 9”)----------, x E (0, l), dpn+1

pn+1|x=0 = Р0((п + 1)т), ——   = 0, dx x=l where sn, 9” and 9n+1 are given the same as for problem (14a) - (14b).

Finally, on { пт < t <  (n + 1)т } we find the approximate velocity of filtration v ⁿ⁺¹ from the reduced Darcy law (6d), in which the functions s ⁿ⁺¹ , 9 ⁿ⁺¹ and p ⁿ⁺¹ found above stand on respective places of the functions s , θ and p . Also, the approximate intensity of the phase transition I ” on the segment { (n — 1)т < t <  пт } for n >  1 is defined from the discrete approximation of equation (6c), more certainly, we have

I ” = р 0 Ф 6 x ,9 )(9 " ⁺¹ — 9”)/т.

After this, we iterate the above constructed procedure with n := n + 1 and find successively 6 ⁿ+² , s ⁿ⁺ 2 , p ⁿ⁺ 2 , and v ⁿ⁺ 2 on the temporal interval { (n + 1)т < t <  (n + 2)т } and I ⁿ⁺¹ on the temporal interval { пт < t <  ( n + 1)т } .

Definition 1. The set of problems (13a) – (13b) , (14a) – (14b) , (15b) – (15c) , (6d) (for v ⁿ⁺¹ ), and (16) to be solved successively for n = 0,1,... , is called the Rothe scheme for Problem B. The set of five functions {6 _T , s _T ,p _T , v _T : Q _T ^ R , I _T : Q _T-T ^ R} ( Q _T-T := Q x (0, T — т ) ) defined via the Rothe scheme by the formulas

Ф _т (x, t ) = Ф ⁿ (x) where Ф := 6, s,p, v for t E (( n — 1)т, пт ], 1 <  n <  M, or

Фт(x, t) = In(x) for t E ((n — 1)т, пт], 1 < n < M — 1, is called an approximate solution to Problem B (according to the Rothe scheme).

• 3.2.    On Justification of Rothe Method

The following result on local (in time) solvability of the Rothe scheme follows directly from the classical theory of ordinary differential equations.

Theorem 1. The Rothe scheme has a unique solution { 6 _T , s _T , p _T , v _T , I _T } belonging to the space Step ( 0, (n + 1)т; C ^2+a [0,l] ) 4 x Step ( 0, пт; C ^2+a [0,l] ) provided that e * <  s _T (x,t) <  1 — e * for all (x, t) E [0, l] x (0, пт ] , where e * = const > 0 , e * < 1 — e * .In particular, due to (6i) 1 , for any т > 0 the Rothe scheme has a unique solution at least on the first temporal segment, i.e., there exists a unique solution { 6 1 , s ¹ , p ¹ , v ¹ , 1 0 } E ( C ^2+a [0,l] ) 5 .

Here by Step (0,t ₀ ; C ^2+a [0,l]) we denote the subspace of U ” (0,t 0 ; C ^2+a [0,l]) consisting of step-functions (i.e., piece-wise constant functions) t ^ ф( - , t) .

Proof. The both assertions of the theorem follow directly from the well-known facts of the theory of linear second-order ordinary differential equations.

Remark 2. Note that due to the structure of the equations of Problem B and the Rothe scheme for it, we should not expect the fulfillment of the classical maximum principle for saturation in the form 0 <  s <  1 (resp., 0 <  s _T <  1 ) in the whole space-time continuum Q _T . In turn, due to representation (5a) 1 , coefficient a remains finite and positive for s = 0 only if k ₀₁ (s)p C (s) —> C ** = const > 0 , and for s = 1 the coefficient a vanishes, due to ^c    s→ +0

which equations (6a) and (14a) become degenerate. For s / [0,1] the functions p c , k ₀₁ and k 02 are not defined at all. From a physical viewpoint, the possibility of violation of the classical maximum principle is due to the very strong restrictions imposed in Section 2 on the given data of Model A. In particular, the classical maximum principle may fail due to the fact that the genuine densities of the components are considered to be constant and unequal to each other but at the same time the ‘ice-water’ phase transition is allowed. This, together with the law of conservation of the total mass, can in principle lead to the fact that the masses of pure phases become (within the framework of the model) negative. Accordingly, Problem B adequately describes real processes in a snowpack not more than for a limited range of values of thermomechanical characteristics. If the values of these characteristics are outside the range of reasonable values, for example, if s < 0 or s > 1 holds on some sets O C Q _T , Problem B fails to describe any real process on O .

Let us carry out an additional regularization of the Rothe scheme in order to obtain a ‘global’ result in time. For the values s- and s+ defined in conditions (6i), we introduce the regularized functions koi, ko2, Q and Pc by the respective formulas koi (s) := koi(s-) for s < s-, koi(s) := koi(s) for s E [s-,s+], koi(s) := koi(s+) for s > s+;

^k 02 (1 - s) := k o2 (1 - s — ) for s <  s — , ^k o2 (1 — s) := k o2 (1 — s) for s E [s — , s + ], ^k 02 (1 — s) := k o2 (1 - s + ) for s > s + ;

Q(s, ф( х , 9)) := Q ( s — , ф ( х , 9)) for s < s — , Q ( s, ф ( х , 9)) := Q ( s, ф ( х , 9)) for s E [s — , s + ], Q ( s^ ( x ,9 )) := Q ( s + ,ф ( x ,9 )) for s > s + ;

P c (s) := P c (s - ) + p ‘ ( s — )(s - s — ) for s < s — , p_c( s ) := p_c( s ) for s E [s — , s ₊ ],

P c ^{(s) :}= P c ^(s + ) + p c' ^(s + ^{)(s - s} + ) ^{for s} > ^s + .

Definition 2. Consider Problem B with function Q on the place of Q and with functions a, b, f, k, K, F , π, and V defined by formulas (5a) – (5d) and (5e) with the functions koi, ko2 and pc on the places of koi, ko2 and Pc, resp. We call this problem the regularized problem B and we say that the Rothe scheme for it is the regularized Rothe scheme.

Due to the form of the regularized functions ^k _oi , ^k _o2 and p c , we have that coefficients a and K in the regularized Rothe scheme are uniformly positive, i.e., there exists a positive constant C * independent of n such that a(s ⁿ ,ф ( х ,9^п )) >  C * , K ( s ⁿ , ф ( х , 9 ⁿ )) >  C * for all x E [0, l] , n = 0,..., M . Due to this and the well-known facts of the theory of linear ordinary differential equations, similarly to Theorem 1 the following result takes place.

Theorem 2. (On global solvability) For any given s _o , p _o and 9 _o satisfying conditions (6i) , for any fixed т >  0 (Т/т E N ) the regularized Rothe scheme has a unique solution { 9 _T , s _T , p _t , v _T , I _T } belonging to the space Step (0,T; C ^2+a [0,l]) ⁴ x Step (0,T — т ; C ^2+a [0,l]) .

Note that, similar to the original Rothe scheme, the possibility of a violation of the estimate 0 <  s _T <  1 on some subsets of O E Q _T is not excluded for the regularized Rothe scheme. On such subsets, the solution to the regularized Rothe scheme cannot be regarded to as some reasonable approximation of a real process in a snowpack.

Justification of the Rothe method, i.e., the passage to the limit as т ^ 0+ , on a rigorous mathematical level for Problem B or for the regularized problem B may not be achievable for the most general input data, due to the complexity of the equations. Therefore, we restrict ourselves to a formal justification for the regularized problem B, based on the introduction of the following condition.

Condition 1. For all sufficiently small т >  0 , there exists a constant C * independent of т such that the solution of the regularized Rothe scheme satisfies the uniform (in τ) estimate

^ ⁹ T IIc(Q?t ) + II s T IIc(Q) T ) + ^ ^P T ||C ( Q t ) + ||^d9 T /^dx^ L ^ ( Q _T ) + ||^ds T /^dx^ L ^ ( Q _T ) + ||^dP T /^dx^ L ^ ( Q T )

⁺ ll ^{d 9 t /dx} II l ² (Q _t ) ⁺ ll ^d s ^{T /dx} II l ² (Q _t ) ⁺ ll ^{d P t /dx} II l ² (Q _t ) ⁺ H ^ ^{T 9 t} H L ^{2 (q} t - t )

+ H A t s t H l 2 (Q t - t ) + H A t ( д9 т /dx) H L 2 (Q t - t ) + H A t ( ds T /dx) H L 2 (Q T-T ) <  C * .         (17)

Remark 3. It is worth noticing that Condition 1 fully corresponds to the case when porosity ф does not depend on 9 (i.e., ф = ф(х) ) and, as a consequence, there exist no phase transitions due to (6c), i.e., I = 0 . Indeed, in this case, from (6b), (6d) and (6g) we easily deduce that v = v ( t ) = f (s ₀ (l, t), ф(1)) . Substituting now v ( t ) and ф(х) into (6a), (6b) and (6e), we find that the system of equations in the formulation of the regularized problem B splits into a sequence of linear equations of uniformly parabolic and elliptic types, which should be solved successively to find first s , and then p and θ . Accordingly, the form of the Rothe regularized scheme is greatly simplified and estimate (17) is constructed for its solutions in a quite standard way, see, for example, [9, Ch. III, §5].

In cases where φ depends on θ , Condition 1 is postulated and the passage to the limit as т ^ 0+ is carried out for arbitrarily given ф satisfying requirement (4).

The following theorem gives a formal justification to the Rothe method.

Theorem 3. Let { θ τ , s τ , p τ , v τ , I τ } τ> 0 be the family of solutions to the regularized Rothe scheme (in the sense of Definitions 1 and 2) satisfying Condition 1; then there exist a sequence { т ^ 0+ } and five limiting functions { 9, s,p,v, I } such that

θτ -→ θ , sτ -→ s, pτ -→ p strongly in L2(0, T; W21(0, l)) П C(Qt), weakly in L2(0, T; W22(0, l)),(18)

AT9T —> d9/dt,  ATsT —> ds/dt,  IT —> I weakly in L2(O,T-5),  V5> 0,

τ→0                τ→0

vT —> v strongly in L2(0,T; W21(0,l)),

τ → 0 ²

and the five limiting functions { θ, s, p, v, I } are a strong generalized solution to the regularized problem B.

Proof. Limiting relations (18) – (20) follow from (17) due to the Alaoglu, Ascoli – Arcel and Kolmogorov – Riesz theorems. In turn, these limiting relations allow to pass to the limit in the regularized Rothe scheme and thus establish that the limit functions serve as a solution to the regularized problem B.

Conclusion

For a highly nonlinear one-dimensional problem of description of thermomechanical processes in a snowpack, the approximate Rothe scheme is constructed. This scheme is uniquely solvable locally in time and its regularized version is uniquely solvable globally in time. Under additional condition, it is established that the family of solutions to the regularized Rothe scheme converges to the solution of the original problem strongly in suitable functional spaces.

Acknowledgment. The work was carried out in accordance with the State Assignment of the Russian Ministry of Science and Higher Education entitled “Modern methods of hydrodynamics for environmental management, industrial systems and polar mechanics” (Govt. contract code: FZMW-2020-0008, 24 January 2020).