Compaction in media with inner boundaries Текст научной статьи по специальности «Математика»

RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 4, NO. 5, PAGES 301-329, OCTOBER 2002

Compaction in media with inner boundaries

A. V. Karakin

VNIIGeosystem

Abstract. A general formulation of compaction problems is described and the classification of these problems in accordance with boundary conditions is presented. A basically new class of compaction boundary problems in a finite region with moving inner boundaries is specified. New models of compaction describe porous media in which qualitative microstructural changes and associated sharp changes in material properties take place. At these time moments, the permeability of the medium drops to zero (at minimum porosity) or the skeleton loses its cohesion and breaks down. In the latter case, the porous medium is transformed into concentrated dispersive mixtures of the suspension, emulsion or sol type. Phenomena associated with jump-like changes in material properties and resembling phase transformations in homogeneous media arise at inner boundaries. Several physicochemical and mathematical aspects of these problems are discussed. Examples of analytical solutions of some problems are presented. In particular, the related wave processes are analyzed. The theory of compaction has numerous applications in geophysics and technological problems of the chemical, petrochemical, food-processing and biochemical industries. This class of compaction models is advantageous to the study of such processes. As an example, the paper presents a review of several geophysical applications including problems of mud volcanism, sedimentation and motion of partially molten mantle rocks.

1. Introduction

General description. Deformation processes in porous saturated media with a viscous skeleton are analyzed. These processes include the deformation of a two-phase medium and the motion of a fluid through the viscous skeleton. For brevity, such mixtures will be referred to as poroviscous media. The dynamics of poroviscous media is called “viscous consolidation” (in Russian literature and databases) and “compaction” (in Western literature). The term compaction is not always appropriate because in modern databases it is often associated with a variety of diverse notions. Moreover, Biot has already introduced the term “consolidation” and there is no reason to change it. Below, I use both terms as equivalent ones. In the narrow sense of this word, compaction will be meant as the compression of skeleton in the

Copyright 2002 by the Russian Journal of Earth Sciences.

Paper number TJE02103.

ISSN: 1681-1208 (online)

The online version of this paper was published 20 November 2002. URL: http://rjes.wdcb.ru/v04/tje02103/tje02103.htm

process of viscous consolidation. Some results of compaction studies widely known in the Western literature had previously been obtained by Russian authors but remained unknown to Western readers. One of the goals of this work is to fill this gap.

As distinct from equations of hydrodynamics, elasticity, poroelasticity and dynamics of suspensions, compaction equations of motion contain the intrinsic length of compaction Hc that has the following physical meaning. I introduce the notions of small, large and intermediate regions of definition whose linear sizes are, respectively, much smaller than, much larger than and on the order of Hc. In the intermediate region, the forces of viscous friction between the melt and skeleton are of the same order as the forces due to viscous deformations of the skeleton. If the linear size of the compaction region are comparable with Hc, the viscous friction forces of the fluid in the Darcy law are of the same order of magnitude as the forces producing volume and shear deformations of the skeleton. If the compaction region size is smaller than the compaction length, the viscous friction forces of the fluid can be neglected due to the smallness of the pore pressure gradient term in equations of motion and the skeleton is deformed in the drained regime (as if no fluid were present). If the characteristic size of deformations

of the porous medium is much larger than the compaction length, the volume strain of the skeleton (associated with the fluid filtration) is much smaller than the shear strain of the two-phase medium. The latter is deformed in this case as a single-phase Newtonian fluid, and the filtration velocity is much smaller compared to the skeleton velocity. In other words, given a sufficiently large characteristic size of the flow, the compaction equation degenerates into ordinary equations of the Stokes hydrodynamics, although the medium still consists of two phases and the filtration process is in progress. Thus, characteristics of physical processes and the type of constitutive equations depend on the scale of processes.

Basic equations of compaction are essentially nonlinear. Volume deformations of the skeleton are accompanied by variations in porosity entering in the continuity and constitutive equations (in the expressions for the permeability and both viscosities). The geometrical nonlinearity caused by a nonlinear dependence of equations on the porosity should be distinguished from the nonlinearity of rheological equations (the physical nonlinearity). In what follows, linear rheologi-cal equations are mostly examined.

Another basic distinction of compaction processes consists in the fact that the static state in a poroviscous medium is only achievable in the absence of external volume forces. Given different densities of the skeleton and fluid, a porovis-cous medium in the gravity field cannot be in a static state. If the external action stops or, to the contrary, becomes overly intense, the poroviscous medium inevitably changes its state after some time interval, separating into fluid and skeleton phases. For this reason, a change in the state of medium is not an external induced process but a natural phenomenon of a deep physical origin inherent in such media. Therefore, the theory of a poroviscous medium should necessarily include processes of state changes. Otherwise, the theory is incomplete and internally inconsistent and is not free of diverse logical paradoxes. This basic aspect is a starting-point of this work.

Brief review. For the first time, formal hydrodynamic equations of a compressible viscous fluid including the shear and volume viscosities, as well as an independent value that can be interpreted as the pore pressure, were derived by Landau and Lifshits [1953] on the basis of general physical principles. The expression for the volume viscosity of a two-phase gas-water system was obtained by Taylor [1954]. Gas bubbles in this system are surrounded by water. The same equations can be interpreted as a limiting case of the viscoelastic equations of Biot. As distinct from the two-velocity model of Biot, the aforementioned gas-liquid mixture is characterized by a single velocity. Thus, the general formulation of the viscous consolidation equations stems from classical works of the mid-20th century. Main difficulties of this theory are related to boundary conditions and its specific effects.

In describing the dynamics of a two-phase (melt-skeleton) two-velocity system, Sleep [1974] wrote out mass and momentum conservation equations for both phases, but he did not account for the volume viscosity and did not attempt to solve these equations. Several authors believed that the skeleton compression is fully controlled by the flow of the melt from the matrix and does not depend on the skele-

ton viscosity [Ahern and Turcotte, 1979; Maal\e and Scheie, 1982; Walker et al., 1978]. Complete equations of compaction describing the dynamics of melting ice and including the volume and shear viscosities of the skeleton, as well as the fluid flow, were derived for the first time in [Karakin, 1974]. In this work, constitutive equations were derived on the basis of thermodynamic principles of irreversible processes, and the notion of the intrinsic length of compaction was introduced. The solution of compaction equations describing the dynamics of partially molten rocks in axial zones of mid-ocean ridges was used for the first time in [Karakin and Lobkovskii, 1979a, 1979b, 1982]. Two of these works, written in Russian, remained unknown to Western readers.

McKenzie [1984] derived once more the general compaction equations and was first to derive nonisothermal equations including phase transitions between the skeleton and melt. He also examined some applications of this model. This and subsequent studies were published in central Western journals and were widely acknowledged in the Western scientific community. Afterward, the compaction equations were widely applied to the modeling of motion of partially molten rocks including filtration of the melt [Cox et al., 1993]. Modeling of nonisothermal compaction incorporating phase transformations significantly complicates a model, impeding its analysis [Karakin, 1988; McKenzie, 1984; Schubert et al., 2001]. For this reason, the motion of partially molten rocks in mantle and volcanoes was mainly studied in terms of isothermal models.

Numerical methods were applied to the evolution of an initially inhomogeneous melt distribution in terms of a 1-D model in an infinite region [Scott, 1988; Scott and Stevenson, 1984]. Larger disturbances were shown to overtake smaller ones. When colliding, these disturbances behave like solitary waves; i.e. they preserve their form after the collision. Periodical disturbances associated with a periodical boundary condition at the lower surface below which partial melting starts were also studied. Scott and Stevenson [1986] examined a 2-D flow and, in particular, solitary waves of the cylindrical shape (magmons). Using numerical and analytical methods, they inspected the stability of a 1-D solution. (The latter was superimposed by a white noise disturbance). The 1-D solution was found to be unstable in two dimensions: it was transformed with time into a periodical 2-D flow that also had a wave-like form. The resulting numerical 2-D solution was consistent with results of an analytical study. The second problem formulated in their work related to the collision of two 2-D disturbances. Similar to the 1-D case, the larger disturbance overtook the smaller one; after their collision, the larger disturbance, preserving its shape, departed from the smaller one. The latter did not preserve its shape but was transformed into an irregular disordered perturbation.

Scott [1988] studied in more detail the motion of a 2-D magmon in an infinite region. As distinct from the previous study, a boundary condition (a rigid impermeable freeslip boundary) was specified at infinity, which gave rise to a curl component of velocity that was absent in the solution obtained in [Scott and Stevenson, 1986]. The solution was shown to preserve its waveform during motion. Similar results were obtained by Spiegelman [1993a, 1993b] and

Spiegelmann and McKenzie [1987]. The latter authors examined, in particular, 2-D disturbances under mid-ocean ridges and island arcs. The movement of a 3-D disturbance in an infinite region was studied by Wiggins and Spiegelman [1995].

Richardson [1998] numerically examined the 2-D motion of a magmon in the presence of preset homogeneous devia-toric strains in the mantle resulting in vertical compression and horizontal extension. He showed that, without the preset deformations, the magmon propagates as a wave preserving its shape. However, with the strain applied in the mantle, the magmon is gradually deformed. Its circular shape stretches upward and becomes similar to a spindle. Moreover, these deformations are asymmetric relative to the horizontal axis. The compression of the magmon in its lower part is much larger than in its upper part. Thus, viscous flows in mantle lead to the instability of magmons and make them elongated in the vertical direction.

The melt concentration effect on the matrix viscosity was investigated in [Khodakovskii et al., 1995; Trubitsyn et al., 1998]. An increase in the melt concentration decreases the matrix viscosity by many orders (from 1018 to 1010 Pa s in the asthenosphere and to a lesser extent in magma chambers). This dramatic drop in the viscosity takes place at a porosity of about 5%. This effect enhances the instability of uniform melt flow, giving rise to a new type of waves that arise even in the absence of a density inversion in a partially melted medium.

Wave solutions in an infinite region were analyzed in [Bar-cilon and Lovera, 1989; Barcilon and Richter, 1986]. The method of dynamic systems was applied by Gordin and Karakin [1990] to the study of wave solutions of 1-D compaction equations.

Problems with inner boundaries. The novel aspect of this work is related to the examination of changes in the state of a poroviscous medium at inner boundaries. This class of boundary problems of compaction was specified in [Karakin, 1986, 1990, 1999; Karakin and Lobkovskii, 1979a, 1979b]. The state changes in question are largely similar to phase transformations in a homogeneous medium. Therefore, a microstructural change in a poroviscous medium is referred to below as a structural-phase transformation as distinct from actual phase transformations in homogeneous media.

Two types of inner boundaries are recognizable. Given a certain minimum value of porosity, pores at boundaries of type A become isolated and the fluid flow stops (or contrariwise, isolated pores become open). Boundary problems with A-type boundaries were numerically solved in [Karakin and Levitan, 1993; Karakin et al., 2001]. As distinct from my previous works, nonisothermal conditions of closure and opening of pores are examined in the present paper.

The skeleton is broken down at boundaries of type B, giving rise to the formation of a concentrated suspension or cluster structures. This type of boundaries was studied in [Karakin, 1999], where boundary conditions depending on physical processes at such boundaries were formulated. Unlike the paper cited, this work admits an unstable state of a two-phase mixture transitional between the states of suspension and compaction. This transitional state allows the presence of cluster structures of various scales that can break down and arise under the action of various factors.

Moreover, somewhat different boundary conditions are formulated.

The incorporation of processes involving inner boundaries substantially widens the class of viscous consolidation models. Addendum A presents the classification of their types and gives various forms of initial equations adapted to specific boundary problems.

Applications. Partially molten mantle rocks have been studied in most detail both theoretically and experimentally. A fundamental problem of the Earth’s evolution is the differentiation of the primordial protoplanetary material into lighter and heavier fractions. Viscous consolidation (and related processes) is one of the main mechanisms of this differentiation. This process is most effective at mid-ocean ridges. The motion of a two-phase mixture in magmatic and mud volcanoes can also be described in terms of the dynamics of a poroviscous medium.

Compaction can also take place in liquefied grounds, silts, running sands and in numerous technological processes involving liquid and paste-like mixtures and colloid solutions, whenever both low- and high-viscosity phases are present. Such mixtures are used in the chemical, petrochemical, biochemical and food-processing industries involving their mixing and processing. Similar processes can arise during the movement of liquid multiphase mixtures through pipelines. All these technologies involve phenomena that are described in terms of the viscous consolidation.

Gigantic tailings ponds form during the development of bituminous sands with the help of solutions containing surfactants. These ponds contain colloid solutions of hydrocarbons that are difficult to be recultivated. Sedimentation and consolidation of colloidal structures in the ponds last for decades. Ecological problems restrict the potential of bitumen extraction with the help of such technologies [MacKinnon, 1989]. Colloidal structures form a coherent skeleton through which water flows. The skeleton has a very low viscosity and a density differing only insignificantly from the water density. Consequently, compaction of a poro-viscous colloidal mixture can last for a very long time. To recultivate such ponds, one should know mechanical properties of this mixture. In this respect, compaction mechanics and the results obtained in this work can prove helpful.

