Shallow Compaction Modeling and Upscaling: A One-Dimensional Analytical Solution and Upscaling

Natural depositional processes frequently give rise to the heterogeneous multilayer system, which is often overlooked but essential for the simulation of a geological process. The sediments undergo the large-strain process in shallow depth and the small-strain process in deep depth. With the transform matrix and Laplace transformation, a new method of solving multilayer small-strain (Terzaghi) and large-strain (Gibson) consolidations is proposed. The results from this work match the numerical results and other analytical solutions well. According to the method of transform matrix which can consider the integral properties of multilayer consolidation, a relevant upscaling method is developed. This method is more effective than the normally used weighted average method. Correspondingly, the upscaling results indicate that the upscaled properties of a multilayer system vary in the consolidation process.


INTRODUCTION
Sediment compaction involves the process of reduction in pore volume of the sediment accompanied by an increase in density (Bjørlykke et al., 2009). The physical characteristics of the sediment change after deposition due to stress from overburden (gravitational), biological, or chemical reactions. To explain the major processes for the sediment compaction and help visualize the relationship of porosity loss with depth, various models were developed over the years to better capture the compaction process. Athy's (1930aAthy's ( , 1930b compaction model illustrates that the decrease in porosity with depth is as a result of expulsion of the interstitial fluid within the pores. Hence, porosity reduction and density increase are directly proportional to the increase of overburden and tectonic stresses. Hedberg (1936) classified the process of sediment compaction into four stages: 1) mechanical rearrangement and dewatering of sediments as porosity reduces from 90 to 75%; 2) loss of adsorbed water as porosity reduces from 75 to 35%; 3) mechanical deformation as the sediment resists further compaction, leading to grain recrystallization with porosity from 35 to 10% to even below 10%; and 4) chemical readjustment stage. Athy's, Hedberg's, and Terzaghi's data were adopted for Weller 's model (1959) which states that clay particles occupy the void spaces as the non-clay particles deform and share mutual contact. In addition, Power's (1967) compaction model is based on changes in clay mineralogy with burial depth and explains clay transformation and changes in the adsorbed water content at different depths. Teodorovich and Chernov (1968) in their model explained that compaction occurs in three stages: 1) fast expulsion of a large volume of fluid takes place with initial porosity loss, from 66 to 40% for clays and sandstones, and from 56 to 40% for siltstones; 2) porosity falls sharply to approximately 20%; and 3) porosity plunges to about 7-8% for shale and sandstones, at 1,400-6,000 m burial depth. Burst's (1969) compaction model resembles the previous models of clay transformation and dehydration, and states that the amount of water in movement should constitute 10-15% of the compacted bulk volume. However, this model has not been substantiated by experimental investigation. Interestingly, Beall's model (1970) is based on his study on core data from offshore Louisiana and from high-pressure experiments on marine muds. Beall's model involves the expulsion of fluid during initial stages of sediment burial with a pore throat angle at approximately six around a depth of 1,006 m. Here, the porosity decreases at a slower rate during the third stage of compaction where the angle is > 1 Å. In Overton and Zanier's (1970) model, the compaction of sediments in four stages resembles Beall's model. This model states that in the Gulf Coast, sands and shales are difficult to distinguish on self-potential electric log at depths less than 3,000 ft due to similarities of water in them. Consequently, Overton and Zanier's model focused on the different water types in four stages at different depths.
Natural deposition processes frequently give rise to layered soil deposits with alternating or random layers, which are characterized by varying properties such as permeability, compressibility, and thickness. The deposited sediments undergo a large-strain process at shallow depth and a smallstrain process in deeper locations. The existence of these processes has been recognized in geology/geotechnical engineering to influence compaction. Small-strain mechanical compaction typically involves minimal deformation of compacted grains due to vertical load and captures the behavior of sediments buried deeply in a basin, whereas largestrain compaction involves large deformation of compacted grains due to its interaction with varying loads at shallow depths. These heterogeneous fine-grained sediments at shallow burial (<1000 m) below the seafloor experience not only large strain but also variable degrees of overpressure in their pore space as a result of disequilibrium dissipation of pore fluid (Mondol et al., 2007). Consequently, the shallow overpressure poses a significant risk to the economics and safety of hydration production and may impact hydrocarbon generation deep in a basin and hydrocarbon migration to traps during basin evolution. In fine-grained sediments, deformations related to mechanical processes are dominant in the very first kilometers of depth (Hedberg, 1936;Maltman, 1994). At greater depths and temperatures, chemically modified consolidation becomes an important porosity-reducing process (Schmid and Mcdonald, 1979;Bjørlykke et al., 1989;Bjørlykke, 1998;Bjørlykke, 1999).
A sketch map for a one-dimensional consolidation model is shown in Figure 1 (left). The surcharge, which is an additional load in the form of a concentrated force or distributed load that acts on a ground surface or inside the soil body, is applied on the top of the sediment with an infinite horizontal width and is surmounted by a certain depth of water on the top of the sediment. The sediment undergoes consolidation processes, in which water flows out from the top and/or bottom boundaries as the sediment height decreases. The top (T) and bottom (B) boundaries may be permeable (P) or impermeable (I) and, hence, can be marked as PTIB (permeable top and impermeable bottom) and PTPB (permeable top and permeable bottom). Figure 1B shows a typical consolidation curve, depicting the consolidation degree versus the time. The consolidation degree is the ratio of the settlement at time t to the final settlement. Hence, the consolidation curve captures the sediment's consolidation characteristics. Some geotechnical engineering studies such as basin modeling require the research object to be discretized into blocks, with sizes in kilometers laterally and hundreds of meters vertically. Each block is then assumed to behave according to a single (homogeneous) compaction and flow relationship, even though the material is typically heterogeneous to variable degrees. The basin modeling ignores the heterogeneity of the sediment, largestrain deformation, and fluid flow conditions that occur at smaller length-and/or time-scales than those at basin scales. This can lead to incorrect predication of shallow compaction and overpressure, and subsequently basin evolution (refer to general software in the basin simulation field, such as PetroMod). Therefore, the effects of intra-block heterogeneity must be taken into account by upscaling, which then substitutes the heterogeneous property region consisting of fine grid cells with an equivalent homogeneous region of a single coarse-grid cell having an effective property value (Jingchen 2015).
Theories have been developed by researchers to describe the large-strain and small-strain consolidation processes. However, the widely adopted theories are the Terzaghi theory (Terzaghi, 1943) for small-strain and the Gibson theory (Gibson et al., 1967) for large-strain consolidations. The Terzaghi consolidation theory is widely adopted for small-strain consolidation due to its convenience and its improved methods which are still widely adopted in geotechnical engineering and other fields (Terzaghi, 1943;ArminKauerauf, 2009). As for the large-strain consolidation, the Gibson consolidation theory is more effective (Gibson et al., 1967;Gibson et al., 1981;Gibson et al., 1982), and the equation solutions are primarily based on numerical solution. However, some analytical solutions have been provided under certain conditions (Xie and Leo, 2004;Morris and Dux, 2010).
As for the multilayer system, analytical solutions such as state space (Ai et al., 2008a), three-dimensional transfer matrix solution (Ai et al., 2008b), and differential quadrature method (Chen et al., 2005) have been developed, and numerical finite difference is also widely adopted. In addition, a great deal of research has been done on multilayer consolidation, considering small-strain and largestrain processes (Schiffman and Stein, 1970;Lee et al., 1992;Xie et al., 1999;Xie et al., 2002;Chen, 2004;Chen et al., 2005;Abbasi et al., 2007;Cai et al., 2007;Ai et al., 2008b;Geng, 2008;Ai et al., 2011). However, there are two main drawbacks associated with the previous research studies. The first is that the solutions are under restricted (some parameters such as permeability, compressibility, and height are the same for different layers), and second, none of those researches focused on upscaling and supplying integral property for the multilayer consolidation system. Our work therefore focuses on overcoming these drawbacks.
As earlier mentioned, the effects of intra-block heterogeneity must be taken into account by upscaling. Since the weighted average method commonly adopted in geological engineering for multilayer systems is presently not supported by theoretical derivation, we then implemented the transform matrix and Laplace transformation to solve the multilayer small-strain (Terzaghi) and large-strain (Gibson) consolidations. According to the method of transform matrix which considers the properties of multilayer consolidation, an upscaling method is developed. Results obtained accurately match the numerical and other analytical solutions. Hence, this method is more effective than the common weighted average method. The upscaling results indicate that the properties of multilayer systems change during the consolidation processes.

