Analytical studies on the tidal wave propagation in a coastal leaky aquifer commonly assume that the low-permeability aquitard is homogeneous. The aquitard is, however, vertically heterogeneous in nature due to varying soil types along the burial depth as can be frequently seen from borehole logs. In this study, an innovative analytical model is developed to explore the vertical heterogeneity in aquitard hydraulic conductivity (K) and specific storage (Ss) on the tidal wave propagation coupling with vertical leakage from the aquitard. The novelty behind the newly derived analytical solution is attributed to the fact that the aquitard along the burial depth can be divided into any number of homogeneous zones, each of which is associated with distinct K and Ss values, so that arbitrary vertical heterogeneity pattern of aquitard K and Ss can be captured. Theoretical analysis results reveal that an interlayer within the aquitard with a smaller K can significantly enhance the amplitude and phase shift of periodic groundwater head fluctuations in the leaky aquifer, while a larger Ss of the interlayer will weaken the amplitude and phase shift provided a relatively small aquitard K. Buried locations of the interlayer also implement nonnegligible effects on the tidal wave propagation. For the scenario of exponentially decaying aquitard K and Ss, which is commonly encountered for a thick aquitard, a larger decay exponent results in smaller amplitude and phase shift. This analytical study highlights the importance of vertical aquitard heterogeneity on tidal wave propagation in a coastal leaky aquifer system.

Sea tide can cause periodic groundwater level fluctuations in coastal leaky aquifers. This phenomenon can be termed tidal wave propagation and has attracted wide attention in the research area of coastal hydrogeology (e.g., [1, 24]). Earlier analytical models for the tidal wave propagation mainly considered the aquifer system to be homogeneous (e.g., [46]). In fact, the heterogeneity is a common place in lithosphere including the coastal leaky aquifer system.

Alluvial plains are often associated with spatial variability of sediments in terms of soil type and grain size. Particularly, alluvial fans formed at the base of mountains usually have coarser sediment at the upstream part of the fan and finer sediment at the downstream part [7, 8]. To deal with such a horizontal heterogeneity in a coastal leaky aquifer system, Trefry [9] developed an analytical solution for a finite composite coastal aquifer, providing prescribed head, prescribed flux, and mixed inland boundary conditions. Trefry and Bekele [10] then applied this analytical solution to handle tidal wave propagation signals monitored at multiple wells placed in a finite-length island aquifer system of Garden Island, Australia. Monachesi and Guarracino [11] and Wang et al. [12] analytically investigated the effects of horizontally linear heterogeneity in aquifer hydraulic conductivity on the tidal wave propagation.

Groundwater leakage from the aquitard has been found to have critical effects on the tidal wave propagation [13]. In nature, the aquitard usually consists of relatively fine sediments associated with lower permeability while larger storativity comparing to sandy aquifers [14]. Consequently, the volume of groundwater released from aquitard storage is a significant portion of groundwater leakage that recharges the aquifer [15]. From this perspective, vertical heterogeneity in aquitard hydraulic properties (i.e., the hydraulic conductivity K and specific storage Ss) can implement nonnegligible effects on tidal wave propagation in the adjacent leaky aquifer. However, currently available analytical models for coastal leaky aquifer system have not dealt with the vertical heterogeneity in aquitard K and Ss.

Vertical heterogeneity in aquitard layers is rather common in the lithosphere and can be directly identified according to lithology analysis data. Generally speaking, the aquitard can be divided into multiple zones along the burial depth profile, and each subunit has a distinct soil type associated with independent K and Ss. Such a type of vertical aquitard heterogeneity can also be identified via interpreting in situ pressure and compaction data [1619].

In addition to vertical heterogeneity of multiple zones, the aquitard can also contain a sand layer, gravel interlayer, and cracks. In other words, there could exist local anomaly of the soil type in the aquitard, leading to local heterogeneity in aquitard K and Ss. Even the small-scale heterogeneity could significantly control groundwater flow in the subsurface. For instance, Meriano and Eyles [20] found that the existence of aquitard interlayers played an important role in supplying the aquifer below. Gerber et al. [21] found that the existence of crack and sand gravels in the aquitard produced abnormal vertical groundwater flow velocity in the aquitard. Harrington et al. [22] studied the effect of sandy interlayer in the aquitard on solute transport and developed a model to explain the radial component of solute diffusion when the cross section of sandy interlayer is circular or elliptical, and the model was successfully applied to explain the source of high chloride ion concentration in the clay-rich area in southern Saskatchewan, Canada.

