Nuclear Magnetic Resonance with Fast Field-Cycling Setup: A Valid Tool for Soil Quality Investigation

Abstract

Nuclear magnetic resonance (NMR) techniques are largely employed in several fields. As an example, NMR spectroscopy is used to provide structural and conformational information on pure systems, while affording quantitative evaluation on the number of nuclei in a given chemical environment. When dealing with relaxation, NMR allows understanding of molecular dynamics, i.e., the time evolution of molecular motions. The analysis of relaxation times conducted on complex liquid–liquid and solid–liquid mixtures is directly related to the nature of the interactions among the components of the mixture. In the present review paper, the peculiarities of low resolution fast field-cycling (FFC) NMR relaxometry in soil science are reported. In particular, the general aspects of the typical FFC NMR relaxometry experiment are firstly provided. Afterwards, a discussion on the main mathematical models to be used to “read” and interpret experimental data on soils is given. Following this, an overview on the main results in soil science is supplied. Finally, new FFC NMR-based hypotheses on nutrient dynamics in soils are described

1. Introduction

The locution “nuclear magnetic resonance”, also indicated with the acronym NMR, refers to a multifaceted technique which can be applied to liquid [1], semisolid [2], solid [3], and gas [4] phases. It can be used to unveil not only the structure of pure chemical compounds [5], but also the chemical complexity of mixtures [6], as well as that of living organisms [7,8].

Notwithstanding the vast areas of expertise covered by the NMR techniques and their different aspects related to sample preparation, instrumentation, and targeted results, all of them share the same principles [9]. In general, an ensemble of nuclei—each with magnetic properties distributed into different energy levels described by the orientation of their magnetic moments with respect to an applied magnetic field (B₀)—are perturbated by an electromagnetic radiation pulse (B₁) applied for a short time interval along a direction perpendicular to B₀. This causes a transfer of energy quanta and an induction of phase coherence in the ensemble. The excited nuclei then return to the ground state via two different relaxation processes. The first one, also referred to as longitudinal relaxation, refers to the regaining of the magnetization equilibrium in the direction of B₀ with the loss of energy quanta. The second process, indicated as transverse relaxation, is related to the loss of phase coherence and to the decay of the excited magnetization in the direction perpendicular to B₀. Both processes contribute to the overall regaining of the magnetization equilibrium. The time constant associated to the longitudinal relaxation is indicated as T₁ (also referred to as the longitudinal relaxation time), while the time constant of the transverse relaxation is indicated with T₂ (also referred to as the transverse relaxation time).

The signals acquired by nuclear magnetic resonance are time dependent and can be handled in many ways. As an example, the application of the Fourier transform generates plots (also named spectra) where resonance signals are reported versus frequency. According to the number of frequency axes, it is possible to distinguish among mono- (1D), bi- (2D) and three-dimensional (3D) NMR spectroscopy [10]. The disposition of the resonance signals in the spectra (named chemical shift) as well as their splitting (due to the spin–spin interactions) are strongly dependent on the chemical environment of each nucleus. Moreover, the intensity of the signals depends on the number of nuclei producing that signal. For this reason, the NMR spectra can be used not only to elucidates structure and three-dimensional conformation of pure compounds, but also to quantify how many resonating nuclei are in each chemical environment. However, when a complex mixture is analyzed, quantification may be limited by the nature of the sample phase, the type of pulse sequence used for the NMR investigation, and hardware constraints [3,11,12].

The time scale of the relaxing nuclei is affected by the molecular size, and by the physical state of the system under investigation. As it will be explained below, the relaxation mechanisms involve local magnetic fields (generated by the electronic and atomic environments) modulated by molecular motions. When the Larmor frequency of the local magnetic fields is in resonance with that of the observed nuclei, relaxation occurs. Once the NMR experiments are set to measure relaxation times (either T₁ or T₂), the time dependent resonance signals can be handled with the inverse Laplace transform to produce relaxograms (either 1D or 2D), where distribution of relaxation times are represented [13,14,15,16,17]. This technique, named time-domain (TD) NMR, is applied to obtain information on molecular dynamics either for simple systems or complex mixtures [18].

Since the end of the seventies of the twentieth century, a new NMR technique has been emerging [19]. It is based on the analysis of the behavior of complex systems when the strength of the applied magnetic field is quickly changed. The technique, named fast field-cycling (FFC) NMR relaxometry, aims at the evaluation of the molecular motions included in the time scale ranging between 10⁻⁸ to 10⁻³ s [20].

Notwithstanding the term “relaxometry” is generally used when both T₁ and T₂ are accounted for [21,22,23], it is noteworthy that it must be referred only to the contexts where T₁ values are evaluated [20]. Therefore, FFC NMR relaxometry is intended as the technique by which the T₁ values are monitored when the strength of the applied magnetic field is varied within a given proton Larmor frequency interval.

Figure 1 shows the time intervals for the molecular motions investigated by applying the corresponding proton Larmor frequencies on the top side of the scale. In particular, molecular tumbling occurs in the time scale < 10⁻⁷ s (corresponding to proton Larmor frequencies as large as 10⁷ Hz), internal motions occur on a time scale comprised between < 10⁻⁹ s and slightly below 10⁻⁵ s (that is proton Larmor frequencies going from around > 10⁹ Hz to ≈ 10⁵ Hz), chemical exchanges take place in the time interval ≈ 10⁻⁵–10⁻³ s (which corresponds to proton Larmor frequencies going from around 10⁵ Hz to around 10³ Hz), and diffusion occurs on a time scale > 10⁻³ s (that is proton Larmor frequency < 10³ Hz). Spanning five orders of magnitude on the time scale reported in Figure 1, fast field-cycling NMR relaxometry can monitor a wide set of molecular motions going from slow tumbling to slow diffusion. When compared to the other NMR techniques, the advantage of FFC NMR relaxometry is its larger sensitivity for the aforementioned molecular motions as a result of the special design of the hardware system [24].

FFC NMR relaxometry versatility makes the technique quite useful in many fields. As an example, it has been applied to monitor food quality [25,26], traceability [27,28,29], conservation [30], and transformation [31,32,33,34,35,36,37]; to evaluate the physical chemical characteristics of systems useful for drug delivery [38]; to design new contrast agents to be used in diagnostic imaging techniques [39,40,41]; to explore tumor development [42]; to monitor environmental quality [18,43,44], and many other different applications that can be found elsewhere [45,46].

The present review paper is intended to explore the suitability in soil science of FFC NMR relaxometry. Firstly, a general overview of the mechanisms of longitudinal relaxation is provided. Secondly, the basic FFC NMR experiment is described. Then, the mathematical models to interpret FFC NMR relaxometry data are indicated. Finally, the reliability of the technique in soil science is described, thereby also suggesting a molecular mechanism in order to explain how water dynamics in soils can influence movement of nutrients towards plant roots.

2. The Meaning of the T₁ Value

T₁ (i.e., the longitudinal relaxation time) quantifies the time needed to recover the longitudinal component of the magnetization along the direction of the applied magnetic field (conventionally the z-axis) [10]. Its value depends on the fluctuating local electric/magnetic fields that are generated by either unpaired electrons, or nuclear dipoles, or charged particles interacting with nuclear quadrupole moments for nuclei with spin number > ½ (e.g., ¹⁴N), or even anisotropy of the chemical shielding tensor, and finally fluctuating scalar coupling interactions and molecular rotations [10]. Among the aforementioned factors, the molecular motions appear to be very important in affecting the fluctuations of the local electro-magnetic fields [10]. Consequently, the evaluation either of the longitudinal relaxation time or the longitudinal relaxation rate (R₁)—which is the inverse of T₁—can provide valuable information on molecular dynamics.

In soils, the longitudinal relaxation time of water hydrogen nuclei (that is the main component of the soil solution) is influenced by the porous boundaries and the presence of quadrupolar systems [47].

Quadrupolar nuclei are those having quantum spin number (S) > ½. Such nuclei show non-spherically symmetric charge distribution, which is represented by a non-zero quadrupolar momentum (μ_Q). The latter affects the effective strength of the applied magnetic field (B₀) felt by the investigated nucleus (that is ¹H for water in soils), thereby shortening its longitudinal relaxation time value. However, the μ_Q values are not only dependent on the size and the charge of the nucleus, but also on the symmetry of the structure of the molecular system where that nucleus lays. Therefore, it must not be a surprise if some of the nuclei that are of interest for soil scientists reveal the quadrupolar effect despite their nuclear quantum spin number being ½ [47]. In fact, while the quadrupolar effect is null in chemical environments with cubic local symmetry (point groups T_d and O_h), most coordination environments in minerals have lower symmetries [47,48].

