On the Validity of the Relationship 1/T = 1/T B+1/T S+1/T D in NMR Techniques With Regards to Permeability Estimation of Natural Porous Media

One of the most relevant feature of geophysical techniques based on nuclear magnetic resonance is their ability to estimate the permeability of natural porous media, since other geophysical techniques, as the use of the formation factor and neutron well-logs, allow to quantify the volume of water in the media. Permeability is conventionally obtained from decay time of the total resonance signal. However, the fluid in the pores of a medium normally has different mobility degree that can be differentiated by the NMR results. Therefore, a detailed estimation of permeability requires decomposing the total resonance signal as a function of the decay times corresponding to the three mechanisms that contribute to the signal: the intergranular free fluid, the surface layer, and the diffusion relaxation mechanism. The relationship currently used to make this decomposition states that the exponential decay rate attributed to the total resonance signal is the sum of the three existing decay rates. We demonstrate that this relationship is not generally applicable in porous media, showing the contradiction with the much more widely accepted relationships as well as computation examples from three typical decay rates in a single pore and from sandstone with bulk and surface relaxation mechanisms. Consequently, we conclude that the assertion whereby the permeability of any porous medium does not depend on the decay time of the free fluid is an overstatement, since it only applies to very small pore sizes.


INTRODUCTION
Besides quantifying the water content, the second historical feature of nuclear magnetic resonance (NMR) techniques is their ability to estimate the permeability of porous media, which is one of the most relevant petrophysical parameters obtained by different methods, both laboratory and field. In Geophysics, NMR techniques have historically demonstrated their validity to estimate permeability from well-logging measurements and, in the last 2 decades, this method has increased its implementation in surface geophysical prospecting (Legchenko et al., 2002;Roy and Lubczynski, 2003;Hertrich, 2008;Yaramanci and Müller-Petke, 2009;Legchenko et al., 2011;Behroozmand et al., 2015;Larsen et al., 2020). Other interesting applications of NMR are fluid typing, estimating of pore size distribution, and characterization of the contribution of pore types to other petrophysical properties of porous media (Mardon et al., 1996;Golsanami et al., 2021), although they are not analyzed in this study.
The key of the ability to estimate permeability resides in analyzing not only the NMR response of the medium but also its relaxation along micro-intervals-times. So, this estimation requires decomposing the total resonance relaxation into the three existing mechanisms: the inter-granular free fluid "B" (termed "bulk" since Korringa et al., 1962), the surface layer "S," and the diffusion relaxation "D." This decomposing is especially important given the different relationship of decay times corresponding to each mechanism with the permeability of the medium analyzed. So, for sandstones the total resonance signal is close to the free fluid response; however, in the case of high clay content as in shales, the relative volume of the surface layer in the total volume is closer to the porosity occupied by "free" fluid, and the surface layer response progressively achieve influence in the total resonance signal. The measured NMR response provide a curve of amplitude versus time which first outcome is a total relaxation time, it is necessary counting with a relationship between this and the different relaxation times of each of these mechanisms. As of today, the criteria adopted for the decomposition of total resonance rate remain a crucial issue for the analysis of the relaxation curve (Lewis et al., 2015;Luo et al., 2015;Petrov and Stapf, 2017).
The evolution suffered in last decades by the theoretical relations of NMR method, has taken to generalize that the exponential decay rate currently attributed to the total resonance response is the sum of the three decay rates, or in relationship form: 1/T 1 1/T 1B +1/T 1S for the longitudinal relaxation and 1/T 2 1/T 2B +1/T 2S +1/T 2D for transverse relaxation.
For transverse relaxation, spot measurements of T 2 at frequencies typically 10-25 MHz, field-cycling measurements over the frequency range 0.01-20 MHz (T 1 ), and 2D correlation T 2 -T 2 measurements probe different dynamics. Among these dynamics, this work focuses on the application in petrophysical studies for the determination of permeability in natural porous media. The time constants in this case may differ considerably from those of applications in biological media (see Figure 1). However, specific experimental frequency range and/ or measurement type are not considered relevant in order to evaluate the validity of the criticized relationship.
The aim of this work is to recommend not using this relationship, which lack is demonstrated showing the contradiction with the much more widely accepted relationships as well as a computation example from three typical decay rates. The discussed issues about permeability focus on applications of one-dimensional time NMR measurements, not considering the increasing applications of twodimensional NMR developments T 1 -T 2 , T 1 -D-T 2 , etc., mainly target on very small pore sizes (Bernin and Topgaard, 2013;Fleury and Romero-Sarmiento, 2016;Faux and McDonald, 2018;Maneval et al., 2019;Li et al., 2020). Hereafter, the influence of temperature on the relationship in question, and some especial conditions as the presence of paramagnetic particles in the pore fluid is not considered.

