aa r X i v : . [ m a t h . A P ] M a r Asymptotic Solutions of Compaction in Porous Media
Xin-She YangDepartment of Applied Mathematics and Department of Fuel and EnergyUniversity of Leeds, LEEDS LS6 9JT, UK
Abstract
Compaction in reactive porous media is modelled as a reaction-diffusion process with a mov-ing boundary. Asymptotic analysis is used to find solutions for the coupled nonlinear compactionequations, and a traveling wave solution is obtained above the reaction zone.Keywords: Reaction-diffusion, Darcy flow, asymptotic analysis, porous media.
Citation detail:
X. S. Yang, Asymptotic solutions of compaction in porous media,
AppliedMathematics Letters , , 765-768 (2001). The accurate modelling of compaction and reactive flow in porous media such as sand and shalesediments is very important in civil engineering and oil industry [1]. The general mathematical modelof compaction and reactive mineral flow considers the fluid-sediment system as a porous mediumconsisting of multiple mineral species. The interstitial volume of the porous solid phase is saturatedwith pore fluids. Due to the action of gravitational overburden loading and the density differencebetween the two phases, the solid phase compacts by reducing the porosity (volume fraction of thepores), thus leading to the expulsion of the pore fluid out of porous media. During compaction andcontinuous burial, the mineral species such as water-rich clay smectite react with pore fluids andare then transformed into a more stable mineral species such as illite, releasing free water into theporous environment[1,2]. In this paper, a reaction-diffusion model together with some asymptoticanalysis is presented.
Let the volume fractions of solid reactant species (smectite) and water be ψ, φ , respectively. Byproper non-dimensionalization and appropriate scalings [3,4,5] in a 1-D basin 0 < z < h ( t ), where h ( t ) is the ocean floor and z = 0 is the basement rock, we can write down the non-dimensionalcompaction equations as ψ t = − e β ( h − z − z ∗ ) ψ − λ − φ [ ψ ( φφ ) m ( φ z − φ )] z . (1) φ t = λ [( φφ ) m ( φ z − φ )] z + a β e β ( h − z − z ∗ ) ψ, m ≥ . (2)The boundary conditions are φ z − φ = 0 , at z = 0 , (3) φ = φ , , ψ = ψ , ˙ h ( t ) = ˙ s + λ − φ ( φφ ) m ( φ z − φ ) at z = h ( t ) , (4)1here λ = O (1) and β ≫ z ∗ and the amount a = O (1) offree water released from the reaction. ˙ s is the rate of new sediment accumulation at the basin top,and thus can be taken as a prescribed constant ( ˙ s = 1). φ and ψ are the initial values of φ and ψ atthe ocean floor z = h ( t ), respectively. In addition, all the variables ( φ, ψ, z, t ) and the parameters( λ, m, φ , ψ , ˙ s ) are non-negative. The reaction term exp[ β ( h − z − z ∗ )] is only dominant within athin region of a width of O (1 /β ) near h − z − z ∗ ≈
0, in other words, the reaction region is locatedat z ≈ θ ∗ defined as θ ∗ = h − z ∗ . (5)Clearly, the present problem is a non-linear diffusion problem with a boundary h ( t ) moving at a speedof ˙ h ( t ), which can be solved numerically by using the predictor/corrector implicit finite-differencemethod. Numerical simulations [3,6] and real data observations [7] imply that the moving boundary h ( t )moving essentially at a nearly constant speed ˙ h = c , although the specific value c depends on theboundary conditions, and is yet to be determined. The constant moving boundary implies thereexist travelling wave solutions. Thus, we define a new variable by ζ = z − θ ∗ = z − h ( t ) + z ∗ , − θ ∗ ≤ ζ ≤ z ∗ , (6)so that the model equations (1) and (2) become − cψ ζ = − e − βζ ψ − λ − φ [ ψ ( φφ ) m ( φ ζ − φ )] ζ . (7) − cφ ζ = λ [( φφ ) m ( φ ζ − φ )] ζ + a β e − βζ ψ, (8)The fact that β ∼ m ≫ β ≫ O (1 /β ) at z = θ ∗ . Above this region ( z > θ ∗ ), we have ζ > − βζ ) → ζ < ψ →
0) and consequently ψ exp( − βζ ) →
0. We shall see that this is true below in equation(24) because ψ → ζ → −∞ (or η → −∞ ). Outer Solutions
In the outer region above the reaction zone ( ζ > cψ ζ = λ − φ [ ψ ( φφ ) m ( φ ζ − φ )] ζ , (9) − cφ ζ = λ [( φφ ) m ( φ ζ − φ )] ζ . (10)The integrations together with top boundary conditions (4) give ψ = ˙ sψ c − λ − φ ( φφ ) m ( φ ζ − φ ) , (11)and cφ + λ ( φφ ) m ( φ ζ − φ ) = cφ + ( c − ˙ s )(1 − φ ) , (12)whose solution can be further written as a quadrature. This solution is only valid in the regionabove the reaction zone ( z = θ ∗ ). On the other hand, the travelling solution will not be appropriate2n the region below the reaction zone ( ζ <
0) because the exponential term ( φ/φ ) m ≪ φ < φ and m ≫
1. However, we can define a typical value of φ ∗ by φ ∗ = φ e − ln mm , (13)so that φ ∼ φ ∗ in the region below the reaction zone ( ζ < φ = φ ∗ e Φ m = φ e Φ − ln mm , (14)then we have ψ t = − φ ∗ λm (1 − φ ) [ ψe Φ ( 1 m Φ z − z . (15) φ ∗ Φ t = λφ ∗ [ e Φ ( 1 m Φ z − z , (16)By using 1 /m ≪
1, the above equations becomes ψ t ≈ , (17)Φ t + λe Φ Φ z = 0 . (18)It is straightforward to write down the solution together with the boundary condition Φ z = m at z = 0. We have Φ = ln( 1 + mz mλt ) . (19)As z = h ( t ), we have Φ ∞ = ln { (1 + mh ) / (1 + mλt ) } → ln( h/λt ) = ln( c/λ ) as t → ∞ or mh ( t ) → ∞ due to dh/dt = c = const and h = 0 at t = 0. However, z = h ( t ) is not usually reached since solution(22) is below the reaction region. When z ≫
1, the above solution shall match the inner solutionsas η → −∞ . Inner Solutions
In the reaction zone, we use the stretched variables defined by η = βζ + ln β, φ = φ ∗ e Φ m = φ e Φ − ln mm , (20)so that equations (7) and (8) become − cψ η = − e − η ψ − Aλφ ∗ β (1 − φ ) [ ψe Φ ( A Φ η − η . (21) − cφ ∗ Φ η = λφ ∗ [ e Φ ( A Φ η − η + a A e − η ψ, (22)where A = βm = O (1). By using 1 /β ≪
1, equation (21) becomes cψ η = e − η ψ, (23)whose solution is ψ = C exp[ − c e − η ] , C = ψ exp[ 1 c e ( βz ∗ − ln β ) ] , (24)where C depends on c , and c will be determined later in (29) by matching. Substituting this solutioninto (22), integrating once from −∞ to η and using Φ → Φ ∞ as η → −∞ , we get cφ ∗ Φ + λφ ∗ e Φ ( A Φ η − − B = − ca A C, (25)where B = cφ ∗ Φ ∞ − λφ ∗ exp(Φ ∞ ). As η → ∞ , we have[ cφ ∗ Φ + λφ ∗ e Φ ( A Φ η − + ∞−∞ = − ca A C, (26)which implies a jump through the reaction region.3 atching
By rewriting solution (25) in terms of φ , we have approximately cφ + λ ( φφ ) m ( φ ζ − φ ) ≈ cφ ∗ Φ ∞ − λφ ∗ e Φ ∞ − C ca A , (27)and matching this to solution (12), we have cφ + ( c − ˙ s )(1 − φ ) = cφ ∗ Φ ∞ − λφ ∗ e Φ ∞ − C ca A , (28)which determines c , leading to c = ˙ s (1 − φ )1 − φ ∗ Φ ∞ + a CA − λφ ∗ e Φ ∞ (1 − φ ∗ Φ ∞ + a CA ) , (29)Since Φ ∞ is a function of c , we now have an implicit equation for c , which depends essentially onthe initial values of φ and ψ .In summary, although it is very difficult to seek directly solutions for the couple nonlinearreaction-diffusion equations, we use a hybrid method to get the asymptotic solutions and travel-ling wave solutions in different regions. The matching of these solutions can thus determine theboundary moving velocity ˙ h ( t ) = c , which shows how the reaction inside the porous media affectthe evolution of its top boundary. References D. Eberl and J. Hower, Kinetics of illite formation,
Geol. Soc. Am. Bull. , (10) 1326-1330(1976). H. J. Abercrombie, et al. , Silica activity and the smectite-illite reaction,
Geology , , 539-542(1994). X. S. Yang,
Mathematical modelling of compaction and diagenesis in sedimentary basins , PhDdissertation, Oxford University (1997). D. M. Audet and A. C. Fowler, A mathematical model for compaction in sedimentary basins,
Geophys. Jour. Int. , (3), 577-590 (1992). A. C. Fowler and X. S. Yang, Fast and slow compaction in sedimentary basins,
SIAM J. Appl.Math. , 1998, (1), 365-385 (1998). X. S. Yang, Nonlinear viscoelastic compaction in sedimentary basins,
Nonlinear Proc. Geophys. , ,1-8 (2000).7