Modelling diagenetic reactions and secondary porosity generation in sandstones controlled by the advection of low‐molecular‐weight organic acids

Higher secondary porosity was observed in the centre of a sandstone unit in the Eocene Shahejie Formation fan delta front sandstones from the Bozhong Depression, Bohai Bay Basin. This differs from past studies showing secondary porosity mainly in the marginal parts of sandstones adjacent to shales. This study utilized reactive transport models involving low‐molecular‐weight organic acids (LMWOA) to discuss potential processes resulting in the contrary distribution of secondary porosity. An interface model simulating LMWOA diffusion from adjacent shales to the sandstone resulted in secondary porosity in sandstones adjacent to shales. In contrast, an advection model simulating advective transport of LMWOA parallel to the sandstone bedding successfully generated higher secondary porosity in the central part. The central part of the sandstone exhibited better grain sorting (greater depositional porosity) and significantly less early carbonate cements compared to the marginal sandstone parts. Consequently, the central part had greater porosity prior to the dissolution through LMWOA. The initially higher porosity in the central part allowed for a higher advective flux of LMWOA‐rich water and associated lower pH, resulting in decreased oligoclase saturation, higher oligoclase dissolution rates and ultimately higher secondary porosity. This study indicates that grain sorting during sediment deposition, early carbonate cementation, LMWOA production in adjacent shales, and advection processes collectively control the diagenetic reactions and the distribution of secondary porosity in sandstones.

LMWOA in sandstone reservoirs can be derived from various sources, including the thermal decomposition of kerogen in shales (Andresen et al., 1994;Barth & Bjørlykke, 1993;Helten et al., 2022), biodegradation of hydrocarbons (Meredith et al., 2000;Wartell et al., 2021) and thermal degradation of hydrocarbons (Fecteau et al., 2018;Yuan et al., 2019).This study focuses on the LMWOA produced from the thermal decomposition of kerogen in shales.After their production, LMWOA are transported by diffusion from shales to adjacent sandstone reservoirs, where they may react with undersaturated minerals, leading to the generation of secondary porosity.Consequently, higher secondary porosity is expected to be observed in the marginal sandstone parts adjacent to shales, as illustrated in Models A, B and C (Figure 1; Ehrenberg, 1991).
In contrast, a model characterized by higher secondary porosity in the central part of a sandstone unit and minor secondary porosity in the marginal parts adjacent to shales is here referred to as the reverse model (Model D, Figure 1) and has been reported earlier (e.g.Liu & Xiong, 2021;Nedkvitne & Bjørlykke, 1992;Sullivan & McBride, 1991;Wang et al., 2019).Nedkvitne and Bjørlykke (1992) and Sullivan and McBride (1991) attributed the reverse model (Model D, Figure 1) to the advection of formation water and meteoric water respectively.However, as the entire sandstone unit can be flushed by water advection, a simple explanation of water advection alone cannot fully account for this heterogeneity.While Liu and Xiong (2021) and Wang et al. (2019) described this phenomenon, they did not provide definitive explanations.Therefore, the mechanism behind Model D remains unclear.
An analogue to Model D was observed in sandstones at depths of ~3194.3-3198.9m at the well B1 from the Eocene Shahejie Formation, Bozhong Depression, Bohai Bay Basin, China (see Section 3.1.2).In the central part of this sandstone interval, the secondary porosity quantified from thin section samples is as high as 6.8 vol.%.In contrast, secondary porosity in the marginal parts is as low as 0.2 vol.%.Therefore, this study employed kinetically constrained reactive transport simulations and incorporated petrographic information from the Bohai Bay Basin.The aim was to examine the potential processes that give rise to the characteristic features of Model D. Fifty core plug samples were collected from a sandstone core sample interval at depths of ~3194.3-3198.9m in the borehole B1 of the Eocene Shahejie Formation, Bozhong Depression, Bohai Bay Basin, China.The Bohai Bay Basin is a Mesozoic-Cenozoic lacustrine petroliferous basin deposited on the Palaeozoic-Precambrian North China craton basement (Hao et al., 2009;Yang et al., 2023).The local depositional environment of the Eocene Shahejie Formation in the study area was interpreted as lacustrine fan delta front (Xu et al., 2017).The sample depth interval is in its historical maximum burial depth and temperature.The present-day burial temperature at the depth interval is ~115°C, which is calculated based on the drilling stem test (DST) temperature gradient (31.6°C/km) and a surface temperature of 14.5°C in the Bozhong Depression (Li et al., 2021).

| Experiments
All 50 sandstone core plug samples were analysed for helium porosity, point-counting volumetric contents of detrital compositions, diagenetic minerals and thin section porosity, cathodoluminescence (CL) observations, X-ray diffraction (XRD) mineralogy and laser grain size.Five thin section F I G U R E 1 Schematic representation illustrating the distribution models of secondary porosity in a sandstone reservoir in contact with shale layers above and below.Models A, B and C were constructed based on the study of Middle Jurassic Garn Formation sandstones from the mid-Norwegian continental shelf (Ehrenberg, 1991).Similar distributions have also been reported in studies by Dou et al. (2017), McMahon et al. (1992), Moncure et al. (1984) and Sullivan and McBride (1991).The absence of higher secondary porosity either in the top or bottom parts of sandstone (Models A and B) is mainly attributed to a lack of significant supply of LMWOA from adjacent shales (Ehrenberg, 1991;Moncure et al., 1984).Model D was constructed based on the study of Middle Jurassic Brent Group sandstones, Huldra field, North Sea (Nedkvitne & Bjørlykke, 1992).Similar distributions have also been reported in studies by Liu and Xiong (2021), Sullivan and McBride (1991) and Wang et al. (2019).The secondary porosity in these cases results from the dissolution of feldspars and rock fragments.The dissolution of cements (e.g.carbonates) is not considered.The sandstone reservoir is assumed to be enclosed by mature shales capable of producing LMWOA (vitrinite reflectance (R o ) >0.5).Young/shallow sandstone reservoirs enclosed by immature shales (R o < 0.5) are not considered.
samples were selected to measure thin section grain size.Fourteen samples were analysed for the homogenization temperatures (Th) of primary aqueous inclusions captured in carbonates and quartz cements.Detailed descriptions for each experimental method are given in Appendix S2.
2.1.3| Classification of early and late carbonate cements Two methods, including petrography (thin section and cathodoluminescence) and Th measurements, were employed to differentiate early and late carbonate cements.Thin section was initially used to identify different occurrences of carbonate cements, including intergranular calcite, secondary pore-filling calcite, ferrocalcite, dolomite and siderite (see Section 3.1.3).Cathodoluminescence was mainly used to aid in distinguishing potential coexisting multiple phases within carbonate cements.Subsequently, Th of carbonate cements were used to determine whether they were early or late carbonate cements.