BACKGROUND
To decompose the spectrum of the longitudinal relaxation time T 1 , Bloembergen et al. (1948) defined a correlation time τ C characteristic of the random motion, thereby simplifying the calculation of the mean value of the correlation function used by these authors. Korringa et al. (1962) later considered the presence of two migration processes, one within the surface layer and the other from the surface layer to the free fluid in the pore. They described the migration within the "surface layer" by a transition time or correlation time τ S , and the transition from the surface layer to the bulk by a time τ. They established an equivalence between the relaxation time in the surface layer T 1S , the relaxation time T 1B corresponding to the bulk mechanism outside the surface layer, and the observable relaxation time T 1 , given by the following relationship, where V S and V B represent the volume of the surface layer and the bulk region, respectively. These terms are related by V S +V B V POR , where V POR V·Ø is the pore volume, V is the total volume, and Ø is the total porosity. In this analysis, the influence of the proton resonance fields of each region (V S and V B ) upon one another is not the only process considered. Indeed, Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 688686 proton migration would have the effect of molecular displacement beyond the sub-atomic instance. Seevers (1966) used the notation r S (T 1S +τ S ) −1 and the inverse (1/T 1 −1/T 1B ) −1 to estimate permeability. In this calculation, assuming a homogeneous system, Seevers (1966) replaced V S with S Sgr ·h·V (1−Ø), where S Sgr is the specific surface of the grains and h is the thickness of the fluid layer influenced by the surface, and he applied V B ≈ V·Ø V POR when V S « V B . This last approximation has been later used in many works, even in cases in which the spectral inverse Laplace decomposition of the relaxation curve shows that the difference between V B and V S is not large enough to omit the term V S in V POR .
Most subsequent publications have adopted the view that τ S is always negligible with regard to T 1S , so that: where S POR is the pore surface area. For some authors, the water must sample each environment during the timescale of the T 1 or T 2 measurement for Eq. 2 to be valid. Seevers (1966) already considered the fact that the temporal amplitude of the total resonance signal A NMR (t) is the sum of the temporal exponentials corresponding to the two mechanisms, weighted by the relative amplitude of each in the total porosity: where Ø N V N /V POR represents the relative porosity of each region (surface and free fluid).
On the other hand, the use of precise instruments for the determination of the transverse relaxation time T 2 could be considered widespread by the final decade of the 20th century. Analysis of these transverse relaxation times included the presence of diffusion relaxation. This is given by an exponential function which decay time, thereinafter T 2D T 2MG , is equal to 3/(D·G 2 ·γ 2 ·t cp 2 ), where D is the diffusion coefficient, G the magnetic gradient in z direction, γ the gyromagnetic ratio of 1 H, and t cp the Carr-Purcell pulse spacing (Carr and Purcell, 1954). T 2MG was found to have much lower typical values than T 2S (Brownstein and Tarr, 1979). The diffusion relaxation's dependence on a high gradient of the static field (and the presence of minerals with high magnetic susceptibility), as determined in that work, illustrates why this phenomenon was established to be significant only for transverse relaxation times.
Regarding the implications derived from the T 2 CPMG measurement system (Carr and Purcell, 1954;Meiboom and Gill, 1958), in 1990 it was stablished (see reference in Straley et al., 1994) that the NMR signal measured at the peaks of the echoes is the product of the exponentials given by the relationship (Kenyon, 1997): Where A NMR (t) is the temporal amplitude of the total resonance signal. As for T 1 , the T 2 of the total NMR signal is obtained by the logarithmic mean, that is equivalent to fitting the total decay curve to a single exponential, A NMR (t) ≈ A (t 0) ·exp (-t/T 2 ). If Equation 4 is considered to adequately represent the CPMG measurement system, it would be correct (by taking logarithms) to accept the relationship for the relaxation rate of the total resonance signal given by: Obviously, if the characteristics of the analyzed medium and the measurement conditions are not appropriate (strong external field and the presence of minerals with enough magnetic properties) to produce the diffusion relaxation, the rate of this resonance mechanism would not appear in Eqs 4, 5.
Many publications in the recent decades state, whether directly or through citations, that the components of NMR, the ones corresponding to the free fluid porosity, to the surface layer, and to the diffusion relaxation, are produced in parallel, and hence, that Eq. 5 can be applied generally (Hürlimann and Venkataramanan, 2002;Al-Mahrooqi et al., 2006;Keating and Knight, 2007;Daigle and Dugan, 2009;Grunewald and Knight, 2011;Behroozmand et al., 2015).
After the above background overview, the authors of this communication agree that the three NMR mechanisms occur in parallel, but disagree that it implies Eq. 5. Therefore, the objective of this paper is to demonstrate that Eq. 5 is not suitable for use in permeability estimations. To achieve this goal, the methods employed have been to show the contradiction with the much more widely accepted relationships and to show a computation example from three typical decay rates, which proves that the sum of the three exponential signals does not comply with Eq. 5.
The influence of Eq. 5 on the permeability estimation lies in that most used relationships to estimate permeability are the Schlumberger Doll Research (SDR) (Kenyon et al., 1988) and the Timur-Coates (T-C) (Timur, 1969;Coates et al., 1999). The SDR is based on the mean decay time of the total resonance signal, while in the T-C the estimation of the free-fluid volume (FFI) is based on the assumption that the producible fluids reside in larger pores (Coates et al., 1999). However, given that the decay time corresponding to the intergranular free fluid is typically ten times greater than that corresponding to the surface layer, if Eq. 5 is considered correct, the free fluid decay rate is currently neglected (Latour et al., 1995;Weller et al., 2010;Grunewald and Knight, 2011;Osterman et al., 2016, among many others). This removal may lead to the assumption that the permeability of a medium does not depend on the free fluid decay time. In addition, this approach contrasts with the fact that relaxation time of the total resonance signal already provides an appropriate value for the estimation of the permeability, generally dependent on the content of free fluid, although it can correspond to surface relaxation for very small pore sizes.

