ISSN 1515-4904

MAT

SERIE A: CONFERENCIAS, SEMINARIOS Y

TRABAJOS DE MATEMATICA

No. 6

SOME PROBLEMS CONCERNING WITH MASS AND HEAT TRANSFER IN A MULTI-COMPONENT SYSTEM

FEDERICO TALAMUCCI Dipartimento di Matematica “U. Dini”,

Università di Firenze,

Viale Morgagni, 67/A, 50134 Firenze, ITALY.

E-mail: talamucci@math.unifi.it Rosario, Enero 2002

Abstract. Continuum mechanics is largely used to approach physics-mathematical problems. In this paper we concentrate our attention on a multi-component system in a non-uniform thermal field and we discuss from a general point of view the conservation equations. A special care is devoted to the possible change of phase of some components of the medium and to the conditions at the material surfaces (possibly free boundaries) which are very often encountered in such problems. As an application, we also present two different physical situations (ground freezing, non-standard filtration in a porous medium) which are modelled from the mathematical point of view according to the general equations we found out.

Resumen. La mecánica del continuo es muy usada en la resolución de problemas de la física-matemática. En este trabajo la atención se concentra sobre un sistema multi-componente en un campo térmico no uniforme y se discuten las ecuaciones de conservación desde un punto de vista general. Una atención particular est dedicada a los posibles cambios de fases de algunas componentes del medio y a las condiciones sobre las superficies materiales (fronteras libres) que se encuentran muy a menudo en tales problemas. Como una aplicación, se presentan dos diferentes situaciones físicas (congelamiento de suelos porosos, filtación no standard con remoción de substancia de la matriz porosa) cuyos modelos matemáticos se escriben siguiendo el esquema general introducido.

Riassunto. La meccanica dei continui rappresenta un approccio molto diffuso nei problemi di fisica matematica. In questo lavoro si focalizza l’attenzione su un sistema composito in un campo termico non uniforme e si discute da un punto di vista generale le equazioni di conservazione. Particolare attenzione viene riposta ai possibili cambiamenti di fase di alcuni componenti del mezzo e alle superfici materiali (frontiere libere) che spesso sono presenti in problemi di questo tipo. Come applicazione, si presentano due differenti problemi fisici (congelamento dei suoli porosi, filtrazione non standard con rimozione di sostanze dalla matrice porosa), i cui corrispondenti modelli matematici sono scritti seguendo lo schema generale introdotto.

Keywords: Continuum Mechanics, Conservation Equations, Free Boundaries, Ground Freezing, Filtration in a Porous Medium.

Palabras claves: Mecánica del Continuo, Ecuaciones de Conservación, Fronteras Libres, Congelamento de Suelos, Filtracioón en un Medio Poroso.

AMS Subject Classifications: 35R35, 76S05, 76T05, 76R50, 80A20.

The present paper is the content of a mini-course given by Dr. F. Talamucci at the Department of Mathematics of FCE-UA, Rosario, February 14-19 2001, in the frame of the Project of Scientific Cooperation “Free Boundary Problems for the Heat Equation”, jointly promoted by CONICET (Argentina)-CNR (Italy). The manuscript has been received and accepted on October 2001.

Se publica en este esemplar el contenido de un minicurso dictado por el Dr. F. Talamucci en el Departamento di Matemática de la FCE-UA, Rosario, del 14 al 19 febrero de 2001, durante su visita realizada en el marco del Pryecto de cooperación Científica “Problemas de Frontera Libre para la Ecuación del Calor” a través del convenio CONICET (Argentina)-CNR (Italia). El manuscrito fue recibido y aceptado en Octubre de 2001.

Some problems concerning with mass and heat transfer in a multi-component system

F. Talamucci ^{∗}

Introduction

The aim of the paper is to give som examples of mathematical models which arise from different physical situations, having however in common several formal aspects. Roughly speaking, we will consider heat and mass transport phenomena in a multi-component medium. The first part (Section 1) is devoted to introduce the conservation equations, which will be discussed from a general point of view. Secondly, we will consider two special problems which can be framed into this formalism. The first one (Section 2) refers to ground freezing: a porous medium saturated with water is subjected to a freezing process and a volume expansion of the medium is observed (frost heave). The second example (Section 3) is a mathematical model for a filtration process under high pressures with removal of small particles away from the porous matrix. The starting point of the formulation of the latter problem has been the understanding of the complex phenomena taking place during the preparation of the expresso coffee.

1 The continuum approach

We consider a medium consisting of several components, say N species. We adopt the continuum approach in order to describe the properties of the medium. The actual multi-components medium is replaced by a fictious continuum: we assign kinematic and dynamic variables and parameters that are continuous functions with respect to the spatial coordinates. Arguing as in [1], we take a volume ΔV centred in some point P of the medium. The species k, k = 1, …, N occupies the volume ΔV_{k}. The typical behaviour of the ratio ΔV_{k}/ΔV is the following:

[figure: ΔV_{k}/ΔV; ΔV; −; 1; ΔV_{0}]

Below a certain value ΔV_{0} the fluctuations in the ratio are large. From the mathematical point of view the ratio tends to 0 or to 1, as ΔV approaches 0. The fluctuations tend to decay in a interval starting from ΔV_{0} and for large values of ΔV the ratio may undergo gradual changes. We call ΔV_{0} the representative elementary volume (REV) which corresponds to the physical (or material) point of the porous medium at the mathematical point P.

∗Dipartimento di Matematica “U. Dini”, Università di Firenze, Viale Morgagni, 67/a, 50134 Firenze, Italy.

1.1 Volumetric fractions, densities and velocities

We define the volumetric content of the species k as

μ_{k} = lim_{ΔV →ΔV_{0}} \frac{ΔV_{k}}{ΔV} (1.1)

Obviously ∑_{k=1}^{N} ΔV_{k} = 1.

If Δm_{k} is the mass of the k-species contained in ΔV_{k} we define

ϱ_{k} = lim_{ΔV →ΔV_{0}} \frac{ΔV m_{k}}{ΔV} bulk density (1.2)

d_{k} = lim_{ΔV →ΔV_{0}} \frac{ΔV m_{k}}{ΔV_{k}} specific density (1.3)

In other words, the bulk density is the mass of the k-species in a unit volume of the mixture, the specific density is the mass of the k-species with respect to a unit volume occupied by the same species. The following relation can be easily obtained from (1.1), (1.2) and (1.3):

ϱ_{k} = d_{k}μ_{k}, k = 1, …, N.

If a species is moving, we call v_{k} the corresponding molecular velocity, averaged in a REV. The volumetric flux or volumetric velocity is defined as q_{k} = μ_{k}v_{k} (1.4)

The physical meaning is the following: q · n is the volume of the k-species passing through a unit surface with normal n in the unit time.

1.2 The general conservation equation

Let us consider any quantity defined in a volume D with concentration (i. e. amount for unit volume) g(x, t) and current j (i. e. j · n is the amount passing through a unit surface with normal n in the unit time). The general conservation principle writes

\frac{∂}{∂t} ∫_{V} g(x, t)dx = − ∫_{∂V} j · ndσ + ∫_{V} I(x, t)dx (1.5)

where V is an arbitrary subset of D and I is the production (or loss) of the quantity g per unit volume and unit time. Equation (1.5) entails

\frac{∂g}{∂t} + divj = I (1.6)

1.2.1 Mass conservation

If one takes g = ϱ_{k}, the corresponding current is j = ϱ_{k}v_{k}. Hence, applying (1.6) we have:

\frac{∂}{∂t}ϱ_{k} + div(ϱ_{k}v_{k}) = I_{k} or \frac{∂}{∂t}(d_{k}μ_{k}) + div(d_{k}q_{k}) = I_{k}, k = 1, …, N (1.7)

where I_{k} is the rate at which the k-mass is produced. By summing up equations (1.7) we find

\frac{∂}{∂t}ϱ + div(ϱv) = 0 (1.8)

where ϱ = ∑_{k=1}^{N} ϱ_{k} and ϱv = ∑_{k=1}^{N} ϱ_{k}v_{k} (composite mass velocity).

Remark 1.1 If it is known that I_{r} + I_{s} = 0 for some index r, s, then it may be convenient to write (1.6) in the form

\frac{∂}{∂t}(ϱ_{r} + ϱ_{s}) + div(ϱ_{r}v_{r} + ϱ_{s}v_{s}) = 0 (1.9)

1.2.2 Heat conservation

We define the heat energy per unit volume as

Q_{V} = ∑_{k=1}^{N} ϱ_{k} (∫_{0}^{T} c_{k}(y)dy + E_{k}(T)) (1.10)

where c_{k}, k = 1, …, N, are the specific heats per unit mass of each species and the quantities E_{k} are related to the possible presence of latent energy at some temperature. In order to explain (1.10), we consider the case of one single material of density ϱ and mass Δm contained in a volume ΔV. If the temperature increases from T to T + dT, we have

\frac{ΔV (T)}{ΔV (T + dT)} = \frac{ϱ(T + dT)}{ϱ(T)} (1.11)

If Q is the heat amount per unit mass, we have by definition

Q = c^{∗}(T)T, c^{∗}(T) = \frac{1}{T} ∫_{0}^{T} c(y)dy (1.12)

The heat contained in Δm passes from ΔQ_{m} = QΔm to ΔQ_{m} + dΔQ_{m}, with

dΔQ_{m} = Δmc(T)dT. (1.13)

Computing now the variation of heat in the volume ΔV, we have from (1.11)-(1.13):

dΔQ_{V} = ΔV (ϱ(T + dT)c^{∗}(T + dT)(T + dT) − ϱ(T)c^{∗}(T)T) (1.14)

Recalling (1.12) and considering the total amount of heat at a some temperature T in a unit volume where several species are present, we finally get (1.10). At this point, we use the conservation equation (1.6) with

g = Q_{V} and j = −k∇T + Q_{c} (1.15)

where k is the thermal conductivity (depending on the properties of the medium and possibly on T) and Q_{c} is the energy flux due to convection:

Q_{c} = ∑_{k=1}^{N} ϱ_{k} (∫_{0}^{T} c_{k}(y)dy + E_{k}) v_{k}. (1.16)

Equation (1.6) applied to energy conservation, i. e.

\frac{∂Q_{V}}{∂t} + div(−k∇T + Q_{c}) = 0

can be written in the following form, by making use of mass balances (1.7):

(C(T) + ∑_{k=1}^{N} ϱ_{k} \frac{∂E_{k}}{∂T} ) \frac{∂T}{∂t} + div (−k∇T) + (1.17)

+ ∑_{k=1}^{N} (c_{k} + \frac{∂E_{k}}{∂T} ) ϱ_{k}v_{k} · ∇T + ∑_{k=1}^{N} I_{k} (∫_{0}^{T} c_{k}(y)dy + E_{k}) = 0

In (1.17) we defined C(T) = ∑_{k=1}^{N} ϱ_{k}(T)c_{k}(T) (equivalent heat capacity).

1.3 Conditions at the interfaces

The presence of material surfaces at which some quantities involved may have a discontinuity corresponds to a very frequent situation in such problems. We model an interface as a moving surface S described by the equation S(x, t) = 0. We consider again the balance (1.7), with I = 0. Assuming that the balance holds on both sides of an open set in R^{4} cut by S, we get the Rankine-Hugoniot condition (see, e. g. , [2])