| Evaluation of mechanical compaction
The compactional porosity loss (COPL) was calculated according to Ehrenberg (1995) where OP is the depositional porosity, and IGV is the pointcounting intergranular volume.
The OP was calculated by Scherer (1987): where S 0 is the Trask grain sorting coefficient.Moreover, vermicular and blocky kaolinites typically have an intergranular microporosity of 15%-60%, as quantified through computer-assisted image analysis using back-scattered electron micrographs (Hurst & Nadeau, 1995).Therefore, an averaged microporosity of 38% was employed to correct the point-counting kaolinite volume.

| Model overview
Reactive transport models were implemented using the Geochemist's Workbench® Pro version 13.0.Two potential mass transport processes occurring in sandstone reservoirs were modelled (Figure 2): (1) A one-dimensional (1D) diffusion-dominated model simulating the reactive transport at the shale-sandstone interface, and (2) A two-dimensional (2D) advection-dominated model simulating transport and geochemical reactions at the fringe of the advecting organic acid plume within the sandstone layer.The mass transport within shales is predominantly controlled by diffusion (Charlet et al., 2017).Therefore, only diffusion was considered in the interface model.Although some studies have suggested that the mass transport in deep sandstones is more likely dominated by diffusion (Bjørlykke, 1993;Giles, 1987), these studies did not consider fault movement.Our earlier study Li, Zhu et al., 2023) has indicated that the active fault movements in the study area can trigger significant fluid flow, as evidenced by the positive correlations between fault displacement rates and fluid accumulation and leakage amounts.This supports the potential feasibility of advection in the study area.Therefore, both advection and diffusion were employed in the 2D model.
The 1D interface model had a length of 5 m (horizontally) and consisted of 50 cells with each cell having a size of 0.1 m (Figure 2).One boundary of the model domain (X position = 0 m) represented the shale-sandstone contact (Figure 2).The 2D advection model had dimensions of 1 m (horizontally; 10 cells) × 5 m (longitudinally; 50 cells) with each cell having a size of 0.1 × 0.1 m (Figure 2) representing the reactive front of the advecting plume.
The mineralogy and porosity for each cell were derived from the analysis of the 50 sandstone core plug samples in well B1.The temperatures of all models were consistently 100°C.This temperature generally corresponds to the highest concentration of dissolved LMWOA in formation waters (Carothers & Kharaka, 1978;Fisher, 1987;Lundegard & Kharaka, 1994;MacGowan & Surdam, 1990).The potential influence factors, including the inflowing water pH, the LMWOA concentration and the advective velocity, were individually tested.Acetic acid and oxalic acid were selected as the input LMWOA species (see Section 2.2.5).The inflowing water pH was varied by choosing pH = 5, pH = 6 and pH = 7.5 to represent acidic, mildly acidic and mildly alkaline conditions at 100°C.The LMWOA concentration ranged from 200 to 1000 mg/L.Advection velocities ranging from 0.01 to 1.0 m/year were tested, which was based on the estimations of both compaction-driven (Giles, 1987) and fault-triggered fluid flows (Pedersen & Bjørlykke, 1994).

| Thermodynamic dataset
The thermodynamic parameters were mainly derived from the LLNL dataset (Delany & Lundeen, 1990), which (2) OP = 20.91 + 22.90 S 0 was mainly based on the SUPCRT dataset (Johnson et al., 1991).The thermodynamic parameters for acetic acid and oxalic acid, as well as the dissociation constants of aqueous acetate and oxalate complex species, were obtained from the MINTEQ dataset (Table S1, Apendix S2; Morrey et al., 1985) and incorporated into the LLNL dataset (Li et al., 2024).The LLNL database utilizes an extended Debye-Hückel equation to calculate the activity coefficient of aqueous species (Bethke, 2008): where γ i is the activity coefficient of aqueous species i, z i is the charge of aqueous species i, I is the ionic strength of the solution, A, B and Ḃ are temperature-dependent coefficients and å i is the ion size of species i.