ANALYTICAL SOLUTION AND UPSCALING FOR MULTILAYER TERZAGHI CONSOLIDATION Governing Equations of Solution and Upscaling for Multilayer Terzaghi Consolidation
The Terzaghi theory for one-dimensional consolidation states that all quantifiable changes in the stress of a soil (compression, deformation, and shear resistance) are a direct result of changes in effective stress. The effective stress σ′ is related to total stress σ and the pore pressure p by the following relationship: The overpressure dissipation is described by the following equation: where C v is the coefficient of consolidation, E s is the modulus of compressibility, k is the hydraulic conductivity, c w is the unit weight of water, and u is the excess pore pressure.
(2) The soil is fully saturated. It is therefore necessary to point out that according to the principle of effective stress, total stress increment is produced by the load applied on the multilayer systems. Resultantly, the effective stress increases with excess pore pressure decrease as time goes on.
According to Eq. 4, the following equation with the variable of effective stress increment can be obtained, which will benefit our solution.
The common weighted average method generates a C v for the whole multilayer system according to The Laplace transform is a widely used integral transform with many applications in physics and engineering. Hence, Stehfest numerical inversion of Laplace transforms is adopted (Stehfest, 1960) for our study, and we utilize the transformation matrix to connect parts within the multilayer system. According to the integral multilayer transformation matrix and transformation matrix between different layers, the distribution of effective stress increments and excess pore pressure can be obtained. Moreover, results are verified with an implicit finite difference numerical solution. Figure 2 shows a schematic plot of the multilayer Terzaghi consolidation.
Consequently, the Laplace transform of Eq. 5 can be expanded to whereσ′(z, s) is the Laplace transform of σ′(z, t).
At the beginning of consolidation, according to the effective stress principle, pore pressure is equal to overburden stress, which is zero initial effective stress; thus Eq. 7 can be obtained.
Then, the general solution of the ordinary differential equation for Eq. 7 is as follows: where c 1 and c 2 are constants, β s cv . Combining Eq 9 错误!未找到引用源。 and its partial derivative about z, the following expression can be derived: When z 0: Combining Eq 10 and Eq 11 yields Eq. 12: And when z i is not zero: Then the relationship between the top surface stress and the bottom stress can be derived. When considering the equation of continuous stress, and flow conservation, between two layers, the relationship between different layers can be derived.
Combining Eq 12 and Eq 13 yields Eq. 16: Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762176 Stress distribution in the same layer and in the interface can be derived, respectively, by using Eq. 11 and Eq. 14. With the equation of each layer combined together, a transform matrix T can be obtained to express the relationship between z 0 and z z n .
Here, this article only considers the situation of the pervious top surface and impervious bottom PTIB for illustration.
The corresponding Laplace transformation: where σ is the pressure on the surface and T 11 is the value of the first column and the first row of T.
With σ′(z n , s), the stress at each upper point σ′(z, s) can be obtained by the transformation matrix. Moreover, the real stress distribution can be derived by the inverse of Laplace transformation.
As for an n-layer consolidation system, the multilayer consolidation transform matrix is The common weighted average method will lead to the following weighted average method transform matrix: When t → ∞, s → 0, β → 0, by applying the following Taylor expansion we obtain e x ∞ n 0 x n n! .
As for a multilayer transform matrix, Eq. 18 will develop into Then, the weighted average method transform matrix given in Eq. 20 will develop into It can therefore be concluded that when conductivities in different layers are nearly the same, the weighted average method can be adopted for the whole multilayer system, which is the situation of homogenous.

