aa r X i v : . [ phy s i c s . g e o - ph ] M a r A Mathematical Model for Voigt Poro-Visco-Plastic Deformation
Xin-She Yang
Department of Mechanical Engineering, University of Wales Swansea, Singleton Park, Swansea SA2 8PP, UK
Abstract
A mathematical model for poro-visco-plastic com-paction and pressure solution in porous sedimentshas been formulated using the Voigt-type rheologicalconstitutive relation as derived from experimentaldata. The governing equations reduce to a non-linear hyperbolic heat conduction equation in thecase of slow deformation where permeability isrelatively high and the pore fluid pressure is nearlyhydrostatic, while travelling wave exists in theopposite limit where over-pressuring occurs and thepore fluid pressure is almost quasi-lithostatic. Fullnumerical simulation using a finite element methodagree well with the approximate analytical solutions.
Citation detail:
X. S. Yang, A mathematicalmodel for Voigt poro-visco-plastic deformation,
Geo-phys. Res. Lett. , (5), 10.1029/2001GL014014(2002). Many physical properties such permeability, viscos-ity, Young’s modulus and thermal conductivity varywith porosity or fluid content in sedements and min-erals. The porosity in turn depends on the de-formation and compaction state and can be calcu-lated from the compaction curves resulting from theproper compaction modelling [Audet and Fowler,1992]. Furthmore, compaction is also related to over-pressuring, mineral deposition, hydrocarbon genera-tion and oil migration in reservoir. Thus the cor-rect modelling of compaction is both of scientificimportance as well as industrial interest. However,due to the nonlinear feature in the compaction pro-cess and difficulty in formulating the accurate andrealistic rheological relationship in sediments androcks, most of existing models use simplified rheol-ogy such as poroelastic, purely viscous or viscoelasticrelations [
Rutter , 1976;
Wangen , 1992;
Holzbercher ,1998;
Revil , 1999;
Yang ,2000]. The rheological prop-erties of realistic granular sediments are usually vis-coelastic or viscoplastic as implied by experiments.
Yang [2000] presents a viscoelastic model of Maxwelltype and comparison of analytical solutions with thenumerical simulations shows very good agreement.However,
Revil ’s [1999] work suggests that it maybemore appropriate to use a poro-visco-plastic modelof Voigt type. We intend to do the analysis similarto
Yang ’s [2000] but using the Voigt-type consitutiverelation as derived by
Revil [1999].Much of the work in this area has been re-viewed by
Rieke and Chilingarian [1974],
Birch-wood and Turcotte [1994] and later by
Fowler andYang [1999]. This paper aims at providing a new ap-proach to compaction and pressure solution by us-ing
Revil ’s visco-poro-plastic relation of Voigt type[
Revil , 1999]. The nonlinear partial differential equa-tions are then analysed by using asymptotic methodsand the obtained analytical solutions are comparedwith numerical simulations. Although the presentwork mainly concerns the 1-D theorectical formula-tion and analytical solution procedure, however, weintend to provide a simplified and yet realistic frame-work for further research in this area and shows howcompaction mechanism is related to rheological re-lationships and material properties of porous sedi-ments, so that more realistic constitutive relation-ships can be formulated and analysed.
The fundamental model of compaction and pressuresolution is essentially similar to the model of soil con-solidation process. The solid sediments act as a com-pressible porous matrix, so mass conservation of porefluid together with Darcy’s law leads to an equationof the general type. Based on earlier work by
Au-det and Fowler [1992],
Revil [1999] and
Yang [2000],we can write down the poro-visco-plastic compactionmodel of Voigt type. In the one-dimensional case, wehave the following governing equations ∂ (1 − φ ) ∂t + ∂∂z [(1 − φ ) u s ] = 0 , (1) ∂φ∂t + ∂ ( φu l ) ∂z = 0 , (2)1 ( u l − u s ) = − k ( φ ) µ [ ∂p∂z + ρ l g ] , (3) − G ∂p e ∂z − ∂p∂z − ρg = 0 , ρ = ρ s (1 − φ ) + ρ l φ, (4)where φ is porosity. u l and u s are the velocities offluid and solid matrix, respectively. k and µ are thematrix permeability and the liquid viscosity, ρ l and ρ s are the densities of fluid and solid matrix, p e is theeffective pressure, p l is the pore pressure, and g is thegravitational acceleration. G = 1 + 4 η / ξ is a con-stant describing the material properties with η and ξ being the shear modulus and bulk viscosity [ Birdet al , 1977]. The first two equations are the conser-vation of mass for the solid phase and liquid phase,respectively. The third equation is the Darcy’s law in1-D form and the last equation is actually the forcebalance in a simplified form whose detailed deriva-tion can be found in [
Fowler and Yang , 1998]. Com-bining equation (3) and (4), we have φ ( u l − u s ) = k ( φ ) µ [ − G ∂p e ∂z − ( ρ s − ρ l )(1 − φ ) g ] , (5)In writing the above equations, we have used an up-ward coordinate z originating from z = 0, whichcorresponds to the bottom of the sedimentary col-umn, so that the ocean floor z = h ( t ) moves as com-paction proceeds. We use such a coordinate systembecause it simplifies the analytical solution proce-dure and also in keeping with the similar lines ofearlier work in this area [ Audet and Fowler , 1992;
Yang , 2000]. However, the conventional depth co-ordinate is simply z − h ( t ), thus the transformationshall be straightforward once the basin thickness h ( t )is known. As we shall see in the later sections, weprovide an explicit formula for h ( t ) as a very goodapproximation.In addition, a rheological compactional relation-ship derived from experimental data [Revil, 1999] isneeded to complete this model in the form p e = − ξ ∂u s ∂z − E Z t ∂u s ∂z dt, (6)where E is the elastic modulus and ξ is the viscositymodulus. There are essentially the same parame-ters as introduced by Revil [1999]. The first term ofthe right hand of the equation is the usual contribu-tion by viscous plastic deformation, while the secondterm corresponds to the poro-elastic deformation.
To write the governing equations in dimensionlessforms, typical length and time scales are required. For a typical sedimentation rate ˙ m s , the correspond-ing typical time scale is d/ ˙ m s where the typicallength scale d can be defined as d = { ξ ˙ m s G ( ρ s − ρ l ) g } / , (7)so that the dimensionless pressure p = Gp e / ( ρ s − ρ l ) gd = O (1). Meanwhile, we scale z with d , u s with ˙ m s , time t with d/ ˙ m s , permeability k with k .By writing k ( φ ) = k k ∗ , z = dz ∗ , ..., and droppingthe asterisks, we thus have − ∂φ∂t + ∂∂z [(1 − φ ) u s ] = 0 , (8) ∂φ∂t + ∂ ( φu l ) ∂z = 0 , (9) φ ( u l − u s ) = λk ( φ )[ − ∂p∂z − (1 − φ )] , (10) p = − ∂u s ∂z − Ξ Z t ∂u s ∂z dt (11)where λ = k ( ρ s − ρ l ) gµ ˙ m s , Ξ = EG ( ρ s − ρ l ) gd . (12)Adding (8) and (9) together and integrating fromthe bottom, we have u s = − φ ( u l − u s ) , (13)where u = φ ( u l − u s ) is the Darcy flow velocity. Now,we have ∂φ∂t = ∂∂z [(1 − φ ) u s ] , (14) u s = λ ( φφ ) m [ − ∂p∂z − (1 − φ )] . (15) p = − ∂u s ∂z − Ξ Z t ∂u s ∂z dt, (16)where we have used the nonlinear constitutive re-lation for permeability k ( φ ) of typical form [ Smith ,1971] k ( φ ) = ( φφ ) m , (17)where φ is the initial depositional porosity. Theexponent m has a typical value of 3 ∼ ∂p∂z − (1 − φ ) = 0 (or equivalently , u s = 0) , at z = 0 , (18) φ = φ , p = 0 , h = ˙ m s + λ ( φφ ) m [ ∂p∂z − (1 − φ )] at z = h ( t ) . (19)It is useful to estimate these parameters by usingvalues taken from observations. By using the typi-cal values of ρ l ∼ kg m − , ρ s ∼ . × kg m − ,k ∼ − − − m , µ ∼ − N s m , ξ ∼ × N s m − , ˙ m s ∼
300 m Ma − = 1 × − m s − , g ∼
10m s − , E ∼ N / m , G ∼ , d ∼ λ ≈ . ∼ ∼
40. We can see that themain parameters λ and Ξ, which govern the evolu-tion of the fluid flow and porosity in sedimentarybasins, are the ratios of permeability to sedimenta-tion rate and material modulus to the typical pres-sure scale. As λ is essentially controlled by the hy-draulic conductivity, so it becomes the dominant pa-rameter controlling the whole compaction and pres-sure solution processes. Since the nondimensional parameter λ ≈ . ∼ λ ≪ λ ≫
1) will have verydifferent features in porosity and flow evolutions. Infact, λ = 1 defines a transition between slow com-paction ( λ ≪
1) and fast compaction ( λ ≫ λ ≪ Fowler and Yang , 1998,1999] to obtainsome analytical asymptotic solutions. λ ≪ ) In the nearly hydrostatic case of λ ≪ z ∼ t ∼ p ∼ u s ≪ ∂φ∂t ≈
0, then φ ≈ φ .We thus have ∂φ∂t ≈ − λ (1 − φ ) ∂ p∂z , (20) u s ≈ λ [ − ∂p∂z − (1 − φ )] , (21)and using (20), we have p ≈ − ∂u s ∂z − Ξ Z t − φ ) ∂φ∂t dt = − − φ ) ( ∂φ∂t + Ξ φ ) , (22) combining these above three equations, we have asingle equation for φ∂φ∂t = λ Ξ ∂ φ∂z + λ ∂ φ∂t∂z . (23)As the Ξ ≫ λ ≪
1, we can use the approxi-mation φ t ≈ λ Ξ φ zz so that we have ∂φ∂t = λ Ξ ∂ φ∂z + 1Ξ ∂ φ∂t . (24)with appropriate boundary conditions ∂φ∂z ≈ − (1 − φ ) Ξ , on z = 0 , (25) φ → φ , z → ∞ , (26)This problem is in fact equivalent to the problemof hyperbolic heat conduction or non-Fourier heatequation which is well-documented in heat transferand laser pulse modelling [ Antaki , 1997]. By usingthe Laplace transform method, we can write the so-lution approximately in terms of Bessel functions J i as φ ≈ (1 − φ ) √ λ Ξ t ierfc( ζ ) − (1 − φ ) Ξ [ √ λ Ξ t + λ ∞ X i =1 J ( zα i ) J ( α i ) α i ] , (27)where ierfc( ζ ) = 1 √ π e − ζ − ζ erfc( ζ ) . (28)and α i is the i th real non-negative root of equation J ( α i ) = 0. We can see that compaction essentiallyoccurs in a boundary layer near the bottom with athickness of the order of √ λ . λ ≫ ) In the case of λ ≫
1, the dependence of permeabilityon porosity ( φ/φ ) m decrease dramatically, so that λ ( φ/φ ) m is only bigger enough when φ > φ ∗ = φ exp[ − (ln λ ) /m ]. Thus, we have ∂φ∂t ≈ (1 − φ ∗ ) ∂u s ∂z , (29) ∂p∂z ≈ (1 − φ ) , (30)and using the equation (29), we get p = − − φ ∗ ) ∂φ∂t − Ξ Z t − φ ∗ ) ∂φ∂t dt = − − φ ∗ ) ( ∂φ∂t + Ξ φ ) . (31)3ombining these equation, we can get a single equa-tion for φ (1 − φ )(1 − φ ∗ ) = ∂ φ∂t∂z + Ξ ∂φ∂z . (32)Now we can seek the traveling wave solution of theform φ = φ ( ζ ) with ζ = z − ct , so that we have(1 − φ )(1 − φ ∗ ) = − cφ ′′ + Ξ φ ′ (33)where φ ′ = dφ/dζ . We can easily write the solutionas φ = 1 − (1 − φ ) exp[ [Ξ − p Ξ + 4 c (1 − φ ∗ )] ζ c ] (34)In fact, the above solution is only valid for the toppart when φ < φ ∗ . Using equation (6) and φ ∼ φ ∗ as ζ → −∞ , the travelling wave implies that cφ + (1 − φ ) u s = cφ + ( ˙ m s − c )(1 − φ ) . (35)so that we have c ≈ ˙ m s ( 1 − φ − φ ∗ ) , h ( t ) ≈ ˙ m s ( 1 − φ − φ ∗ t, (36)which means the basin thickness increases linearlywith time. In order to check the accuracy of the above analysis,we used a finite element method to solve the aboveequations (8)-(11). For simplicity, we only presentthe related results in Fig.1 where t = 5 with differentvalues of the λ = 0 . , . , , φ = 0 . λ ≫ λ case. A mathematical model for poro-visco-plasticcomapction and pressure solution in porous sed-iments has been formulated using the Voigt-typerheological constitutive relation as derived fromexperimental data by
Revil [1999]. After the properscalings, the governing equations reduce to a systemof coupled partial differential equations of mixedtype. In the case of small λ where the sedimentationis fast, permeability is small and the pore fluidpressure is nearly hydrostatic, the pressure solutionprocess reduces to the case of hyperbolic heat con-duction equation with a boundary layer forming at the bottom of the compacting column. On the otherhand, for the large λ ≫ Yang , 2000]using the viscoelastic rheological relation, we seethat the there is no essential difference in the casesmall deformation ( λ ≪ λ ≫
1) ismore complicated. Although both viscoelastic andporo-viscous-plastic cases have a transition at thedepth where φ ≈ φ ∗ , however, the mechanismsabove and below the transition are very different.In the viscoelastic case, the top region is essentiallyporoelastic while the lower region is almost purelyviscous, while in the present poro-visco-plastic case,the mechanism is a complicated combination ofviscous-plastic mechanism and poro-elastic defor-mation process controlled by the proper balanceof the parameters λ and Ξ. Furthermore, in com-parison with the earlier results, this work suggeststhat Voigt-type poro-visco-plastic deformationhas much interesting characteristics due to itstime-dependence feature. Acknowledgments:
The author would like tothank the referee(s) for their instructive and helpfulcomments which have greatly improved the originalmanuscript.
References [1] Antaki, P J, 1997. Analysis of hyperbolic heatconduction in a semi-infinite slab with surfaceconvection,
Int. J. Heat Mass Trans. , , 3247-3250.[2] Audet, D.M. & Fowler, A.C., 1992. A mathe-matical model for compaction in sedimentarybasins, Geophys. J. Int. , , 577-590.[3] Birchwood, R. A. & Turcotte, D. L., 1994.A unified approach to geopressuring, low-permeability zone formation, and secondaryporosity generation in sedimentary basins, J.Geophys. Res. , , 20051-20058.4 D ep t h z − h ( t ) Figure 1: Comparison of numerical simulations(solid) with analytical solutions (dashed curves) fordifferent values of λ = 0 . ∼ SIAMJour. Appl. Math. , , 365-385.[7] Fowler, A. C. and Yang, X. S., 1999. PressureSolution and Viscous Compaction in Sedimen-tary Basins, J. Geophys. Res. , B , 12 989-12997.[8] Revil, A., 1999. Pervasive pressure solutiontransfer: a poro- visco-plastic model,
Geophys.Res. Lett. , , 255-258.[9] Rieke, H.H. & Chilingarian, C.V., 1974. Com-paction of argillaceous sediments, Elsevier, Am-sterdam, 474pp.[10] Rutter, E. H., 1976. The kinetics of rock defor-mation by pressure solution, Philos. Trans. R.Soc. London Ser.
A 283 , 203-219.[11] Smith, J.E., 1971. The dynamics of shale com-paction and evolution in pore-fluid pressures,
Math. Geol. , , 239-263. [12] Wangen, M., 1992. Pressure and temperatureevolution in sedimentary basins, Geophys. J.Int. , , 601-613.[13] Yang, X. S., 2000. Nonlinear viscoelastic com-paction in sedimentary basins, Nonlinear Proc.Geophysics ,7