| Mineralogy
The mineralogy and porosity in all models were determined from the obtained samples (see Section 3.1).The sandstone mineralogy is composed of quartz, K-feldspar, plagioclase, calcite, dolomite, siderite, kaolinite and chlorite (Table S2).Ferrocalcite was selected as the secondary mineral.Detailed mineralogy and porosity for each cell are provided in Appendix S1.The formula and thermodynamic data for chlorite (Mg 2.5 Fe 2.5 Al 2 Si 3 O 10 (OH) 8 ) were rearranged and recalculated based on an ideal binary solid solution between the two end-members (clinochlore and daphnite) (Holland & Powell, 1998).The formula and thermodynamic data for ferrocalcite (Ca .95Fe .05CO 3 ) were calculated based on the ideal binary solid solution of calcite and siderite (Duan et al., 2016).
2.2.4 | Kinetics of mineral dissolution and precipitation According to the transition state theory (TST) (Lasaga, 1984;Lasaga, 2014;Oelkers et al., 1994), the rate law for mineral dissolution and precipitation can be expressed as: where r → k is the reaction rate (mol/s), RSA is the reactive surface area (cm 2 ), k + is the intrinsic rate constant (mol/cm 2 / sec), a j is the activity of promoting or inhibiting species, P j is the reaction order with respect to species a j , Q and K are the activity product and equilibrium constant for the reaction, respectively, and ω and Ω are empirical and dimensionless orders of the rate law.
(3) The P j , ω and Ω define whether the rate law is linear or nonlinear.When P j , ω and Ω are simultaneously set to 1, a linear rate law is employed; otherwise, a nonlinear rate law is employed.The mineral dissolution rate is typically significantly lower at close-to-equilibrium conditions compared to far-from-equilibrium conditions (Black & Haese, 2014;Lasaga, 2014).A linear rate law can lead to a predicted rate that is orders of magnitude faster than the experimentally derived rate at close-to-equilibrium conditions (Black & Haese, 2014;Hellmann & Tisserand, 2006).To address this, the second-order nonlinear rate law (ω = 1 and Ω = 2) was employed for all minerals because it could substantially improve the accuracy of predicted reaction rates at close-to-equilibrium conditions (Yuan et al., 2017).
The temperature-adjusted rate constant, k + , can be calculated using the Arrhenius equation (Lasaga, 1984;Steefel & Lasaga, 1994): where k 25°C (mol/cm 2 /s) is the standard rate constant at 25°C (298.15K) and pH = 0, E a is the activation energy (J/ mol), R is the gas constant and T is the Kelvin temperature of the solution (K).
The H + concentration plays a crucial role in promoting the dissolution rate, as the protonation/hydroxylation of the mineral surface can activate surface sites (Nagy, 1995).Considering H + as the primary promoting species, the rate law can be expanded to include the acid, neutral and base mechanisms (Palandri & Kharaka, 2004): where base are the acid, neutral and base enhanced standard rate constants, respectively, E acid , E neutral and E base are the activation energies under acidic, neutral and alkaline conditions, respectively, a H + is the activity of H + and n and m are reaction orders with respect to H + .
The values of k 25°C , E a , n and m for K-feldspar, quartz, oligoclase, calcite, siderite, dolomite, kaolinite and chlorite were obtained from the compilation by Palandri and Kharaka (2004) (Table S2).Parameters of calcite were used as a proxy for ferrocalcite.A nucleus density was introduced to account for the initial surface area for the precipitation of secondary minerals (Bethke, 2008).The nucleus densities for mineral precipitation were arbitrarily set to 60 cm 2 /cm 3 for nonclay minerals (carbonates and quartz) and to 600 cm 2 /cm 3 for kaolinite (Table S2).
Several methods have been suggested to estimate the reactive surface area, RSA, based on the specific surface area (SSA), including the geometric surface area (GSA) (Steefel & Lasaga, 1994), the BET surface area (Brunauer et al., 1938) or the 2D image perimeter-based surface area (Beckingham et al., 2016).These methods can lead to 3-5 orders of magnitude differences in the estimation of RSA (Beckingham et al., 2016).Beckingham et al. (2016) compared the modelling-derived effluent concentrations using the aforementioned three RSA methods with the effluent concentrations obtained from flow-through dissolution experiments using sandstones from the Pleistocene Haizume Formation at a depth of 1093 m in the Nagaoka pilot CO 2 injection site, Japan.They suggested that the image perimeter-based SSA matched the experimental effluent concentrations.The high-BET surface area (highest BET values derived from literatures summarized in Beckingham et al. ( 2016)) provided the second-best match.Therefore, the available image perimeter SSA values reported by Beckingham et al. (2016) were employed for quartz, plagioclase, K-feldspar and calcite (Table S2).The high-BET surface area values were employed for chlorite and kaolinite.The SSA of dolomite (Luhmann et al., 2014) and siderite (Golubev et al., 2009) were derived from mineral-specific dissolution experiments (Table S2).Furthermore, a scaling factor (SF) is typically introduced to constrain RSA (Beckingham et al., 2017): where ASA is accessible surface area.Beckingham et al. (2017) suggested that the SF was a function of mineral abundance.The SF is 1-2 orders of magnitude for minerals >5 vol.%.The SF is 3-5 orders of magnitude for minerals <5 vol.%.In this study, an SF of 0.01 was introduced for minerals >5 vol.% (quartz, Kfeldspar, oligoclase, calcite and dolomite), and a SF of 0.00001 was introduced for minerals <5 vol.% (siderite, kaolinite and chlorite), as well as for secondary ferrocalcite (Table S2).
Moreover, since LMWOA are generated in shales, in the 2D advection models, the LMWOA concentration in marginal sandstones may be higher than that in middle sandstones.However, a preliminary diffusion model indicates that the LMWOA concentrations will be quickly diluted and reach a constant concentration across the sandstone column in 30-40 years, even if the original LMWOA concentrations in marginal sandstones are significantly higher than those in middle sandstones (Figure S3).Therefore, a constant LMWOA concentration across the 2D domain was used in a single simulation.

| Composition of inflowing water
Subsurface formation waters reach equilibrium with in situ diagenetic minerals at temperatures >70-80°C (Bjørlykke & Egeberg, 1993;Helgeson et al., 1993).Therefore, a reaction modelling was employed by using 1.0 M NaCl solution (pH = 6) to react with the averaged initial mineral assemblage (Table S2).The obtained water (Table 1) is in equilibrium with all diagenetic minerals presented in this study and is used as the initial water in all simulations.
The inflowing water compositions were based on the initial water compositions.H-Acetate and H 2 -Oxalate were added into the inflowing water while keeping the other components consistent with the initial water (Clcharge balanced).Moreover, the reported pH in the produced water from gas-bearing shales in American basins predominantly ranges from 5 to 9 (Figure S4a; Blondes et al., 2018).A pH in the range of 5.7-7.4 was reported earlier for the formation water in sandstones with a temperature ranging between 90 and 110°C in the Bohai Bay Basin (Figure S4b; Li, Haese et al., 2023).Therefore, three water pH values were selected, including pH = 5, pH = 6 and pH = 7.5, to represent acidic, near neutral and alkaline conditions at 100°C (Table 1).Moreover, although no additional CO 2 was introduced into the system, the pCO 2 values in the inflowing waters are 0.119 bar (pH = 5), 0.068 bar (pH = 6) and 0.004 bar (pH = 7.5).Therefore, these inflowing waters can represent different pCO 2 values that may be caused by the thermal decomposition of kerogen (Andresen et al., 1994;Helten et al., 2022).

