Advances in low-field nuclear magnetic resonance (NMR) technologies applied for characterization of pore space inside rocks: a critical review

NMR serves as an important technique for probing rock pore space, such as pore structure characterization, fluid identification, and petrophysical property testing, due to the reusability of cores, convenience in sample processing, and time efficiency in laboratory tests. In practice, NMR signal collection is normally achieved through polarized nuclei relaxation which releases crucial relaxation messages for result interpretation. The impetus of this work is to help engineers and researchers with petroleum background obtain new insights into NMR principals and extend existing methodologies for characterization of unconventional formations. This article first gives a brief description of the development history of relaxation theories and models for porous media. Then, the widely used NMR techniques for characterizing petrophysical properties and pore structures are presented. Meanwhile, limitations and deficiencies of them are summarized. Finally, future work on improving these insufficiencies and approaches of enhancement applicability for NMR technologies are discussed.


Introduction
Rock is a sort of porous media with tight texture and multiple components, which is also the research focus in the field of oil and gas exploration and development. Since cores from deep reservoirs are precious and strongly heterogeneous, the reusability of those cores and acquirement of experiment data are of great significance for experimental design. NMR is such a method that satisfies the demand of characterization of pore space inside rocks as the mechanism of NMR and optical measurement are completely different, and the NMR technology is less limited by experimental conditions (Pape et al. 2009;Srivastava et al. 2018). For petroleum engineers and researchers, NMR studies of rocks include three levels: (1) knowledge of the NMR principle, (2) relaxation theories and models of porous media for signal interpretation, and (3) taking characteristics of lithology and fluids into consideration and analysis of special relaxation mechanisms in rocks. Generally, interactions of atoms and molecules are expressed as bulk relaxation, surface relaxation, and diffusion-enhanced relaxation in an inhomogeneous magnetic field at a macroscopic scale. In pores, hydrogen nuclei ( 1 H) work as "probes" to touch the pore walls (Brown and Fantazzini 1993;Weisskoff et al. 1993;Borgia et al. 1995), and present opposite states before and after touches. Therefore, fluids fill up pore voids and outline the pore structure. Additionally, physical and chemical reactions between fluid and solid molecules, or among fluid molecules make the relaxation response more complicated. Song (2013) has divided the development of NMR technologies in porous media into two stages: the classical period  and modern period (2000-now), with a shift from the research bloom of relaxation theories and models to promotion of NMR methodology interpretation. This is because the target cores show considerable complexity in NMR responses, and challenges appearing in applications are MR signal interpretation and collection. Specifically, signal overlapping due to similar resonance responses from different sources and diffusion coupling in well-connected voids generates ambiguous results, and short relaxation signal reflecting hydrogen-bearing components in nanoscale pores tends to be easily lost (Hiejima et al. 2005;Liu et al. 2018a). To reproduce the relaxation process in experiments, find out dominant mechanisms that affect short relaxation signal components, and induce signal overlapping, inhomogeneous magnetic field gradient distribution that influences nuclei diffusion and NMR responses in rock pores have been simulated in several studies (Audoly et al. 2003;Zhang and Hirasaki 2003). However, oversimplified models are of limited accuracy, and idealized results can significantly deviate from real data. To strengthen the ability of NMR signal separation and interpretation, improvements of numerical simulation and experimental methodologies are needed (Hoop and Prange 2007;Mu et al. 2007;Fleury et al. 2015). Scaled-up numerical simulation for relaxation in pore network provides supplementary data for experimental observations, but its effectiveness relies on the full understanding of specific relaxation mechanisms due to divergent rock compositions and pore fluids, inserting appropriate controlling equations, and refining digital core reconstruction based on CT (computer tomography) images (Zhang and Zhang 2015;Wu et al. 2019). As for laboratory tests, subtraction of identical signal achieved by comparison of multiple pulse sequences can extract hidden information. Meanwhile, for applicability enhancement, NMR tests should be checked and corrected by other experimental methods. For example, imaging measurements can be used to support the characterization of pore geometry and topology, and indirect experimental approaches such as core imbibition, flooding, and excavation can be utilized in a dynamic observation process.
This paper aims at introducing the major challenges and insufficiencies in characterization and providing guidance for future work that can make up for these deficiencies. Section 2 will give a brief introduction for the development history of relaxation theories and models for porous media. Section 3 will review mature applications using NMR techniques. Section 4 will summarize the existing challenges and analyze effects of inaccurate models or methods on test results. Section 5 will discuss efforts on studies of relaxation mechanisms, NMR response simulation, and experimental methodology to promote the applicability of NMR technologies and progress of core analysis.

