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

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


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 0 )-are perturbated by an electromagnetic radiation pulse (B 1 ) applied for a short time interval along a direction perpendicular to B 0 . 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 0 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 0 . Both processes contribute to the overall regaining of the magnetization equilibrium. The time constant associated to the longitudinal relaxation is indicated as T 1 (also referred to as the longitudinal 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.

The Meaning of the T1 Value
T1 (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., 14 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 (R1)-which is the inverse of T1-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 nonspherically symmetric charge distribution, which is represented by a non-zero quadrupolar momentum (μQ). The latter affects the effective strength of the applied magnetic field (B0) felt by the investigated nucleus (that is 1 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 Td and Oh), 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 1 H-1 H interactions increases. Hence, reduction in the 1 H T1 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 1 H-1 H dipolar interactions and 1 H T1 reduction. Conversely, as molecular mobility increases because of pore size enlargement (e.g., in silty and sandy systems), the 1 H-1 H dipolar interactions weaken and longer 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.

The Meaning of the T 1 Value
T 1 (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 > 1 2 (e.g., 14 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 1 )-which is the inverse of T 1 -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) > 1 2 . 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 0 ) felt by the investigated nucleus (that is 1 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 1 2 [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 1 H-1 H interactions increases. Hence, reduction in the 1 H T 1 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 1 H-1 H dipolar interactions and 1 H T 1 reduction. Conversely, as molecular mobility increases because of pore size enlargement (e.g., in silty and sandy systems), the 1 H-1 H dipolar interactions weaken and longer longitudinal relaxation times are measured [44]. In other words, shorter T 1 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 1 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 1 values closer to the shortest T 1 limit. On the other hand, water molecules that move inside the largest pores produce T 1 values closer to the longest T 1 limit. All T 1 values between the two limits are due to water molecules moving inside pores having sizes comprised between the two extremes [49]. 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 T1 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 T1 values closer to the shortest T1 limit. On the other hand, water molecules that move inside the largest pores produce T1 values closer to the longest T1 limit. All T1 values between the two limits are due to water molecules moving inside pores having sizes comprised between the two extremes [49]. 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.

The Inversion Recovery Sequence and the Basic FFC NMR Experiment
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 B0), 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 is the signal intensity at the end of the evolution period with duration τ, and is the signal intensity when τ = 0 [50].
The use of any fitting software where -vs-τ is reported, allows the achievement of both and T1. 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 0 ), 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(τ) is the signal intensity at the end of the evolution period with duration τ, and I 0 is the signal intensity when τ = 0 [50].
The use of any fitting software where I(τ)-vs-τ is reported, allows the achievement of both I 0 and T 1 . 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 1 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 1 , 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 1 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 ω RLX = ω POL 2 , where ω RLX and ω 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]. 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 Agronomy 2020, 10, 1040 6 of 33 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).

How to Obtain the T 1 Value from the FFC NMR Experiment
of the experiment).

= +
(2) In both Equations (2) and (3), is the signal intensity at the selected τ value, a is the offset and bi is the magnetization intensity at the Boltzmann equilibrium of the i-th component of the molecular motion at each fixed BRLX 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 hypercorrected 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 nonpolarized sequence, respectively, can be applied [18].  In both Equations (2) and (3), I(τ) 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].
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 Agronomy 2020, 10, 1040 7 of 33 single model, thereby avoiding arbitrary assumptions about the number of exponentials to be applied in modelling the relaxometric results.

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 1 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 1 values reported in Table 1 are obtained. Table 1. Longitudinal relaxation times for the quartz sand sample added with Milli-Q grade water as described in the caption of Figure 3. ω L is the proton Larmor frequency of the applied B RLX . Here, T f ast 1 may refer to the relaxation time of water molecules in the closest proximity of the solid surface, while T slow 1 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 f ast 1 is due to water molecules moving in pores smaller than those where water molecules generating T slow 1 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).
Agronomy 2020, 10, x FOR PEER REVIEW 7 of 33 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.

