A. V. Karakin
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 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 rheological 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 poroviscous 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  on the basis of general physical principles. The expression for the volume viscosity of a two-phase gas-water system was obtained by Taylor . 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  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 skeleton viscosity [Ahern and Turcotte, 1979; Maaloe 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  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  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  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 free-slip 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 . 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 .
Richardson  numerically examined the 2-D motion of a magmon in the presence of preset homogeneous deviatoric 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 10 18 to 10 10 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 [Barcilon and Lovera, 1989; Barcilon and Richter, 1986]. The method of dynamic systems was applied by Gordin and Karakin  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 poroviscous 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.
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:
The following notation is adopted here: sij and ssij, tensors of stresses in a two-phase medium and stresses averaged over the skeleton phase; pf, pore pressure; ws and wf, velocities of skeleton and fluid; g = (0,0,- g ), gravity; r = (1 - f) rs + frf, where rs and rf are the densities of both phases assumed constant; T, temperature; l, heat conductivity coefficient; l, latent heat of melting; F, 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 rc = rf cf + rs cs, where cf and cs are specific heats of the phases. The vertical x3 axis of the Cartesian coordinate system (x1, 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 W bounded by a surface G. 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 sefij depending on grad ws:
This relation was proposed by Terzaghi and is sometimes named the Terzaghi relation. The constitutive equation of the skeleton is
where z and h are the coefficients of volume and shear viscosity. The Darcy law holds true in the fluid phase,
Here, s = f( wf - ws) is the fluid flow velocity, hf is the fluid viscosity and k is the permeability. For the sake of convenience, I introduce the coefficient of hydraulic resistance d hfk. In the general case, z, h and d are functions 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:
where Dr rs - rf 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
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 poroviscous medium ( q = 0), the continuity equations can also be written in a more convenient form:
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 G int. 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: Gmin (type A) and Gmax (type B). Pores close at a Gmin boundary ( f = fmin ). Given a maximum porosity ( f = fmax ), the skeleton breaks down at a Gmax boundary. The entire two-phase region including inner boundaries of the above types divides into two regions, W = Wc UWv: the region of the compaction proper Wc and the contiguous region Wv. The latter includes a suspension, a poroviscous 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 ( fmin, fmax ) 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 Gmin or Gmax, respectively. Regions adjacent to these boundaries are denoted as Wc(A) and Wc(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:
where Al 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 sefij 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 poroviscous medium and suspension are separated by transitional cluster zones rather than an inner boundary Gmax. 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 integrodifferential 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 spacings 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,
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 Gmax 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
Relations (2.8) and (2.9) hold only at the inner boundary Gmax 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 Wc(B) accounts for a small porosity interval Df of the entire porosity range ( fmin, fmax ). The skeleton in this interval weakens prior to its breakdown at the upper boundary Gmax. The width of the region Wc(B) is characterized in this interval by a small parameter e1,
However, this does not mean that the region Wc(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 Gmax cannot be considered small. However, relations in the Wc(B) vicinity of this boundary can be simplified in a different way. If an inner boundary Gmax 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 Wc(B) and use conditions (2.8) or (2.9). Rejecting higher-order terms of the expansion, I obtain
Below I address two situations. In the first case, the volume viscosity remains finite at l = 0 and k = 1 (condition (2.8)):
In the second case l = 1 and k = 1 (condition (2.9)) and all material functions vanish at the boundary:
Dependences (2.12) and (2.13) of material functions are simplest in the region Wc(B) and qualitatively characterizes its specific properties. Other expressions describing these functions are also conceivable.
General boundary conditions. First I consider boundary conditions at inner boundaries G int 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 si (of molecular, electrostatic or other origin) is applied to the inner boundary. The stress vector at the boundary is then also discontinuous:
where [ j]|+- j+ - j-, with j being boundary values of the quantity j at the inner boundary inside and outside the compaction region; sn si ni and wsm = wsi m(a)i, 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; sij 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
If the porosity is continuous and no concentrated forces are present, these conditions are simplified,
I introduce a curvilinear orthogonal system of coordinates ( y1, y2, y3 ) such that the coordinate surface y3 = Const coincides with the inner boundary considered. The continuity conditions for the tangential stresses and velocities at this boundary yield the relations
The initial condition for the porosity in all nonstationary problems is
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 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:
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,
The compaction region Wc 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 Gmin at a minimum value of porosity. The closure of pores means that this boundary is impermeable. Therefore, (3.3a) yields
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
In accordance with the second relation in (3.2), one obtains
The substitution of the time derivative from (2.6a) gives
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 Gmin 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 G(+)min [Karakin, 1999]. If the leading front G(+)min 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
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 fmin:
The function x(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 fmin value.
If pa < p0, the boundary does not move until the porosity reaches the value x(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 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):
Here sef is the first invariant of the effective stress tensor. The critical parameter s0 plays the same role as p0 and satisfies a relation similar to (3.9).
Closure of pores at the trailing front G(-)min results in the energy release and dissipation in space. Condition (3.6) is satisfied at this boundary at fmin = x(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
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 sef ) at this boundary is first determined. If it exceeds the critical value p0 (or s0 ), the leading front conditions are chosen. If it is smaller than the critical value, the boundary value of porosity f is found. If the porosity is fmin = x(0), the trailing front condition (3.7) is chosen. If f > fmin = x(0), the boundary does not move. Finally, if f < fmin = x(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 fmin = x(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 fmax, and the pore pressure is independent,
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
The effective stress tensor vanishes in the suspension. Therefore, the following relation holds in the poroviscous medium:
According to the above, if surface forces are absent, the inner boundary Gmax exists during the suspension consolidation and skeleton formation, i.e. at the leading front. At this boundary
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: div ws , and div s . 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 Gmax. Therefore, the breakdown processes are subjected to additional conditions at this boundary:
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)):
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 Wc 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 ( Wc ) and outer ( Wv ) 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.
Low porosity approximation. I present examples of analytical solutions in 1-D cases simplifying equations (2.5) and (2.6):
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:
Here h0 is the dimensionlees length of compaction and k = h0-1 is a parameter characterizing the size of the definition region. Small, large and normal regions are defined by the respective conditions k 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)):
Boundary conditions (3.7) and (3.9) in 1-D problem IIIA reduce to (B.21) and (B.22):
The dimensionless expressions for initial condition (3.4) (common to all types of problems) and the total mass of disturbance m are
I derive the simplest nonstationary solutions to problems II and IIIA. In a small region ( k 1 ) equation (4.3) is simplified,
This equation has the first integral
where C1(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
where C2(x) is the second function of integration. Condition (4.7) defines one of these functions:
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,
Then, the porosity and fluid flow velocity are readily found as
This is a general solution depending on two unknown functions that can be determined from respective boundary conditions.
where s0(t) is a given function. Equation (B.18c) and the first relation in (4.12) yield
Substituting expression (4.11) into these boundary conditions, I obtain the solution
Now I consider the motion in the problem IIIA with inner boundaries in a small region ( k 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 integrodifferential equations
In order to obtain this solution in an explicit form, I differentiate in t the expression f(x-(t),t) = fmin:
This readily provides the expression for the trailing front velocity u-
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
Let the initial function j(x) be given as the parabola
Substituting this expression into (4.15) yields the ordinary differential equation of the first order describing the motion of boundary points
The system of two equations (4.14b) and (4.17) in two unknowns defines the movements of the leading and trailing fronts of a disturbance. A direct substitution shows that these equations have at least one solution,
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
In this case, the solution of (4.14a) reduces to the expressions
Differentiating in time the last expression in (4.18), I obtain the first-order equation
that can be readily integrated, implicitly defining the function h(t); more specifically, it gives the inverse function
With t, the function h(t) tends asymptotically to zero as h(t) = 32t.
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
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:
The solution of this boundary problem has the form
Problems with the skeleton breakdown. The porosity in problems of type IIIB is not a small value, but another small parameter ( e1 ) is introduced. Transformation to dimensionless variables and the use of the expressions for material functions (2.12) and (2.13) simplify equation (4.1) for each of these cases:
Equation (4.4) is simplified in a similar way:
In the case of these problems, boundary conditions (3.3) and (3.16) assume the form (see also (B.24)):
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
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 compaction. 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
The integration constant is found from boundary conditions B; in particular, asymptotic conditions (4.5) reduce to a simple expression:
Conditions (4.3) that must hold in problems of type IIIA give the relation
where fmin is a material constant. I set f0 f0 - 1u s0 and s0 = fk0 in the infinite region and f0 = fmin in the finite region. Then expressions (5.1) and (5.2) can be written in the general form
In this case, equation (4.3) reduces to the ordinary differential equation
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 x(p0) and x(0) at the leading and trailing fronts are generally different. With x(p0) < x(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 x(p0) = x(0), which simplifies boundary conditions (4.6):
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, this type includes asymptotic condition.
Wave equation (5.4) is a two-parameter equation because it contains the independent parameters f0 and u. However, the meaning of the parameter f0 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 wavelength is equal to h, and x = h2. 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:
This problem is studied in Addendum C, and its general symmetric solution is given by quadrature (C.14):
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 m and is given by formula (C.18):
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  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 f0 and u entering wave equation (5.4) in the form (C.11):
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
The fluid flow velocity is assumed to be directed in the positive direction of the x axis. Since s < smax and f < 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,
Two boundary conditions are possible. One implies the absence of surface forces:
the second condition with surface forces present is
Equations (5.10) and (5.11) have, respectively, the first integrals
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 (5.14a) alone under condition that the integration constant vanishes,
The solution of this equation reduces to the quadrature
where h2 is the position of the leading or trailing front of the boundary Gmax and
This wave has the length
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.
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 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.
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. [, 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 poroviscous 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 clayey rocks. The feeder walls are then washed away and the mud mixture is saturated with solid rock (clay and aleurolite) fragments.
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 chemical 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 quasi-solid 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.
W = WcUWv, total region of definition;
Wc, compaction region ( Wc(A) and Wc(B) are compaction regions adjoining boundaries of types A and B, respectively);
Df and e1, porosity interval and a small parameter characterizing the width of the region Wc(B);
Wv, regions of the porous medium with isolated pores, suspensions or the transitional cluster zone that are adjacent to the compaction zone;
G, any boundary (exterior or inner) of any region;
G ex, exterior boundary of the region W;
G int, any interior boundary of the region W;
Gmin, inner boundary (between the regions Wc and Wv ) at which pores close or open ( G(+)min and G(-)min are the leading and trailing fronts of this boundary);
Gmax, inner boundary (between Wc and Wv ) at which the skeleton of the poroviscous medium forms or breaks down;
pf, sij, sefij and ssij, pore pressure, tensor of total stresses, tensor of effective stresses and tensor of stresses averaged over the skeleton phase in the two-phase region Wc;
ws, wf and s, velocities of the skeleton, fluid and filtration;
f, fmax and fmin, porosity and its critical maximum and minimum values in the compaction region;
z and h, coefficients of the volume and shear viscosity;
d = hfk, coefficient of hydraulic resistance;
hf, fluid viscosity;
k, permeability coefficient;
Al and Bk, rheological constants in expression (2.7);
dij, Kronecker delta;
g = (0, 0, - g), gravity;
r = (1 -f)rs + frf, density of the two-phase medium;
rs and rf, densities of both phases assumed constant;
Dr = rs - rf, density difference between the skeleton and fluid phases;
c, cf and cs, specific heats of the mixture and individual phases;
l, heat conductivity coefficient;
l, latent heat of melting;
F, 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;
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 Wc;
F ( x), (j( x)) initial value of the porosity (and its normalized value measured from the minimum porosity value);
w sol, solenoidal component of the skeleton velocity;
f, scalar potential of the skeleton velocity;
t solij and P pf - 13sefkk, deviatoric stresses and pressure in a hypothetical incompressible Stokesian fluid, hydrodynamic equations of which are used as one of the representations of compaction equations;
tvij, wv and Pv, deviatoric stresses, velocity and pressure in a real incompressible Stokesian fluid in the region Wv;
Pa pa - 13sefkk and pa (Pav pav - 13sefkk, pav), reduced values of P and pf from which the hydrostatic pressure of the poroviscous medium in Wc (Wv) is subtracted;
Pb pb - 13sefkk and pb (Pbv pbv - 13sefkk, pbv), reduced values from which the hydrostatic pressure of the suspension in the region Wc (Wv) is subtracted;
Hc zd, dimensional compaction length;
h0, dimensionless compaction length;
k = h0-1, parameter related to the dimensionless compaction length;
[ j]|+- j+ - j-, j, boundary values of a quantity j at an inner boundary inside and outside the compaction region;
sef = 13sefkk, first invariant of the effective stress tensor;
p0 and s0, 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 = x(p0), minimum porosity as a function of the pore pressure;
si and fi, external surface forces at inner and exterior boundaries;
wsm = wsi m(a)i, where n and m(a) are the normal and tangential unit vectors (the index a = 1, 2 indicates two mutually perpendicular directions on the boundary surface);
k, vertical unit vector;
smax, maximum velocity of filtration at the inner boundary Gmax;
s and f, asymptotic values of the filtration velocity and porosity corresponding to the asymptotics at x ;
s0 and f0, trivial homogeneous solution (the state of homogeneous fluidization), as well as the asymptotic values of the filtration velocity and porosity at x 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.
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 G ex 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 ( w0 ) or forces ( f ) can be specified at boundaries of type II:
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:
Another typical case is a free surface at which external forces in a two-phase medium and the pore pressure vanish,
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 structural-phase transformations of the microscopic structure that lead to qualitative changes in the state of the medium. Such models include interior ( Wc ) and exterior ( Wv ) 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]:
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 P:
where f is the scalar potential of the skeleton velocity.
The solenoidal component of velocity can be brought into correspondence with the tensor tsolij
describing stresses associated with this component. It is connected with the total and effective tensors of stresses through the relations
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 P plays the role of pressure:
Taking the divergence and curl operator div and rot to equation (A.6b) yields the equations for P and rot wsol:
each of these values obeys the Poisson equation. Expression (2.5) can be rewritten in the form
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 pg of the suspension at a critical value of the fluid phase concentration (porousness) f enabling structural-phase transitions:
In the particular case f = 1, the suspension is transformed into the pure fluid,
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 inhomogeneous 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 P are transformed, respectively, to the forms
Equation (A.8a) is transformed in a similar manner to
Using decomposition (A.9) and (A.10), the following components are recognizable in the total stress tensor:
Summarizing the above analysis, I give equivalent representations of equations (A.4) and (A.6a) reflecting the specific character of canonical forms:
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:
Using (2.6a), the filtration velocity can be eliminated in this equation which can then be written in the equivalent forms
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 sij and sefij 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., div ws and pf. Both equations are interrelated through the P 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 w sol.
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). 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 - f) ws + f wf, continuity equations (2.6) can be written in the form
Similar to wsol, the velocity component v is solenoidal. However, this case does not reduce to transformation (A.6). I transform (2.5) to new variables:
The effective stress tensor is then transformed to the form
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.
Application of the P 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:
In the case of nonlinear rheology, these definitions remain the same with the additional condition: l = 0 and z = A0 = ( z0 + 43 h).
The boundary problem of type IIIA depends on eight independent dimensional constitutive parameters ( t, X a, p0, z, d, Drg and H(0) ) and on the dimensionless parameter Fmax. These parameters involve three independent dimensions: time, length and mass. According to the P 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 P theorem states that the number of independent constitutive parameters is exactly seven.
In accordance with the requirements of the P theorem, the following values are chosen as independent dimensionless constitutive parameters:
Here Hc is the dimensional compaction length. With such a choice of dimensionless variables, the value h0 k-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 f( x, t), filtration velocity s, fluid mass m etc.) through the independent values:
Using these relations, equations of motion in velocities A.5), continuity (2.1b, c) and coupling (A.3) assume the form
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:
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:
The initial condition is the same for problems of any type:
The small-porosity approximation is often used in problems of the first and second types,
which simplifies some expressions, in particular, (B.4a, b, d):
In this work, only a trivial asymptotics of the type
is used in the problem of the first type; s0 and f0 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 the equation
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:
Accordingly, boundary conditions (3.3)-(3.9) in problem IIIA are transformed to the following dimensionless form:
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
The 1-D case. Spatially one-dimensional models are most often encountered in applications. The finite region Wc(t) is a segment with moving boundaries x- = x-(t) < x < x+ = x+(t) of the length h(t) = x+ - x-. Here x+ and x- are boundary functions for the leading ( G+ in ) and trailing ( G- in ) 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:
On the strength of these simplifications, expressions for stresses (B.5a) and (B.6a) admit, after some transformations, the first integrals
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
providing additional constraints on the equation of motion and boundary conditions. For example, initial condition (B.7) assumes the form
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
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:
where the values s and f are connected through the relations similar to (B.10b)
Boundary conditions in problem IIIA are the same at both boundaries:
An additional condition is set at the leading front:
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)
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
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:
In the absence of surface forces, we have condition (3.15) transformed to the form
At l = 0 this relation yields
Problems IIIA. I analyze wave equation (5.4)
in finite and infinite regions. A finite region implies problem IIIA with boundary conditions (5.5):
The problem IIIA involves a constant value of the minimum porosity fmin 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 difference between the extreme and boundary points. The wave mass in problem IIIA is given by the expression
In problem I the mass is defined as
I examine wave equation (C.1) in the particular case l = 0 and k = 3:
This equation has the first integral I (the constant of integration)
In the particular case, when the maximum value of f at the upper point of the wave is equal to unity ( fmax = 1, f0< 1 ), the integration constant is fixed,
The solution can be obtained with the use of the quadrature of this expression:
The necessary condition for the existence of this solution is the inequality
implying that the radicand in (C.8) is positive. I analyze this inequality in more detail. Let f(1) and f(2) be the roots of the equation R(f) = 0 and let D be its discriminant,
In the case of an infinite region, both roots coincide ( D = 0 ), and minimum values of the porosity and fluid flow velocity coincide with their values in the trivial solution:
In this case
and quadrature (C.8) reduces to the expression
calculated in elementary functions:
The wave mass is determined by formula (C.4):
Using (C.7) and (C.10), this expression is transformed to
and can be calculated in elementary functions:
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 f0 = 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
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 the method of dynamic systems ( f, f ). Equation (C.5) is transformed into the dynamic system
Phase portraits in the plane ( f, f ) 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
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
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
this curve is tangent to the abscissa axis at the point f = 32f0. Then equation (C.15) has three real roots fi of which one is negative and two positive roots coincide:
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 f1:
If the solution in a finite region exists, quadrature (C.8) can be calculated and is represented through elliptic functions [Karakin, 1986]:
Problems IIIB. Here I analyze the wave solutions of (5.10)-(5.11),
Similar to the preceding case, equations (C.17) are transformed to the dynamic systems
Moreover these equations have the first integrals I
In the absence of surface forces, boundary condition (B.26) 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):
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 f1 and f2 on the f axis of the phase plane ( f, f ):
under the condition
To discuss these situations in more detail, I address a solution of (C.17a) with a positive wave velocity u and with boundary condition (5.12a) implying the presence of a boundary Gmax free of surface forces. This solution corresponds to the circular trajectory of the phase portrait in Figure 9 that is tangent to the vertical line f = 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 h2 is the position of the leading or trailing front of the boundary Gmax and
The length of this wave is given by the expression
All trajectories in Figure 9 that lie in the upper half-plane above the aforementioned specific trajectory but below the separatrix intersect the line f = 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
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.
The above analysis of phase portraits concerned only such situations when the inner boundary Gmax separates a poroviscous 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.
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 Structure 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.
Maaloe, 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 Syncrude's oil sands plant: 1978-1987, AOSTRA J. Res., 5, 109-113, 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, 1099-1102, 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, Wave-like 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.