A general formulation of compaction equations is given for both isothermal and nonisothermal processes. However, most results are obtained in the isothermal case. The examples presented in this work mostly illustrate analytical solutions of the new class of problems including mobile inner boundaries. A review of geophysical applications is given in Conclusion. Techniques related to transformations of initial equations and cumbersome analytical examination of solutions are presented in the addenda.

2. Formulation of Main Equations of Compaction

General derivation of equations. The first group of compaction equations includes balance relations for momentum,

mass and energy that hold true in any heterogeneous two-phase medium:

Oc

do;.

d f fi* —pf f + div i pf f wf

+div

(1 - fpsw1

= -q

= (1 — f )o|j — fpf Si;

pc

dT

-dr- + (wsV) T + pf cf (wsv) T

dt

(2.1a)

(2.1b)

(2.1c)

(2.1d)

(2.1e)

The following notation is adopted here: o. and o|.,

ten-

: °i/- Sij .

(2.2)

ef

dw1 dw1 2 dw1

„dw1

• = "(a. + a. - 3ai j S« ■ (2.3)

where Z and n are the coefficients of volume and shear viscosity. The Darcy law holds true in the fluid phase,

nf f  = - f Si + p/gi.

(2.4)

convenience, I introduce the coefficient of hydraulic resis-In the general case, Z,n and S are functions

n*

tance S = t~ .

k

of porosity and tensor invariants of stress and strain rate. Given a linear rheology of the phases, these values depend on the porosity alone.

Initial equations (2.1), (2.4) can be transformed to equations of motion expressed through the velocities of the skeleton and fluid flow:

grad

Z + 3 n ) divw1

+ nV2w1

= -Ss - Ap( 1 - f)<

(2.5)

sors of stresses in a two-phase medium and stresses averaged over the skeleton phase; , pore pressure; ws and , velocities of skeleton and fluid; g = (0,0,-g), gravity; p = (1 — f )ps + fpf, where ps and pf are the densities of both phases assumed constant; T, temperature; A, heat conductivity coefficient; l, latent heat of melting; $, density of volume heat sources; q, density of mass sources of the fluid; t, time; and c, specific heat of two-phase mixture determined from the expression pc = pf cf + pscs, where cf and cs are specific heats of the phases. The vertical axis of the Cartesian coordinate system (xi,x2,X3) is directed upward. The indexes s and f refer to the skeleton and fluid, respectively. These equations are valid in a compaction region fi bounded by a surface r. The region of definition can be either finite or infinite. In particular, it can occupy the entire 3-D space.

The second group comprises specific relations (constitutive equations) associated with properties and microscopic structure of the medium, as well as with physical processes in the medium proper and at its boundaries. A heterogeneous two-phase medium is a porous medium if it includes a skeleton through which a fluid flows. Mathematically, this implies that the stress tensor consists of the pore pressure and the effective stress tensor <7^/ depending on gradws :

where Ap = ps — pf is the density difference between the skeleton and fluid.

The right-hand side of (2.5) includes the buoyancy force represented by an inhomogeneous term of the equation of motion. It is this term that precludes a static state of a poroviscous medium in the gravity field.

A typical scheme of transition to a dimensionless form of these equations is presented in Addendum B for various types of boundary problems. The choice of characteristic parameters and scales depends on specific properties of concrete problems. This procedure involves a certain combination that has the dimension of length. It is called the compaction length Hc and is defined by the expression

Z

Hc = W ^- =

Henceforward, the asterisk means the scale of the respective value. As shown in Introduction, the ratio of the region size to the compaction length controls the compaction regime and the type of basic equations.

If true phase transformations are absent in the porovis-cous medium (q = 0), the continuity equations can also be written in a more convenient form:

Ijt + div (s + fws) = 0

divS + divw1 = 0

(2.6a)

(2.6b)

This relation was proposed by Terzaghi and is sometimes named the Terzaghi relation. The constitutive equation of the skeleton is

Here, s = f — wis the fluid flow velocity, n/ is the fluid viscosity and k is the permeability. For the sake of

If phase transformations are present, additional conditions are required in order to determine the source function q. This case was analyzed in [Karakin, 1988; McKenzie, 1984; Ribe, 1985].

The key point of this work is the statement that there exist two critical values of porosity fmin and fmax that are achieved on inner boundaries rint. Structural-phase transformations of the medium (i.e. qualitative changes in its microscopic structure) mentioned in Introduction occur exactly at these inner boundaries. The latter can be of two types: rmin (type A) and rmax (type B). Pores close at a rmin boundary (f = fmin). Given a maximum porosity (f = fmax), the skeleton breaks down at a rmax boundary. The entire two-phase region including inner boundaries of the above types divides into two regions, fi = :

the region of the compaction proper fic and the contiguous region fiv. The latter includes a suspension, a poroviscous

d

l

o

o

medium with isolated pores and a medium having a different microstructure.

A specific type of boundary problems of compaction depends on the part of the interval (/min,/max) within which the porosity varies during motion and on the boundary of the poroviscous medium. Problems of types A and B arise if the compaction region is adjacent to boundaries rmin or rmax, respectively. Regions adjacent to these boundaries are denoted as fic(A) and fic (B). Constitutive relations in problems of type A are usually given by inverse power functions of porosity and in most cases one may assume that the shear viscosity does not depend on porosity:

fic(A): Z +4n = Af, 5 = Bk, n = Const , (2.7)

where and Bk are material constants.

Analysis of skeleton breakdown processes in poroviscous media. Note that the breakdown and recovery of the skeleton proceeds differently depending on the physical conditions at the inner boundary. In particular, the poroviscous medium can be transformed into a concentrated suspension. From the physical standpoint, the skeleton breakdown means that the tensor a/ vanishes. At this moment fmax value is on the order of 20-30%. Formally, a suspension can also be considered as a poroviscous medium in which the effective stress tensor is zero. The relative volume concentration of the fluid in suspension f is called porousness and is denoted by the same symbol as the porosity.

Note that the breakdown and recovery of the skeleton occur under basically different conditions, which naturally leads to different boundary conditions. The skeleton breakdown at the trailing front requires a certain energy to be expended in forming new surfaces. For this reason, the failure of the skeleton on surfaces is preferable to its failure throughout the volume. Consequently, the skeleton breakdown in a poroviscous medium can be accompanied by the formation of clusters, i.e. blocks of the poroviscous medium of various sizes. These blocks can be connected and disconnected in a random manner. If clusters are present, the porovis-cous medium and suspension are separated by transitional cluster zones rather than an inner boundary rmax. These zones consist of material that possesses specific properties and does not meet conditions of a heterogeneous medium including compaction or suspension equations. It possesses a memory and, possibly, thixotropic properties. Clusters are comparable in size with the definition region. Therefore, the constitutive equations cannot be local in space and time and must explicitly depend on time. In other words, the medium should be described by certain kinetic integrodiffer-ential equations of stochastic processes.

The situation is basically different in the case of consolidation of suspension and formation of coherent skeleton. Too close convergence of suspension particles is hindered by repulsive hydrodynamic forces described in terms of the lubricating layer approximation of hydrodynamic equations. They arise whenever the longitudinal size of spac-ings between particles is much greater than their lateral size [Karakin and Shklover, 1990]. If molecular forces of converging particles are sufficiently small, the process of their

convergence is uniform. The skeleton formation starts when all particles come into contact with each other through a thin lubricating layer. The contact area is initially small compared to the overall surface of particles. Therefore, the motion of the fluid in the poroviscous medium initially differs little from its motion in the suspension (at similar values of phase concentrations). However, the interaction of contacting particles qualitatively changes. This means that, at the moment when the skeleton arises, its volume viscosity changes in a jump-like manner from zero to a finite value, but the permeability at this moment is very large (i.e. the hydraulic resistance is close to zero). This regime at the inner boundary coinciding with the leading front is described by the condition imposed on material functions,

Tmax : Z + 4n = 0, n = 0, S = 0, at f = fmax . (2.8)

On the other hand, if molecular forces are sufficiently large and exceed the repulsive hydrodynamic forces, the formation of the skeleton can be accompanied by the formation of the same cluster structures as those arising during its breakdown.

It is also possible that clusters cannot form at any surface due to the action of such physical factors as large stresses and strains, vibratory and acoustic effects, and other physical and physicochemical effects. In particular, surface forces at an inner boundary can produce in its vicinity large (formally, infinite) strains. An inner boundary rmax arises in the absence of cluster structures in the transitional layer between the poroviscous medium and suspension. In some cases, if external forces act near the inner boundary, viscous forces can be neglected near this boundary. Then, condition (2.8) changes to

rmax : Z + 4n = 0, n = 0, S = 0, at f = fmax . (2.9)

Relations (2.8) and (2.9) hold only at the inner boundary rmax between the poroviscous medium and suspension. If a transitional region including cluster structures exists between these states, other conditions should be satisfied at this boundary, but such a situation is not discussed in this work.

It is reasonable to assume that the region fic(B) accounts for a small porosity interval Af of the entire porosity range (fmin,fmax). The skeleton in this interval weakens prior to its breakdown at the upper boundary rmax. The width of the region fic(B) is characterized in this interval by a small parameter e1,

ei = fAf « 1 . (2.10)

fmax

However, this does not mean that the region fic(B) is small. Likewise, the viscous forces controlled by the buoyancy force are not small. The parameter e1 is a material constant and is assumed to be known.

The porosity at the boundary rmax cannot be considered small. However, relations in the fic(B) vicinity of this boundary can be simplified in a different way. If an inner boundary rmax of this type exists, all material functions are

assumed to be analytical up to this boundary and to have one-sided derivatives on it. This means that the expression for viscous forces including first derivatives of the effective stress tensor is physically meaningful up to the very boundary. I expand the material functions in the neighborhood of this critical point in fic(B) and use conditions (2.8) or (2.9). Rejecting higher-order terms of the expansion, I obtain

r int : wm ----------- 0

(3.1a)

(3.1b)

4 Z * I

0,(5): Z +3n — ^^(/max - /)

3 fmax

n —

n

/m