At the regional scale, the porosity of geologic media usually shows a decreasing trend with the burial depth, which can also be indirectly manifested as the depth-decreasing hydraulic properties [2326]. Many studies have reported that the aquifer K exponentially decays with burial depth [15, 2730]. Comparatively speaking, a thick aquitard could also demonstrate exponentially decaying hydraulic properties with the burial depth. Ferris et al. [31] measured the aquitard K values at 15 sites in south-central Saskatchewan, Canada, and further compiled the K data through 21 studies to characterize the hydrogeology of aquitard on a regional scale. Their results indicated that the aquitard K followed an exponentially decaying pattern with depth index, and the variation in aquitard K reached 5 orders of magnitude in the range of 80 m. Previous studies also found that the aquifer Ss also decreases with depth [3234]. Recently, Kuang et al. [35] proposed an empirical model of aquifer Ss exponentially decaying with depth and also concluded that Ss has an exponential decay with depth in the shallow crust. Without losing generality, we believe that the aquitard Ss could also exponentially decay with depth at a regional scale and from the statistical perspective, even though few literatures provided relevant evidence.

The objective of this study is to explore the influences of vertically heterogeneous aquitard K and Ss on the tidal wave propagation in a coastal leaky aquifer system. Firstly, we construct a quasi-2D model for a horizontally extended coastal leaky aquifer system with a finite length and with the consideration of depth-dependent aquitard K and Ss. Secondly, the analytical solution associated with tidal-head inland boundary condition will be derived. Lastly, tidal wave propagation bias of the groundwater level fluctuations caused by the vertical heterogeneity of aquitard K and Ss will be discussed.

2.1. Mathematical Model

Here, we consider a horizontally extended coastal leaky aquifer system as shown in Figure 1. To consider the unconfined aquifer overlain, the aquitard of interest has a large specific yield to effectively damp the tidal effect so that the water table in the unconfined aquifer almost remains constant (e.g., [3, 13]). Besides, the difference in density between seawater and fresh groundwater is neglected. It has been verified that ignoring the density difference when one derives the analytical solution results in an error not greater than 2.5% of the tidal amplitude and the error decreases landward [36].

Let the x-axis be positive in the landward direction and coincide with the contact interface between the aquitard and leaky confined aquifer. Let the z-axis be vertical and positive upward. The origin of the Cartesian coordinate system is located at the coastline (Figure 1). The tidal boundary is located at x=0 and the inland-side boundary is located at x=l. Vertical Cartesian coordinate of the zone interfaces inside the aquitard are indicated by z1,z2,,zn,,zN. The aquitard K is at least two orders of magnitude smaller than the underlying confined aquifer, so that the groundwater flow is essentially vertical within the aquitard (e.g., [2, 37]). In comparison, groundwater flow is horizontal in the confined aquifer.

Under the above assumptions, the governing equations for groundwater flow within the confined aquifer and the aquitard are
(1)Shx,tt=T2hx,tx2+K1h1x,z,tzz=0,(2)Ss,nhnx,z,tt=Kn2hnx,z,tz2,zn1zzn,
where h=hx,t is the aquifer hydraulic head [L] at location xL and time tT, TL2T1 and S are the transmissivity and storage coefficient of the leaky confined aquifer, respectively; hn=hnx,z,t is the aquitard hydraulic head at location x,z and time t corresponding to zone n; KnLT1 and Ss,nL1 are the vertical hydraulic conductivity and specific storage of the aquitard corresponding to zone n, respectively.
Based on the Euler formula, the tidal boundary for the confined aquifer at x=0 as shown in Figure 1 is expressed as
(3)hx,tx=0=A0expiwt+c,
where i=1, A0, w, and c are the amplitude [L], angular frequency [T1], and initial radian phase [-] of the tidal fluctuations, respectively.
We also consider a tidal-head boundary at the inland side (i.e., x=l). Here, we consider different amplitudes and initial phase shift at the inland-side boundary comparing with that at the coastline. Therefore, one has
(4)hx,tx=l=Alexpiθexpiwt+c,
where i=1, hmsll, Al, and θ are the mean sea level [L], amplitude [L], and initial phase difference [-] at the inland-side boundary, respectively. Note that Equation (4) is typically representative for an island responding to dual tides [10, 38, 39], where Al=A0 seems to coincide with the real-world situation. If the hydraulic head does not tidally fluctuate (i.e., Al=0), Equation (4) degenerates to a constant-head inland boundary condition, which could be suitable for a horizontally extensive coastal aquifer system.
In terms of the aquitard, at the interface between zones n and n+1 (i.e., z=xn), the groundwater flow within the aquitard has to meet the continuity requirement of head and flux. As a consequence, one has
(5)hnx,z,tz=zn=hn+1x,zn,tz=zn,Knhnx,z,tzz=zn=Kn+1hn+1x,z,tzz=zn.
In addition, at the interface between the aquitard and the underlying confined aquifer of interest, as well as the constant-head overlying unconfined aquifer, hydraulic head is also continuous; thus, one has
(6)h1x,z,tz=0=hx,t,(7)hNx,z,tz=xN=0.

