Diagnosing of clay distribution in argillaceous sandstone by a rock physics template

The elastic properties of argillaceous sandstones are significantly controlled, however, by the perplexing distribution of dispersed, cemented and matrix (including structural and layered) clay. The corresponding rock physics models are established to investigate the influence of different clay distribution on the elastic properties of sandstone. The rock physics modelling and laboratory experimental results exhibit that the higher the content of the matrix clay, the lower the elastic wave velocity. The dispersed and cemented clay increases the sandstone's velocity by reducing the porosity; the increase of dispersed clay only causes a slight increase in the elastic wave velocity in sandstone. In contrast, a small amount of cemented clay can significantly increase the velocity in sandstone. Based on these understandings, we construct a cross‐plot of P‐wave velocity and clay content as a template to diagnose clay distribution. Based on the rock physics model and empirical knowledge, we divide the template into four zones, namely, matrix, dispersed, cemented and mixed clay distribution zones. Then, we use the published experimental data and numerical modelling data to validate the template. Based on diagnosing the clay distribution, we introduce how to use the diagnostic results of the template to select a suitable rock physics model for argillaceous sandstone in Shengli Oil Field, East China. The selected rock physics model, guided by diagnosing the clay distribution, predicts P‐ and S‐wave velocity very well. The proposed rock physics template for diagnosing clay distribution also has potential application in well‐logging data interpretation, rock physics modelling and reservoir characterization.