| Molecular and mass transport
The movement of formation water and dissolved solutes in porous media is governed by three fundamental processes, including advection, molecular diffusion and mechanical dispersion (Bethke, 2008).The advective flux (q A ) of a chemical component is calculated using the following equation (Bethke, 2008) where q is the specific fluid flow velocity, and C is concentration of the chemical component.
The molecular diffusion and mechanical dispersion are combined and referred to as hydrodynamic dispersion (Bethke, 2008).The hydrodynamic dispersive flux in sediments (q D ) is calculated using Fick's law (Schulz & Zabel, 2006) where Ø is the sediment porosity (%), D is the coefficient of hydrodynamic dispersion (cm 2 /s), C is the volumetric concentration (mol/cm 3 ) and x is distance (cm).
The coefficient of hydrodynamic dispersion (D) is calculated by Bethke (2008): where a L is the longitudinal dynamic dispersivity (cm), v x is the average linear groundwater advection velocity (cm/s) and D* is the molecular diffusion coefficient in sediments (cm 2 /s).
The a L is relevant to the flow distance of the domain and reflects the scale effect (Anderson, 1976) and can be calculated using an empirical equation: a L = 0.83(logL) 2.414 , where L is the length of the flow path (cm) (Xu & Eckstein, 1995).The D* is derived from the diffusion coefficient in free solution, D sw , and accounts for tortuosity, θ, according to D * = D sw 2 (Schulz & Zabel, 2006).Tortuosity was calculated according to Boudreau (1996Boudreau ( , 1997): 2 = 1 − ln(Ø 2 ), where Ø is the sediment porosity (%).
Here, we used the diffusion coefficients in free solution as provided by Schulz and Zabel (2006) for seawater and derived the temperature-adjusted diffusion coefficient according to Park (2014) (Appendix S1).The average value (4.90 × 10 −5 cm 2 /s) of the calculated diffusion coefficients for the available species in free solutions at 100°C was employed in this study.Detailed mass transport parameters for each cell are provided in Appendix S1.

| Lithology and detrital compositions
It is anticipated that the initial porosity plays a prominent role in the development of secondary porosity.Therefore, the sampled sandstone core interval was initially divided into three intervals according to the distribution of core plug porosity (Figure 3a).The top marginal interval A (median porosity: 11.18%; N = 12) and bottom marginal interval C (median porosity: 14.95%; N = 12) exhibited lower core plug porosity, whereas the middle interval B (median porosity: 21.50%; N = 26) had higher core plug porosity (Figures 3a and 4a).In the following sections, the petrographic and model results are presented based on the division of these three intervals.
The sandstone samples in the three intervals exhibited a relatively consistent quartz-feldspar-rock fragments (Q-F-R) lithology, which was mainly composed of lithic arkose (Figure 3b).The grain size ranged from fine sandstone to conglomerate (Figure 3a).The grain texture was generally characterized as subangular to subrounded with poor to moderate grain sorting.The rock fragments mainly consisted of basic volcanic rocks.

| Thin section porosity
The total thin section porosity is composed of primary and secondary thin section porosity.In this study, only the quantification of secondary porosity resulting from the dissolution of feldspars and rock fragments was feasible using thin sections (Figure 5a-c), as the identification and verification of carbonate cement dissolution were challenging.For example, petrographic observations showed the dissolution of feldspar grains while the surrounding calcite remained intact (Figure 5b,c).The total thin section porosity exhibited a linear relationship with helium porosity (Figure S5a).The middle interval B displayed higher median total thin section porosity, higher median primary thin section porosity and higher median secondary thin section porosity compared to the marginal intervals A and C (Figure 4b-d).

| Diagenetic minerals
The identified diagenetic minerals mainly consisted of carbonate cements, quartz overgrowth and kaolinite.Carbonate cements were mainly composed of calcite, dolomite, siderite and ferrocalcite.Calcite precipitation occurred in both the intergranular pore and intragranular dissolution pore (secondary pore-filling calcite) (Figure 5b-e).Ferrocalcite precipitation occurred in both the intergranular and intragranular pores (Figure 5f).Dolomite precipitation mainly occurred in the intergranular pore (Figure 5g).Siderite mainly precipitated in the intergranular pore (Figure 5h).In general, cathodoluminescence analysis indicated that calcite exhibited relatively consistent colouration (Figure 5i,j), and dolomite also displayed relatively consistent colouration (Figure 5j).Quartz overgrowths typically precipitated on the grain surface of detrital quartz (Figure 5k,l).Kaolinite typically exhibited in booklet shapes and precipitated in both the intergranular and intragranular pores (Figure 5k,l).The ( 9) abundance of kaolinite, as an example of diagenetic minerals, was higher in the middle interval B but lower in the marginal intervals A and C (Figure 4e).The distributions of other diagenetic minerals are presented in Figure S6.

| Mechanical compaction
All samples were situated in the 'compaction-dominated' area, where the porosity destroyed by compaction was consistently >50% (Figure S7).This indicated that the mechanical compaction was the predominant mechanism leading to the porosity destruction in the studied samples.Quantitative calculation of COPL indicated that the middle interval B had higher COPL than the marginal intervals A and C (Figure 4f).The typical identification mark of chemical compaction, stylolite or suture contacts, was not observed.

| Homogenization temperature of aqueous inclusions
The observed fluid inclusions in the samples mainly consisted of aqueous and oil inclusions trapped in healed microfractures, aqueous and oil inclusions trapped in diagenetic minerals and oil inclusions trapped in the boundary between detrital quartz and quartz overgrowth (Figure S8).The size of aqueous inclusions trapped in diagenetic minerals typically ranged from 2 to 13 μm.The Th of intergranular calcite ranged from 41.8 to 72.5°C (Figure 6a).The Th of dolomite ranged from 41.2 to 69.9°C (Figure 6b).The Th of siderite ranged from 44.4 to 64.5°C (Figure 6c).The Th of ferrocalcite ranged from 86.7 to 108.7°C (Figure 6d).The Th of quartz overgrowth ranged from 87.2 to 112.3°C (Figure 6e).Detailed measurements refer to Appendix S1.

| Interface diffusion model results
In this section and the following sections, the model results from the case with inflowing water pH = 6 were mainly presented as examples.The model results from the other two cases with inflowing water pH = 5 and pH = 7.5 are provided in Appendix S2.
Upon the injection of inflowing water, the pH was rapidly buffered to a stable value of ~6.75 (Figure 7a).Oligoclase dissolution was interpreted in the following sections to represent the generation of secondary porosity.Oligoclase dissolution preferentially occurred in the sections of the formation adjacent to the source of inflowing LMWOA-rich water (shale-sandstone contact) within 20,000 years (Figure 7b).

| Fluid velocity and pH
The fluid velocity and water pH profiles exhibited significant variations.However, the fluid velocity in the middle interval B was generally higher than those in both the top and bottom marginal intervals A and C (Figure 8a,b).In contrast, the water pH was lower in the middle interval B while higher in both the top and bottom marginal intervals A and C (Figure 8c,d).In comparison, the LMWOA concentrations in the column did not show identifiable variations during the simulations (Figure S11).

| Dissolution of oligoclase and generation of secondary porosity
Oligoclase dissolution is presented here showing the generation of secondary porosity.The dissolution of K-feldspar and chlorite is presented in Appendix S2.In general, lower oligoclase saturation (Q/K) and a correspondingly greater oligoclase dissolution rate were observed in the middle interval B compared to both the top and bottom marginal intervals A and C (Figure 8e-h).As a result, greater oligoclase dissolution was observed in the middle interval B compared to both the top and bottom marginal intervals A and C (Figure 9a-g).

| Dissolution and precipitation of other minerals
Kaolinite is presented here as an example of diagenetic mineral formation, while secondary quartz, calcite, dolomite, siderite and ferrocalcite are provided in Appendix S2.In general, kaolinite demonstrated a preference of precipitation in the middle interval B, whereas lower concentrations of kaolinite were observed in both the top and bottom marginal intervals A and C (Figure 10a-g   and 10h,i).For example, the oligoclase dissolution with LMWOA = 1000 mg/L and fluid velocity = 0.1 m/year at the cell position X = 0.95 m and Y = 0.35 m was 5.43, 0.97 and 0.35 vol.% with an inflow pH of 5, 6 and 7.5 respectively (Figure 9h).Under the same conditions and at same position, the kaolinite precipitation was 3.23, 0.58 and 0.21 vol.% with an inflow pH of 5, 6 and 7.5 respectively (Figure 10h).

| LMWOA concentration
When keeping the inflowing water pH and fluid velocity constant, an increase in LMWOA concentration resulted in a slight decrease in oligoclase dissolution (Figure 9j,k) and a slight decrease in kaolinite precipitation (Figure 10j,k).However, the effect of LMWOA concentration is simultaneously influenced by the inflowing water pH.For example, the difference in oligoclase dissolution at the cell position X = 0.95 m and Y = 3.5 m between 1000 and 200 mg/L was 3.89, 0.23 and 0.01 vol.% after 1000 years with an inflow pH of 5, 6 and 7.5 respectively (Figure 9j, Figure S13A8 and S13B8).

| Fluid flow velocity
The comparisons among results using different fluid flow velocities, ranging from 0.01 to 1.0 m/year, at the position X = 0.95 m are presented as examples to discuss the influence of fluid flow velocity.In general, when keeping the LMWOA concentration and inflowing water pH constant, an increase in fluid flow velocity resulted in a significant increase in oligoclase dissolution (Figure 9l,m) and a significant increase in kaolinite precipitation (Figure 10l,m).However, the effect of fluid flow velocity is also influenced by the inflowing water pH.For example, the difference in oligoclase dissolution at the cell position X = 0.95 m and Y = 3.5 m between 0.1 and 1 m/year was 11.03, 6.19 and 2.58 vol.% after 1000 years with an inflow pH of 5, 6 and 7.5 respectively (Figure 9l, Figures S13A9 and S13B9).Secondary porosity in thin sections is linearly correlated with the volumes of both point-counting kaolinite and quartz overgrowth (Figure S5b,c).Furthermore, petrographic observations indicated that kaolinite tended to precipitate within secondary grain dissolution pores, and quartz overgrowth tended to precipitate surrounding secondary pores (Figure 5k,l).These observations agree with earlier studies (Putnis, 2009) and have been described by Ruiz-Agudo et al. (2014) as follows: "The common observation in both nature and experiment that one solid phase can be replaced by another via coupled dissolutionprecipitation suggests that this is a universal mechanism of re-equilibration during fluid-mineral interaction".In sandstone reservoirs, external sources have been suggested to have a limited role in providing Al and SiO 2 for the precipitation of kaolinite and quartz overgrowth (Bjørlykke & Jahren, 2012;Milliken, 2005).This is mainly due to the low solubilities of Al and SiO 2 in both sandstone and shale formation waters (Blondes et al., 2018;Engle et al., 2016), which significantly restricts the mobilities of Al and SiO 2 (Bjørlykke & Jahren, 2012;Milliken, 2005).Therefore, internal sources are more likely to dominate the supply of Al and SiO 2 in sandstone reservoirs.Based on the statistical correlations between the secondary porosity in thin sections and point-counting kaolinite and quartz overgrowth (Figure S5b,c), as well as petrographic observations (Figure 5k,l), it is reasonable to conclude that kaolinite and quartz overgrowth predominantly originate from grain dissolution processes.Therefore, the Th of quartz overgrowth can be used as an approximate indicator to estimate the timing of the formation of secondary porosity and kaolinite.The Th of quartz overgrowth showed a range of 87.2-112.3°C(Figure 6e).Therefore, it can be inferred that the generation of secondary porosity and the precipitation of kaolinite roughly occur within a same temperature window.Moreover, chemical compaction of detrital quartz has been suggested to occur at a temperature of ~80°C (Walderhaug et al., 2000;Worden & Morad, 2000).Therefore, although no typical identification marks (stylolite or suture contacts) were observed, it is possible that chemical compaction of detrital quartz may play a minor role in providing SiO 2 for quartz overgrowth.

| Mechanical compaction
Mechanical compaction was identified as the primary factor responsible for porosity loss (Figure 4f and Figure S6).This conclusion is consistent with previous studies on compaction in sandstones (Ehrenberg, 1995;Houseknecht, 1987;Li, Yang, et al., 2023;Lundegard, 1992;Pittman & Larese, 1991).Moreover, strongly leached feldspar and rock fragment grains maintained their intact shapes without collapsing (Figure 5a-c).This indicated that the main phase of grain dissolution post-dated major mechanical compaction phase and suggested that any additional mechanical compaction during or after grain dissolution phase was likely to be relatively minor.This hypothesis agrees with the widely accepted concept of progressive mechanical compaction with depth (Lander & Walderhaug, 1999;Paxton et al., 2002;Pittman & Larese, 1991).For example, Paxton et al. (2002) collected 368 uncemented sandstone samples (quartz cement-free and carbonate cements-free) with depths ranging from surface to 6027 m and constructed a compaction curve (represented by intergranular volume) versus depth.The intergranular volume rapidly declined from ~40% to 42% at the surface to ~28% at ~1500 m and continued to decline to only ~26% at ~2500 m.A depth of 1500 m is equivalent to a temperature of ~62°C in the study area (Li et al., 2021), which significantly predates the grain dissolution phase inferred from the precipitation of quartz overgrowth at temperatures between 87.2 and 112.3°C (see Section 4.1.1).Therefore, it is reasonable to conclude that any mechanical compaction during or after the grain dissolution phase is negligible.

| Early and late carbonate cements
Intergranular calcite exhibited a relatively consistent colour in CL images (Figure 5i,j), indicating no distinctive coexisting multiple phases of precipitation.At the shale-sandstone contacts, there were abundant intergranular calcite cements (Figure 5d,e).Some grains were separated by intergranular calcite cements and exhibited a floating distribution.This indicated that the intergranular calcite precipitated during the early stage of compaction, which is supported by the Th of intergranular calcite.The Th of intergranular calcite ranged from 41.8 to 72.5°C (Figure 6a).A temperature of 41.8°C is equivalent to a depth of ~867 m based on the thermal gradient in the study area (Li et al., 2021), where the compaction of sandstones is typically weak.Moreover, intergranular calcite, dolomite and siderite exhibited a consistent precipitation temperature range of ~40-75°C (Figure 6a-c).Ferrocalcite and secondary pore-filling calcite precipitated within grain dissolution pores, indicating that they formed during or after the grain dissolution phase.This is supported by the Th of ferrocalcite (86.7-108.7°C; Figure 6d), which is generally consistent with the grain dissolution phase inferred from the precipitation of quartz overgrowth (Figure 6e; see Section 4.1.1).
Therefore, based on the petrography and Th, the carbonate cements can be roughly divided into two stages.The early stage is mainly composed of intergranular calcite, dolomite and siderite, which precipitated at ~40-75°C and preceded the grain dissolution stage.The late stage is mainly composed of ferrocalcite and secondary pore-filling calcite, which precipitated at ~85-110°C and formed within or after grain dissolution stage.The middle interval B generally exhibited significantly lower early carbonate cements compared to the marginal intervals A and C (Figure 4g). 4.1.4| General paragenetic sequence Based on the aforementioned observations and conclusions regarding diagenetic processes, the paragenetic sequence of the studied sandstones can be reconstructed (Figure 11).Specifically, the precipitations of carbonate cements and quartz overgrowth were based on Th measurements.Grain dissolution phase and the precipitation of kaolinite were suggested to be contemporaneous with the precipitation of quartz overgrowth.Mechanical compaction was assumed to predominantly develop before grain dissolution commenced.Any additional mechanical compaction during or after the grain dissolution phase was considered negligible.Chemical compaction was previously estimated to commence at a temperature of ~80°C (Walderhaug et al., 2000;Worden & Morad, 2000).
Since the additional compaction during or after the grain dissolution phase was considered negligible, the thin section porosity prior to grain dissolution can be estimated using the following equation.The reconstructed thin section porosity prior to grain dissolution is shown in Figure 4i.In general, the middle sandstone interval B exhibited a significantly higher thin section porosity before grain dissolution commenced compared to the marginal intervals A and C. The interface model results indicated that the preferential dissolution of oligoclase occurred in the parts adjacent to the sandstone-shale contact (Figure 7b).This result is consistent with the patterns illustrated in Models A, B and C (Figure 1).
In fact, it is important to note that the interface model employed in this study assumes that the shales have a sufficient supply capacity of LMWOA over 20,000 years (Figure 7).However, due to the thin thicknesses of the top and bottom shales observed in the samples (Figure 3a), the LMWOA yield in these shales is expected to be limited.Barth et al. (1996) determined the yield of LMWOA in shale samples collected from the North Sea continental shelf using hydrous pyrolysis experiments.They utilized a multivariate approach to construct a quantitative equation for predicting the yield of LMWOA.The yield of acetic acid was found to correlate with shale maturity (represented by depth) and total organic carbon (TOC) according to the following equation (Barth et al., 1996).
This equation was used here to provide a rough estimation of LMWOA yield in the studied shales.The stepwise calculation is provided in Appendix S2.The TOC of the Eocene Shahejie Formation shales in the Bozhong Depression typically ranges from 0.7 to 2.5 wt% (Xu et al., 2019).Assuming that the produced LMWOA from both the top and bottom shales can be fully exported to the interlayer sandstones, the dissolved LMWOA concentration in the sandstone pore water would range from 2.7 mg/L (TOC = 0.7 wt%) to 7.7 mg/L (TOC = 2.5 wt%).These estimated LMWOA concentrations are substantially lower than the commonly reported values in sandstone formation waters (Figure S2).It should be noted that the sample location is adjacent to the root area of the fan delta, which results in the low shale thickness (Figure 3a).However, the far end of the fan delta is expected to be in contact with thick shales, where the LMWOA yield should be high (~388-2200 mg/L; see Appendix S2).Therefore, the LMWOA transport by advection from the far end to the root area may explain the commonly reported high LMWOA concentrations in formation waters in areas away from major shale intervals.Moreover, it is expected that the produced LMWOA would initially react with undersaturated minerals within the shales (Barth et al., 1996;Barth & Bjørlykke, 1993), which is supported by the observed dissolution pores within shales (Adeyilola et al., 2022;Baruch et al., 2015).This process may consume LMWOA and result in the near neutral to slightly alkaline pH commonly observed in shale waters (Figure S4a).If the acidity consumption within shales is disregarded and all acidity is exported to adjacent sandstones and reacts with oligoclase within the marginal sandstones (e.g.within 0.1 m thickness near the top and bottom contacts), the produced secondary porosity ranges from 0.33 vol.% (TOC = 0.7 wt%) to 0.95 vol.% (TOC = 2.5 wt%) (Appendix S2).This range is generally consistent with petrographic observations, which show that the secondary thin section porosity in the marginal intervals is 0.23-2.23 vol.% (median: 0.52 vol.%) in the top interval A and 0.43-2.29 vol.% (median: 1.01 vol.%) in the bottom interval C (Figure 4d).In conclusion, the interface model provides a reasonable explanation for secondary porosity in intervals A and C, but not for interval B.

| Validation of the advection model
The 2D advection model results indicated that the middle interval B exhibited higher fluid velocity compared to the marginal intervals A and C (Figure 8a,b).This increased fluid velocity in the middle interval B allowed for more LMWOA-rich water to flow through, resulting in lower pH in the middle interval B compared to the marginal intervals A and C (Figure 8c,d).As a result, the middle interval B had significantly lower oligoclase saturation (Q/K) (Figure 8e,f) and a significantly higher oligoclase dissolution rate (Figure 8g,h).
Moreover, the 2D advection model results demonstrated that the middle interval B exhibited higher kaolinite cements (Figure 10a-g).These modelling results agreed with the petrographic observation, which also indicated higher abundances of kaolinite in the middle interval B (Figure 4e).Moreover, the petrographic observations of quartz cements and late carbonates also agree with the advection model results (see Section 7, Appendix S2).This further supports the validity of the advection model.
The remaining question is why the advection velocity in the middle interval B is higher than in the marginal intervals A and C. The advection velocity is positively correlated with the sediment porosity (Bethke, 2008).Therefore, the higher sediment porosity observed in the middle interval B (Figures 3a and 4a) results in higher advection velocity compared to the marginal intervals A and C (Figure 8a,b).In other words, higher initial porosity in the middle interval B results in the subsequent generation of increased secondary porosity (Figure 9a-g).This conclusion is consistent with the petrographic data, which demonstrates a positive correlation between the reconstructed thin section porosity prior to the dissolution phase and secondary thin section porosity (Figure S5d).
Moreover, the potential factors, including inflowing water pH, LMWOA concentrations and fluid flow velocity, may influence the degree of oligoclase dissolution (Figure 9h-m).However, all these factors did not change the preferential formation of secondary porosity in the middle interval B. 4.2.3 | Mechanisms leading to the higher initial porosity before dissolution phase in the middle interval B Section 4.2.2 emphasizes the importance of the initial porosity prior to dissolution phase in the development of the higher secondary porosity in the middle interval B. Therefore, the question arises as to why the middle interval B had higher initial porosity before grain dissolution commenced (Figure 4i).The additional compactional porosity loss during and after the dissolution phase has been argued to be negligible (see Section 4.1.2).Therefore, based on the reconstructed paragenetic sequence of diagenetic processes in the studied samples (Figure 11), the porosity before dissolution phase can be calculated using the following equation.Among these three parameters, the middle interval B exhibited higher COPL (median: 21.86%) compared to the marginal intervals A (median: 19.43%) and C (median: 21.50%) (Figure 4f).In other words, mechanical compaction had a negative impact on the porosity before dissolution phase in the middle interval B. The middle interval B had higher depositional porosity (median: 33.2%) compared to the marginal intervals A (median: 30.04%) and C (median: 30.82%) (Figure 4h).Moreover, early carbonate cements in the middle interval B (median: 6.25%) were significantly lower than those in the marginal intervals A (median: 10.38%) and C (median: 9.05%) (Figure 4g).Therefore, both the higher depositional porosity and lower early carbonate cements had positive impacts on  2).A smaller grain sorting coefficient (better grain sorting) results in higher depositional porosity.Therefore, the combination of better grain sorting and lower early carbonate cements predominantly resulted in the higher porosity before dissolution phase observed in the middle interval B.