As stated above, the longitudinal relaxation is also dominated by possible strong relaxation sinks present at the pore surface. These are caused by temporary adsorption of water molecules by means of weak interactions [44]. When water tumbling is reduced due to the interactions with the solid surface, the efficiency of the dipolar ¹H–¹H interactions increases. Hence, reduction in the ¹H T₁ value is observed [44]. In particular, the strength of the water–solid surface interactions is influenced by the pore size of the porous system. In fact, the smaller the pore size (such as in clayey soils), the slower the molecular motions due to space restrictions. This leads to the aforementioned more efficient ¹H–¹H dipolar interactions and ¹H T₁ reduction. Conversely, as molecular mobility increases because of pore size enlargement (e.g., in silty and sandy systems), the ¹H–¹H dipolar interactions weaken and longer longitudinal relaxation times are measured [44]. In other words, shorter T₁ values can be associated with the movement of water into small-sized pores. These are typical features of clay-rich systems such as clay primary particles and small soil aggregates. Longer T₁ values may be led back to the movement of water in silt, sand particles, and large aggregates. However, because the size and number of particles in soils is heterogeneous, water molecules perceive a wide variety of differently sized pores, which results in broad continuous distribution of longitudinal relaxation times [49]. Water molecules moving inside the smallest pores provide a range of T₁ values closer to the shortest T₁ limit. On the other hand, water molecules that move inside the largest pores produce T₁ values closer to the longest T₁ limit. All T₁ values between the two limits are due to water molecules moving inside pores having sizes comprised between the two extremes [49].

3. The Fast Field-Cycling NMR Experiment and the Evaluation of the T₁ Value

3.1. The Inversion Recovery Sequence and the Basic FFC NMR Experiment

Figure 2A reports the scheme for the classic 180–τ–90 inversion recovery (IR) experiment designed to measure the value of the longitudinal relaxation time. The three periods of preparation, evolution and acquisition are also indicated.

At the end of the preparation period (which is a waiting time needed to let the magnetization completely regain its alignment along the direction of B₀), an inversion pulse is applied in order to flip the magnetization along -z. Then, the evolution period is set to be varied within a series of progressively increasing time intervals, τ. During each τ, the magnetization evolves towards the equilibrium condition along the direction of the applied magnetic field. Finally, at the end of each evolution time with the duration of τ, a 90°_x pulse is applied in order to generate the observable and acquire the free induction decay (FID) signal [50].

Signal intensity depends on the τ values as described in Equation (1), where $I\left(\tau \right)$ is the signal intensity at the end of the evolution period with duration τ, and $I_0$ is the signal intensity when τ = 0 [50].

$$I\left(\tau \right)=I_0\left(1-2e^{-\frac{\tau }{T_1}}\right)$$

(1)

The use of any fitting software where $I\left(\tau \right)$–vs–τ is reported, allows the achievement of both $I_0$ and T₁.

The inversion recovery sequence is quite simple to be applied. However, it is designed to be used in high-field NMR spectroscopy and, hence, to retrieve information about fast molecular motions in purified molecules. As a general remark, it must be also mentioned that the acknowledgment of the T₁ values via the inversion recovery sequence allows one to understand the 3D conformation of pure systems [50].

When dealing with very complex mixtures (e.g., dissolved organic matter, soils soaked with water solutions, etc.), one cannot apply NMR spectroscopy to recognize the three-dimensional structure of a molecular system. In fact, mixture complexity makes spectroscopy only suitable to gain information on the general chemical nature of the mixture components. In other words, the NMR signals are assigned to general classes of functional moieties such as aromatic groups (resonating in the chemical shift interval 160–110 ppm), carbohydrate units (acquired in the chemical shift interval 110–40 ppm), alkyl systems (resonating in the interval 40–0 ppm), etc. [1]. Therefore, the application of the inversion recovery sequence to unveil 3D structures and conformations is useless. However, the IR experiment can still be applied to understand the molecular dynamics of the aforementioned classes of compounds. As a matter of fact, the direct relationship between T₁, temperature, and correlation time (τ_C, that is, the time needed for a molecule to rotate one rad or to move within a distance corresponding to its length) [51] permits one to understand how strongly the motion of molecular groups is restrained [10].

There are at least two limitations in the use of the inversion recovery sequence to study molecular dynamics of complex mixtures. The first one is the necessity to apply variable temperature experiments in order to obtain information about correlation times, while the second limit is due to the application of static high magnetic fields. Variable temperatures can be detrimental when temperature sensitive organic systems are analyzed. In fact, the latter can decompose, and artifacts can be obtained. Static high magnetic fields needed to acquire high sensitivity spectra (i.e., spectra with a signal-to-noise ratio at least > 10) allow detection of fast molecular motions only (Figure 1). As a consequence, the dynamics of slow motions, such as those occurring on the surface of porous systems, are missed. Therefore, the fast field-cycling sequence reported in Figure 2B must be applied [18]. In fact, it allows the monitoring of the T₁ values of complex systems in a range of proton Larmor frequencies going from 0.01 up to 40 MHz with a single instrument. In addition, the correlation times can be measured at room temperature or in conditions that are able to prevent thermal decomposition of organic matter. Of course, being a low field technique, all the spectroscopic information is lost.

In the FFC NMR sequence of Figure 2B, the preparation time corresponds to the application (for a fixed period of time also referred to as polarization time, T_POL) of a polarization field (B_POL). Its value can be either non-null or null. In the first case, a pre-polarized (PP) sequence is achieved, while in the second case a non-polarized (NP) sequence is retrieved.

The pre-polarization is needed to generate the magnetization that evolves to reach a new equilibrium condition (evolution period) under the action of a relaxation field (B_RLX) applied for a variable period τ. The acquisition time starts at the end of the evolution: the magnetic field intensity is switched to a new value (indicated as B_ACQ), while a 90° pulse is applied to generate the observable magnetization, and the FID is finally acquired.

The use of the PP sequence is recommended whenever the intensity of the relaxation field is very low, therefore enhancement in sensitivity is needed. The crossover field between the PP and the NP sequences might be fixed empirically at ${\omega }_{RLX}=\frac{{\omega }_{POL}}{2}$, where ${\omega }_{RLX}$ and ${\omega }_{POL}$ are the B_RLX and the B_POL proton Larmor frequencies expressed in MHz, respectively.

In Figure 2B, a switching time (SWT) is also indicated, i.e., the time needed to switch between different magnetic field intensity values. SWT is usually set at ca. 3 ms. Shorter SWT values can be only applied providing that peculiar electronic precautions are adopted [52].

3.2. How to Obtain the T₁ Value from the FFC NMR Experiment

Figure 3 shows a typical dataset for the decay and recovery curves obtained by the PP and NP sequences, respectively. The full lines represent the data fitting by using the bi-exponential form of Equation (2) for the PP sequence and Equation (3) for the NP sequence. The dataset depicted in Figure 3 has been acquired by saturating a quartz sand with Milli-Q grade water (see Appendix A for the details of the experiment).

$$I\left(\tau \right)=a+{\displaystyle \sum }_{i=1}^N{b}_i{e}^{-\frac{\tau }{T_{1i}}}$$

(2)

$$I\left(\tau \right)=a+{\displaystyle \sum }_{i=1}^N{b}_i\left(1-e^{-\frac{\tau }{T_{1i}}}\right)$$

(3)

In both Equations (2) and (3), $I\left(\tau \right)$ is the signal intensity at the selected τ value, a is the offset and b_i is the magnetization intensity at the Boltzmann equilibrium of the i-th component of the molecular motion at each fixed B_RLX intensity [18].

The limit of Equations (2) and (3) lies in the need to make assumptions on the number of exponential components to be used for data fitting. From a mathematical point of view, the larger the number of exponential components, the better the quality of the fitting, at least when a hyper-corrected model is used. Therefore, any assumption made on the number of components to be considered in the aforementioned equations is arbitrary, and, hence, questionable, when the system under consideration is unknown and complex (e.g., the case of soil systems). In order to overcome the aforementioned limits, the stretched Equations (4) and (5), for the pre-polarized and non-polarized sequence, respectively, can be applied [18].