EVIDENCES OF THE WEAKNESS OF THE CRITICIZED RELATIONSHIP
The main question raised is whether the distribution over time of the resonance total signal is equal to the sum of exponentials (used for T 1 in older studies) or to the product (assumed for T 2 since 1990). The values of T 1 and T 2 , which most authors consider to be of the same order of magnitude for water saturated media T 1 /T 2 ∼ 1.5 (in hydrocarbon saturated media this ratio ranges from 5 to 15), lead to conclude that T 1 and T 2 distribution must correspond to similar decomposition. From the result of algebraic operation, if the resonance total signal were the sum (in-parallel mechanisms), the total decay time would be close to the longest decay time, and if it were the product (in-series mechanisms), the total decay time would be close to the shortest time. Thus, If it were the product, in water-saturated media, where the wetting layer is present, the decay time of the resonance total signal would be close to the surface decay time (tens of ms), which for the authors' knowledge, is contradicted by the results shown in the literature. The total relaxation time extracted from well logs when there is a significant presence of producible water, is similar to that for the free fluid (few hundred of ms) (see Figure 5A). Therefore, the logic of the decomposition of the relaxation curve A NMR (t) leads to that must be realized by the sum of exponentials, which can be expressed in the form, where I MG , I S and I B address the relative amplitude of each relaxation mechanism in the analyzed fluid volume.
In this communication, we state that Eq. 5 is not adequate for two reasons. The first reason is that, as NMR in porous media is a phenomenon in which the three mechanisms are produced in parallel, we stated this implies that NMR relaxation must be considered as a sum of exponential functions. That is, Eq. 6 is considered to be more appropriate, since essentially maintains the same conceptual decomposition as Eq. 3, and is the current basis for the decomposition in porous media of the total relaxation curve in the main NMR effects [the determination of the distribution of the A NMR (t) relaxation curve is not the subject of this communication]. Thus, if Eq. 6 is correct, which in this work is considered indisputable, it is mathematically impossible to consider that Eq. 4 is also correct, and therefore, to consider Eq. 5 valid.
Both Equations 4, 6 are presented in most works as relationships for a single pore, not specifically reflecting that they apply to the decomposition of the signal corresponding to all the pores of a porous medium.
Regarding Equation 4, for multiple pores of radius sizes R, Kleinberg and Horsfield (1990) considered a Gaussian distribution around a central value R 0 , given by exp (-R/R 0 ) 2 , so that the total NMR signal for multiple pores is given by: In this equation, the exponential function corresponding to the magnetic gradient relaxation was substituted by the expression that relates its decay rate with the molecular diffusion coefficient. For the statement of Eq. 7, Kleinberg and Horsfield (1990) defined ρ 2 h/T 2S (named "relaxivity" since Howard et al., 1993) and considered ρ 2 /h ρ 2 ·(S POR /V POR ) which is equal to 3·ρ 2 /R for spherical grains. In Eq. 7, the exponential functions of bulk and diffusion relaxations do not depend on pore radius. From a theoretical standpoint, the influence of pore size on relaxation rates can be seen in Faux and McDonald (2018); besides remarking the surface relaxation at distance of a few nm and the influence of measurement frequency on T 1 values, the bulk relaxation dependence on pore type and size is also showed.
One aspect that has strongly influenced the meaning of the decomposition of the T 2 distribution, has been its correlation with the pore sizes in the medium (Loren and Robinson, 1970;Vogeley and Moses, 1992;Howard et al., 1993;Straley et al., 1994;Kenyon et al., 1995). Following Coates et al. (1999), the multi-exponential decomposition in water-saturated media reveals the distribution of pore sizes and each T 2 value corresponds to each pore size, some authors consider that the T 2 distribution only reflects the size of each pore. In contrast, we consider that the T 2 distribution directly shows the weight of each exponential; in some cases, this distribution can be correlated with the pore size, but always reflects the existing relaxation mechanisms. Thus, as in the conventional conception for a saturated granular medium with a wetting fluid, e.g., water, the distribution values for times around hundreds of seconds are "proportional" to the fraction of pores with grain sizes larger than 4 μm (silts, sands and gravels), and the distribution values for low times are "proportional" to the amount of small pores. This is the reason why there are local maxima, or at least an inflection . Disregarding magnetic gradient relaxation, these local maxima do not occur because two main pore sizes are common (bimodal grain size distribution), but because there are two main relaxation mechanisms.
In the external reviews of this study, the extension to all the pores of a porous medium starts from A(t) Σ j (A j ·exp (−t/T 2j )), where j denotes each pore. Then, in the absence of magnetic gradient relaxation, if contrary to that supported in this communication, the relationship 1/T 2j 1/T 2Sj +1/T 2Bj is considered valid, it would follow that A j ·exp (−t/T 2j ) A j ·exp (−t/T 1Sj )·exp (−t/T 1Bj ). The latter relationship would indicate that in each pore, the two mechanisms -surface and bulk-are not additive processes, but the relaxation mechanisms occur in series. The general extension of this relationship for multiple pores would be: However, this signal does not match with the Laplace decomposition, whose conventional expression is given for the different decay times by: Moreover, the T 2 distribution addressed in Eq. 8 would reflect only the decomposition into different pore sizes, and from the product of exponentials, T 2B and T 2S could not be extracted (nor T 2MG in the case of this mechanism occurs).
Therefore, when analyzing a porous medium in which there are different pore sizes, the generalization of the decomposition of the total NMR signal into the sum of exponentials is given by: where I B , I S , and I MG are the relative amplitudes stated for Eq. 6, which are typically obtained by inverse Laplace transform, and j denotes the different spaces involved in each pore for each relaxation mechanism. Eq. 10 maintain the same decomposing showed in Equation 2 from Behroozmand et al. (2015) (Eq. 11), although in that work the magnetic field relaxing in the transverse direction with relaxation times T 2j are not specifically related to the different relaxation mechanism.
Using properly the T 2 CPMG system, although the transverse total resonance curve is the result of successive spots, this curve is not the series of juxtaposed values, but the concatenation of the values that reflect the state of relaxation of the resonance at times separated by a specific time interval. This fact is that justifies that the transverse total resonance curve shown in literature cases on sandstones (see Figure 5A) presents a T 2 (time for which the value of the curve drops to 1/e of its initial value) close to that corresponding to bulk relaxation.
Equation 10 does not means that the observed transverse relaxation rate will be given by the weighted sum of bulk, surface and diffusion relaxation rates.
When the three relaxation mechanisms are considered in a single pore, the Laplace spectrum results in the sum of three type Dirac delta functions. Otherwise, if a porous medium with different pore sizes is analyzed (truly more realistic), the part of Laplace spectrum corresponding to each relaxation mechanism will be wider, and the result is the sum of three distributions centered on each relaxation mechanism. An example of the distribution of the three relaxation mechanisms provided by Eq. 10 is shown in Figure 2 (adapted from Allen et al., 2000).
Hence, except for the relationship criticized, the renowned foundations on the decomposition of the total signal, especially the results obtained in practical cases through the inverse Laplace transform, show that this decomposition consists of the sum of the different resonance processes, whether for a single pore size or for multiple sizes.
Unlike Equation 4, the physical significance assumed in this work to consider that the different relaxation mechanisms occur in parallel, is that they are produced by the same external fields. Each resonance mechanism occurs in its corresponding time interval although it may be overlapped, and the relative amplitude of each, by its own character, is weighted in the normalized curve of the total resonance response. However, each NMR decay rate does not significantly depend on the others as it happens with intensity in an electrical circuit when the current meets several resistors in parallel. As a representative example of this disengagement, figure 10 from Straley et al. (1994) of a sandstone sample has been selected (see Figure 3). T 2 distribution of the surface relaxation for the centrifuged sample at 100 psi vs. air (the free fluid is then negligible) results very similar to the part of T 2 distribution corresponding to the surface when the sandstone is fully saturated.
FIGURE 3 | T 2 distribution of a sandstone sample before and after extracting free fluids (adapted from Straley et al., 1994).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 688686 The second reason for disagreement with Eq. 4 is the corresponding Eq. 5, because in most real cases of porous media, none of the "measurable" relaxation times T 1 or T 2 meet Equation 1/T (1,2) (1/T (1,2)N ), where N indicates the three relaxation mechanisms (B, S and MG).
This fact is easily verified by means of any practical example on a single pore, such as that shown in Figure 4. The values used for this example are I MG 0.1 with T 2MG 5 ms for the diffusion relaxation, I S 0.4 with T 2S 30 ms for the surface layer, and I B 0.5 with T 2B 270 ms for the free fluid, which are typical values for sandstones (see the similar values used in Figure 1 from Timur, 1969, and see also Straley et al., 1994;Allen et al., 2000;Mohnke, 2014;Osterman et al., 2016). As it can be seen in Figure 4, the "observable" value of 1/T 2 (the obtained rate when total NMR curve is fitting to an exponential function) is 1/255 0.0039 ms, while the result of 1/T B +1/T S +1/T MG 0.20 + 0.033 + 0.0037 0.237 ms is 100 times higher, demonstrating that it is generally not possible to adopt Eq. 5.
The same would occur if only two relaxation mechanisms were considered or if several relaxation times are considered for each mechanism. Thus, the conclusion is the same for the three exponentials of the single pore model, as for the n B +n S +n MG exponentials of a multiple pore sample. In short, it is mathematically obvious that if a given curve is the sum of more than two exponentials functions with clearly differentiated decay times, the rate showed by the curve is not equal to the sum of rates of each exponential function.
In spite of the above, it should be noted that as the distribution of the V S and V B fractions exhibit a progressively greater difference, the decay of the total relaxation curve would become more similar to the decay time corresponding to the free porosity. It is necessary to exclude the cases in which measurements are carried out within a total time span close to the time of surface relaxation, and where the relative volume of the surface layer in the total volume is close to the porosity of the "free" fluid (which is the case for high clay content). It should also be noticed that although Eq. 5 appears to indicate that longer decay times (T B » T S » T MG ) have less influence on the slope of the average straight line on a logarithmic scale, in practice the mean relaxation time T ML is more similar to the longer decay time.