2.2. Analytical Solution

Based on the separation-of-variable method, which has been widely employed to handle analytical models of tide-induced aquifer head fluctuations (e.g., [4, 39]), the periodic hydraulic head of the confined aquifer, as well as that of the aquitard, can be, respectively, expressed in the following complex form as
(8)hx,t=hpxexpiwt+c,(9)hnx,z,t=hnpx,zexpiwt+c,
where hpx contains the information of amplitude and phase shift of periodic aquifer hydraulic head fluctuations at the horizontal location x, and similarly, hnpx,z contains the information of amplitude and phase shift of periodic aquitard hydraulic head fluctuations at the horizontal location x and vertical location z within the zone n. It is easy to find that the derivation of the analytical solution for hx,t is equivalent to the derivation of hpx.
Substituting Equation (8) into Equation (1) yields
(10)2hpxx2+K1Th1px,zzz=0iwSThpx=0.
Substituting Equation (9) into Equation (2) yields
(11)2hnpx,zz2ξn2hnpx,z=0,
where ξn=iwSs,n/Kn. The general solution of Equation (11) is
(12)hnpx,z=Enxcoshξnz+Fnxsinhξnz.
The matching conditions at the interface between zone n and zone n+1 are
(13)hnpx,zz=zn=hn+1px,zz=zn,(14)Knhnpx,zzz=zn=Kn+1hn+1px,zzx=xn.
Substituting Equations (13) and (14) into Equation (12) yields
(15)EnxFnx=ΨnEn+1xFn+1x,
where
(16)Ψn=coshξnznsinhξnznKnξnsinhξnznKnξncoshξnzn1·coshξn+1znsinhξn+1znKn+1ξn+1sinhξn+1znKn+1ξn+1coshξn+1zn.
Thus, one has
(17)E1xF1x=Θn1EnxFnx,
where
(18)Θn=Πk=1nΨk=Ψ1Ψ2Ψn.
It is noteworthy that if the aquitard is vertically homogeneous, Ψn and Θn will become the identity matrix and one has E1x=Enx and F1x=Fnx. In order to handle the zonational vertical heterogeneity of the aquitard, i.e., to determine Enx and Fnx in Equation (17), the boundary conditions for zones 1 and N (i.e., Equations (6) and (7)) must be incorporated. In terms of zone 1, substituting Equations (8) and (9) into Equation (6) yields
(19)h1px,zz=0=E1x=hpx.
In terms of zone N, substituting Equations (8) and (9) into Equation (7) yields
(20)hNpx,zN=ENxcoshξnzN+FNxsinhξNzN=0.
Substituting Equation (17) into Equations (19) and (20), one can further derive
(21a)ENx=tanhξNzNΘN11,2tanhξNzNΘN11,1hpx,(21b)FNx=1ΘN11,2tanhξNzNΘN11,1hpx,(21c)F1x=ΘN12,2tanhξNzNΘN12,1ΘN11,2tanhξNzNΘN11,1hpx.
After substituting Equations (19) and (21c) into Equation (12), one can further obtain
(22)h1px,zzz=0=ξ1ΘN12,2tanhξNzNΘN12,1ΘN11,2ΘN11,1tanhξNzNhpx.
Substituting Equation (22) into (11) yields
(23)2hpxx2η2hpx=0,
where
(24)η2=iwSTK1ξ1TΘN12,2tanhξNzNΘN12,1ΘN11,2tanhξNzNΘN11,1.
The general solution of Equation (23) is
(25)hpx=Ccoshηx+Dsinhηx,
where C and D are constants that need to be determined by incorporating the boundary conditions at the coastline and inland side. Substituting Equations (3) and (4) into Equation (25) yields
(26)C=A0,D=AlexpiθcschηlA0cothηl.
As a result, the analytical solution for hpx is
(27)hpx=A0coshηx+AlexpiθcschηlA0cothηlsinhηx.
The normalized amplitude α (i.e., the ratio of the amplitude of the periodic head fluctuation to that of the tide) and radian phase shift (τ) of aquifer head fluctuations are computed as
(28)αx=hpxA0,τx=arctanImhpxRehpx,
where hpx is the modulus of hpx that represents the amplitude at the horizontal location x, and Im and Re denote the imaginary and real parts of the expression, respectively.