$$I\left(\tau \right)=a+b\times e^{-{\left(\frac{\tau }{T_1}\right)}^k}$$

(4)

$$I\left(\tau \right)=a+b\times \left(1-e^{-{\left(\frac{\tau }{T_1}\right)}^k}\right)$$

(5)

In Equations (4) and (5), a and b have the same meaning as in Equations (2) and (3), and k is a heterogeneity parameter that is related to the stretching of the decay/recovery processes. The use of Equations (4) and (5) offers the advantage to deal with a large variety of different behaviors within a single model, thereby avoiding arbitrary assumptions about the number of exponentials to be applied in modelling the relaxometric results.

4. From Time Domain to Time Domain

The mathematical models given in Equations (2) and (3) can be considered as a discrete representation of the T₁ distribution of a molecular system. They provide valuable information only when the longitudinal relaxation times in a complex system are very different to each other. As an example, with the assumption of a biexponential decay/recovery for the water in equilibrium with the surface of a quadrupolar-less-nuclei quartz sand used as porous system model (Figure 3), the T₁ values reported in

Table 1. Longitudinal relaxation times for the quartz sand sample added with Milli-Q grade water as described in the caption of Figure 3. ${\omega }_L$ is the proton Larmor frequency of the applied B_RLX.

${\mathit{\omega }}_{\mathit{L}} \left(\mathbf{MHz}\right)$	$$\begin{gathered} {\mathit{T}}_1^{\mathit{s}\mathit{l}\mathit{o}\mathit{w}} \\ \mathbf{ms} \end{gathered}$$	$$\begin{gathered} {\mathit{T}}_1^{\mathit{f}\mathit{a}\mathit{s}\mathit{t}} \\ \mathbf{ms} \end{gathered}$$
0.1	129	34
30	445	96

are obtained.

Here, $T_1^{fast}$ may refer to the relaxation time of water molecules in the closest proximity of the solid surface, while $T_1^{slow}$ may be the relaxation time of the water molecules flowing farther from the surface (Figure 4A). However, according to the mechanisms outlined in Section 2, one can also stand that $T_1^{fast}$ is due to water molecules moving in pores smaller than those where water molecules generating $T_1^{slow}$ lay (Figure 4B). Actually, without the aid of any additional information, there is no reason to place trust in one of the different aforementioned hypotheses. It is also possible that a combination of all the aforesaid dynamics may occur (Figure 4C). However, this latter information has not been accounted for in the biexponential assumption used to apply Equations (2) and (3).

Due to the very good fitting (R² > 0.99) provided by the biexponential assumption (Figure 3), there is no mathematical reason to apply for the three-exponential form of Equations (2) and (3) in order to introduce a hypothetical complexity in the dynamic behavior of water molecules in equilibrium with the solid surface of the quartz sand. Certainly, the approximation given by Equations (4) and (5) does not help in solving the complexity of water dynamics in quartz sand. In fact, those equations provide only one average T₁ value, and a k ≠ 1. It is not yet clear what the relationship between the heterogeneity parameter and the number of the different T₁ components associated with water dynamics is.

When the various components of the molecular dynamics in multi-phase frames are described by longitudinal relaxation rates (the values of which are very close), their T₁ distributions can be more suitably obtained by applying an inverse Laplace transform. It can be expressed as in Equation (6) when pre-polarized experiments are performed and in the form of Equation (7) when the non-polarized experiments are carried out [18].

$$I\left(\tau \right)={{\displaystyle \int }}_{T_1^{min}}^{T_1^{max}}D\left(T_1\right)e^{-\left(\frac{\tau }{T_1}\right)}d\left(T_1\right)+\sigma$$

(6)

$$I\left(\tau \right)={{\displaystyle \int }}_{T_1^{min}}^{T_1^{max}}D\left(T_1\right)\left(1-e^{-\left(\frac{\tau }{T_1}\right)}\right)d\left(T_1\right)+\sigma$$

(7)

In Equations (6) and (7), $T_1^{min}$ and $T_1^{max}$ are the suitable limits within which all the T₁ values range; $D\left(T_1\right)$ is the relevant distribution function that must be determined by solving either Equation (6) or Equation (7); σ is a parameter accounting for a suitable unknown noise component. The most likely distribution of T₁ values may be obtained when some constraints, such as variance of the experimental data or smoothness of the solution, are taken into account.

Two algorithms have been developed to switch from $I\left(\tau \right)$ to $D\left(T_1\right)$. These are the continuous distribution, also referred to as CONTIN, [53,54], and the uniform penalty regularization, also referred to as UPEN [13,14,15]. The two algorithms differ in the smoothing procedure used. Nevertheless, it worth underlying that the two algorithms provide de facto similar T₁ distributions, also referred to as relaxograms (Figure 5), regardless of the procedure used to obtain the most probable distribution of relaxation times.

Figure 5 shows the relaxogram obtained by applying the UPEN algorithm to the data points reported in Figure 3. It is clear, now, that water dynamics in the quartz sand is more complex than the simple data interpretation based on the biexponential form of Equations (2) and (3). In fact, the inspection of Figure 5 reveals that three different T₁ components are distinguishable. They are centered at around 51, 126, and 292 ms in the relaxogram acquired at the ¹H Larmor frequency of 0.1 MHz, while they are displaced at around 80, 185 and 505 ms in the relaxogram obtained by using the proton Larmor frequency of 30 MHz. Both relaxograms in Figure 5 do not show any band at around 2.5 s, which is the typical longitudinal relaxation time of free water [38]. This finding can only mean that all the water is trapped inside the quartz sand. It is, then, conceivable that water can be subjected to a horizontal surface diffusion (Figure 6A); a diffusion in the spaces among quartz sand particles, which can be indicated as interparticle diffusion (Figure 6B); and a diffusion inside the internal pores of each sand particle, also referred to as intraparticle diffusion (Figure 6C) [49].

5. The Nuclear Magnetic Resonance Dispersion (NMRD) Profile and Its Modeling

The T₁ values obtained via Equations (2) to (5) can be plotted versus B_RLX proton Larmor frequency in order to obtain an NMRD profile. Its typical shape is a sigmoidal curve (Figure 7), which is mathematically described by Equation (8) [20].

$$R_1=\frac{1}{T_1}=\frac{2{\tau }_c}{1+{\left({\omega }_L{\tau }_C\right)}^2}$$

(8)

In Equation (8), R₁ is the longitudinal relaxation rate (R₁ = 1/T₁) which represents the distribution of the motion frequencies in a molecular system; ${\omega }_L$ is the B_RLX intensity given as proton Larmor frequency; and ${\tau }_c$ is the relevant correlation time [10]. The latter, indeed, describes the random molecular motions of molecular systems either in solution or in porous media. In particular, the larger the ${\tau }_c$ value, the slower the molecular motions, thereby accounting for the existence of restraints in the motional degrees of freedom. Conversely, as a molecule encompasses faster motions in larger spaces, shorter T₁ values are expected [18] (see also Appendix B for further mathematical details).

Luchinat and Parigi [55] reported that protein dilute aqueous solutions can be investigated by using Equation (9), also referred to as the Bloembergen–Purcell–Pound (BPP) model. The latter has been obtained for aqueous systems where bulk water is dominant [51]. Therefore, protein dynamics can be indirectly revealed by using water molecules as target for the relaxometric investigations. As will be clarified below, Equation (9) as well as its modifications [55] can also be applied to recognize the properties of environmentally relevant porous systems via the indirect observations made on water molecules in equilibrium with the solid surfaces [49,56].

$$R_1=\frac{1}{T_1}=\alpha +\beta \left[\frac{0.2{\tau }_c}{1+{\left({\omega }_L{\tau }_C\right)}^2}+\frac{0.8{\tau }_c}{1+{\left(4{\omega }_L{\tau }_C\right)}^2}\right]$$

(9)

In Equation (9), all the terms in the squared brackets have the same meaning as in Equation (8), while α is the high-field relaxation rate and β is a parameter that relates to the dipolar interactions. The latter contains the Planck constant, the proton quantum-spin number, the gyromagnetic ratio, and the electron–nuclear hyperfine coupling constant (which in turn accounts for the interactions between resonant protons and unpaired electrons). The higher the β value, the stronger the dipolar interactions responsible for the relaxation.