Verification
This analytical solution is verified against the results of Lee et al. (1992), as shown in Figure 3. The parameters used in Lee's model are shown in Table 1, and the model sketch map follows the illustration in Figure 1 (right), showing a simplified sketch map of four layers with a PTIB and without overlaying water.
To compare different upscaling methods, an implicit finite-difference numerical model is developed. This study utilizes a numerical model as a benchmark model and is verifiable with both laboratory tests and a basin model (Jingchen Zhang., 2015). Analytical results are also compared with implicit finite difference numerical solution for verification. The parameters of the three different layers are shown in Table 2. Comparison of the results is then shown in Figure 4. Hence, a conclusion can be reached that this method can be applied to multilayer Terzaghi consolidation compaction. It should be noted that the values of m vl and surcharge ensure the small Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762176 strain, and correspondingly nearly no settlement can be seen in Figure 4. Resultantly, the small changes in height are shown in Figure 5.

Comparisons of Different Upscaling Methods
Here, T represents the properties of the multilayer system; then a new consolidation coefficient C v is needed to represent the whole multilayer system. With the new C v , a new transform matrix for the multilayer consolidation T′ can be obtained. The new C v should be the one that makes the minimum of Eq. 24, to numerically match T and T′.
C v changes with time according to Eq. 25, as shown in Figure 6.
In light of the long-time consolidation, C v is set to be 1.4 × 10 −8 (m 2 /s). A homogeneous layer with the new upscaling C v can then be compared with the three-layer system. Specifically, the transform matrix for a three-layer system is FIGURE 3 | Verification with Lee-excess pore water pressure isochrones and PTIB (q u is the surcharge on the surface and H is the height of the multilayer system).
Meanwhile, the common weighted average method will lead to another C v . Here, C v 3.0093 × 10 −8 (m 2 /s) in the weighted average method is used. With this C v , a new β and a weighted average transform matrix can be derived. Comparisons between the three methods are shown in Figure 7.
According to Eq. 27, the relative error between the weighted average method or the new upscaling method and the three-layer numerical results can be obtained, as shown in Figure 8.
When it comes to 20,000 days, there is nearly no overpressure; hence, the multilayer system's characteristics are studied within 20,000 days. As can be seen from Figure 9, in the first 100 days, the weighted average method is more efficient than the upscaling method. As fluid flows out through the top surface, the whole system is determined by the properties of the first layer before pressure reduction reaches the second layer. Overpressure distributions of one homogeneous layer (same properties with the upper layer) and three nonhomogeneous layers' consolidation are the same. This can explain why the results of changing C v and the weighted average method compact faster than the real situation. The possible explanation is that bigger C v of the upper layer is applied to the whole layers with small C v characteristics. The multilayer consolidation will show the integral properties more accurately after the stimulation reaches the bottom. FIGURE 7 | Overpressure evolution comparison between the two upscaling methods (weighted average and the new method) and numerical results for 10, 100, 1,000, and 10,000 days.  Figure 10 shows the comparison results, which are consistent and proves the effectiveness of this method in solving large-strain multilayer consolidation.