-(/max - /y , '(/max - /)k .

(2.11)

Below I address two situations. In the first case, the volume viscosity remains finite at l = 0 and k = 1 (condition (2.8)):

4

0c(B) : Z + g n — Const , n — Const

/m

(/max - /) .

(2.12)

In the second case l = 1 and k = 1 (condition (2.9)) and all material functions vanish at the boundary:

4

°C(B) : Z + 3n r

3 /m

Z *

(/max - /)

n—

5 :

/m

" (/max - /) , (/max - /) .

(2.13)

(1 - /)PS

(1 - / ) pS-wn

(3.1c)

(3.1d)

where [p]|+ = <p+ — ip-, with ip± being boundary values of the quantity <p at the inner boundary inside and outside the compaction region; sn = Si«i and = w|m(a), with (n and m(a)) being the normal and tangential unit vectors; the index a (a = 1, 2) indicates two mutually perpendicular directions on the boundary surface; nj is the stress vector; and u is the boundary movement velocity relative to the skeleton of the porous medium (directed along the normal to the boundary).

Conditions (3.1) are quite general and can be simplified in specific cases. If no phase transitions (changing the matrix density) are present in the matrix, the velocity of the inner boundary is evidently given by the expression

Ti,

P

0

(3.2)

If the porosity is continuous and no concentrated forces are present, these conditions are simplified,

Dependences (2.12) and (2.13) of material functions are simplest in the region fic(B) and qualitatively characterizes its specific properties. Other expressions describing these functions are also conceivable.

—0

—0

(3.3a)

(3.3b)

+

+

u—

/

+

*

5

a

*

5

+

*

n

s

U — — w

n

/

*

5

+

n

+

s

w

n

3. Inner Boundaries

General boundary conditions. First I consider boundary conditions at inner boundaries rint of any type. Similar to the main equations, boundary conditions at inner boundaries are subdivided into those of general and specific types. General conditions at all inner boundaries are based on conservation laws for mass, energy, momentum and other conservative quantities and reduce to continuity conditions for their values. The tangential velocity must usually be continuous at inner boundaries with structural-phase transformations. If a boundary moves relative to the medium, the conservation law of the mass of phases implies that the porosity is discontinuous at this boundary. Then, one may assume that a normal surface force 7i (of molecular, electrostatic or other origin) is applied to the inner boundary. The stress vector at the boundary is then also discontinuous:

0 or

(3.3c)

I introduce a curvilinear orthogonal system of coordinates (yi,y2,y3) such that the coordinate surface = Const coincides with the inner boundary considered. The continuity conditions for the tangential stresses and velocities at this boundary yield the relations

—0

dys .

— 0 .

The initial condition for the porosity in all nonstationary problems is

/ — f(V)

at t — 0

(3.4)

Specific boundary conditions. The classification of types of boundary conditions and related boundary problems is given in Addendum A. Problems of types I and II are defined in an

+

+

+

f

j'j

1

2

n

n

infinite region or in a finite region with exterior boundaries. They are known in the scientific literature and have been comprehensively studied in both 2-D and 3-D variants.

Problems of type I can arise in a poroviscous region whose characteristic size considerably exceeds the compaction length. An important particular case of this type is the situation when a homogeneous stationary porous fluid flow is specified at infinity and no external force other than gravity is applied. This simplest class of trivial solutions of compaction problems is characterized by a homogeneous time-independent solution. The deviation of dynamic quantities (stress, strain, flow velocity etc.) from the trivial state tends toward zero at infinity:

s = sofc, Bfcso = /0 Ap(1 — /o)g, pa = 0, w1 = 0 .

Here, k is the unit vector directed upward along the vertical axis. If the buoyancy force is absent, the trivial solution is static (s0 = 0). This boundary problem involves asymptotic attenuation of disturbances, with its solution tending to a concrete trivial solution,

tains

->■ 0,

f = fm

0

d/

rmin : "dt — u/,n = 0 where /,n = (ngrad/) .

The substitution of the time derivative from (2.6a) gives

rn

divs + div (fws)

Tn

(3.7)

(3.5)

The compaction region fic in problems of type III is bounded (at least in part) by inner structural-phase boundaries. These problems are less studied and represent the actual object of study in the present work.

In problems of type A, pores close at an inner boundary rmin at a minimum value of porosity. The closure of pores means that this boundary is impermeable. Therefore, (3.3a) yields

Relation (3.7), determining the velocity of the inner boundary, is neither boundary nor additional condition because it is a consequence of the existing equations and boundary conditions. However, a fairly strong disturbance of the poroviscous medium can qualitatively change the behavior of this process. A pore pressure jump at the inner boundary rmin can lead to the failure of partitions in closed pores, with the pressure jump delivering certain work that makes this boundary move.

In this case, boundaries of two types, namely, leading and trailing fronts, can arise. An additional condition determining the energy expended in changing the microstructure is set at the leading front rm¿ [Karakin, 1999]. If the leading front rm+n opens pores due to the fracture of partitions in the pore space, this work is done by the pressure jump. If the amount of energy is insufficient to open pores, the boundary does not move.

The partitions are generally fractured not instantaneously but over a certain time interval during which dislocations and microfractures accumulate. This process is described by a kinetic relation with a characteristic time. If the kinetic processes are comparable in time with the inner boundary movement, the velocity of the latter depends on the excess of the pore pressure over its critical value. This situation is described by the boundary condition

(3.6)

r

(+) .

The closure and opening of pores are controlled by the following factors: phase concentrations, temperature at the phase boundary and the stress difference at this inner boundary. Phase transformations, porosity jump and surface forces are assumed to be absent at this boundary (conditions (3.2) and (3.3)). However, boundary jumps in the pore pressure and effective stress tensor are admissible here. These jumps give rise to infinite interphase forces mutually compensating each other. Any of these forces can deliver work in the related phase during the pore opening process.

A microstructural change corresponding to boundary condition (3.6) can occur with no or very slow movements of the medium. In this case, the main factor controlling structural-phase transformations in a heterogeneous mixture is the phase concentration. For example, partial melting of rocks can lead to the opening and closure of pores due to a change in temperature without any additional effects. In this case, the position and velocity of an inner boundary are governed solely by relations (3.6). Differentiation of the first equation in (3.6) with respect to time yields

df

rmin : -djf + (w sgradf ) = 0 .

In accordance with the second relation in (3.2), one ob-

f = fmin, (pa - Po) 0 (pa - Po) = au a > 0, p0 > 0, pa = - pG ,

(3.8)

where pa is the so-called piezometric (or hydrodynamic) component of the pressure, i.e. the deviation of the pore pressure from the hydrostatic pressure pG in the poroviscous medium; p0 is the critical value of the pore pressure; and a is a material constant determining the kinetics of the pore fracture process.

The critical pore pressure p0 is a function of the minimum porosity at the leading front /min :

fmin = £ (Po)

(3.9)

The function £(p0) is found from experimental data. Its form is determined by physical processes proceeding on a microscopic level at the time moment of opening of pores. This is a decreasing function because the action of pore pressure lowers the opening threshold of pores, i.e. the /min value.

If pa < p0, the boundary does not move until the porosity reaches the value £(0). From then on the movement of this boundary is governed by boundary condition (3.6) and its velocity is given by expression (3.7).

If opening of pores is associated with skeleton expansion (without loss of its coherence), the work at the leading front is done due to a jump in the effective stress tensor. In this

S

s

n

case as well, the fracture does not occur instantaneously but is a result of accumulation of fatigue microcracks. The fluid enters the latter and reduces the ultimate strength due to the Rehbinder effect. With higher stresses exceeding a certain threshold, the decrease in strength becomes more intense and this boundary moves faster. This phenomenon is also characterized by a specific time of kinetic processes. If these processes are comparable in time with the boundary movement, one obtains a condition at the inner boundary analogous to condition (3.8):

г£+П : / = /min, (CTef - сто) в (o-ef - сто) = «

a > 0, ст0 > 0 .

(3.10)

(+)

г

(+)

/ = /m

/ = /m

: p0, p0 > 0 .

CTef = сто,

сто > 0 .

(3.11)

(3.12)

=0

(3.13)

Clusters and an intermediate cluster layer, rather than inner boundaries, can arise between two states (suspension and poroviscous medium). However, clusters can break down in the presence of sufficiently large concentrated surface forces and the poroviscous medium can then come into contact with the suspension. I analyze this situation in more detail. Relation (3.1d) yields

ef

стг/ П

The effective stress tensor vanishes in the suspension. Therefore, the following relation holds in the poroviscous medium:

Here <ref is the first invariant of the effective stress tensor. The critical parameter <r0 plays the same role as p0 and satisfies a relation similar to (3.9).

Closure of pores at the trailing front rj- results in the energy release and dissipation in space. Condition (3.6) is satisfied at this boundary at /min = £(0). The velocity of the front is given by expression (3.7). No additional condition is required.

If the time of kinetic processes is much smaller than the time of inner boundary movement, the velocity of the boundary is independent of the stresses or pore pressure at this boundary. Moreover, the pore pressure and first invariant of the stress tensor in boundary conditions of types (3.8) and (3.10) exactly coincide with their critical values. Accordingly, this yields the boundary conditions

Tn

ef

(3.14)

According to the above, if surface forces are absent, the inner boundary rmax exists during the suspension consolidation and skeleton formation, i.e. at the leading front. At this boundary

/ = /m

ef

ст - UiUj ij г j

=0

(3.15)

In this case, the velocity of the inner boundary movement is found from condition (3.7).

In constructing the algorithm for choosing the type of boundary condition at an inner boundary, the following scheme can be applied. The value of pa (or <ref) at this boundary is first determined. If it exceeds the critical value p0 (or <r0), the leading front conditions are chosen. If it is smaller than the critical value, the boundary value of porosity / is found. If the porosity is /min = £(0), the trailing front condition (3.7) is chosen. If / > /min = £(0), the boundary does not move. Finally, if / < /min = £(0), pores are still interconnecting and the pore closure condition must be specified at the trailing front. The simplest assumption is that the minimum porosity at the trailing front is independent of the pore pressure and other conditions. One may then set /min = £(0) at this front. In the general case, more general conditions depending on the physical formulation of the problem are chosen at the trailing front.

In problems of type B, breakdown or formation of the skeleton takes place at inner boundaries at the maximum porosity value /max, and the pore pressure is independent,

The effective stress tensor is finite at the inner boundary in the presence of surface forces. If the volume viscosity of the skeleton near this boundary tends toward zero (see (2.13)), the strain tensor at this boundary tends toward infinity: divws ^ to, and divs ^ to. It is in this case that conditions favorable for breaking down the skeleton and clusters by the surface forces arise. They break up clusters and create a dense packing of the skeleton at the leading front, and they destroy the skeleton at the trailing front. In any case the surface forces do additional work, similar to the case with opening of pores. This work controls the velocity of the boundary rmax. Therefore, the breakdown processes are subjected to additional conditions at this boundary:

Tn

/ = /m

ef

bu

(3.16)

The sign of the kinetic coefficient b in (3.16) depends on the direction of the surface force and on its effect (i.e. on whether it is breakdown or formation of the skeleton). The breakdown processes governed by condition (3.16) can be accompanied by the formation of a transitional cluster zone. Further specification of this condition requires an additional study beyond the scope of this paper.

At first glance, conditions (3.16) and (3.7) contradict each other because the velocity cannot be arbitrarily given but is controlled by values that are already determined. Actually, no contradiction is present here. If the numerator and denominator in (3.7) tend toward infinity, the velocity of both fronts is an indeterminate value to be found from an additional relation, namely, condition (3.16).

Conditions at lateral boundaries. In addition to mobile boundaries at which adequate boundary conditions are set, there can exist boundaries representing natural restraints (e.g. walls). They are referred to as lateral boundaries and belong to the second category (see Addendum A). An example of a lateral boundary is a rigid impermeable slip-free wall (condition (A.2)):

+

г

ст

max

ст

“

p

+

f

г

p

max

sn = 0

wS = 0

(3.17a)

(3.17b)

Boundary conditions in problems of the third category include a piezometric component of pressure. For this reason, the problem formulation in the presence of the gravity at inner boundaries of the region fic requires a preliminary calculation of the pore pressure with due regard for its piezometric and hydrostatic components. Otherwise, hydrodynamic processes must be analyzed jointly in both the inner (fic) and outer (fiv) regions. In some cases this can be avoided by subtracting the hydrostatic component in the viscous consolidation equations. Note that there are various types of the hydrostatic equilibrium state; these are the states of fluid and skeleton phases and the state of incomplete hydrostatic equilibrium of a suspension in which the skeleton phase is surrounded by a fluid. On the other hand, the hydrostatic equilibrium state cannot take place in the viscous consolidation regime in the presence of gravity. Therefore, the aforementioned subtraction of the hydrostatic component can be effected by various means described in Addendum A. According to the terminology adopted here, reduced canonical forms are compaction equations subjected to the subtraction procedure and preserving the form of hydrodynamic equations.

tion region. Small, large and normal regions are defined by the respective conditions k C 1, k ^ 1 and k ~ 1.

The 1-D asymptotic conditions (3.5) can or cannot coincide in the positive and negative directions of the x axis (see (B.19)):

pa ^ 0, s ^ 0, s ^ so, f ^ f0 at I x I ^ to , (4.5a)

к2SQ = /q , ^s± = /± .

(4.5c)

Boundary conditions (3.7) and (3.9) in 1-D problem IIIA reduce to (B.21) and (B.22):

rm-n : f = £(0), s = 0 at x = x_(i) and x = x+(t) , (4.6a)

гі+n : / = too) ,

dx+ 1 f -і дs \ f i дs .

~dT = -a V дї +poJв l-/ дї - pQ) .

(4.6b)

The dimensionless expressions for initial condition (3.4) (common to all types of problems) and the total mass of disturbance p. are

/(x, 0) = <^(x) + /m

ß = I ^(x)dx .

(4.7)

4. 1-D Case (the Simplest Solutions)

Low porosity approximation. I present examples of analytical solutions in 1-D cases simplifying equations (2.5) and (2.6):

_д_

дx

Ss - (1 - /) Дрд

д f/-і д^ 2,-k -,

TxV sj =к/ s-1

д/ дs ді + дx

I derive the simplest nonstationary solutions to problems

II and IIIA. In a small region (k C 1) equation (4.3) is simplified,

— /-і д£ дx V дx

1.

(4.1)

This equation has the first integral

r дs

dC

/-1 ^ = -x + Ci(t), Ci = -^, Сі(0) = 0

dt

(4.8)

(4.9)

(s + ws)= 0 . (4.2)

The spatial coordinate of 1-D models is denoted as x. Similar simplifications are used for designating other quantities (their indexes are omitted). The procedure of the transformation to a dimensionless 1-D case is described in more detail in Addendum B. For brevity, dimensionless values in this addendum retain the notation used for dimensional quantities. These equations have the following form for problems I, II and IIIA in the low porosity approximation:

where Ci(t) is an arbitrary function. Eliminating s with the use of (4.4) and integrating this equation, I obtain the second integral of the initial equations

/1

I- 1

= - xt + Cj^t) + C2(x) at l = 1

where C2 (x) is the second function of integration. Condition (4.7) defines one of these functions:

/1

I- 1

- xt + C1(t) +

[^(x) + /min]

(4.10)

(4.3)

(4.4)

Expression (4.9) in the initial form also reduces to a quadrature that can be easily calculated in some particular cases. Let l = 0 and let the initial function be constant,

^(x) + fmin = 1 .

Here h0 is the dimensionlees length of compaction and Then, the porosity and fluid flow velocity are readily found k = h—1 is a parameter characterizing the size of the defini- as

г

ex

“

p

0

I

aX3

• « • • • • • і

• • » • • • • •

••V

%v

%v

•Ä

••V ••V ••V ••V • •

I

к • • •

і

•V

•V

>v

îèV

•V •V •V •V • •

!

Solid

phase

Suspension

‘I

X,

Porous medium ^with connected pores

permeable grid

Fluid

V

/ = xt - C^t) + 1 , S :

г : |î = 0.

дx

(4.13)

Substituting expression (4.11) into these boundary conditions, I obtain the solution

x2

f = xt +1 , s =----------2---+ so(t) .

Now I consider the motion in the problem IIIA with inner boundaries in a small region (k C 1). As before, I set l = 0. In this case, a formal solution to equations (4.8) and (4.4) with boundary conditions (4.6) reduces to the integrodiffer-ential equations

1