The clay distribution has significant influence on the elastic properties of argillaceous sandstone (Minear, 1982;Sams and Andrea, 2001;Tyagi et al., 2009;Knackstedt et al., 2010;Andrea et al., 2003;Cosenza et al., 2014;Ali et al., 2016). The previous studies show that the distribution of clay in argillaceous sandstone could be divided into three types (Thomas and Stieber, 1975;Minear, 1982;Han et al., 1986;Chu, 1987;Sams and Andrea, 2001;Ali et al., 2016;Zhao et al., 2020): dispersed clay, cemented clay and matrix clay. The matrix here is rock physical rather than the geological definition, which indicates the solid component of the rock. The dispersed clay is distributed in the pore space discretely. The increase in the dispersed clay will reduce the porosity, and it is often regarded as suspended matter in pore fluid. The cemented clay can be regarded as a special kind of dispersed clay distributed at the contact of the matrix particles and bonds the independent particles. A small amount of cemented clay can significantly increase the elastic wave velocity of argillaceous sandstone (Dvorkin et al., 1994(Dvorkin et al., , 1996. With its elastic behaviour similar to the matrix, matrix clay can bear the overburden pressure of the formation and does not reduce the intergranular porosity. The matrix clay can be subdivided into discretely distributed structural clay and continuously distributed layered clay. The classification of the clay in argillaceous sandstones is shown in Figure 1. Furthermore, the clay distribution is closely related to diagenesis. The diagenetic processes control different distributions of the clay in the rock. The dispersed and cemented clay usually occurs in unconsolidated sandstone with weak diagen-esis (Sams and Andrea, 2001). The rock skeleton is supported by matrix particles in unconsolidated sandstone, as shown in Figure 2(a). It is easy to distinguish the skeleton from the dispersed and cemented clay and determine the coordination number, while the structural clay usually exists in consolidated sandstone with strong diagenesis. For well-consolidated sandstones, due to long-term compaction, cementation and pressure solution, the framework is almost continuous, so it is difficult to distinguish the skeleton from the structural clay and determine the coordination number. The pores and clay can be considered as inclusions in the framework, as shown in Figure 2(b). However, the layered clay exists in all diagenetic stages and is mainly controlled by sedimentation.
In the study of elastic properties of the argillaceous sandstone, Anstey (1991) indicates that the dispersed clay has little influence on the bulk modulus of the matrix but increases the density, leading to a decrease in the P-wave velocity. Minear (1982) used different rock physics models to calculate the effects of structural, layered and dispersed clay on the elastic properties of the argillaceous sandstones. The calculation results of the equivalent medium model show that the structural and layered clay have similar effects on the P-wave velocity of the argillaceous sandstones, and the elastic wave velocity decreases with the increase of the clay content (Minear, 1982;Sams and Andrea, 2001). For dispersed clay, it is generally considered to be a part of the pore fluids. However, the calculation results of the dispersed clay model show that the clay content has a weak influence on the velocity of the argillaceous sandstones (Minear, 1982;Sams and Andrea, 2001). Sams and Andrea (2001) establish different rock physics models to systematically analyse the effects of the dispersed, layered and structural clay on the elastic wave velocity of the argillaceous sandstones. In the study of clay distribution diagnosis, Thomas and Stieber (1975) presented a volume template to diagnose the clay distribution of argillaceous sandstones using the relationship between porosity and clay content. Dejtrakulwong and Mavko (2016) proposed a meshbased fluid substitution method for interbedded shaly sands based on the Thomas-Stieber model. Saxena et al. (2019) proposed a solid and mineral substitution model to investigate the elastic properties of layered, structural and dispersed clay based on the Thomas-Stieber volume model. Besides, some laboratory experimental methods can also be used to diagnose the clay distributions on a micro-scale, such as thin-section analysis, scanning electron microscopy and computerized tomography imaging. However, the classification of the clay distribution often depends on the subjective cognition of the researchers.

Figure 2
The schematic diagram of unconsolidated and consolidated sandstone, where the grey, blue and black represent the quartz skeleton, the pores and clay, respectively. (a) The unconsolidated sandstone with weak diagenesis, in which the cemented clay (connecting two or more grains) and dispersed clay (connecting only one grain) exist; (b) the consolidated sandstone with strong diagenesis, in which the structural clay exists. In this paper, we propose a rock physics template to diagnose the distribution of the clay for argillaceous sandstone. First, we investigate the effect of different types of clay distributions on the elastic properties of argillaceous sandstone based on rock physics models. Second, we propose a rock physics template to diagnose the clay distribution based on rock physics modelling and empirical knowledge. Third, we introduce a case study on using the clay distribution template to guide rock physics modelling. Finally, discussions and conclusions are presented.

Matrix clay
As mentioned above, the matrix clay can be subdivided into structural and layered clay according to its macroscopic distribution. The structural clay is usually deposited as particles or debris (Kurniawan, 2005), often occurs as small sludges of the same size almost as sand particles (Sams and Andrea, 2001). The structural clay is a part of the matrix shown in Figure 3(a), which can also be considered as part of the sandstone matrix replaced by the structural clay; it does not affect the intergranular pores. The non-porous sandstone replaced by the porous clay results in an increase in total porosity (Thomas and Stieber, 1975;Dejtrakulwong and Mavko, 2016;Saxena et al., 2019). However, the pores in clay cannot store fluids. So, we treat the clay as an equivalent medium without considering porosity and pore structure inside the clay. Therefore, the porosity in this paper refers to effective porosity. From a macroscopic point of view, the layered clay can be regarded as a continuous structural clay. The characterization of the layered clay is usually related to its scale. There is an apparent difference among the classification of the layered clay Table 1 Rock physics modelling parameters (Sams and Andrea, 2001;Mavko et al., 2009)

Properties Values
Bulk modulus of solid frame K s (GPa) 37 Shear modulus of solid frame G s (GPa) 44 The density of solid frame ρ s (g/cm 3 ) 2 . 6 5 The bulk modulus of clay K c (GPa) 25 Shear modulus of clay G c (GPa) 9 The density of clay ρ c (GPa) 2.55 The bulk modulus of fluid K fl (GPa) 2.18 The density of fluid ρ fl (g/cm 3 ) 1 . 0 The aspect ratio of the structural clay α 0.05 The aspect ratio of pores α p 0.2 Critical porosity φ 0 (%) 0.36-0.4 Coordination number C 9-12 Hydrostatic pressure P (MPa) 15 on a seismic, logging and laboratory scale. Hence, the scale factor should be considered when it comes to a strict definition of layered clay. The layered clay can lead to anisotropy. This paper only considers the elastic properties of the horizontally layered argillaceous sandstone in the vertical direction, as shown in Figure 3(b).

Figure 4
The relations between the P-and S-wave velocities and the clay content modelled by the structural and layered argillaceous sandstone model. The subscript s and l represent the structural and layered argillaceous sandstone, respectively. The red, blue and green represent porosity of 5%, 15% and 25%, respectively. The solid and dashed lines with different colors represent the P-and S-wave velocity of the structural clay sandstone with different porosities, respectively. The dotted and dotted dashed lines with different colors represent the P-ans S-wave velocity of the layered clay sandstone with different porosities, respectively.

Figure 5
The diagram of (a) the dispersed clay and (b) the cemented clay.
Since the structural clay is distributed in wellconsolidated sandstone, it can be considered as the inclusions in the framework. Therefore, it is adopted to use the effective medium model (EMT) to study the effect of the structural clay on the elastic properties of the argillaceous sandstones. The commonly used EMT theories are Kuster-Toksöz equations (Kuster and Toksöz, 1974), self-consistent approximation (Budiansky, 1965;Berryman, 1981) and differential equivalent medium (DEM) theory Norris, 1985).
The DEM model is applied to rock physics modelling for structural argillaceous sandstone (see Appendix A) because the distribution of structural clay in the rock is discrete and random, which meets the assumptions of DEM. The quartz sand is chosen as the baseline matrix mineral in which the clay fills. It should be emphasized that the elastic parameters of the clay we quote are equivalent properties of wet clay, and hence we do not need to consider bond water and porosity in clay.
The Backus model is adopted to calculate the effective bulk and shear modulus for layered argillaceous sandstone (see Appendix B). Both the layered sand and clay are assumed to be horizontally layered isotropic mediums. The relative volume content of the sand and clay layer is set to 0.5.
We first use the above models to calculate the effective dry rock modulus of sand and clay layer. Then, the Gassmann equation is applied to calculate water-saturated elastic modulus, P-and S-wave velocities. The rock physics modelling parameters are shown in Table 1.

Figure 6
The relations between water-saturated P-wave velocity and cemented and dispersed clay content were calculated by rock physics models. The solid red and green lines represent the cemented and dispersed clay sandstone, respectively. The black arrow indicates the direction of the increase in porosity. Figure 4 shows the effects of the structural and layered clay on the elastic properties of the water-saturated argillaceous sandstones by rock physics modelling. Both the P-and Swave velocities decrease with the increase of the clay content, because the increase of the content of the matrix clay decreases the rigidity of the matrix and softens the rock skeleton. Besides the layered and structural clay have almost the same effects on the elastic properties of argillaceous sandstone. Therefore, it is still challenging to distinguish the layered and structural clay based on the elastic wave velocities.

Dispersed clay and cemented clay
The presence of the dispersed and cemented clay will reduce the porosity. The size and shape of the dispersed clay particle depend on the chemical reaction of minerals in the formation (Eshimokhai et al., 2011). It usually happens at the end stage of the deposition process, during which the dispersed clay has remained in the pore space. The diagram of the dispersed clay is shown in Figure 5(a). The cemented clay exists in the particle contact zone (Dvorkin et al., 1994(Dvorkin et al., , 1996, which bonds two independent particles and increases contact stiffness, as shown in Figure 5(b).
Since the dispersed and cemented clay are distributed in unconsolidated sandstone, they can be considered as the pore fillings (Thomas and Stieber, 1975;Dejtrakulwong and Mavko, 2016;Saxena et al., 2019). There is a quantitative relationship between the content of the dispersed clay and porosity expressed as V sh + φ = φ 0 , where V sh , φ and φ 0 represent the clay content, effective porosity and the critical porosity, respectively (Nur et al., 1995). The value of φ 0 usually ranges from 0.38 to 0.42 depending on sorting and grain shape.
The unconsolidated sandstone model (USM; Dvorkin et al., 1996) is applied to rock physics modelling for dispersed argillaceous sandstone (see Appendix C). In contrast, the contact cement theory (Dvorkin et al., 1994;Guo and Han, 2016;Han et al., 2020) (see Appendix D) is used to calculate the elastic properties of the cemented argillaceous sandstone. The input parameters of rock physics modelling are also shown in Table 1. The influence of dispersed and cemented clay on the elastic properties under water-saturated conditions of sandstone is shown in Figure 6. Figure 6 indicates that the P-wave velocity rises slowly with the increase in the content of the dispersed clay (green line), which is consistent with previous studies (Minear, 1982;Sams and Andrea, 2001). For dispersed clayey sandstone (green curve in Fig. 6), when the clay fills nearly 40% of the pores, its velocity increases rapidly, driven by the Hashin-Shtrikman's lower bound (Hashin and Shtrikman, 1962). For the cemented clay, a small amount of clay increases the P-wave velocity significantly (red curve in Fig. 6). As the clay content increases, the velocity increases slowly. Because the cemented clay is distributed far away from the particle contact and has less influence on the elastic properties of the framework. The effects of the cemented clay and dispersed clay on the elastic wave velocity of argillaceous sandstone are pretty different. Therefore, it is easy to distinguish them.

D I AG N O S I N G O F C L AY D I S T R I B U T I O N B Y A RO C K P H Y S I C S T E M P L AT E
The rock physics template to diagnose the clay distribution can be established based on the above analysis. In Figure 7, the template is the cross-plot of P-wave velocity and clay content. The upper bound (solid red line) is predicted by the contact cement theory (CCT) model, while the lower bound (solid green line) is predicted by the unconsolidated sandstone model (USM) model. The two red dashed lines are the reference lines, which are 1.1 and 0.9 times the values of the solid red line, respectively. Similarly, the two green dashed lines are also the reference lines representing 1.2 and 0.8 times the values of the solid green line, respectively. These reference lines provide constraints to differentiate clay. The values of the reference lines are empirical and chosen according to the actual situation. We can divide the template into four zones corresponding to the cemented clay (marked in blue colour), dispersed clay (marked in green colour), matrix clay (marked in red colour) and mixed clay (marked in black colour), respectively.

Cemented clay zone
In the cemented clay zone, its clay content is less than 5% and its P-wave velocity is near the upper bound. We divide this zone based on the CCT, so that we can accurately describe the elastic properties of cemented argillaceous sandsotne.
We divide this zone based on the CCT, so that we can accurately describe the elastic properties of cemented argillaceous sandstone. The actual cemented argillaceous sandstone formation usually has high velocity, low clay content but high porosity. The typical cemented sandstone data are placed on the template (Dvorkin et al., 1994) and distribute in the cemented clay zone. However, the cemented clay zone is limited by clay content because CCT is not valid for excessive cement content. Therefore, we set the threshold of clay content as 5%. This value is empirical and should be chosen according to the actual situation.

Dispersed clay zone
In the dispersed clay zone, the P-wave velocity is near the USM lines. Therefore, we divide this zone based on the USM, so that we can accurately describe the elastic properties for dispersed argillaceous sandstone. In the actual dispersed argilla- ceous sandstone formation, the dispersed clay usually has low velocity. The typical dispersed argillaceous sandstone data are also placed on the template (Marion et al., 1992;Han et al., 2013) and falls on the dispersed clay zone.

Matrix clay zone
In the matrix clay zone, its P-wave velocities are also distributed near the CCT lines, but its clay content is more than 5%. Unlike the cemented and dispersed clay, the data of the matrix clay zone are not driven by the rock physics model but the empirical understanding of rock physics.
The CCT is unsuitable for the equivalent medium rock physics modelling such as matrix argillaceous sandstone. Even though the CCT can still output a calculation result, these results are numerically close to the P-wave velocity of the matrix argillaceous sandstone. Hence, we consider the data distribute in this zone are related to the matrix clay. It should be noted that the CCT lines in this zone are only used as references and should not be over-interpreted for their physical meaning. The data of typical matrix argillaceous sandstone are placed on the template. As a result, they are mainly distributed in the matrix clay zone. These data include the actual published data (Han et al., 1986) and the numerical model data used in Figure 4.
Besides, we note that there are some data of structural argillaceous sandstones which velocities distribute in the dispersed clay zone. Because, the properties of the frame work may have been transformed from quartz-grain supported to clay supported due to the high clay content (larger than 30%).

Mixed clay zone
The mixed clay zone is surrounded by the above three zones. It indicates that the clay distribution may be a mixture of dispersed clay, cemented clay and matrix clay.

The validated data
We use some published data to validate the template. The data of structural argillaceous sandstone include 1. Han's data (Han et al., 1986) marked in circles are watersaturated structural argillaceous sandstone. 2. Numerical modelling data marked in black points are water-saturated structural argillaceous sandstone.
The data of dispersed argillaceous sandstone include 1. Han's data (Han et al., 2013) marked in symbols '♦' are artificial quartz sands under dry conditions, representing the dispersed clay sandstone.

Figure 8
The cross-plot of P-wave velocity, clay content and porosity. (Marion et al., 1992) marked in symbols '✢' are water-saturated artificial cemented quartz sands, representing the dispersed clay sandstone. 3. Numerical modelling data marked in the green line are water-saturated dispersed argillaceous sandstone.

Marion's data
The data of cemented argillaceous sandstone include 1. Dvorkin's data (Dvorkin et al., 1994) marked in symbols ' × ' are artificial epoxy-cemented glass beads used to verify the CCT model under dry conditions, representing the cemented clay sandstone. 2. Numerical modelling data marked in the red line (clay content less than 5%) are water-saturated cemented argillaceous sandstone.

A CA S E S T U DY O N D I AG N O S I N G O F C L AY D I S T R I B U T I O N S U S I N G T H E RO C K P H Y S I C S T E M P L AT E
In this section, we show how to apply our template to guide rock physics modelling. The target is argillaceous sandstone of Guantao Formation, Shengli oilfield, East China. The relationship between P-wave velocity, clay content and porosity is shown in Figure 8. The Guantao formation consists of loose and unconsolidated sandstone. The buried depth ranges from 1200 to 1800 m, and the porosity is between 25% and 35%. The cross-plot of velocity, porosity, and clay content shows nearly no statistical relations between them. In this case study, we aim to diagnose the distributions of the clay in sandstone to select a suitable rock physics model to predict the P-and S-wave velocity.

Figure 9
The well-logging data and the predicted P-and S-wave velocities using the HM model. (a) The cross-plot of P-and S-wave velocity and clay content, where '+' and 'o' represent the predicted Pwave velocity and well-logging data, while the ' ' and ' ' represent the predicted S-wave velocity and well logging data, respectively; (b) the cross-plot of the predicted and actual logging P-wave velocities; (c) the cross-plot of the predicted and actual logging S-wave velocities.

Figure 10
The diagnosis of the clay distributions of the logging data using our rock physics template.
Usually, the granular medium model (Mavko et al., 2009) is suitable for high porosity and loose sandstone, for example, the Hertz-Mindlin (HM) model (Mindlin, 1949) (see Appendix E). Figure 9 displays the predicted P-and S-wave velocity using the HM model without considering the distribution of the clay in the rock, which shows that the predicted velocities are lower than the well-logging data. Coordination number C in the HM model could be increased to achieve the best fitting of the predicted and well-logging data. However, coordination number C is restricted by porosity. The higher the coordination number is the lower the porosity will be. The porosity of the target layer is 25-35%, the value of coordination number C should be 9-12 (Mavko et al., 2009). If we let the value of the coordination number be 12, in such a case the predicted velocities are still lower than the well-logging data, which means the HM model without taking the clay distributions in the sandstone into account cannot accurately predict the elastic wave velocities of the target sandstones. We illustrate the well-logging data of the rock physics template in Figure 10 to diagnose the clay distribution.
The diagnosed results show that the clay distribution of target formation is divergent, including the cemented, dispersed clay and some mixtures between them. The distribution of well-logging data on the template behaves like the characteristics of the typical weakly cemented sandstone formation; that is, the degree of cementation should be the main factor affecting the P-wave velocity. Therefore, a suitable model was chosen for weakly cemented formation, the constant cement model (CCM) (Avseth et al., 2005(Avseth et al., , 2014) (see Appendix F),

Figure 11
The well-logging data and the predicted P-and S-wave velocities using CCM. (a) The cross-plot of P-and S-wave velocity and clay content, where '+' and 'o' represent the predicted P-wave velocity and well-logging data, while the ' ' and ' ' represent the predicted S-wave velocity and well logging data, respectively; (b) the cross-plot of the predicted and well logging P-wave velocities; (c) the cross-plot of the predicted and well logging S-wave velocities.

Figure 12
The comparisons of the elastic velocities predicted by the CCM and HM model. The black curves are well-logging data, the blue curves are P-and S-wave velocities predicted by CCM and the red curves are P-and S-wave velocities predicted by CCM. The shale zone is not predicted.
which characterizes not only the cemented and dispersed clay but also the mixtures of them. Figure 11 shows the effects of the predicted velocities using CCM. Figure 12 shows the comparisons of P-and S-wave velocities predicted by the CCM and HM model. The first column shows the porosity, the second column shows the clay content and the third and fourth columns represent the P-and S-wave velocity, respectively. The results predicted by CCM are perfectly in agreement with actual well-logging data. Therefore, the diagnosis of the clay distribution using a rock physics template effectively guides the selection of a suitable rock physics model for P-and S-wave velocity predictions. Thomas and Stieber (1975) presented a porosity-clay-content template to diagnose clay distribution, as shown in Figure 13. The model predicts volumetric properties of sand-clay mixture that the properties of clean sands (denoted by 'A') and low-porosity shale (denoted by 'C') are known. The point 'B' represents a dispersed sand point, and it represents all the soft clay that completely fills the pores of the sand. As for point D, all sand grains are replaced by load-bearing clay clast without disturbing the packing, and hence it represents the so-called structural clay point (Saxena et al., 2019). The trend line AB represents the decrease in the porosity as the dispersed clay or cemented clay filling the pores. The trend line AD shows the increase in the total porosity as the structural clay is replacing the matrix. While the trend line AC represents rock whose total porosity and the volume of clay fall on the laminated trend. This template could be used to diagnose the clay distribution qualitatively. For example, point 'E' can be interpreted as a laminar clayey sandstone where the net to gross is 60%, with the dispersed clay filling 40% of the pore space in the sand layers. Furthermore, the laminated trend accurately describes the total thickness of the lamination but carries no information on the precise distribution of sand-shale layers (Saxena et al., 2019).

Figure 13
Porosity versus shale content using the Thomas-Stieber model (Mavko et al., 2009). V_(disp shale) is the volume of the dispersed shale in the sand pore space, and φ_(clean sand) is the porosity of the clean sands. The minimum porosity occurs when the sand porosity is completely filled with shale, that is V_(disp shale) = φ_(clean sand). The N/G (net-to-gross) represents the fraction of sand in the laminated rock.
The premise of applying the Thomas-Stieber template (TST) is to know the porosity in sand and shale. However, to obtain the accurate porosity of shale is a challenging task even in laboratory conditions. Therefore, diagnosing clay distribution with TST has specific limitations. The TST cannot distinguish between the dispersed and cemented clay because they both behave as pore-filling clay to reduce the intergranular porosity. Compared with TST, the proposed template in this paper has fewer limitations in application. What is required to know is only the P-wave velocity and content of shale. Because the elastic wave velocity of shale can be easily obtained compared with the porosity of shale, our template is easier to use. However, our template cannot distinguish between structural and layered clay because their effects on elastic properties are similar. So, we recommend combining of our template and TST to obtain a more assured result of clay distribution.
There may be some potential ambiguities about clay distribution and some other diagenesis, such as quartz cementation. For example, in Figure 7, scattered data of Han's data (Han et al., 1986) are distributed in the cemented clay area. Actually, Han's sandstones have quartz rather than clay cement. It is inappropriate to interpret them as clayey-cemented sandstone. So, it is not sufficient to use only clay distribution to explain all the differences in P-wave velocity. However, most of Han's data are from structural clayey sandstone, consistent with the diagnosed conclusions.
In addition, diagnosing the diagenetic processes of sandstone is exceptionally complex. For example, the over-coating clay may hinder the secondary growth of quartz and inhibit quartz cementation, while structural or layered clay can enhance quartz cementation via stylolite generation. So, inferring these complex geological diageneses according to elastic wave velocity is still a great challenge.

C O N C L U S I O N S
Based on rock physics models, the effects of different types of clay distribution on the elastic properties of argillaceous sandstone were analysed. The structural and layered clay have similar effects on the elastic properties of sandstone, and the higher the content of matrix clay the lower the velocity of the argillaceous sandstone. The dispersed and cemented clay increases the velocity of the sandstone by reducing the porosity. The dispersed clay has a weak influence on the velocity of the sandstone. In contrast, the cemented clay significantly increases the velocity of the sandstone.
Based on these understandings, the rock physics template is proposed to diagnose the clay distribution and validated by the published data. A case study in Shengli oilfield, East China, shows that the rock physics template can effectively guide rock physics modelling for P-and S-wave velocity prediction.

AC K N OW L E D G E M E N T S
We would like to thank the anonymous reviewers for their valuable comments and suggestions to make this paper greatly improved.

A P P E N D I X A D I F F E R E N T I A L E QU I VA L E N T M E D I U M
The structural argillaceous sandstone could be modelled using differential equivalent medium (Norris et al., 1985a;Norris, 1985): In the above formula, K * (0)=K 1 , μ * (0)=μ 1 are used as the initial conditions of the coupled differential equations (A1) and (A2). The K 1 and μ 1 represent the bulk and shear modulus of the background medium (the first phase), respectively. The K 2 and μ 2 represent the bulk and shear modulus of inclusions medium (the second phase), respectively. The y represents the volume fraction of phase 2. The P and Q are the geometric factors used to characterize the geometry of inclusions. The superscript ( * 2) refers to the geometric factors for inclusion material 2 in the main phase background medium with equivalent modulus K * and G * . In this paper, the geometric shape of inclusion is ellipsoid, and the expression as follows: Among them, the tensor T ii j j relates the uniform far-field strain field to the strain in the ellipsoidal inclusion (Wu, 1966). The relevant scalars needed to calculate P and Q are as follows (Mavko et al., 2009): where For oblate spheres, the aspect ratio of inclusions α < 1, then the function of θ is where f = α 2 1 − α 2 (3θ − 2) . (A20)

A P P E N D I X B BAC K U S M O D E L
Assuming that each pure sandstone layer and shale layer is isotropic, then the equivalent medium is also isotropic, and the equivalent medium properties could be described as (Backus, 1962) where λ and μ represent the Lame coefficient and the shear modulus, respectively. Angle brackets represent the weighted average of the corresponding media properties in the brackets according to the volume ratio. The expressions of the P-wave and S-wave velocity in the equivalent anisotropic medium are where V P,h represents the P-wave velocity propagating in the horizontal direction and V SH,h represents the S-wave velocity propagating in the horizontal direction and polarized in the vertical direction. The ρ represents the average density which can be weighted by the volume fraction.
In the above formula, the C represents the coordination number (average contact number of particles). The ϕ 0 represents the critical porosity, and P represents the effective pressure. The G and v represent the shear modulus and Poisson's ratio of the framework mineral particles, respectively. The K eff and G eff represent the effective bulk and shear modulus for unconsolidated sandstone, respectively. The ϕ represents the porosity, while the K s and G s represent the bulk and shear modulus of the solid framework of the argillaceous sandstone, respectively.
In the above formulas, the ρ c is the density of the cement. The V PC and V SC represent the P-wave and S-wave velocity of the cement, respectively. The S n and S τ represent the normal and tangential contact stiffness between particles, respectively. The G and v represent the shear modulus and Poisson's ratio of the framework particles, respectively. The G c and v c represent the shear modulus and Poisson's ratio of the cement, respectively. The α represents the cement content, which could be calculated by the thickness of the contact cement zone a and the radius of the particles R. When it is assumed that the reduction of the pores of sandstone is only related to the cement, which is distributed at the contact area of the two independent particles in a specific way, the relationship between the cement content α and the existing porosity ϕ could be established as follows: where C represents the coordination number, ϕ 0 represents the critical porosity while ϕ represents the existing porosity.

A P P E N D I X E H E RT Z -M I N D L I N M O D E L
The Hertz-Mindlin model (Mindlin, 1949) could describe the elastic properties of loose-packing sandstone and cemented sandstone. The equations are shown as follows: In the above formula, the C represents the coordination number (average contact number of particles). The ϕ 0