It is well recognized that neither Equation (8) nor Equation (9) satisfactorily fit the experimental data points when a stretched dispersion is retrieved. The stretching depends on how complex the re-orientational dynamics within the molecular systems are, as well as on the heterogeneous distribution of intermolecular dipole couplings and proton exchange rates [18]. An example of stretched dispersion is given in Figure 8, where the NMRD profile of one of the water saturated soils analyzed in Reference [57] is reported.

In order to account for the NMRD stretching, Halle and coworkers [58] developed the free model analysis given in Equation (10):

$$R_1={\displaystyle \sum }_{i=1}^N\frac{C_i{\tau }_C^i}{1+{\left({\omega }_L{\tau }_C^i\right)}^2}$$

(10)

where the i index refers to the different components of the motion, ${\tau }_C^i$ is the correlation time of the i-th component, and C_i is a fitting parameter. The sum of the C_i values constitutes the mean square fluctuation bearing the information on the equilibrium structure of the system and is not dependent on its dynamics.

To limit the number of components to be used in the free model analysis, Halle et al. [58] suggested the application of the F-test. In other words, ‘any fit obtained by using N Lorentzian terms to M data points J(ω_i), with errors σ_I, provides a χ²(N) value. The latter must be compared with the χ²(N + 1) value obtained by applying N + 1 Lorentzians. If the fit improves, the ratio F_{(N, N+1)} = χ²(N)/χ²(N + 1) increases [i.e., χ²(N + 1) < χ²(N)]. The procedure must be iterated until F_{(N+m−1, N+m)} > F_{(N+m, N+m+1)}, where m is an integer and ≥0. When the latter condition is satisfied, the acceptable number of Lorentzians to be applied is N + m’ [58].

It is noteworthy that the free model analysis is only a convenient mathematical approach to fit the experimental data such as those shown in Figure 8. The set of $\left\{C_i,{\tau }_C^i\right\}$ parameters from Equation (10) have no physical meaning, unless independent information suggests that the investigated system can be modelled by a fixed number of Lorentzians. In the latter case, a direct physical interpretation of the parameters can be attempted. Nevertheless, according to Halle et al. [58], the aforementioned parameters can be used to calculate a weight-averaged correlation time as in relation (11):

$$\langle {\tau }_C\rangle =\frac{{{\displaystyle \sum }}_{i=1}^N{C}_i{\tau }_C^i}{{{\displaystyle \sum }}_{i=1}^N{C}_i}$$

(11)

The combination between the free model analysis given in Equation (10) and the BPP model reported in Equation (9), provides the BPP free model analysis as in Equation (12) [49,55,56].

$$R_1=\alpha +\beta \left[0.2{\displaystyle \sum }_{i=1}^N\frac{C_i{\tau }_C^i}{1+{\left({\omega }_L{\tau }_C^i\right)}^2}+0.8{\displaystyle \sum }_{i=1}^N\frac{C_i{\tau }_C^i}{1+{\left({\omega }_L{\tau }_C^i\right)}^2}\right]$$

(12)

All the parameters in Equation (12) have been already defined. The set of $\left\{C_i,{\tau }_C^i\right\}$ parameters from the BPP free model analysis are used to obtain the weight-average correlation time given in relationship (11).

Instead of the aforementioned models, the NMRD profiles can be analyzed by applying a very elegant model provided by Korb and coworkers [59]. This model was designed to interpret NMRD data points acquired for systems presenting paramagnetic impurities (Equation (13)).

$$R_1\left({\omega }_I,{\omega }_S\right)=R_0+K{\tau }_m\left[3ln\left(\frac{1+{\omega }_I^2{\tau }_m^2}{{\left(\frac{{\tau }_m}{{\tau }_s}\right)}^2+{\omega }_I^2{\tau }_m^2}\right)+7ln\left(\frac{1+{\omega }_s^2{\tau }_m^2}{{\left(\frac{{\tau }_m}{{\tau }_s}\right)}^2+{\omega }_s^2{\tau }_m^2}\right)\right]$$

(13)

where ${\omega }_I$ is the intensity of B_RLX expressed as ¹H Larmor frequency; ${\omega }_S=659{\omega }_I$ is the Larmor frequency of the unpaired electron in the paramagnetic species; R₀ is the sum of the water bulk relaxation time of the order of 2.5 s and the frequency independent water-bound contribution, which depends on the amount of water-binding sites. The term in squared brackets in Equation (13) is the surface contribution relevant to the surface diffusion of H₂O molecules near the paramagnetic relaxation source (Fe³⁺, for instance). K is provided by the relationship:

$$K=\frac{\pi }{30{\delta }_{water}^3}{\sigma }_S{\rho }_{water}{S}_{P,NMR}{\left({\gamma }_I{\gamma }_S\eta \right)}^2{S}_S\left(S_S+1\right)$$

(14)

where ${\delta }_{water}$ = 0.3 nm is the width of a monomolecular water layer; ${\sigma }_S$ is the density of the paramagnetic ions at the pore surface (i.e., sources of relaxation expressed in a number of paramagnetic spins per cm²), ${\rho }_{water}$ = 1 g cm⁻³ is water density; $S_{P,NMR}$ (given in m² g⁻¹) is the NMR-specific surface area calculated as $S_P\times F$ (where $S_P$ is the surface area achieved by the Brunauer–Emmet–Teller (BET) analysis, while $F<1$ is the solid-to-liquid ratio); ${\gamma }_I$ and ${\gamma }_S=659{\gamma }_I$ are the proton and electron gyromagnetic ratios, respectively (for the most abundant paramagnetic species). Finally, $S_S$ is the spin quantum number of the paramagnetic species. ${\tau }_m$ is the surface diffusion correlation time, while ${\tau }_s$ is indicated as surface residence time. The former represents the “hopping” time of water among surficial binding sites, while the latter (${\tau }_s$) is the time of residence of water on the surface of the porous medium [18]. The $\frac{{\tau }_s}{{\tau }_m}$ ratio is referred to as NMR surface wettability [60,61].

6. How to Choose the Right Model to Correctly Interpret NMRD Profiles in Soil Science

In the previous paragraph, several mathematical models to fit the NMRD profiles have been provided. How should one choose the best one for a given set of experimental data dealing with soils and soil related systems? Unfortunately, there is not a general rule. In order to apply any of them, either the congruence between the results coming out from the models and the general properties of the porous system under investigation must be accounted for or, more simply, when not additional information is available, one choses the model which best fits the experimental data.

6.1. The Application of the BPP Model: The Quadrupolar-Less Nuclei Quartz Sand

Let us consider, for instance, a sample containing no quadrupolar nuclei such as quartz sand, which has been used as a sample model for the application of Equations (2), (3), (6), and (7). As reported above, a three-exponential behavior appeared to explain the decay/recovery curves acquired at different B_RLX intensities. For this reason, the three NMRD dispersion curves reported in Figure 9 have been obtained.

Visual inspection of Figure 9 allows one to differentiate among a fast component (black data points), an intermediate component (red data points), and a slow component (blue data points) of the water molecular motion. The terms “fast”, “intermediate”, and “slow” refer to the values of the longitudinal relaxation rates in the entire range of ¹H Larmor frequencies which was accounted for. In fact, the fast component revealed the fastest R₁ values, the slow component showed the slowest R₁ values, while the R₁ values for the intermediate component fell all between those measured for the other two components. The three NMRD profiles can be assigned to water molecules subjected to intraparticle diffusion (fast component), to water systems subjected to horizontal diffusion (intermediate component), and to water molecules subjected to interparticle diffusion (slow component) (Figure 6). According to this qualitative evaluation and to the meaning of the T₁ (or R₁) values reported in Section 2, one may expect that the correlation times obtainable by the models listed in Section 5, must follow the order: ${\tau }_C^{fast}>{\tau }_C^{intermediate}>{\tau }_C^{slow}$.

Among the different NMRD models, the only one providing the expected results was the BPP model depicted in Equation (9).

Figure 10 reports all the parameters accounted for in Equation (9). The values of the parameter α, representing the high-field relaxation rate, follow the relative position of the three profiles in Figure 9. Therefore, ${\alpha }_{fast}>{\alpha }_{intermediate}>{\alpha }_{slow}$.