s(x) = — (x+ - x) (x - x^ , t

/ = xt - 1 J(x+ + x-)dt/ + <^(x)

+ dx+ 1 h h

u+ = -dT + a U - poJ eU- po

h = x+ - x- .

(4.14a)

(4.14b)

In order to obtain this solution in an explicit form, I differentiate in t the expression /(x_(t),t) = /min :

d^(x-) dx- h tdx-

ds dt 2 + dt

= 0 .

Figure 1. 1-D motion of a poroviscous medium in a vessel with free-slip walls and a permeable bottom. The bottom is made of a rigid net that let the fluid though but is impermeable for the matrix material. Arrows show movement directions of the phases. Modeled is a laboratory experiment in which the porosity at the lower boundary is higher than the critical minimum value /min but lower than /max.

This readily provides the expression for the trailing front velocity u-

dx-

dt

h dx— — t-------

2_______dt

d^(x_)

ds

(4.15)

- y + C2 (t) + C3(t) . (4.11)

This is a general solution depending on two unknown functions that can be determined from respective boundary conditions.

Now I address the compaction problem of a poroviscous medium in a vessel with slip-free walls and permeable bottom (Figure 1). The bottom is an external boundary Г at x = 0. Therefore, this is a problem of category II. The vessel bottom is assumed to be made of mesh that lets pass a fluid but is impermeable to a skeleton substance. Hence, the piezometric component of the pore pressure pa vanishes at the outer boundary. The fluid flowing through this boundary is measured, which makes the fluid flow velocity a known function. This yields the boundary conditions

Note that a similar expression for the leading front velocity u+ cannot be obtained because the boundary is displaced into the region in which the solution has not been found as yet. Therefore, the value u+ is found from condition (4.14b).

Condition (4.14b) implies that two alternative situations are possible: the right-hand boundary either moves or not. At h > 2p0 its velocity is

dx+

dt

(4.16)

Let the initial function ^(x) be given as the parabola

<^(x) = 2 (x+ (0) — x)(x — x- (0)), h(0) > 2p0 .

Substituting this expression into (4.15) yields the ordinary differential equation of the first order describing the motion of boundary points

г : p“

0

So (t)

(4.12)

- (x+ (0) + x- (0)) -x- (t) + t

dx-(t) - h =0 . (4.17) dt 2

where s0(t) is a given function. Equation (B.18c) and the first relation in (4.12) yield

The system of two equations (4.14b) and (4.17) in two unknowns defines the movements of the leading and trailing

u

s

fronts of a disturbance. A direct substitution shows that these equations have at least one solution,

x_ (t) = x_ (0) + t , x+ (t) = x+ (0) + t ,

describing a solitary stationary wave of a unit velocity.

The leading front stops at h < 2p0. For definiteness, I assume x+(0) = 1 and x_(0) = —1. Then

<^(x) = 1 (1 — x2) at t = 0 .

In this case, the solution of (4.14a) reduces to the expressions

s(x) = 1 (h2 — x2) ,

dh

dt

h

2( t + h)

д2s

дx2

s1

Problems with the skeleton breakdown. The porosity in problems of type IIIB is not a small value, but another small parameter (єі) is introduced. Transformation to dimension-less variables and the use of the expressions for material functions (2.12) and (2.13) simplify equation (4.1) for each of these cases:

lx (Ü) = к2(— - /)s -

1

1 - /max

S (<‘- f) D = "Ч1 - f)'- (T-^

Equation (4.4) is simplified in a similar way:

(4.18)

t

/ = xt + 2 J h(t/)dt/ + — (l - x2) ,

Q

t

—ht +2 J h(t )dt + 2 (l — h ) = 0 .

Differentiating in time the last expression in (4.18), I obtain the first-order equation

д/ + (2_ / ) ^

+ V /maxJ дx

дx

0.

(4.19a)

(4.19b)

(4.20)

In the case of these problems, boundary conditions (3.3) and (3.16) assume the form (see also (B.24)):

г : / — / —x — 2 (1 — / )1 (l — / )

г max : / / max , , — k V / / V ./max у

dt

pb

b

0

дx ’ (4.21)

стзз

that can be readily integrated, implicitly defining the function h(t); more specifically, it gives the inverse function

t = 3 (— h

With t ^ to, the function h(t) tends asymptotically to 3

The value smax in (4.21) is not an a priori given constant but is to be determined in the process of solution. Therefore, a correct formulation of the problems requires an additional condition (the second relation in (4.21)) determining this quantity.

Condition (B.26) following from (3.15) in the absence of surface forces is

(4.22)

zero as h(t) = y —.

Now I address linear problems arising if l = 0, k = 0 and k = 1. Equation (4.3) is then simplified,

and is easily integrated.

Problem 1. A fluid flows from below into the half-space x > 0, with a pore pressure at its boundary given by a known function

ds

r: pa = -dx = A(t).

A stationary flow at infinity is given by the trivial solution s = 1, f = f0. Disturbances at infinity decay and the solution tends to the trivial solution. A homogeneous initial solution corresponds to this trivial solution:

s ^ 1 , f ^ f0 at x ^ to ,

<^(x) = fo .

The solution of this boundary problem has the form

= 1 + A(t)e-x , /(xt) = e-x / A(t/)dt/ + /о .

5. Wave Processes

General description. All quantities in 1-D wave solutions are functions of the wave argument x' = x — ut. A particular case of wave processes is traveling stationary waves. As a rule, they have leading and trailing fronts or attenuate at infinity. However, a wave solution can be involved in the description of essentially nonstationary and nonwave processes. For example, regions of the leading and trailing fronts can be described by the wave solution if the movement of these fronts is stationary. The fronts can be linked with the main nonstationary solution.

Traveling solitary waves naturally arise in problems of types I and III. Periodic wave solutions belong to the class of type-II problems if periodical boundary conditions are considered as conditions in a finite region. Below, primes and other symbols distinguishing wave solutions are omitted. The wave process develops from an initial disturbance but with time the wave “forgets” its initial state, resulting in a somewhat different set of constitutive parameters in reducing equations to a dimensionless form. The characteristic initial size H(0) is replaced by the intrinsic length of com-

b

s

s

ст

ax

paction. In other respects, the reasoning remains unchanged. This implies that a renormalization of dimensionless values should be made and k = 1 should be set in all of the above expressions.

Low-porosity approximation. First I address low porosity problems of types I, II and IIIA. The change to the wave argument simplifies equation (4.4) and its integration readily yields

s = uf + Const .

The integration constant is found from boundary conditions B; in particular, asymptotic conditions (4.5) reduce to a simple expression:

' — so = u(/ — /o) .

(5.1)

this type includes asymptotic condition.

Wave equation (5.4) is a two-parameter equation because it contains the independent parameters /0 and u. However, the meaning of the parameter /0 is not the same in problems of different types. It is easy to see that this equation is symmetric with respect to the change of variables x ^ — x. Therefore, it is reasonable to identify symmetric wave solutions and to place their symmetry center in the origin of coordinates. In a finite region (e.g. in a problem IIIA) the

h

wavelength is equal to h, and x± = ± —. These solutions can be interpreted in terms of solitary waves. Asymmetric wave solutions, the physical meaning of which is analyzed below, also exist.

Now I address problem (5.4)-(5.5) at l = 0 and k = 3:

Conditions (4.3) that must hold in problems of type IIIA give the relation

(5.6)

s = u(/ — /min)

(5.2)

where fmin is a material constant.

1

I set /0 = /Q — — sQ and

s0 = /o in the infinite region and /0 = /min in the finite region. Then expressions (5.1) and (5.2) can be written in the general form

s = u(/ - /Q) .

-x /-1-x ) = (/ - /Q)/-k - —

hh

--< x < -

2 2

г

(-)

/ = £(0)

^in : / = £(0), ( ) -/ = —au — pQ

(£(0))

dx

at x = — .

rmin: f=£(0), udx =—au—P0 at x =h. (5.7)

This problem is studied in Addendum C, and its general symmetric solution is given by quadrature (C.14):

(5.3)

x(/) = ±

d/

In this case, equation (4.3) reduces to the ordinary differential equation

(5.8)

(5.4)

Two types of wave equations exist. The first type is a solitary wave in infinite or finite regions. In this case, the wave equation must be complemented by asymptotic condition (4.5a) or boundary conditions (4.6). Note that the respective values £(p0) and £(0) at the leading and trailing fronts are generally different. With £(p0) < £(0), any disturbance in a finite region with moving boundaries infinitely increases with time. Under this condition, a solution in the form of a solitary wave does not exist.

Given small perturbations, one may approximately set £(p0) = £(0), which simplifies boundary conditions (4.6):

1/

©(/) = -/TR(/),

R(/) = 2/2 — (2 — /0)/ + /0 .

u

Its solution in an infinite region with asymptotics (4.5a) is provided by expression (C.17). The velocity of the wave depends on its mass ^ and is given by formula (C.18):

= 3 TW— - /о (2 + 2/q - /q2) .

(5.5a)

(5.5b)

As seen from this formula, the velocity is proportional to the squared mass. This is quite conceivable: a larger mass results in higher fluid flow velocity and therefore the velocity of the center of mass coinciding with the wave velocity. Richter and McKenzie [1984] were the first to obtain this analytical wave solution in an infinite region. Importantly, the solution in an infinite region is matched by a curve in the space of the parameters /0 and u entering wave equation (5.4) in the form (C.11):

8/0

Stationary solitary waves can exist in this particular case.

The second type of wave solutions (in infinite or finite regions) does reduce to stationary waves and usually describes various transient processes at the leading or trailing fronts of more complex nonstationary disturbances. In particular,

(2 - f°)2

on the other hand, wave solutions in a finite region with moving inner boundaries occupy regions of a finite area in the space of these parameters. In other words, solutions in a finite region are much more numerous compared to the solutions in an infinite regions. The latter are rather an exception to the general set of solutions.

Solutions in finite region with boundary conditions (5.7) exist not for all values of parameters. Conditions under

і

which such solutions exist are examined by the method of dynamic systems in Addendum C, where a solitary wave in a finite region is also calculated analytically for equation (5.6).

Problems with the skeleton breakdown. Now I analyze the wave process in problems of type B and try to find their solutions in the form of a traveling solitary wave. Integrating the continuity equation (4.20) under condition (4.21) yields

1 - /max s - smax = u / - 1 .

(5.9)

д2/

1

1x2 = (i - /)(/ -1+/Q) - u

(5.10)

(5.14a) alone under condition that the integration constant vanishes,

2 /2 — 1 (1—/ )s+i (1 — / 0)(1 — / )2—=0.

2

2

The solution of this equation reduces to the quadrature

The fluid flow velocity is assumed to be directed in the positive direction of the x axis. Since s < smax and / < 1 in the physical region, (5.9) implies that u < 0; i.e. the inner boundary movement is opposite in direction to the fluid flow velocity. The substitution of this expression into (4.19) yields, in dimensionless variables,

f3

d/ h h

x(/) = ± ------< x < -

(/ ) J У© 2 2

f

(5.15)

where ±— is the position of the leading or trailing front of the boundary rmax and

e=3 (1 - / )3-(l - / 0)(1 - / )’+2(—

This wave has the length

If(i-/)Ix)= (1 -/)(/-1+/0)- —, (511)

where

/0 = - (1 - /max )

max max

Two boundary conditions are possible. One implies the absence of surface forces:

г : / =0

г max : r-, — 0

ix

(1 - / ) ix

(5.12a)

(5.12b)

the second condition with surface forces present is

h=2

y©

Expression (5.15) provides a solution in the form of a traveling solitary wave and is examined in detail in Addendum C. As shown, the solitary wave can only exist in the presence of surface forces. All other cases are consistent with nonsolitary wave solutions. The latter are usually obtained for the leading or trailing fronts of nonstationary solutions, provided that these fronts move slowly enough.

Conclusions

(i - / ) |x.

(5.13a)

(5.13b)

Equations (5.10) and (5.11) have, respectively, the first

integrals

2- 3Iі - /)3 + + — і— - / o)(i - / )’ -

(5.14a)

Const

(5.14b)

2u

Conditions (5.12) and (5.13) determine the constants of integration. However, the solution formally written in terms of quadratures provides no constraints on its existence. The analysis of the existence of wave solutions by the method of dynamic systems is given in Addendum C. The solution in the form of a traveling wave is shown to exist for expression

A general formulation of inertia-free compaction problems accounting for structural-phase transformations at inner boundaries is presented. The compaction models known in the literature are shown to belong to a particular class of a more general group of compaction models including problems with inner boundaries.

The problem of phase transformations of a special type giving rise to qualitative microstructural changes at inner boundaries of a heterogeneous medium is treated for the first time within he framework of compaction models. Processes of opening and closure of pores (with a minimum porosity), as well as processes of breakdown and formation of the skeleton, can take place at inner boundaries. These processes are analogous to phase transitions in homogeneous media and have many common features. As regards the mechanics of porous media, such phenomena were only considered in the percolation theory of media with a rigid skeleton. A transitional zone with a cluster structure can exist between the compaction and suspension regions and it should be described by kinetic and stochastic constitutive equations. The analysis of the latter is beyond the scope of this work.

Various canonical forms of compaction equations admitting a correct formulation of boundary problems are shown to exist. They are distinguished by boundary conditions for

з

г

0

max

г

max

г

b

max