From Time Domain to Time Domain
The mathematical models given in Equations (2) and (3) can be considered as a discrete representation of the T1 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 T1 values reported in Table 1 are obtained. Table 1. Longitudinal relaxation times for the quartz sand sample added with Milli-Q grade water as described in the caption of Figure 3.
is the proton Larmor frequency of the applied BRLX. Here, may refer to the relaxation time of water molecules in the closest proximity of the solid surface, while 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 is due to water molecules moving in pores smaller than those where water molecules generating 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 2 > 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 . Different dynamics of water molecules in the quadrupolarless-nuclei quartz sand used as a porous system model (white circles). The biexponential assumption used to apply Equations (2) and (3) is valid only for A, and B. In case C, a three-exponential model should be applied.
Due to the very good fitting (R 2 > 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 Agronomy 2020, 10, 1040 8 of 33 provide only one average T 1 value, and a k 1. It is not yet clear what the relationship between the heterogeneity parameter and the number of the different T 1 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 1 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].
In Equations (6) and (7), T min 1 and T max 1 are the suitable limits within which all the T 1 values range; D(T 1 ) 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 1 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(τ) to D(T 1 ). 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 1 distributions, also referred to as relaxograms ( Figure 5), regardless of the procedure used to obtain the most probable distribution of relaxation times. 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 T1 value, and a k ≠ 1. It is not yet clear what the relationship between the heterogeneity parameter and the number of the different T1 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 T1 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 nonpolarized experiments are carried out [18].
In Equations (6) and (7), and are the suitable limits within which all the T1 values range; 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 T1 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 to . 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 T1 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 T1 components are distinguishable. They are centered at around 51, 126, and 292 ms in the relaxogram acquired at the 1 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  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 1 components are distinguishable. They are centered at around 51, 126, and 292 ms in the relaxogram acquired at the 1 H Larmor frequency of 0.1 MHz, Agronomy 2020, 10, 1040 9 of 33 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].
Agronomy 2020, 10, x FOR PEER REVIEW 9 of 33 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].

The Nuclear Magnetic Resonance Dispersion (NMRD) Profile and Its Modeling
The T1 values obtained via Equations (2) to (5) can be plotted versus BRLX 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].

The Nuclear Magnetic Resonance Dispersion (NMRD) Profile and Its Modeling
The T 1 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].
Agronomy 2020, 10, x FOR PEER REVIEW 9 of 33 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].