| Diagenetic evolution model
Based on the aforementioned analyses, the sequence of diagenetic phases can be reconstructed to explain the formation of the high secondary porosity in the middle interval B (Figure 12).During the depositional stage, the sandstones in the middle interval B had better grain sorting compared to the marginal intervals A and C, resulting in higher depositional porosity (Figures 4h  and 12a).During the mechanical compaction phase, the middle interval B experienced slightly stronger compaction.However, there was a significantly lower precipitation of early carbonate cements in the middle interval B. These processes further widened the porosity difference between the middle interval B and the marginal intervals A and C (Figures 4i and 12b).Taking the comparison between intervals A and B as an example, the difference in the depositional porosity was 3.16% (Figure 4h).This difference was increased to 4.38% after compaction and the precipitation of early carbonate cements (Figure 4i).During the dissolution phase, the production of LMWOA in the adjacent top and bottom shales was limited, resulting in minor secondary porosity in the marginal sandstone parts.During the advection phase, the higher porosity in the middle interval B facilitated higher advection flux of LMWOA-rich formation water.This resulted in lower pH, decreased oligoclase saturation, a faster oligoclase dissolution rate and ultimately higher secondary porosity in the middle interval B (Figure 12c).Meanwhile, the dissolution byproducts (e.g.kaolinite) were preferentially precipitated in the middle interval B.