Enlightenment and origin
The accumulation of 1 H relaxation produces the total MR signal detected, and fluid throughout rock pores is the carrier of 1 H, of which the distribution is used to outline geometrical and topological properties of pore networks. Since spins belong to quantum category, the probability density m(r, t) is introduced to describe the state evolution of hydrogen nuclei in pore space. Zimmerman and Brittin (1957) found that transition between free and bounded water triggered by molecular diffusion brought out significantly distinct relaxation behavior. They thought the fact was caused by exchange of spin states and in line with the Markov process. Initially, the Chapman-Kolmogorov equation (Olivares-Robles and García-Colín 1994) firstly gave an expression of 1 H probability density: Supposing that the spins share n kinds of energy state, the subscripts "s," "i," and "f" denote the start (complete polarization), intermediate, and final (complete relaxation) states, respectively. For example, "sf" denotes that a spin experiences the process of "start → final" state. c is the change probability from the present state, while m represents the probability of transition from one state to another. Specially, m if is the stationary distribution, given by Eq. (2): To simplify Eq. (2), only the start and final states (denoted by "f") are considered, and the probability density m sf is obtained by solving Eq. (3): where 1/T s is the relaxation rate of the s-state spin, and E sf is the transition rate between the two states. Meaningfully, this formula compares the relative speed of inherent relaxation rate and molecular diffusion: The whole relaxation is regarded as the weighted summation of spin relaxation when the former is faster than the latter; on the contrary, the density of polarized 1 H in each part is almost the same.
Without sufficient understanding of the relaxation mechanism in porous media, Eq. (3) was solved by the first-order perturbation theory at early time, contributing to the loss of some useful information in calculation (Brownstein and Tarr 1977).