The Nuclear Magnetic Resonance Dispersion (NMRD) Profile and Its Modeling
The T1 values obtained via Equations (2) to (5) can be plotted versus BRLX 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].  In Equation (8), R 1 is the longitudinal relaxation rate (R 1 = 1/T 1 ) which represents the distribution of the motion frequencies in a molecular system; ω L is the B RLX intensity given as proton Larmor frequency; and τ 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 τ 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 1 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].
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.
Agronomy 2020, 10, x FOR PEER REVIEW 10 of 33 In Equation (8), R1 is the longitudinal relaxation rate (R1 = 1/T1) which represents the distribution of the motion frequencies in a molecular system; is the BRLX intensity given as proton Larmor frequency; and 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 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 T1 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].
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 reorientational 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. Owing to the complexity of the dynamics of H2O molecules on the porous surface, the shape of the dispersion is not sigmoidal, such as in Figure 7, but stretched.
In order to account for the NMRD stretching, Halle and coworkers [58] developed the free model analysis given in Equation (10): In order to account for the NMRD stretching, Halle and coworkers [58] developed the free model analysis given in Equation (10): where the i index refers to the different components of the motion, τ i C 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 χ 2 (N) value. The latter must be compared with the χ 2 (N + 1) value obtained by applying N + 1 Lorentzians. If the fit improves, the ratio F (N, N+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 C i , τ i C 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): 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].
All the parameters in Equation (12) have been already defined. The set of C i , τ i C 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)).
where ω I is the intensity of B RLX expressed as 1 H Larmor frequency; ω S = 659ω I is the Larmor frequency of the unpaired electron in the paramagnetic species; R 0 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 2 O molecules near the paramagnetic relaxation source (Fe 3+ , for instance). K is provided by the relationship: where δ water = 0.3 nm is the width of a monomolecular water layer; σ 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 2 ), ρ water = 1 g cm −3 is water density; S P,NMR (given in m 2 g −1 ) is the NMR-specific surface area calculated as S P × 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); γ I and γ S = 659γ 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. τ m is the surface diffusion correlation time, while τ s is indicated as surface residence time. The former represents the "hopping" time of water among surficial binding sites, while the latter (τ s ) is the time of residence of water on the surface of the porous medium [18]. The τ s τ m ratio is referred to as NMR surface wettability [60,61].

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.

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.
Agronomy 2020, 10, x FOR PEER REVIEW 12 of 33 specific surface area calculated as × (where is the surface area achieved by the Brunauer-Emmet-Teller (BET) analysis, while < 1 is the solid-to-liquid ratio); and = 659 are the proton and electron gyromagnetic ratios, respectively (for the most abundant paramagnetic species). Finally, is the spin quantum number of the paramagnetic species. is the surface diffusion correlation time, while is indicated as surface residence time. The former represents the "hopping" time of water among surficial binding sites, while the latter ( ) is the time of residence of water on the surface of the porous medium [18]. The ratio is referred to as NMR surface wettability [60,61].

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.

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 BRLX intensities. For this reason, the three NMRD dispersion curves reported in Figure 9 have been obtained.  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 1 (or R 1 ) values reported in Section 2, one may expect that the correlation times obtainable by the models listed in Section 5, must follow the order: τ f ast 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, α f ast > α intermediate > α slow .
Agronomy 2020, 10, x FOR PEER REVIEW 13 of 33 component) ( Figure 6). According to this qualitative evaluation and to the meaning of the T1 (or R1) values reported in Section 2, one may expect that the correlation times obtainable by the models listed in Section 5, must follow the order: > > . 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, > > . Figure 10. Values of α, β, and for the three components of the water motion in quartz sand. These parameters were obtained by fitting the NMRD profiles in Figure 9 by the BPP model given in Equation (9).
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 T1 and the faster results the R1. For this reason, as expected by the relative position of the three NMRD profiles, the β values follow the order > > .
Finally, the correlation times are in the order: > > , which confirms the expectations reported above.

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 Figure 10. Values of α, β, and τ C for the three components of the water motion in quartz sand. These parameters were obtained by fitting the NMRD profiles in Figure 9 by the BPP model given in Equation (9).
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 1 and the faster results the R 1 . For this reason, as expected by the relative position of the three NMRD profiles, the β values follow the order β f ast > β intermediate > β slow .
Finally, the correlation times are in the order: τ f ast C > τ intermediate C > τ slow C , which confirms the expectations reported above.

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 , τ 1 , τ 2 , τ 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).
Agronomy 2020, 10, x FOR PEER REVIEW 14 of 33 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). Figure 11. Correlation times obtained by applying a six-parameter free model analysis (Equation (10) and (11)) on the NMRD profiles for soils under A. Eucalyptus camaldulensis Dehnh., B. E. occidentalis Endl., C. Pinus halepensis Mill., and D. Cupressus sempervirens. The longest value is due to water moving in clay-type pores. The shortest value is related to water moving in sand-type pores. The intermediate value is attributable to water moving in loam-type pores. The histogram has been built by using the data from Reference [62].
The free model analysis was also applied in Conte et al. [63] to unveil the nature of the waterbiochar 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 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 electronpoor 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.  (10) and (11)) on the NMRD profiles for soils under A. Eucalyptus camaldulensis Dehnh., B. E. occidentalis Endl., C. Pinus halepensis Mill., and D. Cupressus sempervirens. The longest τ C value is due to water moving in clay-type pores. The shortest τ C value is related to water moving in sand-type pores. The intermediate τ C value is attributable to water moving in loam-type pores. The histogram has been built by using the data from Reference [62].
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 τ 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. 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].

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. Contact angle-vs-NMR wettability as obtained from the data reported in Reference [60]. The full line is only for the eyes. CA stands for contact angle. The shorter the contact angle (e.g., CA < 90°), the higher is the surface hydrophilicity. Conversely, as the contact angle increases (e.g., CA > 90°), the more hydrophobic (or the less hydrophilic) the solid surface.

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.

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. [ [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 . Contact angle-vs-NMR wettability as obtained from the data reported in Reference [60]. The full line is only for the eyes. CA stands for contact angle. The shorter the contact angle (e.g., CA < 90 • ), the higher is the surface hydrophilicity. Conversely, as the contact angle increases (e.g., CA > 90 • ), the more hydrophobic (or the less hydrophilic) the solid surface.
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].

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.

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 1 H-1 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: 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.
Agronomy 2020, 10, x FOR PEER REVIEW 17 of 33 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 1 H-1 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: > > which is the inverse of the motion rate, i.e., < < . 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. Figure 13. Bimodal water movement into soils. Surface water molecules are hooked to the soil surface being subjected to horizontal diffusion. When surface water "jumps" to the bulk, one of the bulk molecules replace it. The frequency of the "jumps" and the duration of the horizontal diffusion are affected by the strength of the interactions between water molecules and soil surface.

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 Figure 13. Bimodal water movement into soils. Surface water molecules are hooked to the soil surface being subjected to horizontal diffusion. When surface water "jumps" to the bulk, one of the bulk molecules replace it. The frequency of the "jumps" and the duration of the horizontal diffusion are affected by the strength of the interactions between water molecules and soil surface.
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.

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 137 Cs has been developed and applied worldwide [87]. The attraction for 137 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, 137 Cs is now present on a global scale. Therefore, monitoring its fate provides useful information on soil transformations [87]. However, 137 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(T 1 ), obtained by the integration of any relaxogram (Figure 14), provides two different types of information.
Agronomy 2020, 10, x FOR PEER REVIEW 19 of 33 soil particles. Therefore, we can indicate the aforementioned ratio as the functional connectivity index (FCI).

Figure 14.
The black curve is the relaxogram acquired at the proton Larmor frequency of 30 MHz for one of the soils described in Reference [89]. This relaxogram has never been published, and is only discussed in the aforementioned reference. The blue curve is the empirical cumulative frequency obtained by the integration of the relaxogram. The segments CA, having length TA (expressed in ms), and BD, with the length TB (in ms), represent the amount of fast-and slow-relaxing nuclei, respectively, corresponding both to 1% of the empirical cumulative frequency curve. The ratio TB/TA is the functional connectivity index (FCI). The slope Δ = GF/EF is the structural connectivity index (SCI).
In order to define the structural connectivity, the middle part of the empirical cumulative frequency distribution can be used. In fact, the 0.01 < < 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 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 T2 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. Figure 14. The black curve is the relaxogram acquired at the proton Larmor frequency of 30 MHz for one of the soils described in Reference [89]. This relaxogram has never been published, and is only discussed in the aforementioned reference. The blue curve is the empirical cumulative frequency obtained by the integration of the relaxogram. The segments CA, having length T A (expressed in ms), and BD, with the length T B (in ms), represent the amount of fast-and slow-relaxing nuclei, respectively, corresponding both to 1% of the empirical cumulative frequency curve. The ratio T B /T A is the functional connectivity index (FCI). The slope ∆ = GF/EF is the structural connectivity index (SCI).
Let us arbitrarily divide the empirical cumulative frequency curve in three parts: F(T 1 ) < 0.01, 0.01 < F(T 1 ) < 0.99, and F(T 1 ) > 0.99. The two extremes represent the fast-(F(T 1 ) < 0.01) and slow-(F(T 1 ) > 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 1 H nuclei in the interval 0.01 < F(T 1 ) < 0.99 are considered to belong to water molecules moving in storage pores. The longitudinal relaxation time interval corresponding to F(T 1 ) < 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(T 1 ) > 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(T 1 ) < 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(T 1 ) 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 2 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.

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) 13 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.

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 H2O directly interact with the solute, is made by immobilized molecules. A second hydration shell made by more disordered and mobile H2O 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,

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 2 O directly interact with the solute, is made by immobilized molecules. A second hydration shell made by more disordered and mobile H 2 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, CaCl2, K2CO3, NaNO3 and NH4NO3 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, CaCl2, and K2CO3, whereas NaNO3 and NH4NO3 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 NH4 + ) when concentrations were < 4 M and by the prevalence of the kosmotropic nature of nitrate when salt concentration was > 4 M.

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.
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 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 2 , K 2 CO 3 , NaNO 3 and NH 4 NO 3 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 2 , and K 2 CO 3 , whereas NaNO 3 and NH 4 NO 3 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 4 + ) 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].

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.
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).
Agronomy 2020, 10, x FOR PEER REVIEW 23 of 33 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]. Figure 17. Kosmotropic ion dissolved in water. This is surrounded by a low-density water (LDW) microdomain which moves towards a high-density water (HDW) microdomain as the distance from the ion increases. The green arrows indicate the molecular exchanges between HDW and LDW domains.
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 intakeouttake 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 N2 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]. Figure 17. Kosmotropic ion dissolved in water. This is surrounded by a low-density water (LDW) microdomain which moves towards a high-density water (HDW) microdomain as the distance from the ion increases. The green arrows indicate the molecular exchanges between HDW and LDW domains.
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 2 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 4 + , 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 2 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. 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 NH4 + , 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 H2O 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 ⇄ 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.  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 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. 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 2 O molecules for the interactions with other H 2 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 2 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.

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.
Author Contributions: Both authors equally contributed to the paper. Both authors have read and agreed to the published version of the manuscript.
Funding: This study received no external funding.