| CONCLUSIONS
The fan delta front sandstone bed from the Eocene Shahejie Formation, Bozhong Depression, Bohai Bay Basin, can be subdivided into three intervals.The middle interval B exhibited higher secondary porosity compared to the marginal intervals A and C. Simulations using the interface diffusion model consistently resulted in minor oligoclase dissolution in the central sandstone parts and failed to explain the higher secondary porosity observed in the middle sandstone interval.
In contrast, simulations using the advection model successfully generated higher secondary porosity in the central part of the sandstone.The distribution of diagenetic minerals further supported the validity of the advection model.During the depositional stage, the middle interval B had better grain sorting, resulting in higher depositional porosity.During the mechanical compaction stage, the middle interval B experienced slightly stronger compaction but had significantly lower early carbonate cements.This further widened the porosity difference between the middle interval B and marginal intervals A and C.During the grain dissolution stage, the adjacent top and bottom shales had limited production of LMWOA, resulting in minor secondary porosity in marginal sandstones.During the advection stage, the higher porosity before dissolution phase in the middle interval B facilitated higher advection flux of LMWOA-rich formation water.This led to lower water pH, decreased oligoclase saturation, higher oligoclase dissolution rate and ultimately higher generation of secondary porosity.
The main variables contributing to the dissolution of oligoclase and respective secondary porosity formation are pH, LMWOA concentration and fluid velocity.However, these factors did not change the preferential formation of the higher secondary porosity in the middle sandstone interval.