Previous studies have figured out the effects of horizontal aquifer heterogeneity on the tide-induced aquifer head fluctuations. This study mainly focuses on analyzing the effects of vertical heterogeneity in the aquitard K and Ss. Here, we consider a coastal leaky aquifer system of a long but narrow island with a horizontal extent of 500 m. To the best of our knowledge, the analytical model with the consideration of vertically heterogeneous aquitard K and Ss has not been yet developed for an island leaky aquifer system. By referring to the ranges of hydraulic parameters in leaky confined systems summarized by Li and Jiao [13], the transmissivity and storage coefficient of the confined aquifer are set to be T=1.5×102 m2/d and S=1×104, respectively. The aquitard thickness is set to be 10 m. For the tidal waves, the amplitude is A0=AN=1 m, the initial phase shift is c=0, the phase difference is θ=0, and the angular frequency is w=6.283 d-1 corresponding to the period of 1 day. Two scenarios of vertical aquitard heterogeneity are considered in the following analysis. The first scenario is related to the aquitard interlayer, which is associated with K and Ss anomaly, indicating the local heterogeneity. Another scenario is related to the exponentially decaying aquitard K and Ss, indicating the generally varying heterogeneity at the regional scale.

3.1. Interlayer with K and Ss Anomaly

To explore the effects of an aquitard interlayer on the tidal wave propagation, we first assume that the interlayer is 0.5 m thick and located at the middle of the aquitard. Consequently, the aquitard of interest can be divided into three zones (i.e., N=3) with z1=4.75 m, z2=5.25 m, and z3=10 m. Zones 1 and 3 are prescribed with the same hydraulic parameters that are K1=K3=1.0×102 m/d and Ss,1=Ss,3=1.0×104 m-1. The interlayer (i.e., zone 2) is prescribed with dramatically different K (case 1) or Ss (case 2) values to illustrate local anomaly of aquitard hydraulic parameters.

For case 1, we consider aquitard K anomaly, while the aquitard keeps vertical homogeneity in Ss, i.e., Ss,2=1.0×104 m-1. Two K values are, respectively, prescribed to the interlayer, which are K2=1.0×100 m/d (a high-K interlayer) and K2=1.0×104 m/d (a low-K interlayer). In addition, the situation associated with the homogeneity in aquitard K, i.e., K2=1.0×102 m/d, is also considered for the aim of comparison. A high-K interlayer is generally corresponding to a thin subunit consisting of higher-content sands or lower-content clay in the aquitard. Comparatively, a low-K interlayer is generally corresponding to a thin subunit that consists of higher-content clay. Both of the two situations are frequently seen in practical hydrogeological investigations. The situation of inland constant head boundary was also considered to reveal the generality of the newly developed analytical solution and to demonstrate the difference in the tidal wave propagation between the common coastal leaky aquifer system and the island leaky aquifer system.

Figure 2 shows that a high-K interlayer seems to produce negligible effects on tidal wave propagation from both the perspectives of amplitude and phase shift. However, a low-K interlayer can produce a remarkable bias of tidal wave propagation. To be specific, as shown in Figure 2, the low-K interlayer will significantly increase α, as well as τ (refer to corresponding absolute value hereinafter). This is mainly due to the fact that the low-K interlayer blocks the groundwater leakage into the confined aquifer, and the tidal wave propagation tends to be increasingly dominated by the confined aquifer.