Comparisons of Different Upscaling Methods
To evaluate this upscaling method, comparison with the weighted average method is carried out here, utilizing the same three-layer model in Verification. When it comes to 1,000 years, there are nearly no excess pore pressures according to the weighted average method; hence, we focus within 1,000 years. Upscaling k vo changes with time as shown in Figure 11. In consideration of the long geology process, the k o value is set to be 2.8356 × 10 −10 (m/s), and C v is 7.2337 × 10 −9 (m 2 /s). While for the weighted average method, k o is 7.1867 × 10 −10 (m/s) and C v is 1.8333 × 10 −9 (m 2 /s). The comparisons of the two upscaling methods and numerical solution are shown in Figure 12.
The whole model only shows properties of the first layer before the pressure reduction reaches the second layer. Also, the weighted average k o is closer to the first layer's k o than the upscaled k o . Hence, within the first 10 years, the results obtained through the weighted average method are closer to the numerical results than those obtained through the upscaling method. However, as a whole, the upscaling method is more effective than the weighted average method.
The multilayer system only shows properties of the place affected by the stimulation, and the integral properties will change with the increase in the affected region. This can partly explain the result of a changing upscaled k o . Consequently, this upscaling method is more accurate than the common FIGURE 10 | Overpressure evolution comparison of numerical and analytical results for the three layers with different properties for 1, 10, 100, and 1,000 years.
FIGURE 11 | lot of upscaled permeability k o changes with respect to time.
FIGURE 12 | Overpressure evolution comparison between the two upscaling methods (weighted average and the new method) and numerical results for 1, 10, 100, and 1,000 years.
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762176 weighted average method as a whole for multilayer largestrain consolidation.

CONCLUSION
This work studies the modeling and upscaling of a shallow compaction using a one-dimensional analytical solution approach. The basic law of one-dimensional basin sedimentary simulation is revealed considering the effects of intra-block heterogeneity for both small-strain (deep compaction in the basin) and large-strain (shallow compaction in the basin) consolidations. And the following conclusions are made: 1. Multilayer small-strain (Terzaghi) and large-strain (Gibson) consolidations are solved with the transform matrix and Laplace transformation. The transfer matrix can upscale the heterogeneous multilayer system into one homogenous layer; therefore, it is more convenient and effective in both physical significance and the numerical form than the common weighted average method. 2. The upscaling properties of the whole multilayer system are dynamic, and the multilayer systems only show integral properties of the place affected by the stimulation. The integral properties vary with the increase in the affected region for both small-strain and large-strain consolidations.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.