Highlights 1 .
Higher secondary porosity was observed in the centre of a sandstone unit in the Bohai Bay Basin.2. 2-D advection model successfully generated higher secondary porosity in the sandstone centre.3. Grain sorting, early diagenesis, LMWOA in adjacent shales and advection control the distribution.4. A stepwise evolution model was constructed.
illustrating the shale-sandstone interface diffusion model domain and the intra-sandstone advection model domain.The left domain boundary of the interface model (X position = 0 m) represents the shale-sandstone contact.

F
I G U R E 3 (a) Schematic representation of the sampled sandstone core interval and the distribution of core plug porosity with respect to depth.(b) Ternary diagram illustrating the detrital compositions and lithology of the sandstone samples.Sandstone lithology classification referred to Folk (1974).

F
Core porosity (a), point-counting thin section porosities (b−d), point-counting diagenetic minerals (e and g) and derived diagenetic parameters (f, h and i) of sandstones split into intervals A, B and C. COPL = compactional porosity loss.Early carbonate cements = intergranular calcite + dolomite + siderite.Thin section porosity before dissolution = total thin section porosity -secondary thin section porosity + quartz overgrowth + kaolinite + late carbonate cements.F I G U R E 5 Images of thin section and cathodoluminescence (CL) of the sandstone samples.Pore space is shaded in blue.(a) Plane-polarized light (PPL) view showing the dissolution of feldspars (FD) and rock fragments (RFD).(b) PPL view showed that the feldspar grain was dissolved, while the surrounding calcites (Cal; stained in light red) were not.(c) PPL view showing feldspar dissolution and intergranular calcite.One feldspar grain was significantly dissolved (photo middle left), while the surrounding calcites were not.(d) PPL view showing a sample at the shale-sandstone contact.The sandstone pore space was significantly filled with calcite cements and no obvious grain dissolution.(e) Orthogonal light (XPL) view showing a sample at the shale-sandstone contact, where the sandstone pore space was fully filled with calcite.The fracture may be handmade resulting from sample preparation.(f) PPL view shows that ferrocalcites (FC) filling in both intergranular and secondary pore spaces.(g) XPL view showing intergranular dolomite (Do).(h) XPL view showing intergranular siderite (Sid).(i) CL image showing intergranular and secondary pore-filling calcite (SPC).Calcite exhibited a consistent colour.(j) CL image showing intergranular dolomite and calcite.(k) PPL view showing feldspar dissolution with kaolinite (Kao) filling in feldspar dissolution pore and quartz overgrowth (Qa).(l) PPL view showing kaolinite and quartz overgrowth. ).