Figure 2. Schematic layout of a mud volcano. The volcanic cone and channel are located in the upper and central parts, respectively; the mechanism feeding the volcano from the fractured zone of a waveguide is schematically shown in the lower part of the figure. Under certain conditions, the mud mixture moving through the channel assumes the compaction regime.

the pore pressure and, accordingly, by the procedure of subtraction of its hydrostatic component. Boundary problems of compaction are systematized. The generalization of the well-known compaction theory made by the author to the case of bounded regions with mobile boundaries is analyzed in detail. As an example, analytical solutions are derived for some of the simplest (including wave-like) processes.

To date, most applications of the compaction theory lie in the field of geophysics, where they are used however in a comparatively narrow aspect related to the study of partially melted rocks in the asthenosphere and volcanoes. The new class of models is appropriate for not only the analysis of geodynamic and fluid-dynamic problems in the crust

and mantle but also the investigation of a wider range of problems. The latter includes, in particular, the dynamics of heterogeneous mixtures used in chemical, petrochemical and biochemical technologies, as well as in petroleum and food processing industries.

As mentioned in Introduction, a melt is released by a partially melted mantle material at mid-ocean ridges and hotspots [Cordery and Morgan, 1993]. This process is schematically illustrated in Figure 2 according to the results of Schubert et al. [2001]. The figure shows the 2-D distribution of flow lines in a poroviscous medium and contours of the mantle material melting extent as determined from the well-known Ringwood diagram of phase equilibrium. The melting extent contours were constructed under the assumption of no filtration in the two-phase medium. Filtration must change their pattern. Figures 3b and 3c show the velocity of the melt and porosity in the partial melting zone as a function of the melt viscosity and the velocity of the ascending mantle flow. These dependences are obtained using balance relations.

In order to examine this process in detail, one should solve the compaction problem in a partially molten asthenosphere under a mid-ocean ridge or a hotspot. This problem includes conditions at inner mobile boundaries where pores open and close. The problem cannot be correctly formulated without the incorporation of this factor. The present work opens the way for solving this class of problems. The major difficulty here is the incorporation of the temperature regime and phase transformations. These problems lie beyond the scope of this paper.

The process illustrated in Figure 3 was numerically modeled by Trubitsyn et al. [2003], who analyzed the motion of a partially molten material in the axial zone of a mid-ocean ridge. The axial flow feeds the magma chamber, and the peripheral flow strikes the cold lithosphere and, upon solidifying, accretes it. The axial flow was shown to be rather insensitive to external effects. On the other hand, the peripheral flow is unstable and can accumulate the melt at the lithosphere boundary. However, the actual melt does not accumulate but flows away toward the magma chamber. This unstable regime of the peripheral flow is responsible for the formation of the chemical heterogeneities observed in the oceanic crust of mid-ocean ridges.

Compaction processes can take place in both magmatic and mud volcanoes. If the magma moving through a volcanic channel is a two-phase mixture (e.g. gas-melt or melt-solid), this mixture can become a poroviscous medium under certain conditions and can then be described by a compaction model. The problems arising here are exactly the same as at mid-ocean ridges.

Pelitic silt and solid rock fragments can form a porovis-cous mixture in mud volcanoes. Unlike usual volcanoes, phase transformations between the skeleton and flowing fluid do not occur here. Therefore, an isothermal compaction model of a mud volcano involving the appropriate motion of the mud mixture can be used without reservations and limitations associated with the modeling of usual volcanoes [Karakin et al., 2001]. Figure 2 presents a schematic layout of a mud volcano with its feeding roots. A poroviscous mixture forms if the volcano feeder crosses a thick sequence of

Figure 3. Processes at mid-ocean ridges (after Schubert et al. [2001]). (a) Movement pattern of partially molten mantle material in mid-ocean ridge and hotspot regions. Broken lines show the matrix motion. Solid lines show the base of the lithosphere and contours of the melting extent. The melt is separated from the partially molten mantle material in ascending flows. The melt in the central part of the flow concentrates in a magma chamber, whereas the peripheral melt solidifies to form the lithosphere material. (b) Melt velocity in the partial melting zone as a function of the melt viscosity and ascending mantle flow velocity at the asthenosphere base. (c) Porosity (in percent) in the partial melting zone. At / < 2%, pores are isolated and the porosity coincides with the melting extent. At / > 2%, pores form interconnected channels and the porosity is much smaller than the melting extent due to the upward filtration of the melt. The porosity depends in this case on the melt viscosity and ascending mantle flow velocity.

Figure 4. 1-D motion of the suspension and poroviscous medium in a vessel with rigid impermeable free-slip walls. Arrows show movement directions of the phases. The model imitates sedimentation in sedimentary basins and tailings ponds.

clayey rocks. The feeder walls are then washed away and the mud mixture is saturated with solid rock (clay and aleu-rolite) fragments.

Sedimentation of rocks and colloid structures plays an important role in natural processes of the sedimentary basin formation. The sedimentation is accompanied by the compaction of sediments and the transition of the suspension into the poroviscous state. This process is schematically illustrated in Figure 4. Depending on the sedimentation conditions, the suspension and porous medium can be separated by either an inner boundary or an extended transitional zone with a cluster structure. The results of this study enable the formulation of a mathematical model describing the sedimentation process as a whole, which is helpful for planning laboratory experiments in order to verify the main model assumptions and constrain the compaction parameters and boundary conditions.

Apart from natural sediments, sedimentation processes can take place in tailings ponds involved in various technological processes. Problems of this type usually arise in relation to the recovery of both domestic and industrial waste. For example, the application of surfactants and other chem-

ical reagents to the extraction of bitumen is associated with the formation of gigantic tailings ponds contaminating vast areas. The detrimental environmental impact of these ponds sets limitations to the application of such technologies and calls for other ways for extracting bitumen from their deposits. The main difficulty consists in that the colloidal structures forming in the tailings ponds differ little in density from water. For this reason, the compaction of the quasisolid phase in the ponds develops very slowly, during many tens of years, implying that the recultivation of such pond water areas by natural sedimentation is actually unfeasible. Active treatment of these processes requires the knowledge of laws governing the motion and physicochemical transformations of these mixtures. This is the goal of the present study.

Compaction models involve material functions and constants determined from experimental data. An additional difficulty of the problem consists in the fact that the behavior of poroviscous media is controlled by many parameters that cannot be therefore determined by simple 1-D experiments. A system of experiments admitting a reasonable interpretation is required. For this reason, analytical solutions that can be regarded as a model of laboratory experiments are interesting. This condition is met by analytical solutions of type (4.14). In particular, Figure 1 illustrates the simplest experiment described by this type of solution.

Notation

fi = fic Ufiv, total region of definition; fic, compaction region (fic(A) and fic(B) are compaction regions adjoining boundaries of types A and B, respectively);

A/ and £i, porosity interval and a small parameter characterizing the width of the region fic(B);

fiv, regions of the porous medium with isolated pores, suspensions or the transitional cluster zone that are adjacent to the compaction zone;

r, any boundary (exterior or inner) of any region; rex, exterior boundary of the region fi; rint, any interior boundary of the region fi; rmin, inner boundary (between the regions fic and fiv)

at which pores close or open (Гт+П and Г;— are the leading and trailing fronts of this boundary);

Гmax, inner boundary (between fic and ) at which the skeleton of the poroviscous medium forms or breaks down;

, Oij, a/ and , pore pressure, tensor of total stresses,

tensor of effective stresses and tensor of stresses averaged over the skeleton phase in the two-phase region fic;

w s,wf and S, velocities of the skeleton, fluid and filtration;

/, /max and /min, porosity and its critical maximum and minimum values in the compaction region;

Z and n, coefficients of the volume and shear viscosity;

n f

S = -/, coefficient of hydraulic resistance;

n/ , fluid viscosity;

k, permeability coefficient;

A and Bk, rheological constants in expression (2.7);

g = (0, 0, —g), gravity;

p = (1 — /)ps + /pf, density of the two-phase medium; ps and pf, densities of both phases assumed constant;

Ap = ps — pf, density difference between the skeleton and fluid phases;

T, temperature;

c, cf and cs, specific heats of the mixture and individual phases;

A, heat conductivity coefficient; l, latent heat of melting;

$, density of volume heat sources; q, density of mass sources of the fluid;

(x1,x2,x3), Cartesian system of coordinates (the x3 axis is directed downward);

(y1,y2,y3), auxiliary system of curvilinear orthogonal coordinates; t, time;

x' = x — ut, wave argument;

u, velocity of any inner boundary relative to the skeleton of the porous medium directed along the outward normal to the region fic;

F(x), (^(x)) initial value of the porosity (and its normalized value measured from the minimum porosity value); wso1, solenoidal component of the skeleton velocity;

0, scalar potential of the skeleton velocity;

rf“1 and n = , deviatoric stresses and pressure in

a hypothetical incompressible Stokesian fluid, hydrodynamic equations of which are used as one of the representations of compaction equations;

Tfj , wv and nv, deviatoric stresses, velocity and pressure in a real incompressible Stokesian fluid in the region fiv;

1

П“ = p“ - 3af and p“ (Щ = pS - 3

f“ akk ,p-,)

reduced

values of П and from which the hydrostatic pressure of the poroviscous medium in Qc(Qv) is subtracted;

nb = p3 — 3a

f

1

kk and p (П = pv - 3a

b

pv)

reduced

values from which the hydrostatic pressure of the suspension in the region Оc(Оv ) is subtracted;

Hc = ^ Z-, dimensional compaction length; у о*

hQ, dimensionless compaction length; к = h-1, parameter related to the dimensionless com-paction| length;

[ф] | + = ф+ — ф-, ф±, boundary values of a quantity ф at an inner boundary inside and outside the compaction region;

ef , first invariant of the effective stress tensor;

k, vertical unit vector;

smax, maximum velocity of filtration at the inner boundary rmax;

s± and /±, asymptotic values of the filtration velocity and porosity corresponding to the asymptotics at x ^ ±to;

s0 and /0, trivial homogeneous solution (the state of homogeneous fluidization), as well as the asymptotic values of the filtration velocity and porosity at x ^ ±to in the case when these asymptotics coincide.

An auxiliary notation system is introduced in Addendum B to designate dimensional quantities (using the tilde sign or the capitalization of the respective symbols). This system is not used in other sections of the paper.

Addendum A. General Equations of Compaction

Classification of compaction boundary problems. Special boundary conditions divide into three categories (below referred to as types I, II and III) corresponding to three types of compaction boundary problems. (Below, the latter are referred to as the first, second and third problems of compaction.)

Problems of the first type are formulated in an infinite and boundless region with asymptotic conditions specified at infinity. In the general case, these are nonattenuating disturbances at infinity. An important particular case when the deviation of a dynamic quantity (stress, strain, filtration rate etc.) from the state of static or dynamic equilibrium tends toward zero at infinity. This is the simplest class of compaction problems, best studied in both 2-D and 3-D variants.

Boundaries of type II include exterior boundaries rex in bounded (finite or infinite) regions: a rigid impermeable surface, free surface, moving deformable surface and so on). Problems of the second type are also known in the literature, although to a lesser extent than problems of the first type.

Velocities (W) or forces (/) can be specified at boundaries of type II:

^x : w = w° or aijn¿ = /і ;

(A.1)

filtration conditions must also be set at these boundaries. Usually, if an experiment is carried out in a vessel with rigid walls, boundary velocities are set at zero, and the walls are considered impermeable:

p0 and <70, critical values of the pore pressure and first invariant of the effective stress tensor;

a and b, kinetic constants in boundary conditions at inner boundaries;

fmin = £(p0), minimum porosity as a function of the pore pressure;

ai and fi, external surface forces at inner and exterior boundaries;

wm = where n and mare the normal and

tangential unit vectors (the index a = 1, 2 indicates two mutually perpendicular directions on the boundary surface);

^x : w = 0, s„ =0 .

(A.2)

Another typical case is a free surface at which external forces in a two-phase medium and the pore pressure vanish,

Tex : /і = 0, pf =0 .

(A.3)

This type also includes problems with periodic boundaries, as well as problems of a solid body emplaced in a finite or infinite space.

The third type comprises compaction models involving

oij, Kronecker delta

aef = — a

structural-phase transformations of the microscopic structure that lead to qualitative changes in the state of the medium. Such models include interior (fic) and exterior (fiv ) regions separated by either inner boundaries or a transition zone. The inner boundaries can move relative to the medium. This situation is described by boundary problems of types III.A and III.B. Pores close at a boundary of type III.A, thereby precluding filtration. At a boundary of type III.B, the skeleton breaks down with the subsequent formation of a dispersive medium. If the porosity experiences a jump at such a boundary, the remaining quantities obey relations (3.1). General conditions (3.1)-(3.3), as well as some specific conditions, must hold at these boundaries. Conditions (3.6) and (3.7) (or (3.8)) must be satisfied in problems III.A. Problems III.B require the validity of conditions (3.13) and (3.15) (or (3.16)).

In addition, there exist mixed boundary problems including various types of boundary conditions. For example, a mixture moving through the pipe can encounter inner boundaries, whereas boundary conditions of second kind (e.g. (A.2) and the like) can arise on the walls.