The parameter β is an indication of the strength of the dipolar interactions occurring between water molecules and the surface of the quartz sand. As reported above, the stronger the dipolar interactions, the shorter is the T₁ and the faster results the R₁. For this reason, as expected by the relative position of the three NMRD profiles, the β values follow the order ${\beta }_{fast}>{\beta }_{intermediate}>{\beta }_{slow}$.

Finally, the correlation times are in the order: ${\tau }_C^{fast}>{\tau }_C^{intermediate}>{\tau }_C^{slow}$, which confirms the expectations reported above.

6.2. When the Free Model Analysis is the Best Choice. Two Case Studies: Soil and Biochar

The quartz sand discussed above is a very simple system. When dealing with more complex materials, such as soils, the simple BPP model reported in Equation (9) may not be applied. As an example, Laudicina et al. [62] studied the effects of afforestation with four unmixed plant species (Eucalyptus camaldulensis Dehnh., E. occidentalis Endl., Pinus halepensis Mill., and Cupressus sempervirens) on the soil–water interactions in a semiarid Mediterranean region by applying, among the others, fast field-cycling NMR relaxometry. The six-parameter ($c_1, c_2, c_3, {\tau }_1, {\tau }_2, {\tau }_3$) free model analysis (that is, Equations (10) and (11)) appeared to provide the best fist for the experimental data. The authors found that the soils under the different plant species may trap water in clay-, loam-, and sand-type pores. In particular, the longest correlation time describes the re-orientational behavior of water confined in clay-type pores, the intermediate value refers to water restricted in loam-type pores, and the shortest correlation time is due to water freely moving in sand-type pores (Figure 11).

The free model analysis was also applied in Conte et al. [63] to unveil the nature of the water–biochar interactions. Biochar is a carbonaceous material used to improve soil quality via mechanisms which are still to be clarified [64,65,66,67]. In the aforementioned study, NMRD profiles of a water saturated biochar were acquired at different temperatures. The evaluation of the ${\tau }_C$ temperature dependence suggested that water molecules can be hooked to the biochar surface via different types of interactions. On one hand, the electron-rich oxygen atom in water may interact with the electron-poor metals displaced on the biochar surface. On the other hand, the water electron-deficient hydrogen atoms can interact with the electron-rich aromatic moieties of the biochar backbone. Finally, water may also form H-bonds with the functional groups usually present on the biochar surface [68]. These interactions may justify the exceptionally large biochar affinity for nitrate, phosphate and many other plant nutrients [66,67,69], as well as its role in improving soil structure when it is applied as an amendment [70,71], thereby allowing soil fertility enhancement [65]. In fact, water may act as a “glue” layer between biochar and soil components. As a consequence, improvement in soil structure can be achieved. The same “glue” layer may also permit the adhesion of anions and cations on the soil/biochar surface [23]. Therefore, gradual and constant nutrient supply to plant roots can be obtained.

6.3. The Combination between the BPP Model and the Free Model Analysis: The Saltmasrh Sediment Case Study

The models discussed hereinabove were not appliable to evaluate the characteristics of sediments sampled in a saltmarsh located nearby Trapani (Sicily, Italy) [56]. In this case, the BPP free model analysis given in Equation (12) provided not only the best fit, but also the best numerical results explaining the number of complementary information available on the sediment samples. Indeed, ‘FFC-NMR relaxometry data highlight[ed] that water dynamic properties are linked to the sedimentological features of the saltmarsh sediments as well as with metal analyses and benthic foraminiferal density evaluation’ [56]. As matter of fact, it was possible to ‘subdivide the [sampled sedimentary core] into three time intervals, each characterized by particular biological and geochemical data. In particular, […] the lower part of the studied core corresponds to the interval before the inundation occurred on 1965. During this period, sediments were characterized by low anthropogenic pollution, and the studied area was utilized for the production of halite. The saline pans were abandoned between 1965 and 1990, and agricultural production was intensified. Because of this, runoff increased the deposition of inorganic and organic pollutants into ponds. Level of metals decreased, and benthic foraminiferal density improved after the institution of the natural reserve and with the restarting of the salt production. This latter has permitted to connect again the ponds with the sea’ [56].

6.4. NMRD Evaluation by the Wettability Model: Application to Biochar

The model obtained by Korb coworkers [59] (Equation (13)) can be applied only under the hypothesis that all the parameters reported in relation (14) are known. This model has been designed for the investigations of cement-based materials, plaster pastes and petroleum fluids [61]. However, recently it was also used in a paper by Bubici et al. [60] where biochar water NMR wettability has been investigated.

Water wettability in biochars is a very important parameter to be acquired when the aforementioned carbonaceous material is applied to soils. In fact, it is well recognized that the presence of unsuitable type of biochar or applied amount (due to either uncontrolled forest fires or deliberate application for improvement of soil quality) may cause soil repellency in such a way that affects the ability of soil in absorbing water, inhibiting microbial activity, altering filter, as well as storage, buffer, and transformation functions of the soils [72,73,74]. It is worth highlighting that up to now the only direct macroscopic way to measure wettability of a porous systems is via the evaluation of contact angles [75,76], although a new predictive, but yet to be confirmed, model recently appeared [77].

Contact angle (CA) is the angle between the intersection of the liquid–solid and the liquid–vapor interfaces. It can be geometrically evaluated by considering in the droplet profile the line tangent to the contact point along the liquid–vapor interface (Figure 12A,B). The technique is based on the surface tension of the liquid. Surface molecules in a liquid do not have neighboring molecules in all directions to provide a balanced net force. By contrast, they are pulled inward by the neighboring molecules, thereby creating an internal pressure. Therefore, the liquid minimizes its surface free by contracting its surface area. The more affine to the liquid the solid surface is, the more spread the liquid on the solid surface (Figure 12A). In contrast, by decreasing the affinity between the two phases, the liquid beads the solid (Figure 12B). Hence, high wettability is achieved when small Cas are measured, whereas low wettability is retrieved when Cas are large.

Figure 12C reports the contact angle-vs-NMR wettability as obtained by the data reported in Bubici et al. [60] for three different gasification produced biochars. In particular, the lower the contact angle, the larger the results regarding the dynamical surface affinity of water (i.e., NMR wettability). This possibly suggests that wettability mechanisms from a nano-scale (i.e., FFC NMR relaxometry measurements) up to a macro-scale (i.e., contact angle evaluation) dimension are invariant.

According to these results, it can be highlighted that the importance of the NMR wettability relies on the possibility to monitor water dynamics at a nanometer length scale. This enables one to predict both the possible transformation mechanisms occurring to biochar as it is applied to soils, and nutrient availability when soils are enriched with biochar [60].

6.5. The FFC NMR Modeling via Molecular Dynamics Simulations

The discussion above evidenced that the choice of the right FFC NMR model to apply to soils and soil related materials is left to investigator’s sensibility. However, very recently, a new study dealing with the combination between FFC NMR modeling and molecular dynamics (MD) simulations appeared [78]. The authors showed that NMRD profiles and temperature trends can be calculated with a very good approximation by using MD simulations. Accordingly, the right mathematical model to be applied to NMRD data can be chosen on the basis of the results of the aforementioned simulations. The limit of this approach lays in the fact that the authors studied simple ionic liquids. Therefore, their procedure does not account for the complexity of soil systems. In fact, the presence of many different inorganic and organic moieties in soils prevents the application of molecular dynamics simulations due to the enormous computer time required to calculate the trajectories of thousands different nuclei. Very likely, MD simulations could be applied successfully on soils when quantum computing will be easily accessible.

7. Understanding Soils by Fast Field-Cycling NMR Relaxometry

7.1. The Role of Soil Pores in Water Dynamics

It is well known that the ability of water to infiltrate into soils depends on soil pore sizes [44]. In particular, three different kinds of pores are recognized [79,80]. The “residual pores” (RP) have a size of ≤ 0.5 μm. Here, strong chemical interactions at a molecular level are supposed to occur between pore wall boundaries and water molecules. Water is strongly trapped in the soil system, thereby becoming unavailable for plant nutrition. Pores with size ranging from 0.5 to 50 μm are referred to as “storage pores” (SP). In these pores, water can be retained, but released against gravity. This permits exchanges and diffusion of nutrients within soil pores and, as a consequence, plant nutrition. Finally, pores with sizes ≥ 50 μm are indicated as “transmission pores” (TP). Here, water is more freely moving than in residual and storage pores. This means that it can be easily leached, thereby leading to nutrient loss and soil fertility reduction.