Classical model
Early research (Torrey 1956;Torrey et al. 1959;Korringa et al. 1962) found that 1 H relaxed faster than ever once they reached the pore surface, indicating that pore walls would exert extra relaxation interaction and limit molecular diffusion. The T s extra effect is called "surface relaxation." Brownstein and Tarr (1979)  where M(t) is the magnetization vector, ⃖� ⃗ D is the spin diffusion coefficient (a second-order tensor), 1/T B is the bulk relaxation rate related to spin inherent properties, n denotes the unit outward normal vector at the surface, ρ is the surface sink strength, and V p is the pore volume. In finite space, the expression of the probability density of spins becomes extremely complex , and eigenfunction expansion is adopted to obtain m(r, t): where A n is the proportion of the nth eigenstate, φ n and 1/T n are the corresponding eigenfunction (related to spin location) and eigenvalue (related to relaxation time), respectively. The minimum eigenvalue 1/T 0 is the total longitudinal or transverse relaxation rate. And for transverse relaxation, 1∕T 0 2 (i.e., 1/T 2 ) is the summation of surface relaxation rate 1/T 2S , bulk relaxation rate 1/T 2B , and diffusion-enhanced relaxation rate 1/T 2D (Kleinberg et al. 1994;Stapf et al. 1995;Kimmich and Anoardo 2004) Actually, pores in the BT model are isolated (several microns or even smaller) and spins reach pore walls easily, satisfying T B ≫ T n . Since the inhomogeneity of the magnetic field is neglected and the diffusion is considered to be instantaneous, 1/T 2D is also ignored. (Senturia and Robinson 1970;D'Orazio et al. 1989): where S is the pore surface area, a is the pore characteristic length, 1/T m is the relaxation rate of the rock skeleton, and h is the thickness of the surface relaxation region. The ratio of the time a 2 ∕D , which spins need to travel to the pore boundary, and the spin relaxation time a∕ , defines three kinds of diffusion domain, i.e., fast diffusion ( a∕D < 1 ), slow diffusion ( a∕D > 10 ), and modest diffusion ( 1 < a∕D < 10 ) Mitra et al. 1993). For slow diffusion, if bulk relaxation is neglected, strict restrictions are required: Since the inequality group (Eq. 8) is difficult to be satisfied simultaneously, the model is less practical. The effect of bulk relaxation was highlighted by simulating spin relaxation in pores with sizes ranging from tens to hundreds of microns (Cohen and Mendelson 1982). Results showed the magnetization evolved evenly inside pores and verified the rationality of the fast diffusion model to some extent.
In the fast diffusion domain, the unrelaxed spin probability density is approximately the same in porous media, so the overall relaxation time, a∕ , contains no knowledge of pore shape and connectivity, and the pore characteristic size a is only an averaged value. Instead, in the slow diffusion domain, because the spin relaxation rate is greater than the replenishment rate caused by molecular diffusion, there is a discrepancy between the probability density in pore void and near boundary, namely m(r, t) ≠ const . 1 H is sensitive to obstacle and all eigenvalues release information of geometries. By comparing the two sorts of diffusion regimes, it is found that spin relaxation reflects the detailed microstructure efficiently as long as the diffusion velocity is sufficiently small. However, this will lead to a complicated analytical expression.

Diffusion model
To include the behavior of spin diffusion in a non-uniform magnetic field, the diffusion propagation function G(r, r � |t) (i.e., the diffusion propagator G, a Green's function) is introduced to replace m(r, t) , and Eq. (4) is then rewritten as Eq. (9): where δ 3 is the Dirac function in 3-D space. Analytical methods, like the perturbation theory, eigenfunction expansion, and cumulant expansion, were used to solve the simultaneous equations .
In an infinite and uniform environment, where the diffusion coefficient is constant ( ⃖� ⃗ D = D 0 ), spins diffuse freely, and G(r, r � |t) is the Gaussian distribution, satisfying 2 = 2D 0 t (Tanner 1968(Tanner , 1970: According to the Einstein diffusion equation, the averaged traveling distance of spins within time t, a statistical amount, In rock porous media, the diffusion is confined and the diffusion coefficient is no longer a constant. If the anisotropy of the diffusion coefficient is not considered (only time-dependent), ⃖� ⃗ D is simplified as a time-dependent diffusion coefficient, D(t) (Woessner 1963;Stejskal and Tanner 1965;Banavar and Schwartz 1987;Basser et al. 1994). The diffusion propagator G(r, r � |t) then deviates from the Gaussian distribution, and its expression is rewritten in Eq. (11), where C(t) is normalized and S(r − r � ) indicates the connectivity of pores.
Enlightened by the thought of random walk, Eq. (12)  shows how D(t) works as a geometrical probe to extract features of pore structure, where d is the spatial dimensionality and k is the wave vector: In fact, the initial and long-term states of spin diffusion are the most easily available knowledge, while the intermediate diffusion process is hardly depicted because of the complicated pore structure and mazy diffusion path. Researchers focused on short-time and long-time behavior of the time-dependent diffusion coefficient (Axelrod and Sen 2001). Within the short-time limit, spins are dominated by the free diffusion regime, and only a small part of them is able to contact the boundary, so the pore structure is outlined by the mean square displacement, i.e., ⟨r(t) 2 ⟩ = 6D(t)t . Regardless of rock connectivity and smoothness of the pore surface, the dominant term, (D 0 t) 1∕2 , in the expression of D(t) , will be remained in the regime. And the variation of D(t) is related to V p ∕S . In the long-time limit, D(t) depends on pore types. For isolated pores, the effective diffusion coefficient D ∞ finally comes to 0, while for connected pore network under the assumption of =0 , D(t) is contributed from those survived spins which reach their destinations after bypassing several voids and passages. Therefore, the path which spins have gone through reflects the connectivity and distortion of the pore structure.
Some structural information is excavated from the shorttime and long-time asymptotic behavior; however, it is still insufficient to draw the exact profiles of pores. The mediumtime diffusion process between the two limits includes messages of pore geometry and morphology, but the analytical solutions of models in this domain are complex and impractical. The anisotropy of diffusion coefficient in a biological cell model with particular geometry has been analyzed in some studies (Kinoshita et al. 2008), but it oversimplifies rock pore.

Effect of internal magnetic field
In an unknown magnetic field, spin diffusion is more complicated. Due to the susceptibility difference △χ between fluid and skeleton particles, or para-/ferromagnetic materials, an internal magnetic field gradient g int with an unclear distribution is produced in B 0 . The expression of g int is simply given by Eq. (13): where C is a dimensionless constant that considers local fluctuation of the internal field gradient. To be more specific, spin diffusion regimes are determined by four lengths that appear in spin motion, which are the diffusion length l D , the dephasing length l g , the structural length l s , and the critical length l * (Hürlimann 1998;Johnson and Daigle 2016): The asymptotic behavior of the diffusion regime depends on the minimum length among the four. When l D ≪ l s and l D ≪ l g are satisfied, the diffusion process is equivalent to free diffusion, in which D(t)=D 0 is met. When (D(t)t) 1∕2 is close to the pore size a, as time increases and D(t) decreases, spins reach pore surface and relax gradually. Both the two occasions are described by asymptotic behavior of the free diffusion regime. When l g ≪ l D and l g ≪ l s are satisfied, or l s > l * is valid, spin diffusion is greatly modulated by the local magnetic field gradient g(r) . Therefore, the spin echo decay is described by the localization regime. In the third asymptotic regime, the motional averaging regime, satisfying l s ≪ l g , l D , spins will arrive at the border before its complete dephasing and relaxation (Lewis and Seland 2016).
In fact, since the internal magnetic field distribution can only be obtained by numerical methods (Sen et al. 1999;Fu et al. 2012;Lewis and Seland 2017), it is challengeable to observe the spin relaxation process in rocks. However, these studies are still conducive to understanding diffusion behavior in different kinds of porous media.

Section summary
In general, the development of relaxation model experiences a blooming period from 1950s to 2000s, and those models become more complex and continually get closer to real relaxation processes. However, they are mainly applied for theoretical analysis since they are solved by analytical methods with simplified assumptions and the solution formats are complicated. Even if the internal magnetic field distribution can be obtained by numerical simulation, there is a big difference between reality and digital results.

Low-field NMR methodology
Practically, as mentioned in our previous study (Wang et al. 2020), the NMR signal in cores is seriously affected by the local magnetic field gradient, and the signal-to-noise ratio (SNR) is low. As a result, low-field NMR is more broadly used. The low-field NMR technology has two branches, namely the NMR relaxation spectrums and magnetic resonance imaging (MRI). For relaxation spectrums, analyses of the line shape and trend change of T 1 , T 2 (longitudinal and transverse relaxation time) spectrums and calculation of spectral parameters are the most powerful tools, while the correlation between the diffusion coefficient and time is often used for fluid identification, diffusion-relaxation studies and analyses of pore structure. Since the collected signal is a statistical result, signal overlapping due to similar NMR responses from different sources tends to occur in T 1 , T 2 spectrums. The relaxation spectrum analysis has been developed from one-dimensional spectrums to two-dimensional spectrums (Hürlimann and Venkataramanan 2002;Sun et al. 2004). To some extent, two-dimensional spectrums are capable of distinguishing different signal sources that produce similar NMR responses by adding relaxation physical information into one map. To be more specific, characterization of petrophysical properties, such as porosity, fluid saturation, and permeability, is successfully achieved by the NMR relaxation spectrum analysis based on NMR signal and variation of fluid distribution. Pore structure characterization, including pore size distribution, geometrical, and topological properties, is also related to fluid distribution and 1 H diffusion. Moreover, the two-dimensional relaxation spectrum analysis is a must for fluid identification, and it is also used for obtaining pore size distribution. Development of those applications is summarized in Table 1.
T 1 tests mainly contain two methods: One is the inversion recovery (IR) method and the other is the saturation recovery (SR) method. For IR, the variation range of magnetization is larger, and the transverse relaxation effect is eliminated. SR is less disturbed by noise and is less time-consuming. For T 2 , the Carr-Purcell-Meiboom-Gill (CPMG) sequence eliminates the dephasing-enhanced relaxation in the non-uniform magnetic field by n spin echoes. Meanwhile, the diffusion-enhanced relaxation is basically eliminated by adopting the minimum echo interval τ. And the time difference between T 1 and T 2 tests demonstrates the superiority of the latter. However, to utilize diffusion coefficient information and avoid the inhomogeneous and unknown internal magnetic fields generating uncertainty on calculation of D(t), the pulse field gradient (PFG) and stimulated echo (STE) sequences are applied (Latour et al. 1993;Sorland et al. 2000). Actually, two-dimensional pulse sequences are mainly D-T and T-T and are a combination of the single relaxation time or diffusion coefficient. For D-T,  Obtaining the full distribution of pore size based on Eq. (7) Tarr (1977, 1979) Obtaining microfracture width based on the modifying shape factor T 1 -T 2 spectrum Separation of NMR responses from water and organic matters Fleury and Romero-Sarmiento (2016) D-T 2 is more widely used and is further divided into two categories. The first type records data of D and T 2 in two separate stages, which allows getting a larger range of D, but may lose short relaxation components easily, such as PFG-CPMG and bipolar PFG-CPMG (BPPFG-CPMG) (Komlosh et al. 2007;Zheng and Price 2008). The second kind of sequences, instead, loads constant field gradient, preserves responses of short relaxation components in the first stage by varying spin echo interval, for example, modified CPMG. As for T-T, correspondingly, usage of T 1 -T 2 is more frequent and is drawn by IR-CPMG or SR-CPMG (Snaar and van As 1992;Rondeau-Mouro et al. 2016). Compared with conventional methods in calculation of rock porosity, fluid saturation, and permeability, NMR techniques are less time-consuming, for core flooding is unnecessary, and in most cases, the latter behaves no worse than those conventional methods. In terms of characterizing core wettability, NMR relaxation spectrums show more convenience since qualitative and quantitative analyses of wettability are related to T 2 variation and NMR saturation, respectively (Chen et al. 2006). Acquirement of pore size distribution and observation of pore structure variation before and after special treatments, such as acidizing, both require water-saturated cores. It is worth noting that water saturation operation will alter the pore structure of clay-mineral-rich cores, such as shale. And filling fluid selection is essential. Moreover, sample preprocessing like drying for gas absorption also changes pore structure. To eliminate the human error, NMR cryoporometry is created of which the calculation is based on the signal difference between 1 H in fluid and icing state. Meanwhile, it also avoids interference of para/ferromagnetic impurities on capture of short relaxation components. Still, the lower bound of NMR cryoporometry is confined by signal detected. In two-dimensional relaxation spectrums, difference of 1 H relaxation behavior in various fluids and part of solids is magnified gradually because molecules containing hydrogen atom(s) in fluids and solids show distinct T 1 and T 2 value (Zhang et al. 2016a, b).
In addition, MRI is used to directly observe the fluid distribution inside rock samples, through which pore structure changes are contoured Balcom 2013, 2014;Zhang et al. 2019). It also helps researchers understand seepage laws and fluid-solid physicochemical interactions. Due to its low resolution in images, MRI is mainly used for observation of loose or artificial samples, and on supplementing relaxation spectrum analyses.

Limitation and deficiency
Although low-field NMR technologies are relatively versatile in characterizing pore space inside rocks, the technology applicability in this field will be challenged by its limitations and deficiencies as hot topics of research shift from conventional formations to unconventional formations, with target cores becoming more complex in compositions and structure. Hence, in characterization of petrophysical properties and pore structures, result analysis is challenged by low SNR, signal loss, and a signal overlapping.

Petrophysical property characterization
For an NMR porosity test, all the relaxation responses are regarded as the contribution from the fluids in pores, but the relaxation rate can be accelerated by other components or impurities inside cores (Almagor and Belfort 1978). The types and properties of pores and fluids remained on pore surfaces also have an impact on NMR signal collection, which produces discrepancy between results and real data. Consequently, the choice of proper fluid media and calibration before tests are of vital importance. As for characterization of saturation, permeability, and wettability, the determination of the T 2cutoff is on the top agenda. However, it is tough to gain a precise value as errors produced in signal collection and inversion of T 2 spectrum are inevitable, and determination of T 2cutoff strongly depends on experimenters' experience (Hu et al. 2019).
The low-field T 1 -T 2 spectrum is used for identification of fluids, and differentiating signal from fluids and organic matters in shale samples to some extent ), but it still remains some overlapped sections, and parts of 1 H-contained components' signals are lost (as shown in Fig. 1).
For cores rich in organic matters and nanopores like shale samples, because the transverse relaxation rates of fluids in small pores and solids contained 1 H are both extremely fast, and the corresponding relaxation signal are approaching the sampling limit, there is little divergence on T 2 values, and the relaxation responses among fluids and solids cannot be sufficiently separated . Thus, shale porosity may be overestimated or underestimated, and organic matter analysis by NMR becomes very hard. In fact, since the relaxation rates of fluids and solids are significantly different, and solid NMR signal is more clearly captured in high field, it is nearly impossible to acquire the whole range of effective signal by an identical device. In other words, low-field NMR technologies are very likely to set a "blind zone" for rock composition analysis (Mehana and Elmonier 2016;Kausik and Hürlimann 2016). Recently, several reports have discussed the method of high-field T 1 -T 2 maps for separating signals of solid components (Papaioannou and Kausik 2015; Song and Kausik 2019), but they are isolated from cores. This means that mineral composition analysis on a core from underground by NMR techniques is yet unaccomplished, for signal and noise will be amplified simultaneously.

Pore structure characterization
Pore size distribution models based on the BT theory describe 1 H relaxation in the fast diffusion domain. In this situation, each pore is isolated, or at most two adjacent pores are connected with pore sizes below 1 μm, and all the pores, throats, and fissures in cores are simplified to the same geometry, such as sphere, ellipse, or cylinder (Brownstein and Tarr 1979;Hürlimann et al. 1994). In most case, spherical pore hypothesis is used . Definitely, this simplification fails to represent the actual pore structure. Figure 2 shows the SEM (scanning electron microscope) images of a shale sample with developed clay slits. It is provided that the storage space merely consists of sphere pores, some of which the radius is r and the characteristic length a is 0.33r. Hence, a 10r × r × r slit will be equivalently treated as a sphere with a equaling 0.24r.
In addition, for cores with well-connected pore networks or with pore sizes ranging from tens to hundreds of microns (Lisitza 2002), such as cavernous/fractured carbonates, water molecules shuttle in large and small pores through throats or fissures. Consequently, the relaxation time is averaged by the periods when 1 H is relaxed in these closed voids alone, also known as the diffusion coupling (as shown in Fig. 3). Because of diffusion coupling, the T 2 distribution fails to be converted to pore size distribution and is difficult for reliable explanation. Comparison between fast diffusion and diffusion coupling shows that whether the T 2 distribution can be converted to pore size distribution is determined by surface relaxation, and better correspondence demands more intensive interplay (Palit and Yethiraj 2008). Although NMR cryoporometry offers new access to pore size distribution, a limited characterization range and inconvenient operation restrict the application of this technique.
In aspect of depiction for pore shapes and connectivity, it is almost impossible to obtain the 1 H diffusion path for tortuous pore networks by analytic methods. As mentioned above, mathematical descriptions are only available for extreme moments, such as short-time and long-time diffusion limits, and we have no idea about what will happen in most time between the two conditions (as shown in Fig. 4). Random walk (RW) simulation has been used to study the medium-time diffusion process. However, the basic models are scale-limited and oversimplified, and relaxation mechanisms are not fully considered (Guo et al. 2016a). Moreover, the 1 H probe is insensitive to pore walls when surface The low-field NMR (2 MHz) T 1 -T 2 map for all the components in unconventional shales in a homogenous magnetic field (from  relaxation is not the decisive factor, resulting in unrecognition of pore shapes (Lagkaditi et al. 2015). As for MRI, fluid signal in tiny pores, which is out of sampling range or covered by noise, is easily lost. This means details about these tiny pores and throats are neglected, which is one of the reasons for low resolution. The onedimensional T 2 mapping technology enhances identifiability of cores' main characteristics on images, but local features, such as fissures, will become blurred due to signal overlapping (as shown in Fig. 5). Take a step back, based on the relaxation theory described in Sects. 2.2 and 2.3, and regardless of diffusion coupling, the surface relaxation rate is controlled by the shortest distance that 1 H travels. For fractures of which the width is close to, or is smaller than the characteristic length of pores, relaxation time of 1 H in fractures and pores is overlapped (as shown in Fig. 2). Thus, pores and fissures cannot be completely identified by spectrums obtained from current NMR methodologies. Porous-type carbonate Fractured carbonate Fig. 3 Schematic diagram of diffusion coupling in porous-type and fractured carbonate. 1 H (red curve arrow) wandering in solo pore will not be affected by diffusion coupling, while 1 H (green curve arrow) crossing large and small pores through throats or fissures will be influenced

Improvement direction
To strengthen the applicability of NMR technologies used in core analyses, progresses are indispensable in NMR signal collection and interpretation. Objectively, signal collection is mainly challenged by physical limits of devices; for example, the sensitivity of RF (radio frequency) coil is crucial for pulse emission and signal acquisition, which directly determines whether the relaxation can be fully observed (Song 2013). Additionally, development of novel pulse sequences for sampling more relaxation information is also an important booster. However, since interpretation of inversion data is more affected by manual interference, more attention should be paid to this direction.

Improvement of theories and models
The MR relaxation theory in porous media has experienced a development for more than 60 years, but the development speed of theories and mathematical models are slow or even stagnant in recent years, causing defects in applications. Specifically, for 1 H relaxation behavior beyond the fast The 1 H walking path in medium time for its large time span and complicated process is not available, and it is difficult to couple the medium-time diffusion process into the short-time and long-time diffusion limits

Slices
One-dimensional T2 mapping Fig. 5 Enhancement of main fractures (yellow oval) and loss of a slim fissure (red rectangle) in MRI images based on the one-dimensional T 2 mapping technology diffusion domain, there is no reliable model or method for qualitative explanation and quantitative analyses. Therefore, it is necessary to improve the theories and refine the models with improved knowledge of the basic physics and chemistry characteristics of rocks.
The first step is to enhance the understanding of multiple relaxation mechanisms related to rock components and pore structures. For instance, prominences including composition diversity, storage, and seepage dominated by nanoscale pores and throats, and coexisting multi-scale pores, make relaxation mechanisms in shales more complex (Washburn and Birdwell 2013;Daigle et al. 2014;Jia et al. 2016Jia et al. , 2017Washburn and Cheng 2017): (1) Organic matter affects fluid NMR signal, and dipole-dipole coupling between 1 H atoms is the controlling factor. (2) The distribution and compaction of clay influence NMR responses, and clay's type and components will affect the wettability of pores and adsorption of water molecules. (3) Most of Fe 3+ ions are distributed in pyrite, and a few paramagnetic impurities are dispersed on the pore surface or in organic matter. Therefore, interaction between fluid molecules and para-/ferromagnetic ions may not be the dominant cause of surface relaxation (as shown in Fig. 6). Furthermore, 1 H diffusion in shale belongs to the free diffusion regime with high probability, but diffusionenhanced relaxation and diffusion coupling should not be overlooked if spins diffuse in fissures or multiple pores.
Secondly, in order to find out dominant mechanisms that induce ambiguous interpretation in characterization of petrophysical properties and pore structures, applications of high-field NMR and solid NMR are necessary (Le Doan et al. 2013;Kausik et al. 2017), and scaled-up numerical simulation is needed to reproduce the relaxation process in experiments and supplement data for experimental observation. NMR response simulation consists of two stages: (1) digital core reconstruction and (2) response simulation (Benavides et al. 2020). As shown in Table 2, methods for NMR response simulation in the two stages are targeted at solving part of the problems mentioned in Sect. 4.
Actually, digital core reconstruction is the foundation of 1 H relaxation simulation in pore networks, but it is never an easy task and existing methods fail to build scaled-up or realistic pore network. Process-based methods are believed to restore the core internal structure to the greatest extent because it simulates rock formation processes like sedimentation and diagenesis. However, some of them are idealized and rock compositions are not fully considered, leading to simple or even single skeleton and void components (Zou et al. 2015). In contrast, image-based methods are more acceptable with continuous improvement in computation and imaging measurement. Still, refined core reconstruction depends on image resolution and data holding and processing capabilities of computers. Since there is a balance between resolution and field of view, imaging devices, such as nano-CT that pursues clear images, can only sacrifice the observation range, and make sample preparation more laborious work. On the contrary, artificially reducing the resolution to get a full picture will cause a unit pixel or voxel size exceeding part of pore sizes, leading to pore information loss. In particular, the micro-X ray CT equipped with detector involving a large field of view in Southwest Petroleum University has a maximum field of view of 51 mm and a corresponding maximum pixel size of 54.2 μm.

Experimental methodologies
Besides breakthrough in theories and models, researchers, who are weak in research and development (R & D) of apparatus and pulse sequence, but interested in signal interpretation should be more concentrated on experimental design innovation by using existing tools and methods. To cope with the issues mentioned in Sect. 4, several corresponding ideas and approaches are summarized in Table 3, followed by some examples.
In order to overcome the limitation of T 2cutoff in characterizing tight sandstone permeability, Fan et al. (2018) introduced the dual T 2cutoff method to subdivide the storage space into three sections: fully movable part, partially movable part, and fully bounded part. When the pore throat is less than the first cutoff value (T 2cutoff1 ), fluid is completely immobile, such as bounded water in clay; when the pore throat is greater than the second cutoff value (T 2cutoff2 ), fluid moves freely, such as free water in large pores; when the pore throat is within T 2cutoff1 and T 2cutoff2 , there are both movable and immobile fluids in pores, allowing to reduce error range. Actually, among four widely used methods for determining T 2cutoff as summarized in our previous study (Wang et al. 2020), the most effective and accurate one is parameter calibration. To be specific, irreducible water saturation and fluid volume of a sample are measured in advance; when the integrated area of T 2 spectrum (fully saturated state) is equal to the irreducible fluid volume, the corresponding T 2 value is the exact T 2cutoff (Wang et al. 2019).
For signal overlapping caused by various 1 H-contained components, signal subtraction means highlighting distinct but overlapped NMR responses by subtracting the identical but dominant relaxation signal in different tests (Müller-Huber et al. 2018). For example, binomial edited CPMG (BE-CPMG) shields 1 H signal in solids, and the T 2 spectrum difference from BE-CPMG and CPMG indicates location of fluids in pores and induced fissures (Washburn and Birdwell 2013). The signal difference from IR-spin echo and IR-solid echo tests reflects the distribution of organic matter on pore surfaces, while the signal difference from IR-solid echo and IR-magic echo tests reflects solid organic matter in shale.
In regard to quantitative analysis of pore size, a DDIF tool utilizing higher eigenmodes is developed by the signal difference from reference signal R and test signal E (Song 2003;Cho and Song 2008). Tests combining DDIF with other pulse sequences like CPMG are also used in studies of pore structure characterization and non-uniform field distribution (Zhang et al. 2018;Connolly et al. 2019). For qualitative analysis, a more direct away is combining processes Process-based methods Reservation of arbitrary and irregular void shape makes it possible to separate NMR signal from pores and fissures Yan et al. (2013) Image-based methods Reservation of mineral compositions makes it possible to capture short relaxation signal of different solid components Ogren (2013), Tiwari et al. (2013) Geometric algorithms A kind of simplified method for study NMR response difference between pores and fissures with regular or simplified shape and structure Jia et al. (2007) Statistics algorithms Simplified version of process-based methods with significant difference in geometrical characteristics among simulation results Sheidaei et al. (2013) NMR response simulation Random walk techniques The most efficient way at present to reproduce travel and relaxation of 1 H in the voids Guo et al. (2016b), Guo and Xie (2017) Finite element methods Finely meshed grid may help raise resolution of simulation results and highlight voids details, but the inconvenience is inserting controlling equations to describe complex relaxation mechanisms Zientara and Freed (1980), Nguyen and Mardon (1995) Finite difference methods like core imbibition, flooding, excavation, centrifugation, or freezing-thawing cycles (Liu et al. 2019), with NMR techniques, to reflect fluid transportation laws and distribution variations, and observe the evolution of pore structure. This requires reasonable experimental design and integrated protocol. Di et al. (2017) used water and gels formulated by deuterium oxide to flood through cores alternately, and the signal change on MRI images outlined the flow channels inside cores. Lai et al. (2020) observed wormhole propagation after acid injection by MRI and T 2 spectrum. It is believed that under the support of image processing and reconstruction technology, these voids can be characterized quantitatively. However, limited resolution MRI images for observing evolution process of pore structure should be supported by other imaging measurements, such as SEM and CT (Safari et al. 2016;Livo et al. 2020). Especially for large pores and cracks that are not suitable for fluid preservation, evaporation or loss of NMR signal sources will seriously influence the results (as shown in Fig. 7). If the sample with large voids in a core holder is continuously flooded to prevent fluid loss, signal of fluids in pipelines will be overestimated.

Conclusions
This paper briefly reviews the development of NMR relaxation theories and methodologies for characterization of rock petrophysical properties and pore structures by low-field NMR. Mature techniques have showed some limitations and deficiencies in practice, which are mainly attributed to signal loss, and ambiguous interpretation caused by low SNR and signal overlapping. To tackle these issues, R & D of NMR instruments and pulse sequences are essential, but it is more realistic and efficient for engineers and researchers with petroleum expertise, rather than physical chemistry or radio physics experts, to overcome their problems by raising the abilities of signal interpretation and experimental design. In summary, porous media relaxation theories have supported applications of low-field NMR technologies, but they are insufficient to handle problems like obtaining explicit pore structure, distinguish pores and fractures, and composition analysis of underground cores, especially for carbonates and shales. This requires researchers to learn more about relaxation mechanisms related to characteristics of various rock sorts, and accelerate speed of theoretical advance. Innovations of experimental design and methodology unleash the potential of existing techniques, like usage of high-field and solid NMR, and combination of imaging measurements with different resolution. It can be expected that NMR will be more interdisciplinary, and core analysis will be a "big data" project, where laboratory experiments should be designed with consideration of addressing unique engineering challenges.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Big wormhole
Non-reactive section Fig. 7 A carbonate sample after acid injection, indicating obvious wormholes through surface that cannot retain fluid. The wormhole and nonreactive section are similarly signal-less on MRI images. (from Lai et al. 2020)