Canonical equations of compaction. The original equations (2.1)-(2.4) can be conveniently written in various forms depending on the type of boundary problems and boundary conditions. In particular, a new problem formulation requires the development of appropriate formalism. In this connection, I derive these equations in various representations and introduce several definition. I consider formalism canonical if it reduces the compaction equations to standard hydrodynamic equations of compressible or incompressible fluid. A transformation irreducible to these forms is not canonical. Canonical forms are helpful if the decomposition of hydrodynamic quantities into components in the compaction equations has a physical meaning.

Equations (2.1)-(2.4) can be transformed into a form coinciding with hydrodynamic equations of a compressible Stokes fluid [Landau and Lifshits, 1953]:

iaef

j

ix,'

ipf

ixi

- pgi

ef

aiejf

n

|w|

ix,

+

iw1 2 iw|

ixi

3 ixk

+Ziw|O

+ Z я---Oij

ixk

describing stresses associated with this component. It is connected with the total and effective tensors of stresses through the relations

a, — j — ПО“, + 2n

i 2ф lxidx.

- v2#i.

ef _sol

a ■■ = T- • aij

+ 2n

д’ф тт2фл

ЇЖ - v #ij

(A.5a)

(A.5b)

The substitution of the velocity decomposition into initial equations (A.4) yields the hydrodynamic equations of Stokes for an incompressible fluid in which the value n plays the role of pressure:

|T Iö1 |Tij ix,

|П

ixi

pg,

where divwIö1 = 0

nv’wisö1 = dx“- pgi.

(A.6a)

(A.6b)

Taking the divergence and curl operator div and rot to equation (A.6b) yields the equations for n and rot-W“1 :

2П

д/

-Др5^-

1x3

nv’rot-wIö1 = Дрго^/д ) ;

(A.7a)

(A.7b)

each of these values obeys the Poisson equation. Expression (2.5) can be rewritten in the form

grad

Z + 3n) v2ф

+ nv2wIö1

+

gradpf - pg

(A.4)

grad + nv2wIö1

Z + 3n) v^

+

-Os - (l - /) Дpg

(A.8a)

(A.8b)

Equations of motion can be expressed through velocities (see (2.5)). These equations can also be written in an equivalent form. To do this, I decompose the velocity of a porous medium into potential and solenoidal components and define a new quantity n :

-»і і j- і -» sol

w = gradф + w

divwso1 = 0

п = pf - (c + 3^ v2ф

where 0 is the scalar potential of the skeleton velocity.

The solenoidal component of velocity can be brought into correspondence with the tensor rf,oi

dwf°

ix,

+

|j

ixi

Forms (A.6) and (A.8) are advantageous to various reductions of the initial equations associated with the separation of the pore pressure into hydrostatic and piezometric (or hydrodynamic) components. It is the latter that enters boundary conditions. Note that, in general, hydrostatic components of the fluid, skeleton and mixture are different. Physical conditions of the problem depend on which of the hydrostatic states mentioned above are realized at boundaries or, asymptotically, at infinity. Therefore, the aforementioned separation techniques can also be different.

Accordingly, one can obtain the first, second and third reduced canonical forms of the initial equations. A hydrostatic equilibrium state of the skeleton phase pG is distinguishable in the first canonical form. The second reduced form is characterized by the hydrostatic component of the suspension at a critical value of the fluid phase concentration (porous-

о

öl

n

J

ness) /* enabling structural-phase transitions:

dpG

pf = p“ + pG,

dx,

s

P gf

= pb + pg

dxf=[<1—/>*+/ y I»,

(A.9a)

(A.9b)

p

/

:p + pg

dpg

dx,

f

PJ g,

(A.10)

dp“ x a

= —ds, — Apg ,

dx,

n“ = p“ — (c + 3v) v20

V2n“

. d/

—Apg/

dx3

dp1

dx.

J>

— — ds,

n b = pb — (c + 3v) v20

(A.11c)

V2 nb

. d/

—ApgTT— . dx3

In the particular case /* = 1, the suspension is transformed into the pure fluid,

Equation (A.8a) is transformed in a similar manner to grad

Z + 3v) v20

+ vV2^

= gradp“ + / Apg

The first and second canonical forms correspond to situations when inner boundaries of the respective types A and B arise within the definition region. Below, these forms are mainly examined. The third canonical form should be examined if the poroviscous region is surrounded by a pure fluid.

A nonreduced (initial) form also exists and it is used if the gravity is absent or if the subtraction procedure does not give an advantage. Boundary conditions in some problems are formulated in an explicit form and do not require any preliminary treatment. Examples are rigid permeable or impermeable boundaries, a free surface etc. With these boundaries, the problem can be correctly formulated for the initial form of the compaction equations even in the presence of gravity.

Other media that are not in the state of hydrostatic equilibrium can also exist at boundaries of the viscous consolidation region. Solutions in the compaction region and in the exterior must then be matched. In this case the question of using a specific canonical form remains open and needs further analysis. Note that other, exotic canonical forms are also possible. For example, this is the case if the viscous consolidation region is surrounded by a suspension inhomo-geneous in density and porosity, which gives rise to a very specific canonical form.

On the strength of (A.9), the Darcy law (2.4), relation (A.7a) and the expression for n are transformed, respectively, to the forms

grad

grad

Z + 3v) V20

+ vV2^1

gradp b + Ap(/ — / *) g

Z +3v) V20

+ vV2t3soi

gradp b — Ap(/ — /) g .

(A.12a)

(A.12b)

(A.12c)

Using decomposition (A.9) and (A.10), the following components are recognizable in the total stress tensor:

0,, = 0“ — pGd,j , 0“ = 0f — p“d,, , (A.13a)

0,, = 0,bj — pg , 0bj = 0,jf — pbdij

(A.13b)

Summarizing the above analysis, I give equivalent representations of equations (A.4) and (A.6a) reflecting the specific character of canonical forms:

drso1 ,j

dx,'

dn“

dx,

+ /Apg, , or ;

d0f _ dp“ , . .

---- = -k-----+ / Apg,

dx, dx,

d0“,

S, = /Apg, ,

(A.14a)

dr,?:

dn

dx,

(A.11a)

d0,ef dPb < )

_ dx: + Ap</ — /*)g

dX,

(A.14b)

or ;

b

or

or

|x- = —ds, — Ap<1 — / *)g,

nb = pb — (z + 3v) V20 , V2nb = — Apg^ ,

dx3

drso1 dnb

% = £ — Ap(1 — / )g.

dx,

(A.11b)

ef

dx-b

(A.14c)

or

b

or

A noncanonical transformation can be exemplified as follows. Neglecting phase transitions between the skeleton and fluid (q = 0), material function (2.7) is substituted into equation (A.8b). Using (2.6b), the latter can be written in the form similar to all canonical forms:

grad (A;/-ldivs)—nV2 wso1 =

= Bkf ks + Ap(l — f)g .

(A.15)

div

/kgrad(f (f +div(/w?S))) + ndiv (/kV2wso1) =

+

(A.16)

Bk (~dt + div (/W)) — Apdiv [(1 — /)g/k]

div

/kgrad ( V20

A; 2

— \7 2

/

+ ndiv (/kV2wsoi) =

+

(A.17)

= Bk V20 — Apdiv[(1 — /) g/k]

Density inhomogeneities in compaction mechanics are due to density differences between phases. A significant distinction is related to the first (filtration) terms on the right-hand sides of these equations.

Another example of a noncanonical transformation under the condition q = 0 is given in [Scott, 1988]. Introducing the center of mass velocity v = (1 — /)ws + /wf, continuity equations (2.6) can be written in the form

Using (2.6a), the filtration velocity can be eliminated in this equation which can then be written in the equivalent forms

dt + div (s + /Ws) = 0

divv = 0 .

(A.18)

(A.19)

Similar to wsoi, the velocity component v is solenoidal. However, this case does not reduce to transformation (A.6). I transform (2.5) to new variables:

-grad

Z + 3 n ) divs

+ nV2 (v — s ) =

= — Ss — Ap( 1 — /)g .

(A.20)

The effective stress tensor is then transformed to the form

ef

T J ----- '

The appropriateness of using a concrete representation of the initial equations is dictated by the problem formulation. The viscous stress tensor in representation (A.4) has a “conditioned” form that is advantageous if boundary conditions are expressed through stresses. The tensors of total and effective stresses <7ij and af in representation (A.6) are connected with the viscous stress tensor via relations (A.5) also containing “redundant” terms. Consequently, the viscous stress tensor has an “irregular” form. However, this representation has its own advantage.

The initial equations in representation (A.6) split into two relatively independent groups of equations. Equation (A.6) describes the fluid motion involving the curl component of velocity. Equation (A.6) contains all values describing the filtration and volume expansion, i.e., divws and . Both equations are interrelated through the n value and boundary conditions. As a result, the actual flow can be represented as the superposition of these two motions. On the one hand, a certain fictive incompressible fluid moves in accordance with equation (A.6). On the other hand, it is superimposed by such processes as filtration and volume compression of the skeleton described by equations (A.6) and the Darcy law (2.4). This representation is convenient if the flow region is much larger than the compaction region (such is the situation, for example, with the flow in the mantle). In this case, the problem of incompressible fluid dynamics can first be solved, after which one has to separately solve the compaction problem with a known distribution of wso1.

Equations (A.4) and (A.6) resemble convection equations compressible and incompressible fluids, respectively. Density differences can be produced by either thermal expansion (thermal convection) or solutes (chemical convection).

(A.21)

Although form (A.20)-(A.21) is not equivalent to any hydrodynamic equations, it is convenient in numerical computations if boundary conditions are expressed through velocities. First, the continuity equation (A.18) alone is evolutionary. Second, two, rather than three (as in the case of (A.6)), values v and s are to be found. Moreover, the sought-for values have clear physical meaning and are measurable quantities.

Addendum B. Transition to Dimensionless Variables

Application of the n theorem to compaction problems. Compaction problem are multifactor and mutliparameter problems. Therefore, the transition to dimensionless variables and the small-parameter simplification are nontrivial procedures. They are specific to each concrete problem. As an example, I consider the transition to dimensionless variables for a nonstationary problem IIIA with an initial condition for the porosity. Other problems are transformed in a similar way.

For convenience, I introduce an auxiliary system of notation used in this section alone. To distinguish between dimensional and dimensionless variables, I use upper- and lower-case Latin letters. If the use of upper-case letters is inappropriate or impossible (because they are already in use), the tilde sign is utilized to indicate dimensional values.

Models of partially molten rocks usually incorporate phase

transformations between the skeleton and fluid. It is then reasonable to assume that the density difference between the skeleton and melt can be neglected in the continuity equations, whereas the buoyancy force cannot be neglected in the equations of motion. This assumption is an analogue of the Boussinesq approximation in the equations of thermal convection and somewhat simplifies the dimensionless equations without loss of their generality.

Let H(t) be the vertical size of the definition region that does not depend on time. The solution depends on the size of the initial region H(0), and the maximum porosity at the initial moment Fmax can be taken as a porosity scale factor. Note that some dimensional values enter equations only in certain combinations (i.e. they are not independent). Therefore, I consider these stable combinations as new dimensional values and, moreover, as scale factors. The scale factors will be marked by asterisks:

/* = Fm

C* =

A;

Fm ax

C* =

Bk

F k

ma

X.

н(о): ^

C*po/,

t=I,

t *

l

cH (О/

l

min

FmaxAl

Др^Н ^)A; ’

H 2(О)

~hT

where

f — Fmax/ (x,t) , /max — Fmaxfmax,

siHcFmax wiHcFmax

Si =------------ Wi =---------------

t *

t *

C ( / ) + 3C (/) = C* (Z(/) + 3n(/)

c(/J = C*O(/) ,

F = Fmax^ (x) + /m

M =

M

Pf FmaxH3

^ (x) max = 1 PO = ДP»HcP0 ,

¡Cy = ДpgHca“j - P%■ =

= ДpgHcabj - PgO, = Дp»Hcaij

(B.3)

Ф

фНс

Дp»HcP, /m

/m

Fm

In the case of nonlinear rheology, these definitions remain the same with the additional condition: l = 0 and С * =

Ao=(Co+3 c).

The boundary problem of type IIIA depends on eight independent dimensional constitutive parameters (t, X a, p0, С*, <5*, Дрд and H(0)) and on the dimensionless parameter Fmax. These parameters involve three independent dimensions: time, length and mass. According to the П theorem, these ten dimensional parameters can be expressed through seven independent dimensionless parameters that can be chosen, in general, arbitrarily. Moreover, all other dimensional and dimensionless values are expressed through these seven independent constitutive parameters. The П theorem states that the number of independent constitutive parameters is exactly seven.

In accordance with the requirements of the П theorem, the following values are chosen as independent dimensionless constitutive parameters:

Using these relations, equations of motion in velocities A.5), continuity (2.1b, c) and coupling (A.3) assume the form

= K2Os - (1 - /*/)k

^ + div (s + /ws) = q ,

- -|t + div [(1 - /) w8] = -q

(B.4a)

(B.4b)

ws = gradф + wso1, divwso1 = О , (B.4c)

v2n = -Др»|/

1x3

(B.4d)

