A Unified Approach to Mechanical Compaction, Pressure Solution, Mineral Reactions and the Temperature Distribution in Hydrocarbon Basins
aa r X i v : . [ phy s i c s . g e o - ph ] M a r A Unified Approach to Mechanical Compaction, Pressure Solution, MineralReactions and the Temperature Distribution in Hydrocarbon Basins
Xin-She YangDepartment of Applied Mathematics and Department of Fuel and EnergyUniversity of Leeds, LEEDS LS2 9JT, UK
Abstract
In modelling sediment compaction and mineralreactions, the rheological behaviour of sedimentsis typically considered as poroelastic or purelyviscous. In fact, compaction due to pressuresolution and mechanical processes in porous mediais far more complicated. A generalised modelof viscoelastic compaction and the smectite toillite mineral reaction in hydrocarbon basins ispresented. A one-step dehydration model of themineral reaction is assumed. The obtained non-linear governing equations are solved numericallyand different combinations of physical parametersare used to simulate realistic situations in typicalsedimentary basins. Comparison of numericalsimulations with real data has shown very goodagreement with respect to both the porosity profileand the mineral reaction.
Key words : Compaction, mineral reaction, vis-coelasticity, sedimentary basins, pressure solution.Published in
Tectonophysics , , 141-151 (2001). When a well-bore is being drilled for oil explo-ration, drilling mud is used in the hole to main-tain its integrity and safety. The mud densitymust be sufficient to prevent collapse of the hole,but not so high that hydrofracturing of the sur-rounding rock occurs. Both these effects dependon the pore fluid pressure in the rock, and drillingproblems occur in regions where abnormal porepressure occurs. Such overpressuring can substan-tially affect oil-drilling rates and even cause seriousblowouts during drilling. Therefore, an industriallyimportant objective is to predict overpressuring be-fore drilling and to identify its precursors duringdrilling. Another related objective is to predictreservoir quality and hydrocarbon migration. Anessential step to achieve such objectives is the sci-entific understanding of overpressuring mechanismsand the evolutionary history of post-depositionalsediments such as shales.The extent of compaction is strongly influencedby burial history and the lithology of sediments.The freshly deposited loosely packed sedimentstend to evolve, like an open system, towards aclosely packed grain framework during the initialstages of burial compaction and this is accom-plished by the processes of grain slippage, rota-tion, bending and brittle fracturing. Such reori-entation processes are collectively referred to as mechanical compaction , and generally take placein the first 1 - 2 km of burial. After this ini-tial porosity loss, further porosity reduction is ac-complished by processes of chemical compaction such as pressure solution (Rutter, 1983; Tada andSiever, 1989). Overpressuring may develop froma non-equilibrium compaction environment (Hed-berg, 1936; Gibson et al, 1967; Audet and Fowler,1992; Connolly, 1997; Connolly and Podladchikov,1998; Fowler and Yang, 1998). A rapidly accu-mulating basin is unable to expel pore fluids suf-ficiently rapidly due to the weight of overburdenrock. The development of overpressuring retards1ompaction, resulting in higher porosity, perme-ability and low thermal conductivity than are nor-mal for a given depth. These changes alter thestructural and stratigraphic features of sedimen-tary units and provide the potential for hydrocar-bon migration. Rieke and Chilingarian (1974) gavea detailed review on this topic. Audet and Fowler(1992) also reviewed compaction briefly and Con-nolly and Podladchikov (2000) provided a more re-cent comprehensive review.Pressure solution has been considered as an im-portant process in deformation and porosity changeduring compaction in sedimentary rocks (Bathurst,1971;Tada and Siever, 1989). Pressure solutionrefers to a process by which grains dissolve at in-tergranular contacts under non-hydrostatic stressand reprecipitate in pore spaces, thus resulting incompaction. The solubility of minerals increaseswith increasing effective stress at grain contacts.Pressure dissolution at grain contacts is thereforea compactional response of the sediment duringburial in an attempt to increase the grain contactarea so as to distribute the effective stress over alarger surface. There has been renewed intereston pressure solution mechanisms because of its im-portant role in the diagenesis of sedimentary rocksand its relation to rock deformation. Two com-mon types of pressure solutions occur in nature.The intergranular pressure solution was first stud-ied by Sorby (1863) based while the stylolitiza-tion was first described by Fuchs (1894) followedby many authors such as Stockdale (1922), Dun-nington (1954) and Heald (1955). Tada and Siever(1989) gave a comprehensive review on pressuresolution and some recent literature can be foundin recent papers (Birchwood and Turcotte, 1994;Schneider, 1996; Connolly and Podladchikov, 1998;Fowler and Yang, 1999; Yang, 2000a). Pressuresolution process is typically assumed to viscous(Weyl, 1959; Rutter, 1983; Birchwood and Tur-cotte, 1994; Fowler and Yang, 1999) and it is usu-ally referred to as viscous compaction, viscous creepor pressure solution creep. Its rheological constitu-tive relation (or compaction relation) is often writ-ten as a relationship between effective stress andstrain rate. However, the more natural and real-istic rheology shall be viscoelastic (Connolly andPodladchikov, 1998; Revil, 1999; Yang, 2000c) andthis is the approach we will use in the present pa-per.Mineral reactions are observed worldwide in sedi-mentary hydrocarbon basins. The close spatial cor-relations between smectite disappearance and illiteformation imply the existence of the smectite-illite reaction. This transformation is one of the mostimportant in clastic rocks. The reaction has re-ceived much attention but the nature of both theillite/smectite mixed-layer and the reaction mech-anism are still under discussion (Eberl and Hower1976; Velde and Vasseur, 1992; Abercrombie et al.,1994). Measured rates of this process in the lab-oratory (Eberl and Hower 1976) suggest that atelevated temperature, the transformation proceedsvery fast from the geological point of view. Onthe other hand, observations suggest that this re-action is initiated relatively suddenly at a temper-ature of 90 C , but then takes place gradually overa depth of several hundred meters, which suggestsa time scale of the order of a million years. Thetransformation can be considered a simple dehy-dration reaction. In fact, the mechanism of such amineral reaction is rather more complicated and isnot completely understood. The mineral reactionmay consist of dissolution of montmorillonite in freepore water and subsequent precipitation of silicaas illite. The following dissolution-precipitation re-actions have been proposed by Yang (2000b) andFowler (2000) Smectite dissolution M s (smectite) → [ − Si l − ] + n [ H O ] , (1) K-feldspar dissolution [K-feldspar] → [ K + l ] + [ AlO − l − ] + [ SiO l ] , (2) Illite precipitation [ K + l ] + [ AlO − l − ] + [ − Si l − ] → I s (illite) + [ SiO l ] , (3) Quartz dissolution and precipitation [ SiO l ] ↔ [quartz] , (4)where s, l denote solid and liquid phase. [ − Si − ]is an aqueous silica combination in such forms as[ − ( Si ) O ( OH ) ]. [ AlO − L − ] is only a general no-tation of the combination such as [ Al ( OH ) − ]. Thereaction rates in the above dissolution and repre-cipitation mechanism may be quite different at dif-ferent steps because dissolution is usually enhancedby pressure solution due to the increased pressureat grain contacts, while precipipation at free porespace is less enhanced by pressure (Rimstidt andBarnes, 1980). In the limiting case, the above fourstep reactions can be approximately considered asa one-step (first order) dehydration process. Thefirst-order dehydration model is a good approxima-tion in the sense of describing the extent of progress2f the overall smectite-to-illite transformation with-out much concern for its detailed geochemical fea-tures. Therefore, we represent it schematically assmectite · [ H O ] k r → illite + n H O [ H O ](free water) , (5)which was suggested in the earlier works (Eberl andHower, 1976; Velde and Vasseur, 1992; Abercromieet al., 1994). This one-step model was studiedbriefly by Yang (2000b) and Fowler (2000).Despite the importance of compaction and min-eral reactions for geological problems, the literaturedealing with quantitative modelling is not exten-sive. Although qualitative features have been re-ceived much attention in the literature, the mech-anismd are still poorly understood. A full under-standing of the mechanism requires an interdisci-plinary study involving soil mechanics, geochem-istry, geophysics and geology (Hedberg, 1936; Riekeand Chilingarian, 1974; Tada and Siever, 1989;Fowler and Yang, 1998; Connolly and Podlad-chikov,2000). This paper provides a unified com-paction relation in a form of a visco-poroelastic re-lation of Maxwell type. The nonlinear partial dif-ferential equations are then solved numerically andcompared with real borehole log data. For the convenience of investigating the effect ofsediment compaction and mineral reactions, thesediment is considered as as a compressible porousmatrix. Combining mass conservation of the porefluid with Darcy’s law leads to model equationsof the general type. Consider a matrix consist-ing of four interdispersed media: free pore wa-ter and three solid phase clay minerals, namely,quartz, water-rich montmorillonite and dehydratedillite. Let the volume fractions of the respec-tive media (montmorillonite, illite,quartz, water)be φ m , φ i , φ c , φ , so that φ m + φ i + φ c + φ = 1. Weassume that all the solids move with the same aver-aged velocity u s , while the pore water has velocity u l . The rate at which montmorillonite is trans-formed is denoted by r m , the rate of illite forma-tion is r i , and the rate at which water released fromthe mineral reaction is r w = r m − r i due to massconservation.Conservation of mass ∂φ m ∂t + ∇ · [ φ m u s ] = − r m , (6) ∂φ i ∂t + ∇ · [ φ i u s ] = r i , (7) ∂φ∂t + ∇ · ( φ u l ) = r w , (8) φ + φ c + φ i + φ m = 1 , (9)Conservation of Energy ∂∂t [ ρ l c l φT + ρ s c s (1 − φ ) T ]+ ∇ · { [ ρ s c s (1 − φ ) u s + ρ l c l φ u l ] T } = ∇ · ( K th ∇ T ) , (10)Darcy’s law φ ( u l − u s ) = − kµ ( ∇ p l + ρ l g j ) , (11)Force balance ∇· σ e − ∇ [ p l ] − ρg j = , ρ = (1 − φ ) ρ s + φρ l , (12)where the momentum equation has been written inthe similar form as Fowler (1990) and Audet andFowler (1992). This present form is also equivalentto that given by (McKenzie, 1984) except for thenotation difference. u l and u s are the velocitiesof fluid and solid matrix, k and µ are the matrixpermeability and the liquid viscosity, ρ l and ρ s arethe densities of fluid and solid matrix, σ e is theeffective stress, j is the unit vector pointing verti-cally upwards. p l is the pore pressure, and g is thegravitational acceleration. K th is thermal conduc-tivity, and T is temperature. c s and c l are specificheat of the solid and liquid, respectively. In for-mulating the mathematical model, the heat releasefrom the mineral reaction has been neglected sinceit is usually very small. In addition, a rheologicalconstitutive relation, which is often refered as com-paction law, is needed to complete this model andthis is described in the following section 2.1 in moredetail.Combing equation force balance and Darcy’s law,we have φ ( u l − u s ) = k ( φ ) µ [ ∇ . σ e − ( ρ s − ρ l )(1 − φ ) g j ] , (13)which is a reduced form of Darcy’s law and j is aunit vector pointing upward. Nonlinear compaction models have been formu-lated in two ways. The early and most commonmodels assumed an elastic or poroelastic rheology,3nd the compaction relation is of the Athy’s type p e = p e ( φ ) (Gibson, England & Hussey, 1967;Smith, 1971; Sharp, 1976; Wangen, 1992; Audetand Fowler, 1992; Fowler and Yang, 1998). Morerecent models assumes a viscous rheology with acompaction relation of the form p e = − ξ ( φ ) ∇ . u s where ξ is the bulk viscosity (Augevine and Tur-cotte, 1983; Birchwood and Turcotte, 1994; Fowlerand Yang, 1999). The poroelastic models are validfor the mechanical compaction while the viscousmodels are used to simulate the chemical com-paction such as pressure solution. We can thusgeneralise the above relations to a viscoelastic com-paction law of the 1-D Maxwell type ∂u s ∂z = − K ( φ ) Dp e Dt − ξ ( φ ) p e , DDt = ∂∂t + u s ∂∂z . (14)and K ( φ ) = C s (1 − φ ) φ , ξ = ξ ( φ φ ) n , (15)where C s is a constant or compaction index de-scribing the degree of compaction of sediments. ξ is the value of viscosity at φ = φ (the initialporosity). Equivalently, we would anticipate a vis-coelastic rheology for the medium, involving mate-rial derivatives of tensors, and some care is neededto ensure that the resulting model is frame indif-ferent. That is to say, the rheological relation ofstress-strain should be invariant under the coordi-nate transformation. This is not always guaranteeddue to the complexity of the rheological relations(Bird, Armstrong & Hassager 1977). Fortunately,for one-dimensional flow, which is always irrota-tional , the equation is invariant and all the differentequations in corotional and codeformational framesdegenerate into the same form (Yang, 2000c). Inthe one-dimensional case we will discuss below, wecan take this for granted. For a 1-D basin, the physical domain is from thebasin basement ( z = 0) to the basin top (oceanfloor) and the basin thickness varies with time t .The length-scale d is a typical length which will bedefined p = Gp e ( ρ s − ρ l ) gd , (16)where G = 1 + µ ξ is a constant describing sed-iment properties and p e is the effective pressure. The choice of scaling d is in such a way that the di-mensionless pressure p = O (1). Assigning the typ-ical thermal gradient to be γ , we also require thatΘ = ( T − T ) / ( γd ) = O (1). Meanwhile, we scale z with d , u s with ˙ m s , time t with d/ ˙ m s , perme-ability k with k , and heat conductivity K th with K where ˙ m s , k and K are set to typically ob-served values. By writing k ( φ ) = k k ∗ , z = dz ∗ , K th = K ˆ K , ..., and dropping the asterisks, wethus have ∂φ m ∂t + ∂∂z ( φ m u s ) = −R e β Θ φ m , , (17) ∂φ i ∂t + ∂ ( φ i u s ) ∂z = R (1 − a ) e β Θ φ m , (18) ∂φ∂t + ∂ ( φu l ) ∂z = a R δe β Θ φ m , (19) φ + φ i + φ c + φ m = 1; (20) ∂∂t { [ α (1 − φ ) + φ ]Θ } + ∂∂z { [ α (1 − φ ) u s + φu l ]Θ } = Λ ∂∂z ( ˆ K ∂ Θ ∂z ) , (21) φ ( u l − u s ) = λk ( φ )[ − ∂p∂z − (1 − φ )] , (22) ∂u s ∂z = − C s φ (1 − φ ) DpDt − ( φ/φ ) n Ξ p, (23)where λ = k ( ρ s − ρ l ) gµ ˙ m s , Λ = K ρ l c l ˙ m s d , Ξ = ξ ˙ m s G ( ρ s − ρ l ) gd . (24) R = k r d ˙ m s , a = n w M w M m , β = E a RT ,δ = ρ s ρ l , α = ρ s c s ρ l c l . (25)and k r and k are the reaction rate and permeabil-ity at the basin top, respectively. C s is a constantcoefficient or compaction index in the known func-tion K ( φ ) in equation (14). n w is the molar wa-ter released from one molar smectite; M w and M m are the molar weights of water and smectite. d isthe typical length scale in the basin; E a is the ac-tivation energy of the mineral reaction; T is thetemperature at the ocean floor and R is the gasconstant. In addition, ( φ /φ ) n expresses the de-pendence of the bulk viscosity ξ on porosity, whichhas been determined empirically to be n = 1 . n = 2 fora wide range of cemented rocks (Paterson, 1995).4he nondimensional parameter λ is essentiallythe ratio of hydroconductivity to the sedimenta-tion rate and thus characterises the relative timescales of fluid escape and sedimentation so that λ ≪ λ ≫ ≫
1) so that a nearly equilibrium stateis reached. The parameter R is the relative timescales of mineral reaction to sedimentation, and β is the reduced activation energy. However, the fac-tor R exp( β Θ) always appears as a combination inthe above model equations, and can be rewrittenas R exp( β Θ) = exp[ β (Θ − Θ c )] and Θ c = 1 β ln 1 R , (26)where the new parameter Θ c , which replaces R ,is a dimensionless critical temperature (with refer-ence to the surface temperature). In the followingdiscussions, we will see that the smectite to illitereaction virtually takes place in a region called the reaction window , at a depth of ∼ Θ c , with its thick-ness controlled by β .In the above derivation, we have used the re-quirement of degenerating to the poroelastic case p = ln( φ /φ ) − ( φ − φ ) when Ξ → ∞ . The consti-tutive relation for permeability k ( φ ) is nonlinear,and its typical form (Smith, 1971) is k ( φ ) = ( φφ ) m , (27)where m = 8 is a typical value for shales (Smith,1971; Audet and Fowler, 1992). However, thepresent model can use any value in the range from m = 3 to m = 11 in the numerical simulations. Ini-tially, the basin thickness is zero or h ( t = 0) = 0.During the whole evolutionary process, new sedi-ment deposits at a constant rate ˙ m s with a constantinitial porosity φ . The basement rock is assumedto be impeameable or u l = u s = 0 at z = 0 while aconstant heat flux (- ˆ K∂ Θ /∂z ) is introduced at thebasement rock. The top boundary is at constantpressure which can be taken to be zero. Thus, theboundary conditions become u l = 0 , u s = 0 , ˆ K ∂ Θ ∂z = − , at z = 0 , (28) φ = φ , φ m = φ m , φ c = φ c , p = 0 , ˙ h = ˙ m ( t ) + u s at z = h ( t ) , (29) where ˙ m ( t ) is the nondimensional sedimentationrate. As the thermal conductivity ˆ K ( φ ) varies with φ (Nield,1991), a constant heat flux will leads to avarying temperature gradient. However, we will as-sume ˆ K = const (in fitting the borehole data) be-cause of its weak dependence on porosity and givesa good approximation in reality (Schneider, 1996).It is useful to estimate these parameters byusing values taken from observations. By us-ing the typical 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 − , c s ∼
500 J Kg − K − , c l ∼ − K − , K ∼ . m − K − , T ∼
280 K, G ∼ , d ∼ , M m ∼ , E a ∼ − , C s ∼ . λ ≈ . − . − ≈ β = 2 . .
0. An initial porosity of φ = 0 . . . . λ and Λare due to the fact that the permeability k variesgreatly over a wide range where the smaller valuesof k correspond to a lower permeability and largervalues correpond to higher peameability. This widerange can thus similate different types of sedimentsin different sedimentation environments (Rieke andChilingarian, 1974). In order to solve the highly coupled non-linearequations, an implicit numerical difference methodis used. Substituting the expression for darcy flowinto the other equations, the essential equations de-scribing for porosities φ and φ m become the stan-dard non-linear parabolic form (Meek and Norbury,1982) φ t = F ( z, t, φ ) φ zz + g ( z, t, φ, φ z ) . (30)The first step gives φ n +1 / as a solution of the fol-lowing equation 2∆ t ( φ n +1 / i − φ ni )= ( 1∆ z ) F ( z i , t n +1 / , φ ni ) δ z φ n +1 / i + g ( z i , t n +1 / , φ ni , z δ z φ ni ) , (31)where δ z φ i = ( φ i +1 − φ i + φ i − ) and δ z φ i =(1 / φ i +1 − φ i − ). The second stage gives φ n +1 i
5s a solution of the following equation1∆ t ( φ n +1 i − φ ni )= ( 12(∆ z ) ) F ( z i , t n +1 / , φ n +1 / i ) δ z ( φ n +1 i + φ ni )+ g ( z i , t n +1 / , z δ z φ n +1 / i ) . (32)We used a normalized grid parameterized by thefixed domain variable Z = z/h ( t ). This will makeit possible to compare in a fixed frame the results atdifferent times and depths and for different valuesof the dimensionless parameters. This transforma-tion maps the basement of the basin to Z = 0 andthe basin top to Z = 1. Numerical results are pre-sented and explained in the following section, andare compared with the real data in the rest of thepaper. The main aim of the numerical model is to showhow the model equations behave, to test the va-lidity of the model, and to predict or simulatereal world situations. In principle, a good modelshould be able to simulate many realistic featureswhen its physical parameters are appropriately cho-sen, but in reality, the accuracy of the simulationis comprised by a lack of information about thereal system. By changing the different parameters( λ, Λ , β, Ξ , m, n and time t ) as well as the boundaryconditions ( φ and heat flux), we can in principleget a best fit to the real data although such a choiceof such parameters is not unique because differentcombinations of parameters in the parameter spacemay correspond to fit the real data equally wellwith a same range of deviations. So the choice ofthe parameter is still rather arbitrary although theknowledge of the geology is used in attempt to geta more realistic grouping of papameters.The borehole log data (designed as LOG I) fromthe DQ-151 borehole in South China Sea Basin inChina was used in the present simulations. DQ-151 has a depth of 6100 m and its porosity, volumefractions of smectite and illite, permeability andtemperature distribution were surveyed at nearly 5meter intervals. The data used below, however, hasbeen smoothed by averaging over intervals of 250m, so that only large scale features are revealed. However the profiles are still detailed enough totest the model presented in this paper.The results of the numerical simulations are com-pared with the data in Figure 1. The rescaledheight Z = z/h ( t ) varies from 0 to 1 and corre-sponds to a variation of 6100 m from the basementto the ocean floor. In simulating LOG I, the val-ues of the parameters which best fit the data were λ = 250, Ξ = 2 . m = 8 . n = 1 . β = 2 . t = 7 .
3, and Θ c = 3 .
4. By using the typical lengthscale d = 970 m, k ≈ . × − m , we can es-timate the sedimentation rate at deposition to be˙ m s ≈ . × − m/s ≈
230 m Ma − . The scal-ing time is about 3 . G = 1 . ξ ≈ × N s m − . Wecan see that the porosity near the top of the basindecreases nearly exponentially with depth. It canbe approximately expressed as φ = φ e − C s ( h − z ) , (33)where C s ≈ .
93 is the compaction index. Poros-ity reduces rapidly from 0.4 to 0.1. This expo-nential profile is consistent with the equilibriumprofile given by Fowler and Yang (1998) found us-ing asymptotic analysis of the nonlinear reaction-diffusion equation for fast compaction in sedimen-tary basins. In the dynamic balance of sedimentionand compaction, the porosity has reached equilib-rium in LOG I in the top 2100 m of the basin.Fowler and Yang (1998) also predicted that equi-librium even for a compacting basin may reach adepth of Π = d ln λm , (34)which is about 700 m in the present case. The nearequilibrium at the top means that (a) compactionin LOG I is mainly poroelastic and may reach theequilibrium state at depths much greater than thatpredicted by asymptotic analysis, and (b) no over-pressuring occurs in the shallow region above about2000 m, which is consistent with the drilling logdata. However, in much deeper regions, the poros-ity profile is no longer exponential. A transitionoccurs at a depth of about 2400 m below which theporosity becomes nearly uniform. This means thatcompaction is in a non-equilibrium state and over-pressuring starts to build up. Compaction gradu-ally becomes viscous due to the mechanism of pres-sure solution which operates at higher pressuresand temperatures than those prevailing at depthsof less than 2000 m.Figure 2 shows the computed temperature pro-file and the corresponding data. The values used6n this figure are Λ = 150 , γ = 0 .
031 (degree/m)with the other values being the same as in Figure1. It is clearly seen that the temperature profile isessentially governed by conduction and evolves ona faster time scale than porosity reduction. Thisis because ˆ K ∂ Θ ∂z ∼ − ≫ K ≈ h − z ).Figure 3 compares the computed volume frac-tions of different mineral species and with datafrom the borehole DQ-151. We can see that thevolume fractions of smectite remains virtually con-stant from the ocean floor to a depth of about 2100m, and suddenly decrease at a depth of about 2400m, and becomes negligible at a depth of about5100 m. Assuming a thermal gradient of nearly30.7 C/km, the mineral reaction commenced atnearly 83 ◦ C, and it ended at nearly 160 ◦ C. Assmectite starts to decrease, the volume fraction ofillite starts to increase and it reaches its maximumat the aforementioned depth of about 5100 m. Thereaction region has a thickness of about 2500 m anda temperature variation of nearly 77 ◦ C. The resultsin this paper are quite consistent with the generalconclusions drawn from other data (Abercrombieet al, 1994) from oceanic and sedimentary basins.For a linear temperature distribution with depth(Yang, 2000b) has analysed the nonlinear equationasymptotically and concluded that the convectionterm ∂ ( φ m u s ) /∂z ≪
1, so that we have ∂φ m ∂t ≈ − e β ( h − z − z ∗ ) φ m , (35)and thus an approximate solution for the volumefraction of smectite is φ m = φ m exp[ − β e − β ( h − z − z ∗ ) ] , (36)where z ∗ is the location of the centre of the reactionregion, and h is the total depth of the basin. Thissolution is remarkably accurate in describing thefeatures within the reaction window. The centre ofthe reaction window is z ∗ = RT E a ln ˙ m s k r d . (37)This clearly means that the higher the sedimenta-tion rate ˙ m s , the higher the critical temperature T ∗ of the mineral reaction, the deeper the reactiveregion z ∗ , and vice versa. A change of 2 ordersin sedimentation rate ( ˙ m s ) will cause a shift of Θ c by 2 (equivalently ∼ C) (with other parame-ters unchanged). In addition, the thickness of the reaction region is the order of (5 /β ) d . A typicalvalue of β ≈ . d = 970 m), orequivalently a temperature range of ∼ C, whichis quite close to the range observed in real data.
The present model of mechanical compaction, pres-sure solution, and the smectite to illite mineral re-action in hydrocarbon basins uses a rheological re-lation which incorporates both poroelastic and vis-cous effects and considers mineral reaction as a one-step dehydration mechanism in a 1-D compactingframe. The nondimensional model equations aremainly controlled by four parameters λ , the viscousparameter Ξ, the thermal conduction parameter Λ,and the reaction parameter R . The highly nonlin-ear coupled partial differential equations have beensolved using a two stage implicit method which isquite robust for the present governing equations.The numerical simulations have shown thatporosity-depth profile is near exponential at shal-low depths and then undergoes a transition to anearly uniform porosity. This is because of thelarge exponent m in the permeability law k =( φ/φ ) m , so that even if λ ≫
1, the product λk maybecome small at sufficiently large depths so thatthe hydroconductivity becomes very small and thuscompaction proceeds extremely slow (Fowler andYang, 1998). Meanwhile, since the bulk viscositydepends on the porosity in a form of ξ = ξ ( φ /φ ) n (typically n ≈ Z (Log I)250 6000 m03000 m Figure 1:scale of the sedimentation rate. Observations sug-gest that the temperature profile is purely conduc-tive and evolves on a faster time scale than porosityreduction, that is to say, Λ ≫ ∂ Θ ∂z ≈ − γ when ˆ K = const . This is consistent with the lineartemperature profile generated by our model.The smectite to illite mineral reaction is charac-terised by the reaction parameter R , which may bedefined in terms of a critical temperature Θ c . Thisstudy reveals that mechanical compaction, which iscontrolled by the strata permeability and sedimen-tation rate, is the most important geological fac-tor in altering the porosity evolutions, while chem-ical compaction, controlled by mineral reaction andpressure solution, causes only small changes in theporosity. The first-order dehydration model is agood approximation since it describes the extentof progress of the overall smectite-to-illite transfor-mation and still reproduces many essential featuresof the smectite-to-illite process if the appropriatereaction rate laws are used based on the knownphysics and chemistry from experimental studies.Naturally, more work is needed on more realisticformulation 2-D and 3-D compaction and mineralreactions in sedimentary environments. Acknowledgements : The author would like tothank the referees, especially Dr. J A D Connollyand Prof. R A Birchwood, for their very helpfulcomments and very instructive suggestions. I alsowould like to thank Prof. Andrew C Fowler for hisvery helpful direction on mathematical modelling. Z Figure 2: Z smectiteillite 6000 m03000 m6000 m03000 m Figure 3:8 igure captions
Figure 1. Comparison of computed porosity (solidcurve) with borehole data (depicted as circles)taken from the South China Sea basin. Z = z/h ( t )is the scaled height. Compaction in the top regionof the basin is essentially at equilibrium (for λ ≈ , Ξ = 2 . , m = 8 . , n = 1 . , t = 7 . λ = 250 , t = 7 . , Θ c = 3 . , β = 2 . REFERENCES
Abercrombie, H. J, Hutcheon, I. E., Bloch, J. D.& Caritat, P., 1994, Silica activity and thesmectite-illite reaction, Geology, v.22, p.539-542.Audet, D.M. & Fowler, A.C., 1992. A mathe-matical model for compaction in sedimentarybasins,
Geophys. J. Int. , , 577-590.Angevine, C. L. & Turcotte, D. L., 1983. Porosityreduction by pressure solution: A theoreticalmodel for quartz arenites, Geol. Soc. Am.Bull. , , 1129-1134.Bathurst, R. G. C., 1971. Carbonate sedimentsand their diagenesis , Elsevier, Amsterdam.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.Bird, R.B., Armstrong, R.C. & Hassager, O.,1977. Dynamics of polymeric liquids, Vol.1,John Willy & Son press.Connolly, J.A.D. and Podladchikov, Y. Y.,2000. Temperature-dependent viscoelasticcompaction and compartmentalization in sed-imentary basins, Tectonophysics , (in press) Connolly, J. A. D., 1997. Devolatilization-generated fluid pressure and deformation-propagated fluid flow during regional meta-morephism.,
J. Geophys. Res. , :18149-18173.Connolly, J.A.D. and Podladchikov, Y. Y., 1998.Compaction-driven fluid flow in viscoelasticrock. Geodinamica Acta , :55-85.Dunnington, H. V., 1954. Stylolite develop-ment post-dates rock induration, J. Sedimen.Petrol. , :27-49Eberl, D. and Hower, J., 1976, Kinetics of illiteformation, Geol. Soc. Am. Bull. , v. 87, 1326-1330.Fowler, A.C., 1990. A compaction model for melttransport in the Earth’s asthenosphere. PartI: the basic model, in
Magma Transport andStorage , ed. Ryan, M.P., John Wiley, pp. 3-14.Fowler, A. C. and Yang, X. S., 1998. Fastand Slow Compaction in Sedimentary Basins,
SIAM Jour. Appl. Math. , , 365-385.Fowler, A. C. and Yang, X. S., 1999. PressureSolution and Viscous Compaction in Sedimen-tary Basins, J. Geophys. Res. , B , 12 989-12 997.Fowler, A. C., 2000. Compaction and diagenesis,Proceeding of IMA, (in press)Fuchs, T., 1894. ¨Uber die Natur und Entstehungder Stylolithen,
Sitzungsber. Akad. Wiss.Wien. Math-Naturwiss , :928-41.Gibson, R. E., England, G. L. & Hussey, M. J. L.,1967. The theory of one-dimensional consol-idation of saturated clays, I. finite non-linearconsolidation of thin homogeneous layers, Can.Geotech. J. , , 261-273.Heald, M. T., 1955. Stylolites in sandstones, J.Geol. , :101-114Hedberg, H.D., 1936. Gravitational compactionof clays and shales, Am. J. Sci. , , 241-287.Lerche, I. 1990. Basin analysis: quantitativemethods, Vol. I, Academic Press, San Diego,California.McKenzie, D. P., 1984. The generation and com-paction of partial melts, J. Petrol. , , 713-765.9eek, P.C. & Norbury, J., 1982. Two-stage,two level finite difference schemes for no-linearparabolic equations, IMA J. Num. Anal. , ,335-356.Nield, D. A, 1991. Estimation of the stagnantthermal-conductivity of the saturated porousmedia, Int. J. Heat. Mass Trans , :1575-1576Patterson, M. S., 1973. Nonhydrostatic thermo-dynamics and its geological applications, Rev.Geophys. Space Phys. , , 355-389.Revil, A. 1999. Pervasive pressure-solution trans-fer: a poro-visco-plastic model, Geophys. Res.Lett. , :255-258item Rimstidt, J. D. & Barnes, H. L.,1980. The kinetics of silica-water reactions, Geochim. Cosmochim. Acta , , 1683-1699.Rieke, H.H. & Chilingarian, C.V., 1974. Com-paction of argillaceous sediments, Elsevier,Armsterdam, 474pp.Rutter, E. H., 1983. Pressure solution in nature,theory and experiment, J. Geol. Soc. London , :725-740Schneider, F., Potdevin, J.L., Wolf, S. andFaille, I. 1996. Mechanical and chemical com-paction model for sedimentary basin simula-tors, Tectonophysics , :307-317Sharp, J. M., 1976. Momentum and energybalance equations for compacting sediments, Math. Geol. , , 305-332.Smith, J.E., 1971. The dynamics of shale com-paction and evolution in pore-fluid pressures, Math. Geol. , , 239-263.Sorby, H. C., 1863. On the direct correlation ofmechanical and chemical forces, Proc. R. Soc.London , :583-600Stockdale, P. B., 1922. Stylolites: their natureand origin, Indiana Uni. Stud. , :1-97Tada,R. and Siever, R., 1989. Pressure solutionduring diagenesis, Ann. Rev. Earth. Planet.Sci. , :89-118.Velde, B. and Vasseur, G., 1992, Estimation ofthe diagenetic smectite illite transformation intime-temperature space, Amer. Mineral., v.77, 967-976. Wangen, M., 1992. Pressure and temperature evo-lution in sedimentary basins, Geophys. J. Int. , , 601-613.Weyl, P. K., 1959. Pressure solution and the forceof crystallization–a phenomenological theory, J. Geophys. Res. , :2001-2025Yang, X. S., 2000a. Pressure solution in sedi-mentary basins: effect of temperature gradi-ent, Earth. Planet. Sci. Lett. , :233-243Yang, X. S., 2000b. Modeling mineral reactionsin compacting sedimentary basins, Geophys.Res. Lett. , :1307-1310Yang, X. S., 2000c. Nonlinear viscoelastic com-paction in sedimentary basins, Nonlinear Pro-cess Geophys. ,7