[[g]] \frac{∂S}{∂t} + [[j]] · ∇_{x}S = 0 (1.18)

where [[ ]] denotes the jump across the interface at any point x_{0} ∈ S:

[[f]] = lim_{x→x_{0}}^{+}f(x) − lim_{x→x_{0}}^{−}f(x)

The symbols ^{+} and ^{−} denote the two regions separated by the interface. Equation (1.18) can be written in the following form: [[g]] V · n = [[j]] · n (1.19)

where V is the velocity of a point on the interface and n = ±∇_{x}S/|∇_{x}S| is the unit normal.

1.3.1 Mass balance at any interface

Applying (1.19) to equation (1.8) we have

[[ϱ]] V · n = [[ϱv]] · n. (1.20)

Equation (1.20) also holds for any species r such that I_{r} = 0. In the same case of remark 1.1, we can write

[[ϱ_{k} + ϱ_{h}]] V · n = [[ϱ_{k}v_{k} + ϱ_{h}v_{h}]] · n (1.21)

1.3.2 Heat balance at any interface

The energy jump condition comes from equations (1.19) and (1.15):

∑_{k=1}^{N} (∫_{0}^{T} c_{k}(y)dy + E_{k})^{+} ([[ϱ_{k}]] V · n − [[ϱ_{k}v_{k}]] · n) = [[−k∇T · n]] + (1.22)

+ ∑_{k=1}^{N} ϱ_{k}^{−} (v_{k} · n − V · n)^{−} [[∫_{0}^{T} c_{k}(y)dy + E_{k}]] .

Assuming, as it is natural, that [[T]] = 0 across the interface, we can reduce (1.22) to

∑_{k=1}^{N} (∫_{0}^{T} c_{k}(y)dy + E_{k}) ([[ϱ_{k}]] V · n − [[ϱ_{k}v_{k}]] · n) = [[−k∇T · n]] (1.23)

2 Freezing in a porous medium

The first example we are going to consider deals with a very well known example of coupled heat-mass transfer process. When a moist soil saturated with water is subjected to a freezing process, a volume expansion is generally observed. The volume increase is not only due to the different specific volumes of water and ice (in situ water freezing), but mainly to freezing of transported water from the unfrozen part of the soil towards a freezing front. The phenomenon, known as frost heave, is the effect of strongly coupled heat and mass transfer processes occurring in the porous soil. Frost heave damages produced on road pavements, pipelines, building foundations in regions with periodical or permanent freezing are well known. On the other hand, the artificial freezing technique for consolidating soils during tunnelling and underground constructions has become more and more used in recent years: predicting the swelling of the soil under specified thermal and mechanics conditions is of crucial importance in such practises. The two aspects of the problem we mentioned are sufficient to explain the attention that industrial and technological research has been dedicating to ground freezing in recent years and the relevant amount of models that have been proposed.

2.1 The physical process

We consider a sample of fine grained soil saturated with water and subjected to a freezing process. The most evident phenomenon occuring during the process is the volume expansion of the medium (frost heave). The heave amount of a frozen specimen can be very considerable: in the tests of ([17]), for istance, the final height of a specimen of Kanto Loam is more than three times the initial one. The most relevant factors controlling the process of frost heave are the freezing rate, the property of the soil (grain size, dry density, permeability, ...), the overburden pressure. The physical problem can be outlined as in the following picture:

[figure: frozen soil; unfrozen soil; freezing region; soil grains and ice; soil grains and water; cooling temperature; heave of soil; change of phase; base of soil; top of soil; POROUS; MEDIUM; water migration]

fig. 2.1: the physical situation.

The freezing process is induced by a cooling temperature acting on the top of the moist soil. The base of the soil is permeable to transfer of liquid from the underground. Both the advancement of the frozen soil on the unfrozen part (heat transfer effects) and the migration of water from the base towards the intermediate region where the change os phase takes place (mass transfer effects) determine the volume expansion. As a first approximation, we may believe that frozen and unfrozen soils are separated by a sharp interface: in this case, the change of phase occurs only at that surface. However, a more careful observation shows the existence of a thin transitional zone, called frozen fringe (region FF in figure 2.2), where water and ice coexist in the porous space and where change of phase takes place. It is comprised between the 0^{o}C isotherm (freezing front) and the lower boundary of the totally frozen soil. As the frozen part advances, if particular thermal and mechanicical conditions are verified, the formation of successive pure ice layers (denoted by L in figure 2.2) may occur. The phenomenon is known as ice segregation or ice lensing. The final configuration of the sample of soil consists, in this case, in a sequence of pure ice layers alternating with regions of completely frozen soil. We may distinguish two different phases in the process: whenever a segregated pure ice layer (region N L in figure 2.2a) is growing immediately over the frozen fringe, we refer to a lens formation process. On the other hand, in case of rapid freezing or remarkable overburden pressure acting on the soil, the lens stops growing and water in the porous space freezes and ice is kept within the porous matrix, even though the porosity of the soil may be affected by the passage of the freezing front. We refer to such a process as frost penetration. The region F 1 immediately over the frozen fringe (figure 2.2b) consists in this case of completely frozen soil.

U grains + water FF grains + water + ice N L forming ice lens F grains + ice L ice F grains + ice L ice F grains + ice L ice F grains + ice

2.2a Lens formation

U grains + water FF grains + water + ice F1 grains + ice L ice F grains + ice L ice F grains + ice L ice F grains + ice L ice F grains + ice

2.2b Frost penetration

fig. 2.2: the two different phases of the freezing process

2.2 The mathematical model: conservation equations

We will consider only one dimensional models, according to the fact that the quantities involved in the process vary appreciabily only along the vertical direction z. The porous medium is formed by three species: water, ice and soil grains: recalling the notations introduced in section 1, we refer to index 1, 2, 3 respectively to those three components and we will use the index w, i, s. For the sake of convenience, we introduce the quantities ν and ε such that (cfr. (1.1))

μ_{1} = ε(1 − ν), μ_{2} = εν, μ_{3} = 1 − ε. (2.1)

The quantities ν ∈ [0, 1] and ε ∈ [0, 1] correspond to the volume fraction of ice (w. r. t. the porous space) and to the porosity of the soil, respectively. Moreover, we call ρ_{w} = d_{1}, ρ_{i} = d_{2}, ρ_{s} = d_{3} (see (1.3)), q_{w} = μ_{1}v_{1}, q_{i} = μ_{2}v_{2}, q_{s} = μ_{3}v_{3} (see (1.4)). We assume that no chemical reaction are present within the soil. Hence, we have simply: I_{1} + I_{2} = 0, I_{3} = 0. (2.2)

Mass conservation, following (1.7) and (1.9), is obtained by taking into account of (2.2):

(2.3)  \frac{∂}{∂t} ((1 − ν)ερ_{w} + ενρ_{i}) + \frac{∂}{∂z} (ρ_{w}q_{w} + ρ_{i}q_{i}) = 0 (water + ice)
(2.4)  \frac{∂}{∂t} ((1 − ε)ρ_{s}) + \frac{∂}{∂z} (ρ_{s}q_{s}) = 0 (soil grains)

Energy conservation comes from (1.17), by taking E_{1}(T) − E_{2}(T) = L, with L latent heat of water per unit mass, E_{3} = 0:

C(ν, T)\frac{∂T}{∂t} + \frac{∂}{∂z} (−k \frac{∂T}{∂z} ) + (ρ_{w}q_{w}c_{w} + ρ_{i}q_{i}c_{i}+ (2.5)

+ρ_{s}q_{s}c_{s}) \frac{∂T}{∂z} + I_{1} (∫_{0}^{T}(c_{w}(y) − c_{i}(y)dy + L) = 0

where C(ν, T) (equivalent heat capacity) is given by

C(ν, T) = (1 − ν)ερ_{w}(T)c_{w}(T) + ενρ_{i}(T)c_{i}(T) + (1 − ε)ρ_{s}(T)c_{s}(T) (2.6)

and

I_{1} = − (\frac{∂}{∂t} (ενρ_{i}) + \frac{∂}{∂z} (ρ_{i}q_{i})) = (\frac{∂}{∂t} (ε(1 − ν)ρ_{w}) + \frac{∂}{∂z} (ρ_{w}q_{w})) .

We remark in (2.5) the dependence on the quantities E_{w} and E_{i} by only their difference L.

2.2.1 Conditions at the interfaces

At any interface z = ẑ mass balance is (see (1.20) and (1.21)):

[[ρ_{w}(1 − ν)ε + ρ_{i}νε]] ˙ ẑ = [[ρ_{w}q_{w} + ρ_{i}q_{i}]] (2.7)

[[(1 − ε)ρ_{s}q_{s}]] ˙ ẑ = [[ρ_{s}q_{s}]] (2.8)

Energy balance at z = ẑ is (1.23), which reduces, by virtue of (2.2):

(∫_{0}^{T} c_{w}(y) − c_{i}(y)dy + L) ([[ερ_{w}(1 − ν)ε]] ˙ ẑ − [[ρ_{w}q_{w}]]) = [[−k \frac{∂T}{∂z} ]] (2.9)

or, equivalently:

(∫_{0}^{T} c_{w}(y) − c_{i}(y)dy − L) ([[−ρ_{i}νε]] ˙ ẑ − [[ρ_{i}q_{i}]]) = [[k \frac{∂T}{∂z} ]] (2.10)

2.3 The equations in each region of the soil

We are going now to write the equations in the different parts of the soil, namely the unfrozen soil (ν = 0), the frozen fringe (0 < ν < 1) and the totally frozen soil (ν = 1), which may include layers of pure ice (ε = 1). The base of the soil corresponds to z = 0, the top is denoted by z_{T}. Moreover:

z_{F} = sup_{0≤z≤z_{T}} {z : ν(z, t) = 0}, z_{S} = inf_{0≤z≤z_{T}} {z : ν(z, t) = 1}.

The two boundaries contain the frozen fringe. If only ice is present immediately over z_{S}, then a lens formation process is occurring, otherwise frost penetration is taking place. Let us examine the two different cases of lens formation and frost penetration separately.

2.3.1 Lens formation

In the region U = {z | 0 < z < z_{F} } (unfrozen soil, ν = 0) it is generally assumed that the porous matrix is undeformable, hence the porosity depends only on z (by eq. (2.4)) and we may assume that ε = ε_{0}, constant. In that case, the water flux q_{w} = q_{w}^{u} depends only on t, owing to (2.3). Energy conservation (2.5) reduces to

C(0, T)\frac{∂T (z, t)}{∂t} + \frac{∂}{∂z} (−k_{u}(T)\frac{∂T (z, t)}{∂z}) + ρ_{w}c_{w}q_{w}^{u}(t)\frac{∂T (z, t)}{∂z} = 0 (2.11)

If even the porous matrix in the frozen fringe FF is rigid, then we conclude, by using (1.20) with ẑ = z_{F} that ε = ε_{0} in the frozen fringe. Equations (2.3) and (2.5) evaluated in FF take the form

ε_{0}(ρ_{i} − ρ_{w})\frac{∂}{∂t}ν(z, t) + \frac{∂}{∂z} (ρ_{w}q_{w}(z, t) + ρ_{i}q_{i}(z, t)) = 0 (2.12)

C(ν, T)\frac{∂T (z, t)}{∂t} + \frac{∂}{∂z} (−k_{ff}(T)\frac{∂T (z, t)}{∂z}) + (ρ_{w}c_{w}q_{w}(z, t) + ρ_{i}c_{i}q_{i}(z, t))\frac{∂T (z, t)}{∂z} + (2.13)

+ρ_{i} (ε_{0} \frac{∂}{∂t}ν(z, t) + \frac{∂}{∂z} q_{i}) ((c_{i} − c_{w})T − L) = 0

Let us consider now the completely frozen region ν = 1, which is made by layers of porous soil saturated with ice (regions F in figure 2.a) alternating with layers of pure ice (regions L). It can be assumed that in the regions F there is not relative motion of ice with respect to the soil grains (i. e. v_{i} = v_{s}). Hence, (2.5) in L reduces to

ρ_{i}c_{i} \frac{∂T}{∂t} + \frac{∂}{∂z} (−k_{i} \frac{∂T}{∂z} ) + ρ_{i}c_{i}ż_{T} \frac{∂T}{∂z} = 0 (2.14)

while in F mass and heat conservation are respectively

\frac{∂ε}{∂t} + ż_{T} \frac{∂ε}{∂z} = 0 (2.15)

C(1, T)\frac{∂T}{∂t} + \frac{∂}{∂z} (−k_{f} \frac{∂T}{∂z} ) + C(1, T)ż_{T} \frac{∂T}{∂z} = 0 (2.16)

An useful boundary condition can be obtained by expressing q_{i}^{+} in (2.7) (evaluated for ẑ = z_{S}) in terms of the heave rate ż_{T}. This can be done by integrating (w. r. t. z) (2.3) and (2.4) in the regions L and F and taking account of the boundary balances (2.7) and (2.8). We get the following formula

ρ_{i}ż_{T} − ρ_{i}q_{i}(z_{S}, t)^{−} − ρ_{w}q_{w}(z_{S}, t)^{−} − ε_{0}(1 − ν(z_{S}, t)^{−})(ρ_{i} − ρ_{w})ż_{S} = 0 (2.17)

which allows us to evaluate the heave rate independently of the quantities concerning with z_{S}^{+}, in particular no matter is the porosity ε(z_{S}, t)^{+}. Boundary conditions (2.7)-(2.9) write, in case of lens formation,

ρ_{w}(q_{w}(z_{F} , t)^{+} − q_{w}^{u}(t)) + ρ_{i}q_{i}(z_{F} , t)^{+} − ε_{0}(ρ_{i} − ρ_{w})ν(z_{F} , t)^{+}ż_{F} = 0 (2.18)

ρ_{i} ( ∫_{0}^{T (z_{F} ,t)}(c_{i}(y) − c_{w}(y))dy − L) (ε_{0}ν(z_{F} , t)^{+}ż_{F} − q_{i}(z_{F} , t)^{+}) = [[−k(z_{F} , t)\frac{∂T}{∂z} (z_{F} , t)]]_{−}^{+} (2.19)

ż_{S} = 0, ρ_{i}ż_{T} = ρ_{i}q_{i}(z_{S}, t)^{−} + ρ_{w}q_{w}(z_{S}, t)^{−} (2.20)

ρ_{w} ( ∫_{0}^{T (z_{S},t)}(c_{i}(y) − c_{w}(y))dy + L) q_{w}(z_{S}, t)^{−} = [[−k(z_{S}, t)\frac{∂T}{∂z} (z_{S}, t)]]_{−}^{+} (2.21)

Moreover, passing from F to L and vice versa we have

[[−k(T (ẑ, t))\frac{∂T}{∂z} (ẑ, t)]]_{−}^{+} = 0 (2.22)

and at any interface the temperature is continuous:

[[T]]_{−}^{+} = 0 (2.23)

We incidentally remark that (2.18) implies that the water flux q_{w} is continuous w. r. t. z at z_{F} only if ν(z_{F} , t)^{+} = 0. We also stress the fact that during the formation of a lens the base of the lens is at rest (condition (2.20).

2.3.2 Frost penetration

Equations (2.11) in region U, (2.12) and (2.13) in FF, (2.14) in L, (2.15) and (2.16) in F are still the same also in the case of frost penetration. In region F1, which lies immediately over FF, some modeler assumes that there is a relative motion of ice with respect to the porous matrix, due to regelation (v_{i} ≠ v_{s}). In that case, mass and heat conservation in F1 are:

\frac{∂ε}{∂t} + \frac{∂q_{i}}{∂z} = −\frac{∂ε}{∂t} + \frac{∂q_{s}}{∂z} = 0

C(ν, T)\frac{∂T}{∂t} + \frac{∂}{∂z} (−k_{f} \frac{∂T}{∂z} ) + (ρ_{i}c_{i}q_{i} + ρ_{s}c_{s}q_{s}) \frac{∂T}{∂z} = 0

Boundary conditions (2.18), (2.19), (2.22) and (2.23) are still the same, while (2.20), (2.21) are replaced by

ρ_{i}ż_{T} − ρ_{i}q_{i}(z_{S}, t)^{−} − ρ_{w}q_{w}(z_{S}, t)^{−} − ε_{0}(1 − ν(z_{S}, t)^{−})(ρ_{i} − ρ_{w})ż_{S} = 0 (2.24)

ρ_{w} ( ∫_{0}^{T (z_{S},t)}(c_{i}(y) − c_{w}(y))dy + L) (q_{w}(z_{S}, t)^{−}(1 − ν(z_{S}, t)^{−})ε_{0}ż_{S}) = [[−k(z_{S}, t)\frac{∂T}{∂z} (z_{S}, t)]]_{−}^{+} (2.25)

Remark 2.1 If one assumes that v_{i} = v_{s} also in F1, then in that region heat conservation is obviously (2.13). In this case the porosity ε has a jump at z = z_{S}, namely (see (2.17))

ε(z_{S}, t)^{+} = \frac{ż_{T} − ε_{0}ż_{S}}{ż_{T} − ż_{S}}

On the contrary, if one assumes that v_{S}(z_{S}, t)^{+} = 0, then ε(z_{S}, t)^{+} = ε_{0}, which means that the porous matrix is not deformed by the passage of the front. In that case ice has to move through the porous space at a speed v_{i}(z_{S}, t)^{+} = ż_{T} /ε_{0}.

2.4 Physical assumptions

In order to close the system of equations we discussed above it is necessary to add constitutive equations which specify the behaviour of each species. In other words, we have to add equations involving the velocities v (or the fluxes q) and the volumetric content ν. The most relevant physical and chemical processes occur within the frozen fringe. Mass and heat transfer processes are strictly related to the interfacial effects between the tree constituents (water, ice and soil grains). More than one hypotheses have been formulated. The lack of a common understanding, due to the experimental difficulties of measuring the quantities involved in the process, motivates the large number of different models in the literature. For a detailed review of models in literaure, we refer to [18]. It is generally accepted that water pressure in the porous space is not uniform and it is related to the movement of water by means of

q_{w} = −K_{1}(p_{w}, T)\frac{∂p_{w}(z, t)}{∂z} − K_{2}(p_{w}, T)\frac{∂T (z, t)}{∂z} (2.26)

where K_{1} and K_{2} are empirical functions. In the literature, two different position are stated: (i) K_{2} ≡ 0 which means essentially a generalized Darcy’s law ([14], (ii) K_{2} > 0 ([15], [12]). A second question concerns with the possible dependence of p_{w} on the temperature T (equilibrium models, [14]). On the contrary, in the model we will discuss T and p_{w} are independent and the Clausius-Clapeyron equation, which relates temperature with pressure, holds only at the phase equilibrium (non-equilibrium models, [10], [5]). A third point consists in discussing the interfacial effects among the three components. Essentially, if water-ice effects are assumed to be dominant, one gets the empirical law p_{i}−p_{w} = ϕ(ν), (capillary models, [7], [9]), if on the contrary water-soil effects prevail, one gets the constitutive equation ν = ν(T) ([22], [11]). The latter one is the point of view we will follow. A further very delicate point is to formulate the mechanism of ice segregation. Main theories are Terzaghi’s stress partition idea ([21]) (it is assumed that P = σ_{n} + σ_{e}, where P is the overburden pressure P, σ_{e} is the effective stress and σ_{n} = χp_{w} + (1 − χ)P_{i} is the neutral stress, the break of the soil is ascribed to the vanishing of σ_{e}, [21]), the existence of a segregation temperature T_{s} (the permeability K_{1} is so small at T = T_{s} that water tends to accumulate somewhere, [8]). In other models ([13]), the stress balance on the top of the frozen fringe is simply p_{w} = P and lens starts to form simply when the speed of the upper boundary of the frozen fringe vanishes. The last question to be faced in formulating a frost heave model concerns with the movement of ice. We already discussed some aspects in the previous section (regelation phenomena). In the model we will discuss ice is at rest with respect to the porous matrix.

2.5 A non-equilibrium model: the constitutive equations

The model we investigated from the mathematical point of view is based on the experimental observations and data performed in ([10], [11]). The statements of the model, according to the questions we introduced in paragraph 2.5, are the following:

• T and p_{w} are independent in the frozen fringe. The value of p_{w} at the top of the frozen fringe is a given positive constant σ which is essentially the overburden pressure.

• The empirical functions K_{1} and K_{2} appearing in the water flux law depend only on T, they are strictly positive, increasing functions and vanishing for very low temperatures. Moreover, K_{1}(T) ≡ K_{0} positive constant if T ≥ 0.

• Interfacial effects are described by the function ν = ν(T), 0 ≤ ν ≤ 1, ∂ν(T)/∂T < 0, ν(0) = 0. No capillarity effects are present.

• Pore ice in the frozen fringe is at rest with respect to the porous matrix (no regelation) and the frozen soil is simply shifted up during thre growth of the lens, without any appreciable deformation.

• The quantities k_{u}, k_{f}, k_{i} (thermal conductivities), ρ_{w}, ρ_{i}, ρ_{s} (specific densities) and L (latent heat per unit volume) are constant. On the other hand, the thermal conductivity in the frozen fringe k_{ff} depends on the temperature T. Empirically, k_{ff} is determined as a function of ν (hence of T). Moreover, k_{ff}(0) = k_{u}.

2.6 The mathematical problem

The conservation and constitutive equations of the model are summarized below (we omit details, see also the statements of the model in section 3.1.):

C(0, T)\frac{∂T (z, t)}{∂t} − k_{u} \frac{∂}{∂z} (\frac{∂T (z, t)}{∂z}) + ρ_{w}c_{w}q_{w}(t)\frac{∂T (z, t)}{∂z} = 0 0 < z < z_{F} , t > 0 (2.27)

q_{w}(t) = −K_{0} \frac{∂p_{w}}{∂z} (z, t), 0 < z < z_{F} , t > 0 (2.28)

ε_{0}(ρ_{i} − ρ_{w})\frac{∂}{∂t}ν(T) + \frac{∂}{∂z} (ρ_{w}q_{w}(z, t)) = 0, z_{F} < z < z_{S}, t > 0 (2.29)

C(ν, T)\frac{∂T (z, t)}{∂t} + \frac{∂}{∂z} (−k_{ff}(T)\frac{∂T (z, t)}{∂z}) + ρ_{w}c_{w}q_{w}(z, t)\frac{∂T (z, t)}{∂z} + (2.30)

+ρ_{i} (ε_{0} \frac{∂}{∂t}ν(T)) ((c_{i} − c_{w})T − L) = 0, z_{F} < z < z_{S}, t > 0

q_{w}(z, t) = −K_{1}(T)\frac{∂p_{w}}{∂z} (z, t) − K_{2}(T)\frac{∂T}{∂z} (z, t) (2.31)

ρ_{i}c_{i} \frac{∂T}{∂t} − k_{i} \frac{∂}{∂z} (\frac{∂T}{∂z} ) + ρ_{i}c_{i}ż_{T} \frac{∂T}{∂z} = 0, z_{S} < z < z_{T} , ε = 1 (ice lens), t > 0 (2.32)

C(1, T)\frac{∂T}{∂t} + \frac{∂}{∂z} − k_{f} (\frac{∂T}{∂z} ) + C(1, T)ż_{T} \frac{∂T}{∂z} = 0, z_{S} < z < z_{T} , ε < 1 (grains + ice), t > 0 (2.33)

\frac{∂ε}{∂t} + ż_{T} \frac{∂ε}{∂z} = 0, z_{S} < z < z_{T} , t > 0 (2.34)

together with the initial and boundary conditions

T (0, t) = h > 0 (2.35)

[[\frac{∂T}{∂z} (z_{F} , t)]] = 0 (2.36)

[[q_{w}(z_{F} , t)]] = 0 (2.37) ρ_{i}ż_{T} = ρ_{w}q_{w}(z_{S}, t) + ε_{0}(1 − ν(T (z_{S}, t)))(ρ_{i} − ρ_{w})ż_{S} (2.38)

ρ_{w} (T (z_{S}, t)(c_{i} − c_{w}) + L) q_{w}(z_{S}, t) = [[−k \frac{∂T}{∂z} (z_{S}, t)]] (2.39)

p_{w}(z_{S}(t), t) = σ (2.40) T (z_{T} , t) = g < 0 (2.41)

Moreover, (2.23) must hold at any interface and (2.22) at any internal interface in the frozen part of the soil. The boundary conditions (2.35) and (2.41) can also be replaced by the assignment of the thermal flux:

\frac{∂T}{∂z} = −α_{0}, \frac{∂T}{∂z} = −α_{1}. (2.42)

with α_{1} and α_{2} non-negative. For the sake of simplicity, in the energy conservation equations we neglect the convective contribution. The unknown quantities of problem (2.27)-(2.41) are the temperature T, the water flux q_{w}, the water pressure p_{w}, the free boundaries z_{F}, z_{S} and z_{T}. As a matter of facts, also the interfaces between the regions L and F are unknown. However, a simpler but reasonable way to treat the problem is to consider just the forming (or just formed lens), that is the lowest one. The base of that lens is fixed (owing to (2.20)) and we call z = z_{L} the upper boundary of the lens. The rest of the frozen soil can be considered as homogeneous from the thermal point of view and equation (2.33) can be assumed to hold in the region z_{S} < z < z_{T}. Such a simplification yields to consider only three regions (unfrozen soil, frozen fringe, frozen soil), as indicated in next picture.

[figure: z = 0 base of soil; z = z_{F} 0^{o}C- isotherm; z = z_{S} freezing front; z = z_{T} top of soil; unfrozen soil; frozen fringe; frozen soil; z]

fig. 3: the different parts of the soil in the one-dimensional case

Even in this case, however, problem (2.27)-(2.41) exhibits non trivial difficulties, mainly for the fact that there is a switch from one process to the other whenever ż_{s} vanishes. In next paragraphs we will focus on two simplified situations, namely the case of lens formation and the quasi-stady approach.

2.7 A free boundary problem for the process of lens formation

We are going to discuss system (2.27)-(2.41) in the case ż_{S} ≡ 0. In a easier but realistic situation, we may assume that the temperature in the unfrozen soil 0 < z < z_{F} and in completely frozen soil z_{S} < z < z_{T} is linear. This fact is consistent with experimental data ([10]) and allows us to concentrate the problem in the frozen fringe z_{F} < z < z_{S}, where the most relevant physical and chemical processes occur. Under such an assumption, the problem of lens formation in z_{F} < z < b ≡ z_{S}(0) reduces to

(C(ν(T)) − Lεν′(T)) \frac{∂T}{∂t} (z, t) = \frac{∂}{∂z} (k_{ff}(T)\frac{∂T}{∂z} (z, t)) (2.43)

ϱ_{w} \frac{∂q}{∂z} (z, t) + (ϱ_{i} − ϱ_{w})εν′(T)\frac{∂T}{∂t} (z, t) = 0 (2.44)

q_{w}(z, t) = −K_{1}(T)\frac{∂p_{w}}{∂z} (z, t) − K_{2}(T)\frac{∂T}{∂z} (z, t) (2.45)

together with the initial and boudary conditions (we consider case (2.42))

T (z, 0) = T_{0}(z) z_{F} (0) = a ≤ z ≤ b (2.46) T (z_{F} (t), t) = 0 t > 0 (2.47)

\frac{∂T}{∂z} (z_{F} (t), t) = −α_{0} t > 0 (2.48)

q_{w}(z_{F} (t), t)z_{F} (t) = −K_{0}p_{w}(z_{F} (t), t) t > 0 (2.49) p_{w}(b, t) = σ t > 0 (2.50)

ϱ_{w}Lq_{w}(b, t) = k_{f}α_{1} + k_{ff}(T (b, t))\frac{∂}{∂z} T (b, t) t > 0 (2.51)

Problem (2.43)-(2.51) is a strictly parabolic problem for the temperature T in a free domain. We indeed remark that Φ(T) = C(ν(T)) − Lεν′(T) is strictly positive, owing to the assumptions we made above. The well posedness of the problem is shown in [6]. We summarize here the main steps of the mathematical analysis.

First of all, we solve the problem isolating the equations for T. Eliminating ν by comparing (2.43) and (2.44) we get

\frac{∂q}{∂z} (z, t) = \frac{ϱ_{i} − ϱ_{w}}{Lϱ_{i}ϱ_{w}} \frac{∂}{∂z} (k_{ff}(T)\frac{∂T}{∂z} (z, t))

Then, by integrating in (z, b), using (2.51) we find

q(z, t) = q(b, t) (1 + ψ(T (b, t))ϱ_{w}L) − ∫_{z}^{b} ψ′(T (ξ, t))k_{ff}(T (ξ, t))\frac{∂}{∂z} T^{2}(ξ, t)dξ −

− ψ(T (b, t))k_{f}α_{1} − ψ(T)k_{ff}(T)\frac{∂}{∂z} T (z, t)

where ψ(T) = (ϱ_{i} − ϱ_{w})εν′(T)/ϱ_{w}Φ(T) At this point, we integrate (2.45) in z_{F} < z < b:

σ + \frac{q(z_{F} , t)}{K_{0}}z_{F} = − ∫_{z_{F}}^{b} \frac{q(ξ, t)}{K_{1}(T (ξ, t))}dξ − ∫_{0}^{T (b,t)}\frac{K_{2}(η)}{K_{1}(η)}dη

and we sobstitute the expression for q(z, t) achieved above:

(2.52) q(b, t) = \frac{1}{(1 + ψ(T_{b})ϱ_{w}L) (\frac{z_{F}}{K_{0}} + ∫_{z_{F}}^{b} \frac{dz}{K_{1}(T)})} (− ∫_{0}^{T_{b}} \frac{K_{2}(η) − ψ(η)k_{ff}(η)}{K_{1}(η)}dη − σ+

+ψ(T_{b})k_{f}α_{1} (\frac{z_{F}}{K_{0}} + ∫_{z_{F}}^{b} \frac{ψ(T_{b})k_{f}α_{1}}{K_{1}(T)} dz) − \frac{z_{F}}{K_{0}} (ψ(0)k_{u}α_{0})+

+\frac{z_{F}}{K_{0}} ∫_{z_{f}}^{b} ψ′(T)k_{ff}(T) (\frac{∂T}{∂z} )^{2} dz + ∫_{z_{F}}^{b} \frac{1}{K_{1}(T)} ∫_{z}^{b} ψ′(T)k_{ff}(T) (\frac{∂T}{∂z} )^{2} dξdz)

where T_{b} = T (b, t). Hence, we find out that the temperature T must verify the following problem T in the frozen fringe z_{F} < z < b:

Φ(T)\frac{∂T}{∂t} (z, t) = \frac{∂}{∂z} (k_{ff}(T)\frac{∂T}{∂z} (z, t))

T (z, 0) = T_{0}(z) z_{F} (0) = a T (z_{F} (t), t) = 0

\frac{∂T}{∂z} (z_{F} (t), t) = −α_{0}

ϱ_{w}Lq(b, t) = k_{f}α_{1} + k_{ff}(T (b, t))\frac{∂}{∂z} T (b, t) (2.53)

In (2.53) q(b, t) is the functional F(T, z_{F}) defined in (2.52). The mathematical problem is easier if we assume that ϱ_{i} ∼ ϱ_{w}. Actually, in this case the water flux depends only on t and the functional (2.52) reduces to

q(t) = \frac{− ∫_{0}^{T_{b}} \frac{K_{2}(η)}{K_{1}(η)}dη − σ}{\frac{z_{F}}{K_{0}} + ∫_{z_{F}}^{b} \frac{dz}{K_{1}(T)}} .

By means of the Kirchoff transformation u(x, t) = ∫_{0}^{T (b−x,t)}k_{ff}(η)dη. one finds that problem T is trans-

formed into the following problem U:

\frac{∂}{∂t}u(x, t) = A(u)\frac{∂^{2}}{∂x^{2}} u(x, t) 0 < x < s(t), t > 0

u(x, 0) = u_{0}(x) = ∫_{0}^{T_{0}(b−x)}k_{ff}(η)dη 0 ≤ x ≤ b − a

s(0) = b − a > 0 u(s(t), t) = 0 t > 0

\frac{∂}{∂x}u(s(t), t) = k_{u}α_{0} t > 0

\frac{∂}{∂x}u(0, t) = k_{f}α_{1} − ϱ_{w}Lq(t) t > 0

where we set A(u) = k_{ff}(T (u))/Φ(T (u)). Problem U is a non-standard free boundary problem, with Cauchy data on the unknown boundary s and a functional-type dependence of the data on the fixed boundary x = 0. Local existence for problem U can be proved in the following way. First of all, we write the problem in the rectangular domain Q_{τ} = (0, 1)×(0, τ), with τ fixed positive value, via the transformation y = x/s(t) and call v(y, t) = u(s(t)y, t). The problem is now

\frac{∂}{∂t}v(y, t) = \frac{Â(v)}{s^{2}} \frac{∂^{2}}{∂y^{2}} v(y, t) + \frac{ṡ}{s}y \frac{∂}{∂y} v(y, t) (y, t) ∈ Q_{τ}

v(y, 0) = v_{0}(y) = u_{0}(ay) 0 ≤ y ≤ 1 v(1, t) = 0 t > 0

\frac{∂}{∂y} v(0, t) = s(t)(k_{f}α_{1} − ϱ_{w}Lq(t)) t > 0

\frac{∂}{∂y} v(1, t) = k_{u}α_{0}s(t) t > 0.

Then, consider the set of functions

S = {(v, s) |v(y, t) ∈ C^{1,1/2}([0, 1]) × [0, τ],
s ∈ C^{1}([0, τ]), v(y, 0) = v_{0}(y) for y ∈ [0, 1], s(0) = b − a}.

For a pair (v, s) ∈ S we solve the following linear parabolic problem L of Neumann type for the unknown function ṽ:

\frac{∂}{∂t}ṽ(y, t) = \frac{Â(v)}{s^{2}} \frac{∂^{2}}{∂y^{2}} ṽ(y, t) + \frac{ṡ}{s}y \frac{∂}{∂y} ṽ(y, t) 0 < y < 1, t > 0

ṽ(y, 0) = v_{0}(y) 0 ≤ y ≤ 1

\frac{∂}{∂y} ṽ(1, t) = k_{u}α_{0}s(t) t > 0

\frac{∂}{∂y} ṽ(0, t) = s(t)(k_{f}α_{1} − ϱ_{w}Lq(t)) t > 0.

The function q(t) depends on v and s in the following way:

q(t) = \frac{− ∫_{0}^{h(v(0,t))} \frac{K_{2}(η)}{K_{1}(η)}dη − σ}{\frac{b − s}{K_{0}} + s ∫_{0}^{1} \frac{dz}{K̂_{1}(v(y, t))}}

where h(w) is the inverse function of the Kirkhoff transformation (i. e. h(v(y, t)) = T (b − s(t)y, t)) and K̂_{1}(v) = K_{1}(h(v)). Finally, we define a new guess (ṽ, s̃) such that ṽ is the solution of problem L and

˙s̃(t) = \frac{Â(0)}{k_{u}α_{0}s^{2}(t)} \frac{∂^{2}v}{∂y^{2}} (1, t)

In [6] we proved that the map (v, s) → (ṽ, s̃) is a contraction. In the same paper the general case ρ_{w} ≠ ρ_{i} is discussed.

2.7.1 Frost penetration

The case ż_{S} < 0 presents a nontrivial complication: actually, one further condition must be added in the model. The physical understanding of the problem has not been solved yet. The following two possibilities can be considered:

(1) the thickness of the frozen fringe remains constant, that is ż_{s} = ż_{F}. This assumption, giving travelling waves solutions,

(2) the temperature T (z_{S}(t), t) corresponds to a known critical value T_{C}, at which all water freezes: ν ≡ 1 for T ≤ T_{C}.

2.8 A quasi-steady approach

We are going now to simplify the mathematical problem in order to obtain more qualitative information of the freezing process. Actually, in practical situations it is important to predict which process will take place (lens formation or frost penetration) and which is the heave rate. We add the hypothesis that the variation with respect to time of the involved quantities is negligible in comparison to the variation with respect to the spatial coordinate. In that case the water flux q_{w} depends only on time, as we see from (2.29). The set of equations corresponding to the model when the boundary temperatures h(t) = T (0, t), g(t) = T (z_{T} (t), t), h(t) > 0, g(t) < 0 are assigned, is

∫_{0}^{T_{S}(t)}k_{ff}(η)dη = \frac{z_{F} (t) − z_{S}(t)}{z_{F} (t)}k_{u}h(t) (2.54)

ϱ_{i}ż_{T} (t) = ϱ_{w}q_{w}(t) + ε(1 − ν_{S})(ϱ_{i} − ϱ_{w})ż_{S}(t) (2.55)

Lϱ_{w}q_{w}(t) = (1 − ν_{S})εϱ_{w}Lż_{S}(t) − k_{f} \frac{g(t) − T_{S}(t)}{z_{T} (t) − z_{S}(t)} − \frac{k_{u}h(t)}{z_{F} (t)} (2.56)

σ + ∫_{0}^{T_{S}(t)}\frac{K_{2}(η) − \frac{q_{w}(t)z_{F} (t)}{k_{u}h(t)}k_{ff}(η)}{K_{1}(η)}dη = −\frac{q_{w}(t)}{K_{0}}z_{F} (t) (2.57)

q_{w}(t) = −K_{1}(T_{S}(t))\frac{∂p_{w}(z_{S}(t), t)}{∂z} + K_{2}(T_{S}(t))\frac{k_{u}h(t)}{z_{F} (t)k_{ff}(T_{S}(t))} (2.58)

z_{S}(0) = b, z_{T} (0) = H > b (2.59)

\frac{∂p_{w}(z_{S}(t), t)}{∂z}ż_{S}(t) = 0 (2.60)

ż_{S}(t) ≤ 0 (2.61)

\frac{∂p_{w}(z_{S}, t)}{∂z} ≥ 0 (2.62)

0 ≤ z_{F} (t) ≤ z_{S}(t) ≤ z_{T} (t) (2.63) q_{w}(t) ≥ 0 (2.64)

The unknown quantities of (2.54)-(2.64) are the water flux q_{w}(t), the boundaries z_{F} (t), z_{S}(t) and z_{T} (t), the freezing temperature T_{S}(t) = T (z_{S}(t), t) and the water pressure gradient ∂p_{w}(z_{S}(t), t)/∂z. Equation (2.54) comes from integrating the energy balance (2.30) in the frozen fringe under the assumption of quasi-steady process. Equations (2.55) and (2.56) (where we defined ν_{S} = ν(T_{S})) are respectively (2.38) and (2.39) (with c_{w} ∼ c_{i}), while (2.57) and (2.58) come from rearranging (2.30), (2.31), (2.37), (2.40).

Equations (2.60)-(2.62) are the most peculiar of the model. In particular, (2.60) is a consequence of the quasi-steady assumption and of (2.40), (2.62) means that no ice lens is forming within the frozen fringe (see [18] for details). Whenever ż_{S}(t) vanishes, a lens is growing at the height z = z_{S}, whenever ż_{S}(t) < 0, the freezing front is advancing (frost penetration); in the latter case, the water pressure gradient at the freezing front must vanish. Condition (2.64) guarantees that the water is sucked up towards the region where the change of phase occurs. We refer to [18], [19] and [20] for the mathematical analysis of the problem. We only summarize the mathematical results, which allow us to predict the kind of freezing process, according to the given initial data. We first define the following functions, which are known once the properties of the soil and the boundary temperatures are given:

φ_{1}(s, t) = \frac{ϱ_{i}b}{ϱ_{w}K_{0}} \frac{k_{u}h(t) − K_{0} ∫_{0}^{s} \frac{k_{ff}(η)}{K_{1}(η)}dη}{k_{u}h(t) − ∫_{0}^{s} k_{ff}(η)dη} (2.65)

φ_{2}(s) = ∫_{0}^{s} \frac{K_{2}(η)}{K_{1}(η)}dη + σ (2.66)

φ_{3}(s, t) = \frac{1}{Lϱ_{i}b} (∫_{0}^{s} k_{ff}(η)dη − k_{u}h(t)) . (2.67)

In [] it is shown that if φ_{2}(s) > 0 for any s ≤ 0, then system (2.54)-(2.64) has no solutions. On the other hand, if φ_{2}(T_{σ}) = 0 for some T_{σ} < 0, then there exists a unique temperature T_{p}(h(t)) satisfying

φ_{2}(T_{p}(h(t))) = ϱ_{w}L\frac{K_{2}(T_{p}(h(t)))}{k_{ff}(T_{p}(h(t)))}φ_{1}(T_{p}(h(t)), t)φ_{3}(T_{p}(h(t)), t).

Thus, admitting that the equilibrium temperature T_{σ} exists, we have the following cases:

(i) if k_{f}g(0) > k_{f}T_{σ} + (H − b)Lρ_{i}φ_{3}(T_{σ}, 0), then a melting process will occur (q_{w} < 0),

(ii) if

k_{f}T_{p}(h(0)) + (H − b)Lρ_{i}φ_{3}(T_{p}(h(0)), 0) (1 + ρ_{w}L\frac{K_{2}(T_{p}(h(0)))}{k_{ff}(T_{p}(h(0)))}) ≤ k_{f}g(0) ≤

≤ k_{f}T_{σ} + (H − b)Lρ_{i}φ_{3}(T_{σ}, 0)

then a lens formation process will take place,

(iii) if k_{f}g(0) < k_{f}T_{p} + (H − b)Lρ_{i}φ_{3}(T_{p}, 0) (1 + ρ_{w}L\frac{K_{2}(T_{p})}{k_{ff}(T_{p})}) , then a process of frost penetration will

take place.

In [19] several qualitative properties of the solution are investigated. In particular, the attention is focussed on the possibility of switching from un process to the other. We end this section with the foolowing remark. A more sofisticated (and realistic) way to switch from one process to the other consists in the possibility that the water pressure attains the critical value σ in some point of the frozen fringe. Actually, the new lens will start growing within the transition region. In the quasi-steady assumption the water pressure profile is non decreasing in the frozen fringe and attains its maximum value at z = z_{S} (see picture ). Therefore, in order to encompass such a possibility the evolutive model presented in section is appropriate. This aspect of the problem is the subject of present mathematical investigations by us.

3 A non-standard problem of filtration in a porous medium

We are going now to present a second application of the general conservation equations we discussed in section . The problem consists in a filtration of a multi-species fluid through a composite porous medium. The main feature of the process is the strong interaction between the liquid and the components of the porous medium, caused by mechanical and chemical effects of the flow on the porous matrix and vice versa. The context of such a research program, promoted by illicaffè of Trieste, Italy, was the understanding of the complex phenomena which take place during the filtration of water through a layer of ground coffee when the rxpresso machine is operating.

3.1 The physical problem

The physical situation can be outlined in the following way. The water filters through the medium made of a compact coffee bed by applying a pressure on the top of the layer (typically: 9 bars). The temperature of water is quite high in a standard coffee machine (about 90 ^{o}C). At the opposite surface, the liquid outflows from the medium. The process is sketched in figure 3.1.

[figure: section; incoming water flow; removal process; compact layer; outcoming flow; liquid; particles bound to the porous matrix; porous matrix; removed particles transported by the flow; compact layer; flow; (given pressure); (atmospheric pressure)]

fig. 3.1: a simplified scheme of the ground coffee layer: the cylinder is the coffee cake, and the lower rectangle is a transversal section of it.

The physical processes developing during the experiment are both of mechanical type (filtration) and of chemical type (extraction of substances). The two processes are strongly coupled: the extraction changes the properties of the liquid (density, viscosity) and of the solid (porosity, hydraulic permeability); on the other hand, the filtration velocity affects the extraction kinetics. If we let down the temperature to few degrees (4^{o} C), (as in the experiments performed in [16])the extraction process is drastically reduced and the mechanical effects are in evidence. However, even in this case the process is complex the same. Actually, two experimental observations point out some non standard aspects, which cannot be explained just modifying Darcy’s law of filtration in a porus medium in one of the ways suggested by the literature ([1]):

(i) flow is not constant for constant injection pressure, but displays an initial transient in which a sharp maximum is reached and then decreases in time towards an apparently asymptotic value (dependent on temperature and pressure),

(ii) even the asymptotic value is not proportional to the injection pressure (as we would expect from Darcy’s law), but it increases until a certain value, and then remains constant or even decreases for larger values of pressure, thus exhibiting a non-monotone behaviour.

Figure 3.2 shows this non standard behaviour of the liquid discharge.

[figure: 2; 4; 6; 8; 10; 12; 5; 10; 15; 20; 25; time (s); discharge (ml/s); 3 bar; 5 bar; 7 bar]

fig. 3.2: pressure dependence on flow during percolation.

A second series of tests puts in evidence further surprising effects. Such tests are performed by percolating water through an over-turnable extraction chamber. The direct flow has the described behaviour and once the asymptotic discharge has been reached the cake seems to behave like an ordinary porous medium. However, as the percolation chamber is rotated, the flow surprisingly goes through the exponential decay once more (see fig. 3.3).

[figure: 24; 22; 20; 18; 16; 14; 12; 10; 8; 6; 4; 2; 2; 4; 6; 8; 10; 12; 14; 16; discharge (ml/s); time (s)]

fig. 3.3: direct/inverse discharge curve.

This effect suggests that the ground coffee particles may exhibit a progressive rearrangement under the action of the flow. Actually, electronic microscopy shows that ground coffee contains a component of fine particles, which can be removed and transported by the flow: after the flow inversion the particles which accumulated first in a lower section of the cake (causing an increase of the hydraulic resistance) now counter migrate, with an initial increase of the hydraulic conductivity in the opposite direction and subsequently leading the system to the previous steady state. In such a setting, the fines are expected to migrate and eventually to accumulate in some part of the bed, typically at the outflow surface.

3.2 The mathematical model

We refer to [3] and [4] for the detailed description of the model. Our conceptual model consists in considering a rigid porous skeleton, where n families of substances are bounded on it. Indexes k = 1, …, n_{1} < n correspond to fine particles species, while k = n_{1} + 1, …, n correspond to dissolvable species. The percolation of water gives rise both to removal of fine particles and to the dissolution of other substances (typically as solutes). Referring to the notations of section 1, we find it convenient to consider N = 2n + 2 species forming the medium, namely

(i) k = 1, …, n_{1} are n_{1} families of fined solid particles, bound to the porous matrix, with concentrations b_{k}, k = 1, …, n_{1} (corresponding to ϱ_{k} of section 1),

(ii) k = n_{1} + 1, …, n are n − n_{1} families of substances bound to the porous matrix which may dissolve in the liquid, with concentrations b_{k}, k = n_{1} + 1, …, n,

(iii) k = n + 1, …, n + n_{1} are the fine particles of point (i) when transported by the flow, with concetrations m_{k}, k = 1, …, n_{1},

(iv) k = n + n_{1} + 1, …, 2n are the substances of point (ii) when dissolved in the multi=species flow as solutes or droplets, with concetrations m_{k}, k = n_{1} + 1, …, n,

(v) k = 2n + 1 is water, with concentration m_{w},

(vi) k = 2n + 2 corresponds to the rigid porous skeleton, which is made by large particles and which is expected to be undeformable.

According to the notations of Section 1, we will denote the volumetric fractions of the species as follows (see (1.1)):

μ_{k} = ϑ_{k}, , k = 1, …, nμ_{n+k} = η_{k}, , k = 1, …, nμ_{2n+1} = η_{w}, μ_{2n+2} = ϑ_{0}.

We assume that the medium is saturated and we define the porosity ε as the volume fraction occupied by water and by the moving dissolved substances:

ε = η^{(n_{1})} + η_{w} where η^{(n_{1})} = ∑_{k=n_{1}+1}^{n} η_{k}. (3.1)

The complementary volume 1 − ε is made by the fine particles at rest and in motion, the porous skeleton and the dissolvable substances still bound to the porous matrix:

1 − ε = η_{(n_{1})} + θ_{0} + θ^{(n_{1})} + θ_{(n_{1})}, (3.2)

where η_{(n_{1})} = ∑_{k=1}^{k} η_{k}, θ_{(n_{1})} = ∑_{k=1}^{k} θ_{k}, θ^{(n_{1})} = ∑_{k=k+1}^{n} θ_{k}.

We will consider a one dimensional model. In a space-time plan, we expect a situation like the one drawn in figure 3.4: D_{T} = {(x, t) : 0 < x < s(t), 0 < t < T } is the region where the removal-dissolution process (species k = 1, …, n_{1}) takes place. The complementary region R_{T} = {(x, t) : s(t) < x < 1, 0 < t < T } is the compact layer where fine particles accumulate giving rise to a region with high hydraulic resistance. The lines x = 0 and x = 1 are respectively the inflow and outflow surfaces; the line s(t) separating D_{T} from R_{T} is a free boundary. The dissolving process (species k = n_{1} + 1, …, n) may take place also in R_{T}.

[figure: D_{T}; R_{T}; t = T; x = s(t); 1; x; t; compact layer; “active” region; inflow surface; outflow surface]

figure 3.4: the regions D_{T} and R_{T}

Since we will focus our attention on the mechanical processes, we will neglect the thermal aspect of the problem, so that we assume that the temperature is nearly constant during the filtration process. Only the species m_{k}, k = 1, … n are moving. We call V_{m_{k}} their molecular velocities. We introduce the volumetric compound velocity q (volume of fluid crossing a unit surface of the medium in a unit time) and the average velocity V:

q = η_{w}V_{w} + ∑_{k=n_{1}+1}^{n} η_{k}V_{m_{k}}, V = q/ε.

We assume that the fine particles (species 1, …, k) are transported convectively by the flux, while the dissolvable species diffuse in the fluid with a diffusive molecular flux J_{i} = −D_{i}∇(m_{i}/ε)(Fick’s law). Hence, we have:

V_{m_{i}} = V, i = 1, …, k J_{i} = m_{i}(V_{m_{i}} − V)/ε = ∇(m_{i}/ε), i = k + 1, …, n (3.3)

As a consequence, we make use of equations (1.7) for the moving species k = n+1, …, 2n, by considering that I_{k+n} = −∂b_{k}/∂t:

\frac{∂m_{k}}{∂t} + \frac{∂}{∂x} (m_{k} \frac{q}{ε}) = −\frac{∂b_{k}}{∂t} , k = 1, …, n_{1}, (x, t) ∈ D_{T} (3.4)

\frac{∂m_{k}}{∂t} + \frac{∂}{∂x} (−D_{k}ε\frac{∂}{∂x} \frac{m_{k}}{ε} ) + \frac{∂}{∂x} (m_{k} \frac{q}{ε}) = −\frac{∂b_{k}}{∂t} , (3.5)

k = n_{1} + 1, …, n (x, t) ∈ D_{T} ∪ R_{T}

The kinetics of fine particles and dissolvable species are regulated by means of the following empirical laws:

\frac{∂b_{k}}{∂t} = −F_{k}(q, b)G_{k}[b_{i} − β_{k}(q, b)]^{+}, k = 1, …, n_{1}, (x, t) ∈ D_{T} (3.6)

\frac{∂b_{k}}{∂t} = −H_{k}(q, b), k = n_{1} + 1, …, n, (x, t) ∈ D_{T} ∪ R_{T} (3.7)

In (3.6), (3.7) F_{i}, G_{i}, H_{i} are positive functions. The threshold concentrations β_{k} play a crucial role in the model. They decrease if the flux increases, since a more intense flux operates in a larger population of bound particles. The dependence of the functions F_{k}, G_{k}, H_{k} on the whole vector b is explained by possible interactions between different species (for instance, a coating of fat substance can interfere with the dissolution of the underlying components). The following equations can be obtained by combining the conservation equations (3.23),(3.24) with (3.1) and (3.2):

\frac{∂ε}{∂t} + \frac{∂q}{∂x} = −\frac{∂}{∂t}θ^{(n_{1})}, (x, t) ∈ D_{T} ∪ R_{T} (3.8)

\frac{∂q}{∂x} + \frac{∂}{∂x} (η_{(n_{1})} \frac{q}{ε}) = 0, (x, t) ∈ D_{T} ∪ R_{T} (3.9)

They are useful from the mathematical point of view and physically expressive. Both of them express saturation and incompressibility. The first one says that the variation in porosity is due to the dissolution of substances from the porous matrix and to the incoming or outcoming flux. We remark that the quantities corresponding to the fine particles are not present in the evolutive equation for the porosity, since they are computed in the complementary volume both if they are bound to the porous medium and if they are in motion. The second equation is the conservation of the total flux of all the moving components: actually, the second term can be seen as the volumetric flux of the fine particles. Finally, we assume that the liquid flux is given in terms of the pressure gradient:

q = −K(q, b, m, ε)\frac{∂p}{∂x}, (x, t) ∈ D_{T} ∪ R_{T} (3.10)

The hydraulic conductivity K is constant in the compact layer R_{T}. The initial conditions of the problem prescribe the concentrations and the porosity for t = 0:

m_{k}(x, 0) = m_{k,0}(x), k = 1, …, n, x ∈ [0, 1] (3.11) ε(x, 0) = ε_{0}(x), x ∈ [0, 1] (3.12) b_{k}(x, 0) = b_{i,0}(x), k = 1, …, n, x ∈ [0, 1] (3.13)

The boundary conditions are one of the most delicate point in the model. At x = 0 (inflow surface), we assume that there is no immission of fine particles and that the velocity of the diffusing species is zero, that is (see (1.21))

m_{k}(0, t) = 0, k = 1, …, n_{1}, 0 ≤ t ≤ T (3.14)

D_{k}ε(0, t)\frac{∂m_{k}}{∂x} (0, t) = \frac{q(0, t)m_{k}(0, t)}{ε(0, t)}, k = n_{1} + 1, …, n, 0 < t < T (3.15)

At the same boundary we prescribe the inflow pressure:

p(0, t) = p_{0}(t) > 0 (3.16)

At the free boundary s(t) we assume, as it is natural, that the liquid pressure is continuous:

[[p]] = 0, x = s(t) (3.17)

Moreover, assuming that the concentrations of the dissolvable species relative to the liquid (m_{i}/ε, i = k + 1 …, n) is continuous at x = s(t), that is

[[\frac{m_{k}}{ε} ]] = 0 k = n_{1} + 1, …, n (3.18)

we get from (1.21)

[[D_{k}ε\frac{∂}{∂x} \frac{m_{k}}{ε} ]] = 0, x = s(t), k = n_{1} + 1, …, n. (3.19)

As to the fine particles, we assume that they are completely at rest in the compact layer. We call M_{k} the concentration of the fine particles in that region, k = 1, …, n_{1}. The structure of the compact layer, which corresponds to a sequence M_{1}, …, M_{n_{1}}, depends on the history of the process. To model the formation of the compact layer means to prescribe a constraint for the concentrations M_{k}, that is to define a packing configuration. A simple but reasonable assumption is to assume that the medium is compact when the incoming particles occupy the maximum volume at their disposal, which is the quantity Θ. By using (1.21) for the fine particles, we get the free boundary condition

(Θ − (η_{(n_{1})} + θ_{(n_{1})})) ṡ = −\frac{q}{ε} ∑_{k=1}^{n_{1}} α_{k}η_{k}, x = s(t), 0 ≤ t ≤ T. (3.20)

Finally, we assume that the outflow pressure is the atmosferic one

p(1, t) = 0 (3.21)

and that there is no diffusion at the outflow surface x = 1:

\frac{∂}{∂x} \frac{m_{k}}{ε} (1, t) = 0, k = n_{1} + 1, …, n (3.22)

3.3 The mathematical result

We are going now to summarize the set of equations of the model:

\frac{∂m_{k}}{∂t} + \frac{∂}{∂x} (α_{k}m_{k} \frac{q}{ε}) = −\frac{∂b_{k}}{∂t} , k = 1, …, n_{1}, (x, t) ∈ D_{T} (3.23)

\frac{∂m_{k}}{∂t} + \frac{∂}{∂x} (−D_{k}ε\frac{∂}{∂x} \frac{m_{k}}{ε} ) + \frac{∂}{∂x} (m_{k} \frac{q}{ε}) = −\frac{∂b_{k}}{∂t} , (3.24)

k = n_{1} + 1, …, n (x, t) ∈ D_{T} ∪ R_{T}

q = −K(q, b, m, ε)\frac{∂p}{∂x}, (x, t) ∈ D_{T} ∪ R_{T} (3.25)

\frac{∂b_{k}}{∂t} = −F_{k}(q, b)G_{k}[b_{k} − β_{k}(q, b)]^{+}, k = 1, …, n_{1}, (x, t) ∈ D_{T} (3.26)

\frac{∂b_{k}}{∂t} = −H_{k}(q, b), k = n_{1} + 1, …, n, (x, t) ∈ D_{T} ∪ R_{T} (3.27)

\frac{∂ε}{∂t} + \frac{∂q}{∂x} = −\frac{∂}{∂t}θ^{(n_{1})}, (x, t) ∈ D_{T} ∪ R_{T} (3.28)

\frac{∂q}{∂x} + \frac{∂}{∂x} ∑_{k=1}^{n_{1}} (α_{k}η_{k} \frac{q}{ε}) = 0, (x, t) ∈ D_{T} ∪ R_{T} (3.29)

m_{k}(x, 0) = m_{k,0}(x), k = 1, …, n, x ∈ [0, 1] (3.30) ε(x, 0) = ε_{0}(x), x ∈ [0, 1] (3.31) b_{k}(x, 0) = b_{k,0}(x), k = 1, …, n, x ∈ [0, 1] (3.32) m_{k}(0, t) = 0, k = 1, …, n_{1}, 0 ≤ t ≤ T (3.33)

D_{k}ε(0, t)\frac{∂m_{k}}{∂x} (0, t) = \frac{q(0, t)m_{k}(0, t)}{ε(0, t)}, i = n_{1} + 1, …, n, 0 < t < T (3.34)

[[p]] = [[\frac{m_{k}}{ε} ]] = [[D_{k}ε\frac{∂}{∂x} \frac{m_{k}}{ε} ]] = 0, x = s(t), k = n_{1} + 1, …, n (3.35)

(Θ − (η_{(n_{1})} + θ_{(n_{1})})) ṡ = −\frac{q}{ε} ∑_{k=1}^{n_{1}} α_{k}η_{k}, x = s(t), 0 ≤ t ≤ T (3.36)

\frac{∂}{∂x} \frac{m_{k}}{ε} (1, t) = 0, k = n_{1} + 1, …, n (3.37)

p(0, t) = p_{0}(t) > 0, p(1, t) = 0, 0 ≤ t ≤ T, s(0) = 1 (3.38)

The 2n + 4 unknown quantities are the concentrations b_{k} e m_{k} (or the specific volumes θ_{k} e η_{k}), the flux q, la pressure p, the porosity ε and the interface s. In [3] we proved that system (3.23)-(3.38) has exactly one solution. We are not going to reproduce here all the details of the proof: we will summarize the main steps of the procedure. The assumptions on the given functions are the following:

1. K(b, m, ε), Φ_{k} = F_{k}(q, b)G_{k}(η), β_{k}(q, b), k = 1, …, n_{1}, H_{k}(q, b), k = n_{1} + 1, …, n are non-negative C^{1} functions such that:

0 < K_{m} ≤ K(b, m, ε) ≤ K_{M} ∥Φ∥ < ∞, ∥H∥ < ∞ |S_{k}(q_{1}, b) − S_{k}(q_{2}, b)| ≤ L_{S_{k}}^{q}|q_{1} − q_{2}|

|S_{k}(q, b^{(1)}) − S_{k}(q, b^{(2)})| ≤ L_{S_{k}}^{b} max_{1≤j≤n} |b_{j}^{(1)} − b_{j}^{(2)}|

|\frac{∂S_{k}}{∂q} (q_{1}, b) − \frac{∂S_{i}}{∂q} (q_{2}, b)| ≤ L_{S_{k}_{q}}^{q}|q_{1} − q_{2}|

|\frac{∂S_{k}(q, b^{(1)})}{∂b_{j}} − \frac{∂S_{k}(q, b^{(2)})}{∂b_{j}} | ≤ L_{S_{k}_{b_{j}}}^{b} max_{1≤r≤n} |b_{r}^{(1)} − b_{r}^{(2)}|

where S = (Φ_{1}, …, Φ_{n_{1}}, H_{n_{1}+1}, …, H_{n}), Φ = ∑_{k=1}^{n_{1}} Φ_{k}, H = ∑_{i=n_{1}+1}^{n} H_{k}

and ∥Ψ∥ = sup_{q≥0,0≤b_{k}≤b_{k,0}} |Ψ(q, b)|.

2. ε_{0}(x), m_{k,0}(x), b_{k,0}(x), k = 1, …, n and p_{0}(t) are C^{1}-functions such that

0 < p_{0}^{m} ≤ p_{0}(t) ≤ p_{0}^{M}
m_{k,0}(x) ≢ 0, b_{k,0}(x) ≢ 0, x ∈ [0, 1]
η_{(n_{1}),0}(x) + θ_{(n_{1}),0}(x) < Θ < 1 − (θ_{0}(x) + θ_{0}^{(n_{1})}(x)), x ∈ [0, 1].
1 − (Θ + θ_{0}(x) + θ_{0}^{(n_{1})}(x)) < ε_{0}(x) < 1 − θ_{0}(x), x ∈ [0, 1].

Under such an assumption it is not difficult to find out that 0 < ε_{0}^{m} < ε_{0}(x) < ε_{0}^{M} < 1, x ∈ [0, 1]

where ε_{0}^{m} = min_{x∈[0,1]} (1 − (Θ + θ_{0}(x) + θ_{0}^{(n_{1})}(x))) , ε_{0}^{M} = max_{x∈[0,1]} (1 − θ_{0}(x))

The mathematical procedure is based on an iterative scheme, which consists in fixing a pair (q(x, t), s(t)) in a suitable space (specified below) and following the steps:

• calculate b_{k}(x, t), k = 1, …, n, by (3.26), (3.27) and (3.32);

• solve the O. D. E. (3.28), (3.31) to find ε(x, t);

• calculate m_{k}(x, t), k = 1, …, n_{1} from the hyperbolic problems in R_{T} given by (3.23), (3.30) (k = 1, …, n_{1}), (3.33) and m_{k}(x, t), k = n_{1} + 1, …, n from the parabolic problems of diffraction type in R_{T} ∪ D_{T} given by (3.24), (3.30) (k = n_{1} + 1, …, n), (3.32), (3.34), (3.35), (3.37);

• compute p(x, t) by (3.25), (3.38).

Equations (3.29) and (3.36) are used in order to find a new guess q̃ and s̃. Formally, we denote by B_{T} the rectangle (0, 1) × (0, T) and we define the set

E_{T} (u_{1}, u_{2}, A_{y}, A_{t}, M_{y}, M_{t}, s_{0}, A_{s}, M_{s}) =

{(u(y, t), s(t)) | u ∈ C^{1,1}(B̄_{T}), u(y, 0) = q_{0}(y), 0 ≤ y ≤ 1, 0 < u_{1} ≤ u(y, t) ≤ u_{2},

|\frac{∂u}{∂y}(y, t)| ≤ A_{y}, |\frac{∂u}{∂t}(y, t)| ≤ A_{t}, (y, t) ∈ B̄_{T} ,

|\frac{∂}{∂y} u(y_{1}, t) − \frac{∂}{∂y} u(y_{2}, t)| ≤ M_{y} |y_{1} − y_{2}| , ∀ y_{1}, y_{2} ∈ [0, 1],

|\frac{∂}{∂t}u(y_{1}, t) − \frac{∂}{∂t}u(y_{2}, t)| ≤ M_{t} |y_{1} − y_{2}| , ∀ y_{1}, y_{2} ∈ [0, 1],

s ∈ C^{1}[0, T], s(0) = 1, 0 < s_{0} ≤ s(t) ≤ 1, −A_{s} ≤ ṡ(t) ≤ 0, 0 ≤ t ≤ 1, |ṡ(t_{1}) − ṡ(t_{2})| ≤ M_{s} |t_{1} − t_{2}| ∀ t_{1}, t_{2} ∈ [0, T]} .

For a pair (u, s) ∈ E_{T}, set

q(x, t) = u(x/s(t), t), (x, t) ∈ D̄_{T}

and define the map F(u, s) = (ũ, s̃)

ũ(y, t) = f(t)\frac{1}{1 + l(s(t)y, t)}, (y, t) ∈ B_{T} ,

s̃(t) = 1 − ∫_{0}^{t} \frac{l(s(τ), τ)q(s(τ), τ)}{Θ − (η_{(n_{1})}(s(τ), τ) + θ_{(n_{1})}(s(τ), τ))}dτ, t ∈ [0, T]

where l = η_{(n_{1})}/ε and

f(t) = q(0, t) = \frac{p_{0}(t)}{∫_{0}^{s(t)} \frac{1}{1 + l(ξ, t)} \frac{1}{K(ξ, t)}dξ + \frac{1 − s(t)}{K_{0}}}

Then define q̃(x, t) = ũ(x/s̃(t), t) in the set {(x, t) : x ∈ [0, s̃(t)], t ∈ [0, T]}. If (ū, s̄) is a fixed point of F in E_{T}, then q̄ = ū(x/s̄, t), s̄ together with (b̄, m̄, ε̄, p̄) are solutions of (3.23)- (3.38).

Remark 3.1 : the initial flux q_{0} is not given, but it has to be calculated:

q_{0}(x) = \frac{p_{0}(0)}{∫_{0}^{1} \frac{1}{1 + l_{0}(ξ)} \frac{1}{K(b_{0}(ξ), m_{0}(ξ), ε_{0}(ξ))}dξ} \frac{1}{1 + l_{0}(x)}

with l_{0} = \frac{η_{(n_{1}),0}}{ε_{0}} given function.

We also remark that ε is not necessarily differentiable w. r. t. x. Using (3.29), we get

\frac{∂ε}{∂x} = ε (\frac{1}{q} \frac{∂q}{∂x} + \frac{1}{ε + η_{(n_{1})}} \frac{∂}{∂x}(ε + η_{(n_{1})}))

Moreover, we replace ε + η_{(n_{1})} with −θ′_{0}(x) + ∂(θ^{(n_{1})} + θ_{(n_{1})})/∂x, which is computed by means of

\frac{∂ζ_{k}}{∂t} = −\frac{∂Ŝ_{k}}{∂q} \frac{∂q}{∂x} − ∇_{θ}Ŝ_{k} · ζ, k = 1, …, n

ζ_{k}(x, 0) = ζ_{0}(x) = θ′_{k}(x), k = 1, …, n.

Then, replace q/ε with E(x, t) such that \frac{∂E}{∂x} = \frac{q}{ε} \frac{1}{ε + η_{(n_{1})}} (ζ_{0} + ∑_{k=1}^{n_{1}} ζ_{k}) .

Problem (3.23) (k = 1, …, n_{1}), (3.30), (3.33) for m_{k}, k = 1, …, n_{1} (fine particles) is replaced by

\frac{∂η_{k}}{∂t} + η_{k} \frac{∂E}{∂x} + \frac{q}{ε} \frac{∂η_{k}}{∂x} = − − \frac{∂θ_{k}}{∂t} , k = 1, …, n_{1}, (x, t) ∈ D_{T}

η_{k}(x, 0) = m_{k,0}(x)/ρ_{k}, k = 1, …, n_{1}, x ∈ [0, 1] η_{i}(0, t) = 0, k = 1, …, n_{1}, 0 ≤ t ≤ T.

On the other hand, problem (3.24) for the diffusive species k = n_{1} + 1, …, n is modified by replacing ∂ε/∂x with

Ẽ(x, t) = \frac{ε}{q} (\frac{∂q}{∂x} − ε\frac{∂E}{∂x} )

Parabolic problem w. r. t. c_{k} = m_{k}/ε, k = n_{1} + 1, …, n in the domain D_{T} ∪ R_{T} becomes

ε\frac{∂c_{k}}{∂t} + ε\frac{∂}{∂x} (−D_{k} \frac{∂c_{k}}{∂x} ) + (q − D_{k}Ẽ) \frac{∂c_{k}}{∂x} − c_{k} \frac{∂θ^{(n_{1})}}{∂t} = H_{k}, (x, t) ∈ D_{T} ∪ R_{T} ,

c_{k}(x, 0) = m_{k,0}(x)/ε_{0}(x), x ∈ [0, 1],

\frac{∂c_{k}}{∂x} (0, t) = \frac{1}{ε(0, t)} (\frac{q(0, t)}{D_{k}ε(0, t)} − Ẽ(0, t)) c_{k}(0, t), 0 < t < T,

[[c_{k}]] = 0, [[D_{k}ε\frac{∂c_{k}}{∂x} ]] = 0 x = s(t)

\frac{∂c_{k}}{∂x} (1, t) = 0, 0 < t < T.

We first obtain a local result, which is based on the proof (given in [3]) of the two following lemmas:

Lemma 3.1 Define

ϵ_{1} = ε_{0}^{m} − A_{y}T/s_{0}, ϵ_{2} = ε_{0}^{M} + ε_{0}^{m} + ∥θ_{0}^{(n_{1})}∥,
C_{1} = ∥η_{(n_{1}),0}∥/ε_{0}^{m} + ∥Φ̂∥u_{2}/ϵ_{2}, Φ̂ = ∑_{k=1}^{n_{1}} Φ_{k}(q, ρ_{1}θ_{1}, …, ρ_{n}θ_{n})/ρ_{k}.

Under the assumptions

(i) ϵ_{2}C_{1} + ∥θ_{(n_{1}),0}∥ < Θ

(ii) \frac{∥Φ̂∥}{ϵ_{1}}K_{M}p_{0}^{M} < 1, p_{0}^{M}(1 + C_{1})^{2} |\frac{∂K}{∂q} | < 1

(iii) the Lipschitz constants w. r. t. x of ε_{0}, η_{(n_{1}),0}, η_{0}^{(n_{1})}, ε′_{0}, η′_{(n_{1}),0} and L_{S_{i}}^{q}, L^{b}S_{i}, i = 1, …, n are sufficiently small

then it exists a set of positive values { T, u_{1}, u_{2}, A_{y}, A_{t}, M_{y}, M_{t}, s_{0}, A_{s}, M_{s}} such that if (u, s) ∈ E_{T}, then (ũ, s̃) ∈ E_{T}.

Hence the map F operating on the closed, convex and bounded set E_{T} has the property F(E_{T}) ⊆ E_{T} for a suitable choice of the constants.

Lemma 3.2 Under the same assumptions, together to (iv) L_{S_{i}_{q}}^{q}, L_{S_{i}_{b_{j}}}^{b}, i = 1, …, n sufficiently small,

the map F is continuous in E_{T} and F(E_{T}) is precompact.

By applying the Fixed Point Theorem of Schauder, we get at least a solution of (3.23)-(3.38) for a suitable T > 0. The well posedness of the problem can be extended far any time t > 0.

For this purpose, define β = \frac{∥Φ̂∥K_{M}p_{0}^{M}}{ε_{0}^{m}}, ϵ_{2} = ε_{0}^{M} + ε_{0}^{m} + ∥θ_{0}^{(n_{1})}∥ and choose C̄_{1} s. t.

(\frac{∥η_{(n_{1}),0}∥}{ε_{0}^{m}} + β) (1 − β)^{−1} < C̄_{1} < \frac{∥θ_{(n_{1}),0}∥}{ϵ_{2}}.

Proposition 3.1 Under assumptions

(i)a \frac{ϵ_{2}/ε_{0}^{m}}{1 − β} ∥η_{(n_{1}),0}∥ + ∥θ_{(n_{1}),0}∥ + ϵ_{2} \frac{1}{1 − β} < Θ, β < 1,

(ii)a p_{0}^{M}(1 + C̄_{1})^{2} |\frac{∂K}{∂q} | < 1,

(iii) of lemma and (iv) of lemma , system (3.23)-(3.38) has exactly one solution (m_{1}, …, m_{n}, b_{1}, …, b_{n}, ε, q, p, s) for any time t ≥ 0.

The proof of Proposition is based on the following arguments:

• (i)a, (ii)a are more restrictive of (i) e (ii), hence it exists a local solution in [0, T]

• The solution verifies the following properties:

(1) s(t) ≥ s_{min} = 1 − \frac{m̄_{0}}{Θ} , m̄_{0} = ∫_{0}^{1}(η_{(n_{1}),0} + θ_{(n_{1}),0})dx

(2) the concentration of fine particles is below the critical value: η_{(n_{1})} + θ_{(n_{1})} < Θ, (x, t) ∈ D_{T};

(3) ε_{0}^{m} < ε(x, T) < ε_{0}^{M}, x ∈ [0, s(T)].

• apply the local result for t ∈ [T, 2T].

• By iteration, find a solution global in time.

3.4 Non-constant temperature

Let us call c_{k}, k = 1, …, N + 2, the specific heats of each species. We can assume that c_{k} = c_{k+n}, k = 1, …, n, as it is natural. Energy balance (1.17) writes (assume no change of phase occurs)

∑_{k+1}^{n}(m_{k} + b_{k})c_{k} \frac{∂T}{∂t} + \frac{∂}{∂x} (−k(m, b, ε, q)\frac{∂T}{∂x} ) + (3.39)

+ ∑_{k=1}^{n_{1}} c_{k}m_{k}V \frac{∂T}{∂x} + ∑_{k=n_{1}}^{N} c_{k} (m_{k}V − εD_{k} \frac{∂}{∂x} (\frac{m_{k}}{ε} )) \frac{∂T}{∂x} + c_{w}ρ_{w} (q − ∑_{k=n_{1}+1} n\frac{m_{k}}{ρ_{k}} V_{k}) \frac{∂T}{∂x} = 0.

In (3.39) we denoted with ρ_{w} and c_{w} specific density and specific heat of water, respectively. By virtue of tha mass balance conditions (3.35), the condition on the free boundary (1.22) is simply

[[−k \frac{∂T}{∂x} ]] = 0. (3.40)

The problem, which is more complicated, because the concentrations c_{k} and m_{k} (and possibly all the coefficient of the equations) depend on T, is an open question by the moment.

References

[1] Bear, J., Dynamics of fluids in porous media, Elsevier, 1972

[2] Evans, L. C., Partial Differential Equations, Graduate Studies in Mathematics, vol. 19, American Mathematical Society, 1998

[3] Fasano, A., Talamucci, F., A comprehensive mathematical model for a multi-species flow through ground coffee, submitted to SIAM J. Math. Anal.

[4] Fasano, A., Petracco, M., Talamucci, F., The Espresso Coffee Problem, in Complex Flows in Industrial Processes, A. Fasano Editor, Series Modeling and Simulation in Science, Engineering and Technology, Birkhauser 2000

[5] Førland, K. S., Førland, T., Ratkje S., K., Irreversible thermodynamics theory and applications, John Wiley & Sons, 1988

[6] Gianni, R., Talamucci, F., Well posedness of the lens formation problem in secondary frost heave, to appear

[7] Gilpin, R. R. , A model of the “liquid-like” layer between ice and a substrate with applications to wire regelation and particle migration, J. Colloid Interface Sci. , 68, 235-251, 1979

[8] Konrad, J. M. , Morgensterns, N. R. , The segregation potential of a freezing soil, Water Resour. Res. 29, 370-376, 1988

[9] Miller, R. , D. , Frost heaving in non-colloidal soils, Proc. of the 3rd Int. Conf. on Permafrost, Edmonton, Canada, vol. 1, 707-713, 1978

[10] Nakano, Y., ”On the stable growth of segregated ice in freezing soil under negligible overburden pressure”, Adv. Water Resour. vol.7, 1986, pp. 93-102

[11] Nakano, Y. , Quasi-steady problems in freezing soils: I. analysis on the steady growth of an ice layer, Cold Regions Sci. Technol. , 17, 207-226, 1990

[12] Nakano, Y. , Transport of water through frozen soils, Proc. of the 6th International Symposium on Ground Freezing, Beijing, P. R. of China, Y. Xiang, W. Changsheng Eds. , vol. 1, 65-70, 1991

[13] Nakano, Y. , Primicerio M. , Application of quasi-steady solutions of soil freezing for geotechnical engineering, Proc. of Calculus of variation, applications and computations, Pont-a-Mousson, Pitman Research Notes in Mathematics,147-166, 1995

[14] O’Neill, K. , Miller, R. D. , Exploration of a rigid ice model of frost heave, Water Res. 21, 281-296, 1985

[15] Oliphant, J. L., Tice, A. R., Nakano, Y., Water migration due to a temperature gradient in frozen soil, Proc. of the 4th Permafrost Conference, Fairbanks, Alaska, 951-956, 1983

[16] Petracco, M., Espresso coffee brewing dynamics: development of mathematical and computational models, 15éme Colloque Scient. Internat. sur le Café, Associat. Scientif. Internat. du Café, Paris 1993

[17] Takeda, K., Nakano, Y., Quasi-steady problems in freezing soils, II. Experiments on the steady growth of an ice layer, Cold Regions Science and Technology, 18, 225-247, 1990

[18] Talamucci, F. , Analysis of coupled heat-mass transport in freezing saturated soils, Surveys on Mathematics for Industry, vol. 7 pp. 93-139, Springer-Verlag, 1997

[19] Talamucci, F., Freezing processes in saturated soils, Mathematical Models and Methods in Applied Sciences, Vol. 8, No. 1, pp. 107-138, World Scientific Publishing Company, Singapore, 1998

[20] Talamucci, F., Some recent mathematical results on the problem of soil freezing. In: Knutsson, S. (ed. ), Proceeding of the 7th International Symposium on Ground Freezing and Frost Action in Soils, Lulea, Svezia, 15-17 April 1997, pp. 179-186, Balkema, Rotterdam - 1997 [21] Terzaghi, K. , The shearing resistance of saturated soils and the angle between the planes of shear,

Proc. of the 1st Int. Conf. on Soil Mechanics, 1, 54-56, 1936

[22] Tice, A. R. , Burrous, C. M. , Anderson, D. M. , Determination of unfrozen water in frozen soil by pulsed nuclear magnetic resonance, Proc. of the 3rd Int. Conf. on Permafrost, Edmonton, Canada, vol. 1, 149-155, 1978