In the last years, the dynamics of water in the aforementioned soil pores has been validated by NMR relaxometry, either at fixed magnetic field strength or by fast field-cycling setup. In particular, Maccotta et al. [56], Laudicina et al. [62], Pohlmeier et al. [81], Bayer et al. [82], Stingaciu et al. [83], Conte et al. [84], and Haber-Pohlmeier et al. [85,86] reported about the direct relationship between soil pore size and relaxation time values (a detailed mathematical approach is reported in Appendix B). All the aforementioned authors showed that the smaller the soil pore size, the faster the NMR relaxation rate, due to the better efficiency of the ¹H–¹H dipolar interactions between water molecules and pore wall boundaries. Moreover, it has been possible to state that the relaxation rate of water in soil pores follows the order: $R{P}_{H_{2}O}>S{P}_{H_{2}O}>T{P}_{H_{2}O}$ which is the inverse of the motion rate, i.e., $R{P}_{H_{2}O}<S{P}_{H_{2}O}<T{P}_{H_{2}O}$. However, soil pore system complexity, where residual, storage and transmission pores coexist in different relative amounts according to soil nature [44], allowed all the aforementioned authors to describe the motion of soil water as bimodal. In other words, water moves in the smallest sized pores by 2D diffusion which can be described as the horizontal motion of water—and nutrients dissolved therein—towards plant roots (Figure 13). Diffusion rate is mediated by the interactions between water molecules and the surface of the pore wall boundaries. The stronger the interactions, the slower the diffusion. For this reason, it is possible to argue that among all the different types of pores, the storage ones are mainly involved in the horizontal transport of water and nutrients towards plant roots.

Water molecules “hooked” on the soil surface can be replaced by bulk water molecules. For this reason, while the former “jump” from the surface to move to the bulk, the latter move on the soil surface, thereby contributing to water diffusion (Figure 13). The frequency of the 3D jump—also referred to as vertical motion of water and dissolved nutrients—is affected by the strength of the interactions between water and pore wall boundaries. The stronger the interactions, the lower the 3D jump frequency. According to this mechanism, it is possible to state that the 3D jumps are common in soils richer in transmission pores, thereby allowing a reduction in soil fertility as a consequence of water and nutrient leaching towards the lowest soil horizons.

7.2. FFC NMR to Quantify Soil Erosion

Soil erosion is a serious problem threatening both land management sustainability and water resource development. It is due to the expansion of agricultural areas because of an increasing need of food supply.

Despite the recognition of soil erosion environmental hazard, there is still a lack of homogeneity in the methodologies used to quantify soil loss associated with the development of agriculture and reduction in forest areas.

Traditionally, the typical techniques applied to monitor soil erosion are short-duration field monitoring programs employing erosion pins and erosion plots, while longer-term experimental studies use soil loss plots, the interpretation of repeated aerial photographs, and the use of satellite imagery [87]. All of the aforementioned techniques are limited by problems related to spatial and temporal variability, operational difficulties, costs and questionable reliability of the resulting data [87]. For this reasons, the more reliable technique based on the artificial radionuclide ¹³⁷Cs has been developed and applied worldwide [87]. The attraction for ¹³⁷Cs investigations lies in the considerations that it has been introduced worldwide in the atmosphere as a consequence of the development of nuclear weapons between the fifties and the seventies of the 20th century. Due to both precipitation phenomena and its strong affinity for soil and sediment particles, ¹³⁷Cs is now present on a global scale. Therefore, monitoring its fate provides useful information on soil transformations [87]. However, ¹³⁷Cs analyses require special equipped labs as well as special precautions because of its radioactive properties. For this reason, the question on the possibility to use FFC NMR relaxometry to quantify soil erosion has risen [57,88,89]. The advantage of the technique is its low cost, the handling of non-radioactive nuclei, and its worldwide diffusion. In fact, most of the NMR labs are also equipped with FFC NMR benchtop machines. However, the little number of NMR labs dealing with environmental problems, and especially soil science, may be the bottleneck for the development of FFC NMR in monitoring soil erosion.

In order to understand how FFC NMR relaxometry can be applied to account for soil erosion, let us consider that the vertical motion described in Section 7.1 is not only related to water and nutrient leaching towards the lowest soil horizons, but also to the surface-water macroscopic flows. Both aforementioned water movements can be described by the concept of “connectivity” [57,88,89]. The latter refers to all the ‘processes involving a transfer of matter, energy, and/or organisms within or between elements of a system such as landscapes, basins and soils’ [90]. This definition implies that a transport vector, such as water, must be involved in moving materials over a range of space and time scales, thus permitting investigations of the effects of heterogeneities of complex and inhomogeneous physical systems [91,92,93].

Three different types of connectivity have been recognized [94]: (i) landscape connectivity, related to the physical coupling of landforms; (ii) hydrological connectivity, referred to the water transfer through the basin; (iii) sedimentological connectivity, related to the sediment transfer. Recently, the landscape connectivity has been also referred to as structural connectivity [95]. This represents the extent to which landscape units are contiguous or physically linked to each other [95]. Moreover, a functional connectivity has also been introduced [95]. It accounts for how the interactions between structural characteristics may affect geomorphological, hydrological and ecological processes [95].

Both structural and functional connectivities have been qualitatively recognized in the relaxograms acquired by FFC NMR on water saturated soil samples [88]. In fact, the shape of the relaxograms was associated to the structural connectivity, while the position of the bands in a relaxogram was related to the functional connectivity [88]. The quantification of both structural and functional connectivity can be achieved by a mathematical empirical model developed by Conte and Ferro [57,89]. In particular, the non-exceeding empirical cumulative frequency curve, $F\left(T_1\right)$, obtained by the integration of any relaxogram (Figure 14), provides two different types of information.

Let us arbitrarily divide the empirical cumulative frequency curve in three parts: $F\left(T_1\right)$ < 0.01, 0.01 < $F\left(T_1\right)$ < 0.99, and $F\left(T_1\right)$ > 0.99. The two extremes represent the fast- ($F\left(T_1\right)$ < 0.01) and slow- ($F\left(T_1\right)$ > 0.99) relaxing nuclei, respectively. The former belong to water molecules constrained in residual pores, while the latter are from water molecules moving in transmission pores. All the ¹H nuclei in the interval 0.01 < $F\left(T_1\right)$ < 0.99 are considered to belong to water molecules moving in storage pores. The longitudinal relaxation time interval corresponding to $F\left(T_1\right)$ < 0.01 is represented in Figure 14 by the segment CA having length T_A, while the segment BD with length T_B corresponds to $F\left(T_1\right)$ > 0.99 (Figure 14). The T_B/T_A ratio can be considered as the relative amount of freely moving water over the water strongly bound to the wall boundaries of the smallest pores in soil. The larger the pore size, the higher the T_B/T_A ratio as a consequence of the huger amount of unconstrained water. Conversely, as pore size decreases, water molecules become more restrained and the T_B/T_A ratio decreases. According to this, it is possible to state that the T_B/T_A ratio can be considered as a measure of the functional connectivity because its value is related to the way water molecules interact with soil particles. Therefore, we can indicate the aforementioned ratio as the functional connectivity index (FCI).

In order to define the structural connectivity, the middle part of the empirical cumulative frequency distribution can be used. In fact, the 0.01 < $F\left(T_1\right)$ < 0.99 interval, which corresponds to water molecules moving in storage pores, represents the main structure of the soil, so that it may suitably represent the spatial pattern inside the soil.

According to Laudicina et al. [88], the width of the relaxogram is correlated to the total amount of pores in the system, whereas the position of the maxima refers to the contribution of residual, storage, and transmission pores. Therefore, it can be assumed that the $F\left(T_1\right)$ S-shaped curve weights the contribution of the storage pores over the other two types of pores. Hence, it is possible to define a structural connectivity index (SCI) as the ratio between the coefficient of variation of the relaxation times corresponding to the empirical cumulative frequency values in the range 0.01 to 0.99, that is, the slope Δ = GF/EF in Figure 14.