Acknowledgments:
The authors acknowledge the help of Gianni Ferrante, Rebecca Steele, and Moreno Pasin from Stelar Srl for their help in providing the FFC NMR relaxometry analyses of the quartz sand described in the present review paper. Anna Micalizzi from University of Palermo is acknowledged for her technical assistance in preparing samples and performing experiments on some of the soil samples reported here.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In Section 3, Section 4, and Section 6, experimental data concerning a quartz sand saturated with Milli-Q grade water have been discussed. These data have never been published. Therefore, here the details of the experiments performed on this model system are reported. The quadrupolar-less nuclei quartz sand (CAS N. 14808-60-7, 50-70 mesh particle size) was purchased from Sigma (Milan, Italy) and used without any pre-treatment. The Milli-Q grade water (resistivity of 18.2 MΩ cm) was produced by using a Millipore system (Millipore Corporation, Billerica, MA, USA). The sample for the FFC NMR relaxometry experiment was prepared as a slurry according to Reference [56]. 1 H NMRD profiles (i.e., longitudinal relaxation rates R 1 or 1/T 1 vs proton Larmor frequencies) were acquired on a Stelar Spinmaster FFC 2000 relaxometer (Stelar s.r.l., Mede, PV, Italy) at a constant temperature of 25 • C. The field-switching time was 3 ms, while the spectrometer dead time was 15 µs. The proton spins were polarized at a polarization field (B POL ) corresponding to a proton Larmor frequency (ω L ) of 29 MHz for a period of polarization (T POL ) of 0.2 s. A recycle delay of 2 s was always applied. The longitudinal magnetization evolution was recorded at values of a relaxation magnetic field (B RLX ) corresponding to ω L comprised in the range of 0.015 to 35 MHz. The NMR signal was acquired with four scans for a period of time (τ) arrayed with 32 values, chosen in an exponential progression for the covering of the entire relaxation curve of interest. Finally, a 1 H 90 • pulse was used at the start of the acquisition period contemporarily to an acquisition magnetic field (B ACQ ) corresponding to a ω L of 16.2 MHz. The observable magnetization was revealed as free induction decay (FID) with a time domain of 100 µs sampled with 1000 points.