nV’rot-wso1 = —/*rot (/k)

(B.1)

where k is the unit vector directed vertically upward. Equations of motion in velocities (B.4) form a complete system that can be integrated under appropriate boundary conditions. If boundary conditions include stresses, the problem should be complemented by the Terzaghi and rheological equations:

F

Fm

Др»н (О)!

Hc =

Z*k*

П/ *

(B.2)

ef f

aii = ai/ - pf °ii

Here Hc is the dimensional compaction length. With such a choice of dimensionless variables, the value h0 = к-1 is the dimensionless length of compaction. The boundary problem is fully described in terms of the independent dimensionless constitutive parameters. The remaining values are expressed through both dimensional and dimensionless constitutive parameters. Now I express the sought-for dimensionless quantities (porosity / (x, t), filtration velocity S, fluid mass м etc.) through the independent values:

ef

T J ------ '

iw.

ix

so1 iwSo1\ ( 4 \ 2

+ (Z + 3n) V2^- .

+

ix.

The procedure of transition to dimensionless variables is quite similar for boundary problems of other types. I present the reduced expressions for the first and second canonical forms, each associated with the specific type of material functions:

*

x=

a=

2

к=

la“.

Ос : ^ = -/*/k, a“,. = a?/ - p“^

gradp“ = — + k ,

п“ = p“ - (z + 3 n) v^

V2 “

п

-/* /

/ ix3

Z + 3n = /i, O = /k ;

i *

О-: if=/*(1 -^

^ = aie/ - Pb°«,

gradp* = —к2Os + (1 — /*) fe

(B.5a)

(B.5b)

(B.5c)

(B.6a)

the equation

к2so = /ok, /o >/m

(B.Wb)

Relations (B.9c) and (B.10a) imply that the solenoidal component vanishes,

Ос : п“

О

and equation (B.9a) is simplified, assuming the form corresponding to the particular case of the potential field of skeleton velocities:

Ос : grad (/ 1divs) = к2/ ks — k .

(B.11)

Accordingly, boundary conditions (3.3)-(3.9) in problem IIIA are transformed to the following dimensionless form:

г

(-) .

Sn = О, / = £(О)

+

— О, ф,п — О ,

(B.12a)

п* = P* - (z + 4n) v2ф

V’п* = -/* / ,

ix3

4

Z + 3n = (1 - /)l , n = (i - /)l, O = (i - /)k .

(B.6b)

(B.6c)

i/ , • д/ t s

—+divs = q, -— + divw = -q it it

V2n = 0, nV2rotwso1 = 0 .

In this work, only a trivial asymptotics of the type

2

p“ ^ О, к2.? = ks0, ws ^ О / ^ /0, at |x| ^ то

гП+П : / = too) ,

(p“ - po) в (p“ - po) = au, a > О, po > О

rex : Sn — О, s — О, ф,п — О .

(B.12b)

(B.12c)

The initial condition is the same for problems of any type:

/ = p (X) + /min at t = 0 . (B.7)

The small-porosity approximation is often used in prob-

lems of the first and second types,

/* « 1 , (B.8)

which simplifies some expressions, in particular, (B.4a, b, d):

Boundary conditions (3.15) or (3.16) in a problem IIIB are transformed in a similar way. For example, condition (3.16) is now written as

rmax : / = /m

ef

i,

О

bu .

(B.13a)

(B.13b)

(B.9a)

(B.9b)

The 1-D case. Spatially one-dimensional models are most often encountered in applications. The finite region Oc(i) is a segment with moving boundaries x- = x_(i) < x < x+ = x+ (t) of the length h(t) = x+ — x— Here x+ and x- are boundary functions for the leading (r+,) and trailing (r- ) fronts. In this case, the above boundary conditions are further simplified. Setting s = s3 and x = x3, equations of motion and continuity reduce to expressions similar to all canonical forms:

(B.9c)

¥ = - - і)

д/ is „

— А------—О .

it ix

(B.14)

(B.15)

(B.Ma)

is used in the problem of the first type; so and /0 in (B.10a) represent the trivial solution of the initial system of equations. This solution is called homogeneous fluidization in the Western literature. These values are interrelated through

On the strength of these simplifications, expressions for stresses (B.5a) and (B.6a) admit, after some transformations, the first integrals

033 = 03f — p3 = Const, 033 = = Const .

In the 1-D case, the medium is in the hydrostatic equilibrium state outside the compaction region. For this reason,

+

the boundary conditions for stresses (B.12a) and (B.13a) lead to the expressions

P“ = af,

a3e3f - p* = a3*

(B.16)

providing additional constraints on the equation of motion and boundary conditions. For example, initial condition (B.7) assumes the form

/(x, 0) = p(x) + /m

(B.17)

The equations of motion for problems I, II and IIIA and the expressions for the total stress tensor and piezometric component of the pore pressure are

Ос : |x(/-1?

ix V ix

2 -k к / s - 1

a33 = — psgx + Const

P“

-i li ix .

(B.18a)

(B.18b)

(B.18c)

(B.19a)

д/ is i/ws

— +-------1----- = 0

it ix ix

is iws о ix + ix

(B.23b)

(B.23c)

and boundary conditions (3.3) and (3.16) implying the presence of surface forces are used. Substituting the expressions for the volume viscosity (2.11) into the effective stress tensor in (3.16) yields

Г max

dx 1 ( ) is ¡

dt = b (2 - ^ix^

/ = /m

. l is ix ,

* ^ * *

* О, a3*3 = a*

(B.24)

The Smax value in this boundary condition is not an a priori given constant but is determined in the process of solution. Therefore, a correct formulation of the problem requires an additional condition determining this value. In some problems of type IIIB, it is convenient to normalize the maximum porosity value:

Asymptotic values in an infinite region of problem I can coincide (and in this case correspond to (B.10a)) or differ along opposite directions of the x axis:

/max = 1 .

(B.25)

In the absence of surface forces, we have condition (3.15) transformed to the form

(l - / )l ix—o

(B.26)

(B.19b)

where the values s± and f± are connected through the relations similar to (B.10b)

At l = 0 this relation yields

|i =o.

ix

(B.27)

к^± = /± .

(B.20)

Boundary conditions in problem IIIA are the same at both boundaries:

Addendum C. Transition to Dimensionless Variables

г(+) , г(-) : mi^ min :

/ = /min ,

s3 = 0 at x = x-(t) and x = x+(t) .

(B.21)

An additional condition is set at the leading front:

Г

(+)

dx+

dt

if /-; is _i_ A a ( /-; is

-в У ix + pv в l-/ ix - po

(B.22)

Problems IIIA. I analyze wave equation (5.4)

dx (/-i dx )=(/ - /o) /-k - u

dx

(C.1)

in finite and infinite regions. A finite region implies problem IIIA with boundary conditions (5.5):

г(+) , г(-) : min , min :

f=fm

£(o)

(C.2a)

Now I address the 1-D formulation of problem IIIB in the dimensionless form incorporating material functions (2.11). The porosity in this problem cannot be considered negligibly small, although it cannot exceed a critical value of about 20-30%. Therefore, the buoyancy term in the equation of motion depends on the porosity; also, the initial continuity equations (B.4b)

Ос

ix (і—- /)ia =

= к2 (1 - /)k s - (1 - /*/)

(B.23a)

Г

( + ^ u/min / = -eu - P0 .

(C.2b)

The problem IIIA involves a constant value of the minimum porosity /min at which pores close and open. If the porosity at the initial moment satisfies the condition f (x, 0) > fmin all along the x axis, a disturbance can assume with time the form of a solitary wave in an infinite region. In this case, the asymptotic conditions of type (B.19a) are set.

An important parameter is the wave mass defined as the

s

s

max

difference between the extreme and boundary points. The wave mass in problem IIIA is given by the expression

M = 2 / (f — fmin)

f f

v7©/)

In problem I the mass is defined as

In this case

2

R(f ) = - (f — fo)2 , f = (2 + fo) fo

(C.3)

fo

2fo

2 — f0

f

o

2fo 2 + fo

and quadrature (C.8) reduces to the expression

(C.10)

M = 2J (f — fo)

fo

f 'df

v7©/) .

(C.4)

f = ©(f ) — 21 .

where

1—f

R(f )

R(/) = / — (2 — /0) / + /0.

The solution can be obtained with the use of the quadrature of this expression:

<f ) = ±

df

f

V7©/) '

R(f ) > 0

f(1) = -(2 — fo — VÂ) , f(2) = -(2 — fo + VA)

o

A = 0, f(1) = f(2) = 4(2 — f°) = fo

— f o)2 — f = 0. -

x(f ) = ±\/ô

fdf

2V (f — fo) V(1 — f ) f

I examine wave equation (C.1) in the particular case l = 0 and k = 3:

calculated in elementary functions:

±

(C.5)

This equation has the first integral I (the constant of integration)

2v/1—/ +

/1 — fo + / — f •n/1 — fo V1 — fo — V1 — f

fo

(C.11)

ln

(C.6)

The wave mass is determined by formula (C.4):

In the particular case, when the maximum value of f at the upper point of the wave is equal to unity (fmax = 1, fo< 1), the integration constant is fixed,

M = 2 J (f — fo)

fo

df

v7©/) '

(C.7)

Using (C.7) and (C.10), this expression is transformed to

fdf

fo

V(1 — f)

and can be calculated in elementary functions:

= 3/T2-/1 — fo (2 + 2fo — fo) -

(C.12)

(C.8)

The necessary condition for the existence of this solution is the inequality

(C.9)

Expressions (C.11) and (C.12) provide exhaustive information on the solution in an infinite region. As seen from (C.12), the wave mass is always positive and the wave velocity monotonically increases with the mass.

In this connection, I examine the solution in a finite region of definition, when the roots f(1) and f(2) are different. Given fo = fmin < 1, the boundary condition at the leading front (C.2b) and the expressions for the roots f(1) and f(2) take the form

implying that the radicand in (C.8) is positive. I analyze this inequality in more detail. Let /and /(2) be the roots of the equation R(/) = 0 and let A be its discriminant,

(1 — fmin)

fm

min I------+ 1 — 1

f

= I® + «

8fm

A= (2 — fo)2 — f .

In the case of an infinite region, both roots coincide (A = 0), and minimum values of the porosity and fluid flow velocity coincide with their values in the trivial solution:

A — (2 — fmin) —

f

f(1) = f (2 — fmin — VA) ,

f(2) = f (2 — fmin + VA) .

(C.13)

The above analysis is incomplete. Expressions (C.11) and (C.13) imply that wave solution in the respective regions exist. However, the conditions of the existence of solutions in finite and infinite regions are somewhat unclear.

For this purpose, it is convenient to qualitatively examine the solutions of wave equations in the phase plane by

1

f

1

1

1

1

2

M

1

2

the method of dynamic systems (/, / ). Equation (C.5) is transformed into the dynamic system

df

G(/)

f

(C.14)

Phase portraits in the plane (/, / ) are constructed using a modified method of dynamic systems based on the analysis of stationary points of the first integral I. For this purpose, expression (C.6) is transformed to the form

I