Noticeably, both FCI and SCI appear independent of the applied magnetic field [57]. Therefore, the measurement of the aforementioned indexes can be achieved by using the magnetic field strength which ensures the best NMR experimental sensitivity. Moreover, the FCI and SCI indexes can be also obtained by the evaluation of T₂ relaxograms. However, this is yet to be verified.

Having introduced the FCI and SCI indexes, let us now reveal why they can be a measure of soil erosion.

Erosion involves changes in soil structure and, hence, in soil pore distribution [96]. For this reason, water molecules change their motion rate according to the size of the pores. As a consequence, the FCI and SCI values also change, thereby allowing the monitoring of soil erosion (ongoing study) and permitting one to apply what is needed to restore soil productivity and functions.

7.3. The Behaviour of Dissolved Organic Matter (DOM)

Dissolved organic matter (DOM) plays a very important role in many biogeochemical mechanisms and responses to changes in ecological processes. As an example, it is involved in the solubilization and transport of both inorganic [97] and organic molecules and colloids [98,99], micronutrient availability [100], rock weathering [101], pedogenesis of topsoils and subsoils [102,103], soil water repellence [73], and soil texture [104].

Due to its environmental importance, DOM chemical–physical properties have been studied for a long time. As an example, many models have been suggested to explain the conformational behavior of dissolved organic matter, such as the linear macromolecular polyelectrolyte hypothesis [105], the supramolecular assemblies of molecules stabilized by weak interactions [106], the heterogeneous Donnan gel phases [107], or the mixture of supra and macromolecules [108].

The complexity and flexibility of DOM are considered to be responsible for its ability to enhance solubility of hydrophobic organic compounds, its capacity to decrease water surface tension, or its capability to trap and transport nutrients across the space [109].

Liquid water can be considered as a three-dimensional network of molecules held together by transient hydrogen bonds [110]. Once organic matter is dissolved, the H-bond network is altered, thereby leading to a new H-bond network which is stabilized by the presence of the dissolved organic matter.

When heating/cooling cycles occur (that is a rule in nature), DOM structure is altered [111] and changes in DOM sorption/desorption capacities can be observed [112]. The temperature dependence of the aforementioned changes were attributed to the dynamic re-arrangement of the water clusters around DOM [113]. This latter hypothesis was verified by NMR relaxometry with fast field-cycling setup [109]. It came out that a delay (also referred to as hysteresis) in re-establishing water network around DOM occurs as a result of the temperature fluctuations. Hysteresis appeared to be more pronounced when more hydrophilic dissolved organic matter was accounted for. The relaxometry results also supported the view that soil DOM consists of a hydrophobic rigid core surrounded by amphiphilic and polar molecules progressively assembled, and forming an elastic structure able to mediate the reactivity of the whole dissolved organic matter [109].

The relaxometry with the FFC setup was also helpful in monitoring the changes in the chemical nature of DOM upon interaction with montmorillonite (M) and kaolinite (K) [84]. In particular, the typical spectroscopic investigations (i.e., FT–IR and cross polarization (CP) magic angle spinning (MAS) ¹³C NMR spectroscopy) indicated that montmorillonite revealed a higher affinity for more hydrophilic –OH containing DOM moieties, whereas kaolinite adsorbed preferentially hydrophobic DOM components. The relaxometric results showed a huge simplification of the chemical nature of DOM after adsorption on montmorillonite as compared to the adsorption on kaolinite. In particular, the number of different DOM relaxation components after adsorption was in the order K > M (Figure 15). In order to explain the relaxometric behavior reported in Figure 15, it must be reminded that dissolved organic matter is a mixture of molecules that result from the degradation of humus, biomasses, plant materials, and root exudates, all of them involved in superstructures with a wide distribution of different molecular sizes [114]. After adsorption on montmorillonite and kaolinite, the residual DOM organic moieties re-aggregate to form new superstructures. Those obtained after interaction with montmorillonite are the simplest (Figure 15B) due to the inter-layer adsorption ability of montmorillonite which enabled the retention of a larger number of organic matter than kaolinite.

7.4. Water Behavior in the Presence of Inorganic Ions: The Dynamics of Soil Solution

Water is crucial to the movement of dissolved species in the environment. Indeed, it is involved in most biogeochemical processes and affects the properties of bio-systems from microscale to macroscale [44]. In order to retrieve information on the role played by water in nutrient dynamics towards plant roots, relaxometry with fast field-cycling setup has been applied on water solutions containing dissolved salts at different concentrations.

Dissolution of a solute in a solvent is indicated as solvation process which is referred to as hydration when the solvent is water. Hydration leads to formation of water shells around the solute. In particular, two different hydration shells are usually recognized [115,116]. The innermost shell, where H₂O directly interact with the solute, is made by immobilized molecules. A second hydration shell made by more disordered and mobile H₂O molecules surrounds the former layer. As solutions become progressively diluted, a third outermost water shell, referred to as bulk, can also be recognized [117]. The frequency of diffusional motion in bulk water is larger than in the first and second hydration shells. The size of each hydration shell depends on the physico-chemical features of the solutes [115,116,117,118]. In particular, the larger the charge density of the ions dissolved in water, the thinner is the hydration layer due to the strong effects of the electrical field generated by the ions [115,116,117,118]. In other words, water molecules are in the closest proximity of the ions. Upon increment of salt concentration, the space available for water in the bulk progressively reduces, thereby restricting the presence of free moving water. At the saturation point, no free moving bulk water is present.

Water molecules in the innermost shells can be exchanged with the ones in the outermost shell. The rate of exchange is affected by the strength of the water–solute interactions. In more detail, the stronger the interactions (as in the case of large solute concentrations or ions with high charge density), the lower the water mobility. Conversely, as the interaction strength weakens (low concentrated solutes, or low charge density ions), water mobility increases, and the exchange rate raises [110]. Water molecules strongly interacting with solutes are arranged in an ice-like shell where they are not flexible enough to accomplish the H-bonds within the first hydration shell. For this reason, strengthening of the H-bonds between the first and the second hydration shell occurs. By weakening the interactions between ions and water, water molecules in the first hydration shell become more flexible. As a consequence, the first shell intra-layer H-bonds strengthen, while the first–second shell inter-layer ones weaken [117]. Ions weaken the first shell intra-layer interactions, thereby allowing stronger first–second shell inter-layer interactions, which are referred to as structure makers or kosmotropes (Figure 16A). Water molecules around kosmotropic ions form a low-density ice-like layer [110,119,120,121,122,123]. Ions which are incapable of weakening the first shell intra-layer interactions, thus leading to weaker first–second shell inter-layer interactions, are indicated as structure breakers or chaotropes (Figure 16B). Water molecules generate high-density microdomains around chaotropic ions [110,119,120,121,122,123].

Fast field-cycling NMR relaxometry has been able to reveal the properties of water molecules in the hydration shells of several salts chosen for their natural presence in the soil solution. In particular, solutions containing different concentrations of KCl, NaCl, CaCl₂, K₂CO₃, NaNO₃ and NH₄NO₃ have been investigated.

Among the aforementioned salts, potassium chloride revealed a chaotropic nature, in the range of concentrations from very diluted up to saturation. The opposite behavior was shown in the same concentration range by NaCl, CaCl₂, and K₂CO₃, whereas NaNO₃ and NH₄NO₃ showed either a chaotropic or a kosmotropic nature depending on their concentration. More in detail, they behaved as chaotropes up to the concentration of 4 M. Conversely, they showed a kosmotropic behavior above the aforementioned concentration. From a qualitative standpoint, the different effects could be explained considering the prevalence of the chaotropic nature of the cations (either Na⁺ or NH₄⁺) when concentrations were < 4 M and by the prevalence of the kosmotropic nature of nitrate when salt concentration was > 4 M.

The model discussed in this sub-section has been confirmed in Payne et al. [41], Yadav et al. [124], Farashishiko et al. [125], and Sharma and Chandra [126].

7.5. A FFC NMR-Based Model for Nutrient Dynamics in Soils

According to the water behavior depicted above, an FFC NMR-based model for nutrient dynamics towards plant roots in soil solution can be suggested.

As already stated, structure maker solutes are surrounded by low-density water (LDW) microdomains, whereas structure breaker solutes generate high-density water (HDW) microdomains [110,119,120,121,122,123].