Appendix B
The first explanation of the relationship between porosity of a porous medium and longitudinal relaxation was given by Brownstein and Tarr [140]. In particular, for wet micro-or nano-porous materials the relaxation rates R 1 can be related to the fraction of water molecules interacting with the pore surface ( f s ) by the equation [38]: Here, R w is the relaxation rate of the bulk water. It is a constant (i.e., ca. 0.4 s −1 ) in the whole proton Larmor frequency range (0.01-40 MHz) usually applied when the typical fast field-cycling NMR relaxometry equipment is used [20]. R s is the intrinsic relaxation rate for the water moving on the porous surface.
R 1 can lead back to the porosimetric parameters via the following equations [38]: In Equations (A2) to (A4), ρ 1 is the longitudinal surface relaxivity (in cm s −1 ). It characterizes the ability of the porous surface to facilitate the longitudinal relaxation in the nearby fluid molecules. Its value is affected by surficial quadrupolar systems. S = mA s where m is the sample mass (in g) and A s is the specific surface area (in cm 2 g −1 ). V is the total pore volume of the saturated porous system, α is a pore shape parameter, D is the average pore diameter, and λ is the thickness of the water layer on pore walls (usually set as large as 0.3 nm). By combining Equations (A1) to (A4), it can be observed that D = λ α f s and f s = λ S V . The inverse relationship between R 1 and pore size (D) given in Equation (A4) explains what is reported in the discussion above.