RESULTS ON PERMEABILITY ESTIMATION
According to the last paragraphs of Background section, it is concluded that it is too bold to state that the permeability of a porous medium is independent of the free fluid decay time, because 1/T ≠ 1/T S . That statement is only true if the area corresponding to the free fluid under the Laplace spectrum of the total resonance signal is negligible compared to the area corresponding to the surface layer (very small pores).
In order to evaluate the error that the criticized equation can produce in the permeability estimation, the T 2 distribution of the work of Kenyon et al. (1995), p. 3, has been taken as an example (see Figure 5). In Figure 5A) the signal amplitude of the total resonance signal together with fitting curve to the digitizing data are shown. Considering the CMR porosity Ø for the total resonance signal (equal to the area under the T 2 distribution) and the corresponging porosity Ø B and Ø S for the bulk and surface mechanisms (see Figure 5B), the values Ø 0.29, Ø B 0.225, and Ø S 0.064 are obtained. Fitting the total resonance signal to a monoexponential function T 2 0.170 s, and fitting the signals corresponging to the bulk and surface mechanisms, the decay rates obtained are T 2B 0.180 s and T 2S 0.020 s. Using the powers established by Kenyon et al. (1988) to estimate the permeability from logarithmic-mean T 2 , k constant·Ø 4 ·T 2 2 , and taking constant 50, the resulting value is k 1.01·10 -2 darcy. However, if only the T 2S component is considered because 1/T 2B «1/T 2S , the resulting permeability value would be k 3.1·10 -7 darcy. The fact that the ratio between these values (independently of the factor 50 adopted) reaches four orders of magnitude, demonstrates the importance of the present study for permeability estimation with NMR techniques.
The results presented in Figure 5A) are a verification that the decay time of the total signal is similar to the bulk relaxation time, contrary to would occur if criticized equation was correct.