(/ — /"0 (/ — /“) -

Stationary points of the integral I in the phase plane (f, f ) are the points on the abscissa axis that are the roots fi of the equation

G' = 0

f3 — ff + ffo = 0 .

(C.15)

D :

(fo

4

f

27

Figure 5. Plot of the function G'(f).

Stationary points of the integral I coincide with the singular points of dynamic system (C.14). To constrain the type of stationary points is easier compared to singular points of a dynamic system. A phase portrait in this plane is constructed depending on type of these points. The type of stationary points depends on the sign of the second derivative G"(fi) at a given point. Given G"(fi) > 0, this point is a center, and it is a saddle at G"(fi) < 0.

Equation (C.15) has the determinant D

Figure 6. Phase portrait of dynamic system (C.14) under condition that singular points coincide and form a cusp-type singular point.

(C.15), I consider the plot G'(/) shown in Figure 5. Under the condition

D = 0 or f

27 (fo

(C.16)

4

Depending on the sign of D, various situations are possible. In order to determine the properties of roots of equation

this curve is tangent to the abscissa axis at the point 3

/ = —/ . Then equation (C.15) has three real roots / of which one is negative and two positive roots coincide:

3

/1 = —3/0, /2 = /3 = -/0 .

The positive roots /2 and /3 are cusps in the phase portrait (Figure 6).

D > 0 or f <

27 (fo 4

This curve lies off the abscissa axis under the condition In this case, equation (C.15) in the physical region does not have roots but has only one negative nonphysical root /1 :

fi

ff

2

+ '

(fo

f

27

+

+

2

fo

4

f

27

<0

If D < 0, the curve G7(/) intersects the abscissa axis, and equation (C.15) has three real roots / one of which /1 is negative and the other two are positive, 0 < /2 <

or

2

2

2

f

2

2

2

f

4

2

2

f

2/min /.

F(Yo,r) — F(Y,r)

Large waves are above the points and have a complex profile including two inflection points.

Problems IIIB. Here I analyze the wave solutions of (5.10)-(5.11),

d 2f

1

df = (1 — f)(f — 1 + f o)—f

I=f(‘—/) f =

dx

(1 — f )(f — 1 + fo ) — -

(C.17a)

(C.17b)

Similar to the preceding case, equations (C.17) are transformed to the dynamic systems

Figure 7. Phase portrait of dynamic system (C.14). Two different singular points, saddle (to the left) and center (to the left), lie on the abscissa axis. “Large” and “small” trajectories are distinguishable. Their respective areas are differently shaded. Dotted shading indicates periodical solutions. The separatrix coming to and going away from the saddle corresponds to the solution in an infinite region.

3

-/0 < /3. Figure 5 shows that G"(/3) > 0 and G"(/2) < 0.

Consequently, the point (/(3), 0) is a center, and the point (/(2), 0) is a saddle. The phase portrait for this case is shown in Figure 7.

If the solution in a finite region exists, quadrature (C.8) can be calculated and is represented through elliptic functions [Karakin, 1986]:

1

df

df

^ — f)(f — 1 + fo) — -

(C.18a)

df

df

(1 — f )(f —1+fo)—1

f

1—f

(C.18b)

Moreover these equations have the first integrals I

i=1 /2 — 3 (1—f )3 +

3

+ 2 (1—/ o)(1—/ )2j1f

(C.19a)

2f2—4 (1—f )2+^o (1—f )—2f.

(C.19b)

In the absence of surface forces, boundary condition (B.26)

A

F(Y,r)— 2E(y,t) — — F(yo, r) — 2E(yo, r) —

A

2 ^Yo4 = 2B2Yo2 + A4 + 2^Y4 + 2B2Y2 + A4

Yo (A2 + Y02)

Y (A2 + Y2)

where

Y=

1—f

Y0=

1 — fm

A

VA2 — B2 V2A ,

fm2 in min

qfmin — pfm

Yo =

+1

fmin (fminp 2)

2 (qfmin — Pfmin + 1) Y2 - A2

ccos Y^^

Y :

1

P = (2 — /min) q, q = 2/minW .

The profile of the wave solution is shown in Figure 8. Two types of traveling waves, small and large, exist. Small waves are below both singular points and have a simple profile.

Figure 8. Profiles of solitary waves corresponding to the phase portrait in Figure 7. “Large” and “small” trajectories are shown. The former have a more complex profile.

I

1

f

B2 =

r=

Yo2 - A2

Yo2 + A2

determines the first integral constant: I = 0. This is sufficient for constraining the single solution. If surface forces are present, it is convenient to use boundary condition (B.24) in its explicit form (5.13):

(C.20)

(C.20b)

The porosity derivative is finite in the first case and tends to infinity in the second case.

Both dynamic systems (C.18) have two singular points /1 and /2 on the / axis of the phase plane (/, / ):

/o

2

f2

fo

1 — T +

Figure 10. Phase portrait characterizing the situation described by equation (C.17b), with the inner boundary rmax separating the poroviscous region and suspension at / = /max. The wave velocity is positive. Two singular points (saddle and center) exist in this case. The circular trajectory tangent to the right-hand boundary rmax characterizes the case when no external surface forces are present at this boundary. Trajectories filling the space between the circular and limiting trajectories (boundary points of the latter are shown by the broken line) also correspond to the absence of surface forces at the boundary. Trajectories filling the space between the limiting trajectory and separatrices correspond to the case when surface forces are present at the boundary r

max •

under the condition

2

o

2

1

> -f

f < 4 (1 — fi

22 max) Sm

Figure 9. Phase portrait corresponding to the situation described by equation (C.17a), with the inner boundary rmax separating the poroviscous region and suspension at / = /max. The wave velocity is positive. Two singular points (saddle and center) exist in this case. The circular trajectory tangent to the right-hand boundary rmax corresponds to the case when no surface forces are present at this boundary. Trajectories filling the space between this circular trajectory and separatrices correspond to the case when external surface forces are present at the boundary rmax.

It is easy to show that at f > 0 these points lie in the physical region of the right-hand half-plane (1 > f2 > f1 > 0) in both cases, with f1 being a saddle and f2, a center. Phase portraits are shown in Figures 9 and 10 for both cases. In the second case, phase trajectories converge to the vertical axis f = 1. As seen from these figures, boundary condition (B.26) (valid in the absence of surface forces as well) can only be met for equation (C.17a). In the second case there is no solution compatible with (B.26).

To discuss these situations in more detail, I address a solution of (C.17a) with a positive wave velocity f and with boundary condition (5.12a) implying the presence of a boundary rmax free of surface forces. This solution corresponds to the circular trajectory of the phase portrait in

or

The length of this wave is given by the expression

+

+ 2 (1 — fo )(1 — f )2—

(1 — f )

(C.21)

=0.

f3

h .

(C.22)

f3

h=2

f

V©

All trajectories in Figure 9 that lie in the upper half-plane above the aforementioned specific trajectory but below the separatrix intersect the line / = 1 at a finite angle. This corresponds to boundary condition (B.24). All these trajectories represent the solution in the form a traveling solitary wave. Trajectories lying above the upper separatrix and under the lower one represent certain wave solutions in a different form involving a single front, either leading or trailing. These waveforms are not the solutions proper but can be constituents of more complex nonstationary solutions. Trajectories lying within the aforementioned specific trajectory have no physical meaning and cannot represent or be a part of solution because they do not reach the inner boundary.

Now I consider trajectories provided by expression (C.19b). They also include a specific circular trajectory described by the algebraic equation

1

1

(C.23)

Figure 11. Phase portrait characteristic of the situation with a negative velocity of the wave described by equation (C.17a). Only one singular point of the saddle type exists in this case. Wave solutions exist only in the presence of external surface forces.

Figure 9 that is tangent to the vertical line / = 1. The substitution of (C.19) into boundary condition (B.26) yields I = 0 or

The solution of this equation gives the wave process in the form of a traveling wave and reduces to the quadrature

where ±— is the position of the leading or trailing front of the boundary rmax and

©=3 (1 — f )3—(1 — f o)(1 — f )2+2

(1 — f )

Figure 12. Phase portrait characteristic of the situation with a negative velocity of the wave described by equation (C.17b). This case is also characterized by a single singular point. Wave solutions exist only in the presence of external surface forces.

This trajectory intersects the line f = 1 at a finite angle. Trajectories lying between the specific trajectory and separatrices represent solutions in the form of traveling waves with boundary conditions (5.13) including concentrated surface forces. Trajectories lying within the specific trajectory are nonphysical, and those outside the separatrices can be, as in the preceding case, parts of nonstationary solutions.

Now I address the model variant with a negative wave velocity (u < 0). It is easy to show that in this case f2 > 1. Two possibilities are possible here. With (1 — fmax) smax — u > 1, only one root fi is present in the physical region (1 > f 1 > 0) and it is a saddle in both cases considered. The phase portraits of both dynamic systems (C.18a, b) are shown in Figures 11 and 12. In the second case, trajectories also concentrate toward the vertical line f = 1 and both situations do not admit solutions meeting boundary conditions (B.26). Thus, if no concentrated forces are present on the inner boundary rmax, the latter cannot move in the negative direction.

The above analysis of phase portraits concerned only such situations when the inner boundary rmax separates a poro-viscous medium from a suspension. The presence of a transitional zone having a cluster structure between these two states requires a special investigation beyond the scope of this paper.

References

Ahern, J. L., and D. L. Turcotte, Magma migration beneath an

ocean ridge, Earth Planet. Sci. Lett., 45, 115—122, 1979.

Barcilon, V., and F. M. Richter, Nonlinear waves in compacting media, J. Fluid Mech., 164, 429—448, 1986.

Barcilon, V., and O. M. Lovera, Solitary waves in magma dynamics, J. Fluid Mech., 204, 121-133, 1989.

Cordery, M. J., and J. P. Morgan, Convection and melting at mid-ocean ridges, J. Geophys. Res., 98, (B11) 19,477-19,503, 1993.

Cox, K. G., D. P. McKenzie, and R. S. White (Eds.), Melting and Melt Movement in the Earth: Proc. R. Soc. Discussion Meeting Held on 3 and 4 March 1992, Univ. Press, Oxford, 1993.

Gordin, V. A., and A. V. Karakin, Solitary wave solutions of viscous consolidation problems (in Russian), Mat. Mode., (2), 98-116, 1990.

Karakin, A. V., On main equations of the melting ice mechanics, in Ice Physics and Techniques (in Russian), pp. 87-97, SO AN SSSR, Yakutsk, 1974.

Karakin, A. V., The Study of Geodynamic Processes and Structure of Heterogeneous Media in the Upper Mantle. Doctoral Dissertation (in Russian), IFZ AN SSSR, Moscow, 1986.

Karakin, A. V., A thermomechanical model of the magma filtration in an asthenospheric layer (in Russian), Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, (1), 21-27, 1988.

Karakin, A. V., Fluid dynamics models of the crust with an inelastic skeleton (in Russian), Fiz. Zemli, (2), 3-15, 1990.

Karakin, A. V., General theory of the low-porosity compaction (in Russian), Fiz. Zemli, (12), 13-26, 1999.

Karakin, A. V., and L. I. Lobkovskii, Mechanics of a two-phase porous viscous medium and some geophysical applications (in Russian), Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, (6), 53-63, 1979.

Karakin, A. V., and L. I. Lobkovskii, Mechanics of porous two-phase visco-deformed medium and its geophysical applications, Lett. Appl. Sci., 17, 797-805, 1979.

Karakin, A. V., and L. I. Lobkovskii, Hydrodynamics and Struc-

ture of a two-phase asthenosphere (in Russian), Dokl. Akad. Nauk SSSR, 268, 324-329, 1982.

Karakin, A. V., and V. E. Shklover, Averaging a granular medium with a ball packing (in Russian), Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, (3), 74-81, 1990.

Karakin, A. V., and S. Yu. Levitan, Modeling of fluid-dynamic processes in rocks with a viscous skeleton, in Mathematical Modeling of Geological Processes (in Russian), pp. 17-32, VNIIGeosystems, Moscow, 1993.

Karakin, A. V., S. A. Karakin, and G. N. Kambarova, Movement of a mud mixture through a mud volcano channel, Izvestiya, Phys. Solid Earth, 37, 812-824, 2001.

Khodakovskii, G., M. Rabinovich, G. Ceuleneer, and V. Trubitsyn, Melt percolation in a partially molten mantle mush: Effect of a variable viscosity, Earth Planet. Sci. Lett., 134, 267-281, 1995.

Landau, L. D. and E. N. Lifshits, Mechanics of Continua. Part 1, Hydrodynamics (in Russian), Nauka, Moscow, 1953.

Maal\e, S., and A. Scheie, The permeability controlled accumulation of primary magma, Contr. Miner. Petrol., 81, 350-357, 1982.

MacKinnon, M. D., Development of the tailings pond at Syn-

crude’s oil sands plant: 1978-1987, AOSTRA J. Res., 5, 109113, 1989.

McKenzie, D., The generation and compaction of partially molten rock, J. Petrol., 25, 713-765, 1984.

Ribe, N. M., The deformation and compaction of partial melt zone, J. R. Astron. Soc., 83, 487-501, 1985.

Richardson, C. N., Two-dimensional flow model for the process simulation of complex shape composite laminates: Melt flow in a variable viscosity matrix, Geophys. Res. Lett., 25, 10991102, 1998.

Richter, F. M., and D. McKenzie, Dynamical models for melt segregation from deformable matrix, J. Geol., 92, 729-740, 1984.

Schubert, G., D. T. Turcotte, and P. Olson, Mantle convection in the Earth and Planets, Univ. Press, Cambridge, 2001.

Scott, D., The competition between percolation and circulation in a deformable porous matrix, J. Geophys. Res., 93, 6451-6462, 1988.

Scott, D., and D. Stevenson, Magma solutions, Geophys. Res. Lett., 11, 1161-1164, 1984.

Scott, D., and D. Stevenson, Magma ascent by porous flow, J. Geophys. Res., 91, 9283-9286, 1986.

Sleep, N. H., Sensitivity of heat flow and gravity to the mechanism of sea floor spreading, J. Geophys. Res., 74, 542-549, 1974.

Spiegelman, M., Flow in deformable porous media, Part I, Simple analysis, Part II, Numerical analysis - the relationship between shock waves and solitary waves, J. Fluid Mech., 247, 17-63, 1993a.

Spiegelman, M., Physics of melt extraction: Theory, implications and applications, Trans. Roy. Soc. London., 342, 23-41, 1993b.

Spiegelman, M., and D. P. McKenzie, Simple 2-D models for melt extraction at mid ocean ridges and island areas, Earth. Planet. Sci. Lett., 89, 137-152, 1987.

Taylor, G. I., The coefficients of viscosity for an incompressible fluid containing air bubbles, Proc. R. Soc. London, A226, 34-39, 1954.

Trubitsyn, V. P., G. I. Khodakovskii, and M. Rabinovich, Wavelike migration of a melt in a porous constant-viscosity medium (in Russian), Fiz. Zemli, (10), 33-39, 1998.

Trubitsyn, V. P., A. V. Karakin, and S. A. Karakin, The wave regime in a magma floating up under mid-ocean ridges, Izvestiya, Phys. Solid Earth, (3), 2003.

Walker, D., E. M. Stolper, and J. F. Hays, A numerical treatment of melt/solid segregation: Size of eucrite parent body and stability of terrestrial low velocity zone, J. Geophys. Res., 83, 1533-1535, 1978.

Wiggins, C., and M. Spiegelman, Magma migration and magmatic solitary waves in 3-D, Geophys. Res. Lett., 10, 1289-1292, 1995.

(Received 17 October 2002)