The LDW microdomains are due to the low flexibility of water molecules around kosmotropes. For this reason, water molecules are forced towards an ice-like behavior, with each molecule occupying a fixed position in order to fulfill the requirements necessary for the formation of interlayer H-bond. By increasing the distance from the solute, water flexibility increases as well. This means that the LDW domain progressively switches to a denser and more disordered domain, with water molecules laying closer to each other (Figure 17). The gradual switch from the LDW to the HDW domains is affected by the intensity of the local solute-generated electric field. For this reason, it can be assumed that the size of the LDW micro-domain increases with the intensity of the electric field. Noteworthy, water molecules in the LDW domain can exchange with those in the HDW domain, while the separation among the domains is maintained (Figure 17).

As water surrounds a chaotrope, its molecular organization can be conceivably considered more disordered because of the its larger flexibility in the first hydration shell. As a consequence, the number of water molecules that can be packed together in the same volume unit is larger than in the case of the structure maker solutes. Hence, a high-density water domain is generated. While a progressive transition LDW → HDW has been described for the structure maker solute, no such a transition can be hypothesized around the structure breaker solute. In other words, water molecules behave as a bulk also in the nearest proximity of the structure breaker solute [119,120,121,122,123]. When a mixture of both kosmotropes and chaotropes is accounted for, an exchange of water molecules among LDW and HDW domains occurs. However, since the motion of water molecules cannot be independent from that of the dissolved solutes, the latter randomly move together with water in order to let the solution reach a water average density homogeneity that may ensure a maximum entropy value, thereby leading also to a solute concentration homogenization throughout the solution. A complete discussion about the thermodynamics of the aforementioned water behavior can be found in Holtzer’s papers [127,128].

Let us now consider how water can be organized in the rhizosphere, which is ‘the zone that includes the soil influenced by the root along with the root tissues colonized by microorganisms’ [129]. Here, root exudation, respiration, redox reactions, and nutrient supply appear to be the main factors affecting pH changes [130,131]. In particular, the amount of H⁺ or OH⁻ in the rhizosphere strongly depend on both the ions crossing the plasma membrane towards the soil and the ions taken up by the plants from the soil. In fact, according to their nature, all the ions coming from the intake–outtake flux must be counterbalanced by either H⁺ or OH⁻ [131]. As an example, it has been well established that plants absorbing nitrogen in the form of nitrate tend to raise the pH in the rhizosphere [132,133]. Conversely, plants absorbing nitrogen in the form of ammonium or N₂ tend to lower the rhizosphere pH [132,133]. However, it must be also pointed out that the degree of the pH changes depends on the type of plant [130,134].

Having acquainted that the composition of rhizosphere solution is very rich in ions, it is possible to state that, in the rhizosphere, water molecules tend to assume the high-density liquid form [135]. In fact, computer simulation studies revealed that the ability to organize water molecules in ordered clusters does not depend only on the charge density of the ions, but also on the solution ionic strength [135]. In other words, an increase in the charge concentration generates kosmotropic conditions, thereby leading water molecules to behave like structured ice [135].

The low-density water domains in the rhizosphere may exchange with the high-density water domains outside the rhizosphere in order to drive water density towards homogeneity, thus ensuring entropy maximization [127,128]. As a consequence, movement of solutes (such as Na⁺, K⁺ and NH₄⁺, that are well-recognized plant nutrients) from/to plant roots can be achieved.

Let us concentrate, now, on the organization of water on soil surface. Due to the presence of hydrogen bond donors/acceptors as well as electron-poor moieties, water molecules hooked on soil surface form a thin ice-like film where LDW can be recognized [136]. The kosmotropic effect of the solid soil surface reduces with the distance, and so activates the 3D “jumps” described in Figure 13. The jump frequency depends on soil wettability which, in turn, is affected by chemical nature of natural organic matter (NOM) [84,109]. For this reason, the higher the wettability, the stronger the surface affinity of H₂O molecules, thereby lowering the number of the possible 3D jumps. At the same time, the 2D diffusion depicted in Figure 13 predominates. The latter is mediated by the formation/demolition of soil–water and water–water interaction dynamics. As reported above, water molecules do not move “alone”. They also transport solvated solutes. Therefore, the dynamics of the solutes from/to the solid soil surface can be explained by the increasing of the solvation effect (see sub-section below) as water moves from the LDW to the HDW domains [110,119,120,121,122,123,136].

Shrinking together the dynamics for the rhizosphere and the soil surface reported above, it is possible to summarize the mechanism for nutrient transport as reported in Figure 18.

Three different water density zones can be recognized (Figure 18). Low-density water layers occur on both soil and root surfaces, while a high-density water layer is identifiable in the soil solution that is the farthest from the two surfaces. Owing to the different densities, the LDW $\rightleftarrows$ HDW equilibrium occurs. This also allows the movement of the inorganic/organic systems adsorbed on both plant roots and soil surfaces. Following the aforementioned equilibrium, the nutrients reaching the rhizosphere solution can then be taken up by plant roots to satisfy plant metabolism.

7.6. The Mechanism of Soil Ionic Exchange Capacity: The FFC NMR Relaxometry Point of View

Although it has not explicated in the sub-sections above, the mechanism suggested for nutrient transport can also explain soil ionic exchange capacity. In fact, dissolution of a solute in water occurs as a consequence of hydration. This should be considered as a true reaction where the solute has to compete with H₂O molecules for the interactions with other H₂O molecules. The weaker the H-bonds between water molecules, the easier it is for the solute to build hydration shells and, hence, to dissolve. However, water domains in which H-bonds are weak are also denser than those where H-bonds are strong. Therefore, it is conceivable that solutes preferentially react with H₂O molecules in the HDW domains [110,119,120,121,122,123]. Due to this, it is conceivable that solutes adsorbed on the solid soil surfaces and included in a surficial LDW domain tend to move away from the surface in order to react with the water molecules belonging to the farthest HDW domains. In addition, the solutes (e.g., cations and anions) removed from the surface are replaced by other solute molecules moving from the HDW domain.

8. Conclusions and Perspectives

We are well aware that fast field-cycling NMR relaxometry is a technique that remains accessible to a small number of specialists. A knowledge of deep quantum mechanics is needed together with the ability to handle electronics, circuits, and programming languages. FFC NMR results are not as easily readable as those from spectroscopy, thereby making the former less popular in soil science than the latter, which is very often misused in the field [3,11,12,137].

Within the present review, we have summarized the main and most recent soil science findings retrieved by FFC NMR relaxometry. The focus point is surely the application of the right model for data interpretation. Although the plethora of models described here gives the idea of an exit-less labyrinth, we have clearly shown that the combination between FFC NMR results and the data from other soil science investigations (e.g., soil texture, pH, soil organic matter nature and content, electron spin resonance results, microscopy etc.) strongly facilitates the choice of the right mathematical model to be applied.

FFC NMR results revealed a paramount importance to unveil molecular mechanisms occurring in soil solution which, in turn, can be used to reinterpret the phenomena arising on the surface of plant roots. In particular, the main results showed two main water movements on the solid surfaces: a 2D diffusion and a 3D out-of-surface motion. While the former is responsible for the distribution of solutes and water on the surfaces, the latter is involved in the transport phenomena toward plant roots and nutrient leaching. Transport phenomena can be efficaciously explained by considering the chaotropic or kosmotropic nature of solutes which allow the identification of an HDW and a LDW domain. The LDW domain is generated by ice-like water which is due to (i) water interacting with kosmotropic solutes, (ii) water in kosmotropic conditions (i.e., large solute concentrations), or (iii) water interacting with solid surfaces. The HDW domains are generated by both chaotropic solutes and chaotropic conditions (i.e., low solute concentrations). Solute movement is affected by water solvation efficiency, which, in turn, depends on the LDW and the HDW domains. The former reduce water solvation efficiency as compared to HDW domains. Consequently, solute transport in soils is affected not only by random water movements but also by the different solubility of the solute in the different water density domains.

The model of water dynamics suggested above can be a subject of concern, as it is based on the combination of the interpretation of FFC NMR relaxometry data and those suggested by Wiggins [119,120,121,122,123], Gallo et al. [135], and Henderson [136]. However, there are few papers where a negligible effect of solutes on the 3D H-bond network of water molecules is reported [138,139]. Therefore, further research is needed to better clarify the role of low-density and high-density domains in soil and rhizosphere water and regarding nutrient mobility in soil systems.