As for the coastal aquifer system, the pattern of head fluctuations in the confined aquifer is similar to that of the bilaterally fluctuating island aquifer system, but the amplitude variation is small. Since the inland is a constant head boundary, the fluctuation amplitude gradually decreases inland until it is 0.

For case 2, we consider the Ss anomaly, while the aquitard keeps vertical homogeneity in K (i.e., K2=1.0×102 m/d). Two Ss values are, respectively, prescribed to the interlayer, which are Ss,2=1.0×103 m-1 (a high-Ss interlayer) and Ss,2=1.0×105 m-1 (a low-Ss interlayer). In addition, the situation associated with the homogeneity in Ss (i.e., Ss,2=1.0×104 m-1) is also considered for the aim of comparison. Generally speaking, the specific storage values should be also highly correlated to the clay content for an aquitard. In this study, a high-Ss interlayer may correspond to the deposits with more soft clay particles showing higher compressibility, and vice versa.

Figure 3 shows that a low-Ss interlayer seems to produce negligible effects on tidal wave propagation, especially the phase shift. However, a high-Ss interlayer will cause a decrease in α and τ. The main reason for this phenomenon is that the higher the Ss of the interlayer, the more water it can release or store. As a consequence, relatively more leakage into the confined aquifer is expected. In other words, the tidal wave propagation tends to be decreasingly dominated by the confined aquifer.

Now, we analyze the effects of the locations of the low-K interlayer and the high-Ss interlayer on the tidal wave propagation. Here, we consider the K values of zones 1 and 3 to be K1=K3=1.0×101 m/d, while zone 2 is associated with K2=1.0×103 m/d (a low-K interlayer). For comparison analysis, we consider the interlayer locating in the upper aquitard associated with z1=7.25 m and z2=7.75 m and the interlayer locating in the lower aquitard associated with z1=2.25 m and z2=2.75 m.

Figure 4 demonstrates the effects of the low-K interlayer location on the tidal wave propagation. Specifically, when the low-K interlayer locates closer to the confined aquifer, the amplitude αwill be significantly increased, while the radian phase shift τ will be also increased. This is mainly due to the fact that the closer of the low-K interlayer to the confined aquifer, the weaker the hydraulic connection between the confined aquifer and the overlying unconfined aquifer and the smaller the amount of aquitard depletion, so that the vertical leakage from the aquitard is far more decreased, making the tidal effect tend be dominated by confined aquifer rather than be influenced by the leakage. On the contrary, when the low-K interlayer locates in the upper aquitard, the amplitude α will decrease and the radian phase shift τ will increase. This is mainly because the low-K interlayer that locates further to the confined aquifer can make the aquitard release a certain amount of groundwater from storage and leake to the confined aquifer to mitigate the tidal wave propagation.

Figure 5 demonstrates the effects of high-Ss interlayer location on the tidal wave propagation. Specifically, when the high-Ss interlayer locates closer to confined aquifer, the decrease in the amplitude α is almost negligible and the radian phase shift τ will be visibly increased. The cause for this phenomenon is that the high-Ss interlayer locates closer to the confined aquifer, and the vertical leakage will be enhanced by more amount of groundwater released from aquitard storage, inhibiting the groundwater level fluctuations in the confined aquifer. However, for a relatively larger aquitard K (i.e., 1.0×101 m/d considered here), suggesting stronger hydraulic connection between the confined aquifer and the overlying unconfined aquifer, the amount of groundwater released from aquitard storage is insignificant compared to the vertical leakage from the unconfined aquifer, so the amplitude α is almost unchanged with the location of high-Ss interlayer. However, we note that when the aquitard K is relatively small (i.e., 1.0×102 m/d considered before), the presence of a high-Ss interlayer decreases the radian phase shift τ (i.e., Figure 3(b)), while when the aquitard K is relatively large (i.e., 1.0×101 m/d considered here), a high-Ss interlayer increases the radian phase shift τ (i.e., Figure 5(b)).

Overall, either a low-K or a high-Ss interlayer in the aquitard can result in significant bias of tidal wave propagation in the confined aquifer, implying that the two aquitard heterogeneity scenarios cannot be ignored and this is particularly true when the low-K or high-Ss interlayer locates closer to the confined aquifer.

3.2. Exponentially Decaying K and Ss