3. 4 |
Sensitivity of parameters in the 2D advection model3.4.1 | pH in inflowing waterThe pH of the inflowing formation water played a prominent role in controlling mineral dissolution and precipitation.In general, acidic inflowing water (pH = 5) resulted in significantly greater oligoclase dissolution and kaolinite precipitation compared to near neutral (pH = 6) and mildly alkaline (pH = 7.5) inflowing waters (Figures9h,i

F
I G U R E 6 (a−e) Histograms showing the homogenization temperature (Th) of primary aqueous inclusions captured in diagenetic minerals.

F
Changes in pH (a) and oligoclase dissolution volume (b) over 5-m flow distance and over 20,000 years.Positive values shown in panel (b) indicate mineral dissolution.The left boundary of the model domain (X position = 0 m) was treated as the shale-sandstone contact.The inflowing water pH was 6.The input LMWOA was composed of 900 mg/L H-acetate and 100 mg/L H 2 -oxalate.The advection velocity was 0 m/year and only diffusion was employed.EAGE LI et al.

F
Line diagrams showing the variations of fluid velocity (a−b), water pH (c−d), oligoclase saturation (e−f) and oligoclase dissolution rate (g−h) at the X position = 0.95 m at different simulation steps.Box histograms showing the respective values in three intervals.The inflowing water pH was 6.The input LMWOA was composed of 900 mg/L H-acetate and 100 mg/L H 2 -oxalate.The advection velocity was 0.1 m/year.F I G U R E 9 Map diagrams showing the reacted volumes of oligoclase across the domain at different simulation steps (a−e).Line diagrams showing the reacted volumes of oligoclase at the X position = 0.95 m (f) and the influences of sensitive parameters (h, j and l).Box histograms showing the reacted volumes of oligoclase in three intervals (g, i, k and m).Positive values indicate mineral dissolution.F I G U R E 1 0 Map diagrams showing the reacted volumes of kaolinite across the domain at different simulation steps (a−e).Line diagrams showing the reacted volumes of kaolinite at the X position = 0.95 m (f) and the influences of sensitive parameters (h, j and l).Box histograms showing the reacted volumes of kaolinite in three intervals (g, i, k and m).Negative values indicate mineral precipitation.
Porosity before dissolution phase = Depositional porosity − COPL − Early carbonate cements the porosity before dissolution phase in the middle interval B. The depositional porosity is controlled by grain sorting (Equation

F
I G U R E 1 2 Schematic representation of the reconstructed diagenetic evolution model (a−c) showing the stepwise mechanisms resulting in the higher secondary porosity in the central part of a sandstone unit.COPL, compactional porosity loss; LMWOA, lowmolecular-weight organic acids.