Phenomenological Modeling of Catastrophic Dilute Gravity Flows

Résumé — Modélisation phénoménologique des écoulements gravitaires catastrophiques dilués — Les sédiments déposés en eaux profondes sont à l’origine de nombreux réservoirs pétroliers. C’est pourquoi les processus régissant leur mise en place suscitent actuellement un intérêt grandissant, en particulier dans le domaine de la modélisation numérique.


INTRODUCTION
The purpose of the present paper is a general exposition of a simple model for the equations of the motion of gravity flows.Special emphasis is placed on seeking analytical solutions for long times.The study shows their evolution when taking new physical phenomena into account.

What is a Gravity Flow?
Turbidites (deposits associated with turbidity flows) represent a significant fraction of the deep marine sedimentary record, and play an important role in hydrocarbon reservoir formation.This is why processes governing their settings have given rise to a growing interest, especially in the field of numerical simulation.
The development of such models would allow for: -a better comprehension of the process dynamics that controls deposits; -a reduction of hazard during the explorations by predicting the: • three-dimensional deposit geometry (mode of organization, lateral extension); • internal facies architecture; • connection between reservoir bodies; -quantification of the impact of external parameters (eustasy and subsidence, nature and amount of sediment supply) on the internal geometry.
In geology, the essential knowledge about natural flows is in fact based on subsequent deposits themselves [1][2][3].The conceptual basis for such an interpretation has therefore often been limited by the three-dimensional, unsteady and tels que la nature et la quantité des apports en sédiments et la géométrie du bassin -sur la dynamique des écoulements gravitaires et l'organisation des séquences de dépôt résultantes (turbidites).
Abstract -Phenomenological Modeling of Catastrophic Dilute Gravity Flows -Numerous hydrocarbon reservoirs originate from sediments deposited in deep water.This is why processes governing their settings have given rise to a growing interest, especially in the field of numerical simulation.Taking the problematic and constraints imposed by the geological data into account, the method that we have adopted to obtain a mathematical model for diluted and turbulent finite gravity flows is the following: the two-dimensional movement created by the instantaneous release of a finite volume of heavy fluid (suspension of sediment particles) into a lighter one (water), on variable slopes, is theoretically studied as a model for gravity flows.These flows develop a characteristic longitudinal structure, comparable to a deformable semi-lens, where the height is small relative to the length.This geometry is imposed on the gravity flow.Using the Boussinesq's approximation, the flow dynamics is supposed to be governed by a balance between gravity driving forces, inertia and turbulent friction.The study of the internal longitudinal flow velocity field allows a law of variation for the spreading velocity to be formulated.An equation, including effects of water incorporation at the suspension-ambient fluid interface, quantifies the variation of the total volume of the flow.Finally, a transport equation for the particles volume concentration is proposed assuming that: -turbulence creates a uniform density distribution in the flow; -particles are advected at the mean flow velocity; -particles fall out or are eroded in the viscous sublayer of the flow.The coupled system of the non-linear differential equations obtained is solved numerically.The model is then validated by experimental small-scale models realized by Laval (1988).The comparison between theoretical predictions and experimental results shows good agreement.An analytical study of the system, by local analysis methods for long times, shows the evolution of solutions when taking new physical phenomena into account and the consistency of the obtained numerical solutions.The obtained model leads to very low computational time on a microcomputer.Simple and complete, it constitutes a first step towards quantitative comprehension of the impact of external parameters-such as the nature and the amount of sediment supply and the geometry of the basin-on gravity flow dynamics and on the organization of subsequent depositional sequences (turbidites).
non-uniform nature of the flows responsible for these deposits.
In mechanics, diluted solid and fluid mixtures are named suspensions.Their study is of great interest for sedimentologists, since they include processes, mainly the deposit of solids in suspension, that do not occur in density flows 1 .Practical examples of such flows are numerous, one can mention: dust storms, fronts of marine breezes [4], snow avalanches [5][6][7][8], ash clouds caused by volcanic eruptions [9,10], water plumes loaded in mud at river mouths [11][12][13], and of course turbidity currents [14][15][16][17][18][19].Two types of turbidity flows can be distinguished: -continuous flows, with relatively long life (being able to be steady and/or uniform), initiated by: • sediments in suspension brought into the sea by rivers in flood or by the glaciers melting (hyperpycnal flows) [13,20]; • storms stirring unconsolidated sediment layers which constitute a deep marine layer with concentrated mud in suspension (bottom nepheloid layer); -unsteady flows, with relatively short life, originating from: • transformation due to mixing with seawater, for example during a hydraulic jump, of slides, slumps or debris flows; all the result of a catastrophic and instantaneous rupture of the slope provoked by an earthquake and/or slope instability [21][22][23][24][25][26][27][28]; • mass instability of sands brought by flows to the underwater canyon head [12,29].Their movement can be divided into three phases: -Phase I includes the initiation of the flow and its acceleration.Even though the flow is not necessarily a turbidity current since its concentration can be high and its behavior can be viscous or even visco-plastic [30], this first phase is characterized by an amount of eroded sediment probably superior to that deposited.In some cases a very great amount of the deposited mass comes from the erosion along its path (nearly 80% during the Nice airport collapse in 1979 [27]); -Phase II can be a period of constant velocity if the slope is constant.During this stage, the amounts of eroded and deposited sediments can be equal.The flow dilutes and sediments are maintained in suspension by turbulence; -Phase III is a period of deceleration, where the sedimentation dominates the erosion.When the flow reaches the slopes of the abyssal plain, it stops, due to the loss of excess density by incorporation of ambient water and the settling of particles; its spreading; a dissipation of the turbulence by friction at the interfaces or at internal density stratifications [15,[31][32][33].
(1) Density flows are characterized by the movement of a heavy fluid into a lighter one (for example: flows of heavy industrial gases and masses of air flowing to the surface of the ground).

Existing Modeling Methods
Since Phase II can be considered, on constant slopes, as a period of constant velocity, with neither erosion nor sedimentation, the first theoretical and experimental researches have focused on it [34].For example, in [35,36], considering that a river is a form of density flow, in which the density contrast between the flow and the ambient fluid (air) is great, and in which there is very little friction resistance at the upper interface, the author determines the speed of a density flow by (Appendix A for notations): ( This formula, analogous to that found in [37] for fluviatile flows, is obtained simply by writing the balance between the apparent weight, g'φhs, and the steady friction resistance . The conception of a mechanical theoretical model for particle and fluid mixtures can be carried out by one of the following three approaches [34]: -microscopic: the resolution is accomplished directly on the motion equations of each sediment particle (ideally rigid or elastic spheres).This discrete method is however limited in so far as it does not permit the simulation of great numbers of particles.It provides nevertheless very precious information on particular phenomena such as sedimentation [38]; -continuous: in this case, material properties of the mixture can be represented by continuous functions and the media can be divided indefinitely without losing its properties.Indeed, the height of the flow is sufficiently great relative to the particle diameter so that the behavior of the granular suspension is locally isotrope and continuous, contrary to the individual particle flow.In this case, the notion of discrete particles or granules is no longer preserved.This method [39][40][41][42][43][44], based on Navier-Stokes equations, which include varying degrees of simplification (steady state, two-dimensionality, boundary layer theory, Boussinesq's hypothesis, neglected sediment flux to the bed), is again limited by high computational times, and therefore does not allow the simulation of series of flows.However, studies that have analyzed the longitudinal flow evolution have shown that processes such as sediment flux to the bed, turbulent (or viscous) friction, marine bottom slope, and turbulent entrainment of ambient seawater, control changes in the state of the flow [44]; -empirical: the modeling carried out by this approach is based not entirely on theory but also on analogical experiments.Dimensional analysis allows the most influent parameters to be determined.Obtained formulae contain constants that are adjusted on experiments.The best known models for geological applications are based on a steady balance between gravity and friction forces acting on the head [45], or on a semi-infinite flow [35,36,46,47] (Eq.( 1)).

A PHENOMENOLOGICAL MODEL IN THIN LAYER THEORY FOR GRAVITY FLOWS
None of the previous approaches seems to be entirely appropriate for our purpose.So, taking the problematic and constraints imposed by the geological data into account, the method that we have adopted to obtain a mathematical model for diluted and turbulent gravity flows is the following.
To characterize the temporal evolution of four characteristic quantities linked to the flow (Fig. 1): -U, the mean motion velocity of the flow object; l, the length (of evolution speed U l ); h, the height; -φ, the particles volume concentration (i.e. the particles mass), we have considered that the gravity flow of a finished volume of particle mix reflects the balance between: -the driving gravity force; -the turbulent friction; -the spreading due to pressure forces; -the incorporation of the ambient seawater; -the deposit and/or the erosion of sediment particles.
To present as clearly as possible the theoretical model in two dimensions, major hypotheses as well as physical and geometrical quantities definitions are stated in the following section.

Hypotheses and Definitions
Observations of atmospheric and pyroclastic density flows [7,10] as well as experiments on submarine density and suspension plumes [14,15,[17][18][19]26] have shown that turbidity flows develop a characteristic longitudinal structure, comparable to a deformable semi-lens.
Geometry and variables of the flow object represented here in two dimensions by a semi-ellipse.

Hypothesis 1
The flow is a deformable object with an imposed geometrical form parameterized by two variables: h, its height and l, its length (Fig. 1).

Hypothesis 2
The flow height is small relative to the flow length: Under these assumptions, the volume of the flow per unit of width is defined by: where C is a geometrical factor.For example, in the case of a semi-ellipse, C = π/2.
The flow-substratum contact surface, per unit of width, is given by: The flow-ambient fluid interface surface, per unit of width, is approximated by: The flow being thin, it is clear that the mean motion, U, will be mainly parallel to the ground.Indeed, Equation ( 2) between characteristic dimensions of the flow implies, via the equation of continuity, that the characteristic velocity scale normal to the wall, w, is smaller than the characteristic velocity scale parallel to this wall, u (u >> w).

Hypothesis 3
The fluid, within the flow, is turbulent.Thus, the Reynolds number Re is high.
It is likely that very rapidly, just after the releasing of the flow (Phase I), friction forces, first viscous or visco-plastic, become turbulent.Classically, turbulent friction is proportional to the square of a characteristic velocity of the flow.In a very simple manner, it is possible to write:

Hypothesis 4
The suspension is sufficiently diluted to justify the use of Boussinesq's approximation: Then it is possible to neglect φ∆ρ in inertial terms, while conserving it in gravity terms.
The velocity of the solid phase is assumed to be equal to the sum of the flow velocity and settling velocity:

Obtained Model
Under the previous hypotheses, the evolution of a gravity flow is determined by the following non-linear coupled differential system, which will be explained in the next section: -object mean motion: (3a) -object spreading: (3b) -ambient water incorporation: For clarity, geometrical constants have been omitted.The spreading velocity U l is defined by: ; the front velocity, U f , by: .

PHYSICAL DERIVATION OF THE MODEL
We give here a short explanation of System (3) in terms of physics.

Motion of the Object: U
The variation of momentum of the flow object, in two dimensions, under the hypotheses in Section 1.1, is modeled by Newton's second law: (rate of change of momentum = longitudinal component of the effective gravity + turbulent friction force along the flow surface) where A is the volume of the flow; U, the mean velocity of the flow, φ, its volume concentration in suspended sediments; g', the gravitational effective acceleration; θ, the slope angle; k t1 , a turbulent friction coefficient; ∂A, the surface of the flow.
The effect of choosing a symmetrical geometry for the flow object eliminates the pressure gradients.One of the implications of such a model for the variation of momentum is that there cannot be sustained mean motion on perfectly horizontal surfaces (θ = 0°).
The velocity, U, determines the position of the object at each time t by: , where X is the curvilinear abscissa.

Geometrical Evolution: l and h
To characterize the geometrical evolution, we consider two effects separately: spreading of the object and incorporation of the ambient water at the upper interface.

Spreading of the Object
The study of the internal velocity field shows that the spreading velocity U l (Fig. 2) follows the equation of evolution: (rate of change of momentum = longitudinal component of the pressure forces + turbulent friction force along the flow surface) where k t2 is a turbulent friction coefficient.The flow lengthens to the velocity U l , therefore .

Incorporation of Ambient Seawater
During the decelerating phase (Phase III), but not exclusively, the flow incorporates the ambient fluid along its Modeling of the spreading phenomenon of the flow object, where U l is the spreading velocity.
Modeling of the ambient fluid incorporation phenomenon.
upper interface (Fig. 3) [16,17].The variation of the flow volume is determined by: (rate of variation of the volume = rate of incorporation of ambient fluid along the upper interface) where E f is the rate of seawater incorporation; ∂A s , the upper surface of the flow.This equation is a simplification of the mass balance equation (continuity equation) under assumption (4).
The dimensional analysis leads to: Considering the mean velocity of the flow, U, as characteristic of the phenomenon, the rate of seawater incorporation, E f , can be written: * is a dimensionless incorporation rate, function of the Richardson number Ri.This dimensionless number is defined by: (4) It is the ratio of gravity to inertia.It characterizes the stability of the flow interface.
The law determining the dimensionless incorporation rate of ambient fluid by the flow, E f * , is an important property.It distinguishes gravity flows from density and fluviatile flows.The water incorporation increases the volume of the flow and reduces the concentration in sediment particles in the flow.The following empirical formula to determine this incorporation rate is given in [43]: (5) where Ri is the Richardson number given by Equation (4), C 1 = 0.00153 and C 2 = 0.0204.
According to Equation ( 5), when Ri approaches 0, a constant value of 0.075 is appropriated for E f * .When Ri tends to infinity, gravity efforts stabilize the interface; the water incorporation diminishes, then E f * decreases as Ri -1 (Fig. 4).

Sediment Mass Evolution: φ
Here again we consider two effects: sedimentation and erosion.However, they will be included in a same equation.Indeed, the variation of the mass of particles can be quantified by a balance equation of the type: where Q is a mass flux, and ∂A, a surface.We will therefore write: (rate of change of particle mass = entrainment and sedimentation rate along the flow bottom interface) where S p is the sedimentation rate; E p , the sediment entrainment rate; ∂A i , the lower interface of the flow.A dimensional analysis of the mass flux Q, of Equation ( 6), shows that: (7) Therefore, for each phenomenon of mass transfer, we will have to find a characteristic volume mass and a characteristic velocity.

Sedimentation
Sediments in suspension fall out of the flow by gravitational effect within the viscous sublayer (Fig. 5).Taking Equation (7) into account, the sedimentation rate S p intervening in the particle mass balance is assumed to be a function of: -the volume concentration of suspended sediments immediately above the bed [42,43]; -the settling velocity V s (function or not of the presence of other particles, i.e. of the other particles volume concentration).The flux S p can therefore be modeled by: Modeling of the sedimentation phenomenon.
Modeling of the erosion phenomenon.
where b is a distance very close to the bed such that: b << h.The terms are indeed evaluated just slightly above the bed to avoid singular behavior due to having neglected the molecular diffusion at the bottom of the flow.It can be assumed that is of the following form: C 3 φ with C 3 ≈ 1.56 [43].

Erosion
As we have underlined, a flow, essentially in its acceleration phase (Phase I), is able to erode, along its lower surface ∂A i , sediments on its path (Fig. 6).Based on dimensional argument (7), the sediment entrainment rate, E p , can be decomposed in the following manner: * is a dimensionless entrainment rate, function of the particles Reynolds number Re p [42,43].This dimensionless number is defined by:

Discussion
Equations (3a) and (3c), with k t1 and k t2 equal to 0, are similar to the ones considered in [48]; Equations (3a), (3c) and (3d) with E p * = 0, to those considered in [49].By characterizing the evolution of the different quantities linked to the flow, we have introduced new unknown parameters, such as the settling velocity, V s , the dimensionless entrainment rates, E p * and E f * , the bottom volume concentration of suspended sediments, , and the turbulent friction coefficients, k t1 and k t2 .Numerous empirical laws describe their evolution.A review of the settling velocity and suspension viscosity laws can be found in [30].

ANALOGICAL EXPERIMENTS
A numerical code in C++ language has been developed to solve non-linear differential system (3) by the fourth-order Runge-Kutta method [50,51] with some computational options (Appendix B).
The objective of the experiments found in [19] on analogical small-scaled models was to study density surges 2 and suspensions.They were carried out in an inclined constant slope channel (0.2 m wide, 0.35 m deep and 4 m long) immersed in a full water reservoir.
Front velocity data allow a validation of the theoretical laws (System ( 3)).Numerical predictions on front velocity, with default values of the models found in the literature [42][43][44] (Appendix B), show a good agreement with experimental data (Fig. 7).

ANALYTICAL STUDY
We present here the analytical study of the above-described model (System ( 3)).Classic methods of local analysis allow analytical solutions for long times to be found [52].The study is accomplished progressively and shows the evolution of solutions when taking new phenomena into account.

Undeformable Sliding
The system is then composed of Equation (3a), which can be cast in the form: .with G > 0 and K > 0. Critical points are: The first critical point, U s1 , is the same steady velocity as in Equation (1) [35,36].The study of the linearized equation shows that U s1 is a stable node, and U s2 is an unstable node.An exact analytical solution is found on constant slopes, for U 0 > U s2 : (8) The numerical model exhibits the same behavior (Fig. 8).As slopes (sin θ) vary as X γ , seeking solutions of the form (analogous to Frobenius series [52]): (9) by local analysis methods for long times (i.e. as → ∞), shows that the first-order term (dominant term) of the velocity is: whatever the initial conditions are.On decreasing slopes (γ < 0), the mean flow velocity, U, decreases.On constant slopes (γ = 0), the leading behavior of the solution is independent of time, it is a constant.This constant is U s1 as show Equation (8) and Figure 8.

Taking the Spreading into Account
This case corresponds physically to a non-miscible density flow.The system is then composed of Equations (3a) and (3b) that can be cast in the respective forms: completed with the volume conservation: The critical point is: (U, l, h) = (0, l ∞ , 0).Again, seeking solutions of form (9) on slopes varying as X γ by local analysis methods shows that, whatever the initial conditions are, the first-order terms are: On decreasing and constant slopes (γ ≤ 0), the flow slows down while thinning down and spreading regardless of γ.The numerical model predicts the same behavior (Fig. 9).

Taking the Water Incorporation into Account
This case corresponds physically to a suspension without sedimentation (for example, on slopes that are too steep for sediments to remain in place) or a density surge.The system  Numerical solutions and leading behaviors as a function of time for the velocity U of an undeformable sliding with quadratic friction law on constant slopes.Numerical and analytical solutions for U 0 = 0 are identical.
is then composed of the three equations (3a), (3b) and (3c) that can be simplified to: The mass of sediment particles remains constant in the flow as expressed by: The critical point of the system is: (U, l, h, φ) = (0, l ∞ , h ∞ , 0).Seeking solutions of form (9) for long times shows that, on slopes varying as X γ , whatever the initial conditions are: On decreasing and constant slopes, the slows down while expanding and diluting, and this regardless of γ.On constant slopes (γ = 0), it expands with a constant aspect ratio.The numerical simulation reproduces these solutions well (Fig. 10).

Taking the Sedimentation into Account
This case corresponds physically to a suspension in sedimentation.The system is then composed of the four equations (3a), (3b), (3c) and (3d) that can be cast in the respective forms: The critical point is: (U, l, h, φ) = (0, l ∞ , h ∞ , 0).Due to Equation (3d), all solutions are no longer of form ( 9) for long times.It is however possible to show that, defining m as: first-order terms are: the Richardson number, Ri, defined by Equation ( 4), being very small.
The leading behaviors of the solutions are also found numerically (Fig. 11).As for the case of Section 4.3, the flow slows down while expanding and diluting, but here, because of the exponential decrease of the volume concentration, the leading behaviors are independent of the slope.We remark that the numerical model predicts a transitional phase (~ 10s), where the flow first thins down and accelerates with a non-constant aspect ratio.After this, it slows down while expanding.This transitional time is of great interest for geoscientists since most of the deposition takes place during this period.The leading behaviors of U, h and l are only reached after the volume concentration has become very small and all the particles have settled out of the flow.

CONCLUSIONS
The model we have developed for gravity flows is a hybrid model, in the sense that: -the flow, assimilated to a geometrical object, can be considered as a discrete deformable particle; -internal variables of the object are continuous functions.
Laws governing the evolution of these functions are deduced directly by the integration of mechanics of continuous media balance laws; -closures are obtained by dimensional analysis and empirical laws.Despite some limitations, due to the geometrical simplification of the flow and the use of empirical laws from fluviatile flow studies, which restrict its validity, the model has the main following advantages: -it conserves a physical description for dynamics and mass transfers; -the flow evolution is governed by a of ordinary differential equations and not by a set of partial differential equations, more complex to implement numerically; -the flow is treated as an object.
Therefore, this complete and original model has a fundamental advantage: it is a simple case.It presents however the disadvantage, at present time, to be limited to the particular cases of turbulent diluted flows.Nevertheless, it can be a good approximation of gravity phenomena.With such a model, computational times are very short, therefore it is foreseeable to simulate series of events-even several events in parallel-and to form multi-event depositional sequences.

Figure 4 Dimensionless
Figure 4Dimensionless seawater incorporation rate, E f * , as a function of the Richardson number Ri.

Figure 7 Comparison
Figure 7 Comparison of front velocity, U f , as a function of the distance X, predicted by numerical and analogical experiments.Default computational values (Appendix B).

( 2 )
Release of a finite volume of dense fluid into a lighter one.

Figure 9
Figure 9Numerical solutions and leading behaviors as a function of time for a non-miscible density flow on constant slopes (γ = 0).Plain lines are leading behaviors, dotted lines are numerical solutions.

Figure 10 Numerical
Figure 10Numerical solutions and leading behaviors as a function of time for a suspension without sedimentation on constant slopes (γ = 0).Plain lines are leading behaviors, dotted lines are numerical solutions.

Figure 11 Numerical
Figure 11Numerical solutions and leading behaviors as a function of time for a suspension in sedimentation as k t1 C 2 /C 1 → 0, i.e. as m → 1/3.Plain lines are leading behaviors, dotted lines are numerical solutions.