DISCUSSION AND CONCLUSION
The concepts and relationships established in Brownstein and Tarr (1979) for the decomposition of the NMR response have been used in some publications; however, in this study it is considered that they do not apply to porous geological media. Brownstein and Tarr (1979) established that the amplitude of the different modes of surface relaxation is differentiated in fast, intermediate and slow regions, considering the occurrence of molecular diffusion, the absence of both, the volume-like sinks and the gradient of the external field. These modes are given as a function of the value M·a/D, where M is the average surface magnetization on the active surface, a the characteristic length of each geometry (flat, spherical and cylindrical), and D the molecular diffusion coefficient. Thus, in the so-called fast diffusion region, for which M·a/D«1, the amplitude of the lowest mode corresponding to the highest decay time T 0 is much higher than the others, so that the total signal of NMR is completely dominated by the lowest mode. In the slowdiffusion region, the majority of the relative intensity still corresponds to the lowest mode and the higher modes contribute a few tens of percent following Brownstein and Tarr (1979). With this model, for which details of volume and surface-like "sinks" are unimportant (Brownstein and Tarr, 1979), the obtained relaxation times are in the order of several ms. However, it must be noted that the magnetic moment per unit volume is given in that work by multiplying a series of constants by the orthogonal spatial eigenfunctions. The final result for the total nuclear magnetization given in Equation 10 of Brownstein & Tarr, produces the same Laplace spectrum as that obtained by Eq. 10 in this work, so both solutions are equivalent.
The results of Results On Permeability Estimation section, it is verified that to extract more information from the T 2 distribution, and then ruling out the bulk relaxation because 1/T 2 ≈ 1/T 2S , provides permeability values very different than the value directly obtained from measured T 2 . In practice of the most works, especially in industry, the permeability values obtained from this T 2 providing satisfactory results. Although in this study, the example presented is considered very representative of typical values in sandstones, more validation of the proposed claims could be done by laboratory works.
In conclusion, the steps in the theoretical development of the expression of the NMR relaxation amplitude with time are unimportant if the exponential decay rate attributed to the total NMR signal is ultimately considered as the sum of the existing decay rates, the values of which are clearly different. Moreover, the characteristics of the measurement system do not matter if the dependence of the echo time, the wait time and the record time in the results are disregarded. If these parameters are adequately selected, the relationship 1/T 1/T B +1/T S +1/T MG is not complied.
The criticized relationship has become irrefutable due to its widespread use, but no mathematical demonstration has been found to support that the total transverse NMR is the product of the three relaxation mechanisms. However, two very different evidences demonstrate that considering the rate of the total signal NMR equal to the sum of rates from each mechanism is not adequate.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.   Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 688686