To consider the scenario of vertical aquitard heterogeneity in K and Ss that can be empirically expressed as an exponential function of depth (e.g., [26, 40, 41]). In this study, the function fz=f0expλz is applied to express the exponential variation pattern, where f0 is the corresponding value at the contact interface between the aquitard and confined aquifer and λ is the decay exponent expressing the exponential relationship. Note that a positive λ produces a depth-decaying exponential relationship of K and Ss versus depth. Particularly, λ=0 refers to the homogeneous situation.

It should be noted that the exact analytical solution with the incorporation of exponential relationship mentioned above is beyond the scope of this study. To serve our analysis, we seek to uniformly divide the entire soil depth profile of the aquitard into enough zones, each of which is assumed to be homogeneous in K and Ss. After that, specific values are assigned to corresponding zones to approximate the exponential relationship. For the 10 m thick aquitard of interest, one can see from Figure 6 that a uniform discretization of 20 zones is accurate enough to capture the exponential relationship of K and Ss versus depth. To consider the situations of slightly, moderately, and highly heterogeneous aquitard K and Ss, we set λ=0.2 m-1, 0.3 m-1, and 0.4 m-1, respectively. Meanwhile, the homogeneous situation, i.e., λ=0.0 m-1, is also included for the aim of comparison.

Figure 7 shows that the tidal wave propagation in the confined aquifer is greatly affected by the exponentially decaying aquitard K. In particular, a larger λ results in a smaller α, as well as τ. Such a phenomenon is easy to understand. A larger λ means the groundwater fluctuations in the confined aquifer can penetrate to a larger distance from the aquitard bottom. Consequently, the aquitard can supply more leakage for the confined aquifer to mitigate periodic head variations.

Similarly, in terms of the exponentially decaying Ss, a larger λ means the aquitard can release or store more water provided the same groundwater level fluctuations inside. Thus, more leakage from the aquitard can be supplied for the confined aquifer and smaller α and τ will be produced (Figure 8).

Overall, we can see from Figures 28 that the vertical heterogeneity in aquitard K and Ss has significant effects on the tidal wave propagation in the coastal leaky aquifer system, implying the importance of practical hydrogeological investigations of vertical aquitard heterogeneity. In addition, the newly proposed analytical model can be applied to determine hydraulic parameters of a leaky aquifer system via interpreting in situ hydraulic response data of surface water-groundwater interaction.

In this study, we derive an analytical solution for the tidal wave propagation in a coastal leaky aquifer system with the consideration of vertical heterogeneity in aquitard K and Ss. Two scenarios of vertical aquitard heterogeneity in K and Ss (interlayer with K and Ss anomaly and exponentially decaying K and Ss) are considered, and their effects on the tidal wave propagation are analyzed. Main conclusions are summarized as follows:

  • (1)

    A low-Kaquitard interlayer makes the amplitude and the radian phase shift of the confined aquifer increase obviously, while the high-Kaquitard interlayer has negligible effects on the tidal wave propagation. A low-Ssaquitard interlayer results in negligible effects on the tidal wave propagation, especially the phase shift. However, a high-Ssaquitard interlayer will cause a decrease in the amplitude and phase shift

  • (2)

    If the low-Kaquitard interlayer locates closer to the confined aquifer, it could significantly increase the amplitude while decrease the phase shift. On the contrary, if the low-Kaquitard interlayer locates further to confined aquifer, the amplitude decreases while the phase shift increases. When the high-Ssaquitard interlayer locates closer to the confined aquifer, the amplitude decreases and the phase shift increases. On the contrary, if the high-Ssaquitard interlayer locates further to the confined aquifer, the amplitude increases while the phase shift decreases

  • (3)

    Exponentially decaying aquitard K and Ss is expected and has significant effects on the tidal wave propagation. A large decay exponent, corresponding to the highly heterogeneous aquitard K and Ss, yields a remarkable decrease in the amplitude and the phase shift. Practical hydrogeological investigations of vertical aquitard heterogeneity are highly recommended

The source code data of the analytical solutions used to support the findings of this study are available from the corresponding author upon request.

The authors declared that they have no conflicts of interest to this work. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

This research was partially supported by the Key Research and Development Program of National Natural Science Foundation of China (2019YFC1804301), the National Natural Science Foundation of China (41902244), the Basic Research Program of Jiangsu Province (Natural Science Foundation) (BK20190497), and the China Postdoctoral Science Foundation (Grant No. 2021M690866).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).