Three-dimensional magnetic resonance imaging for groundwater

The surface nuclear magnetic resonance method (SNMR) is an established geophysical tool routinely used for investigating one-dimensional (1D) and sometimes 2D subsurface water-saturated formations. We have expanded the tool by developing a 3D application. 3D-SNMR is a large-scale method that allows magnetic resonance imaging of groundwater down to about 80 m. Similar to most surface geophysical methods, 3D-SNMR has limited resolution, but it is effective for investigating water-saturated geological formations larger than several tens of meters. Because the performance of the method depends on variable survey conditions, we cannot estimate it in general. For demonstration purposes, we present an example of numerical modeling under fixed conditions. Results show that under certain conditions it is possible to detect a water volume as small as 500 m3 and the detection threshold depends on the ambient electromagnetic noise magnitude and on the location of the target volume relative to the SNMR loops. The 3D-SNMR method was used to investigate accumulated water within the Tête Rousse glacier (French Alps). Inversion of the field measurements made it possible to locate the principal reservoir in the central part of the glacier and estimate the volume of accumulated water. These results were verified by 20 boreholes installed after the 3D-SNMR results were obtained and by pumping water out of the glacier. Very good correspondence between the 3D-SNMR and borehole results was observed.

Abstract. The surface nuclear magnetic resonance method (SNMR) is an established geophysical tool routinely used for investigating one-dimensional (1D) and sometimes 2D subsurface water-saturated formations. We have expanded the tool by developing a 3D application. 3D-SNMR is a largescale method that allows magnetic resonance imaging of groundwater down to about 80 m. Similar to most surface geophysical methods, 3D-SNMR has limited resolution, but it is effective for investigating water-saturated geological formations larger than several tens of meters. Because the performance of the method depends on variable survey conditions, we cannot estimate it in general. For demonstration purposes, we present an example of numerical modeling under fixed conditions. Results show that under certain conditions it is possible to detect a water volume as small as 500 m 3 and the detection threshold depends on the ambient electromagnetic noise magnitude and on the location of the target volume relative to the SNMR loops. The 3D-SNMR method was used to investigate accumulated water within the Tête Rousse glacier (French Alps). Inversion of the field measurements made it possible to locate the principal reservoir in the central part of the glacier and estimate the volume of accumulated 6 Author to whom any correspondence should be addressed.

Introduction
Groundwater plays an important role in modern society. Knowledge of the location and dynamics of groundwater is necessary for water supply purposes and pollution control. Under certain conditions, groundwater may present a natural hazard and must be monitored. We live in a period of global climate change that affects groundwater processes worldwide. General increases in temperature and changes in the amount and distribution of rain and snow cause changes in groundwater circulation. As a result, there is an increasing need for rapid, and in some cases, non-invasive methods for groundwater investigation.
Surface geophysical methods are commonly used for groundwater exploration. However, most geophysical methods are sensitive to the physical properties of the subsurface (for example, electrical resistivity) but not to groundwater itself, so these methods investigate groundwater only indirectly through the changes that the groundwater causes to the physical properties of the rocks.
The surface nuclear magnetic resonance (SNMR) method, also called magnetic resonance sounding (MRS), is a geophysical technique specifically developed for hydrogeological investigations [1]- [3]. The main advantage of SNMR over other geophysical methods is its selective sensitivity to groundwater. In most applications, SNMR is used in a sounding mode (one-dimensional (1D)-SNMR, also called MRS). In this mode, investigated aquifers are assumed to be horizontally stratified (1D formations). It has been shown that even 1D-SNMR can be used to detect bulk water in such complicated geological structures as karst caverns [4]. However, the lateral resolution of 1D-SNMR is limited and may be insufficient for identifying relatively small water-saturated formations for a drilling program. For example, to investigate and address problems related to natural hazards in a karst environment, the location of karst caverns must be known.
Reported applications of 2D-SNMR to the investigation of water-filled karst [5]- [7] show that the resolution of this method may be a significant improvement over 1D-SNMR. The 3 resolution of shallow 2D targets may be further improved by investing more time and labor in the 2D-SNMR tomography approach [8,9]. However, the 2D investigation still does not enable us to resolve 3D subsurface formations.
In this paper, we present a newly developed 3D implementation of the SNMR method. We conducted a numerical study of the sensitivity and resolution available with 3D-SNMR and we present an example of a practical application of the method to the investigation of water accumulated in a glacier. We show that using 3D-SNMR makes it possible to correctly locate water in ice caverns.

The three-dimensional surface nuclear magnetic resonance (3D-SNMR) method
All of the individual SNMR soundings are incorporated into the 3D-SNMR data set. Each SNMR field setup consists of a wire loop laid out on the ground, usually in a square with a side that ranges between 20 and 100 m [10]. The loop is then energized by a pulse of alternating current i(t) = I 0 cos(ω 0 t). The frequency of the current is set equal to the Larmor frequency of the protons ω 0 in the geomagnetic field B 0 (ω 0 = γ B 0 with γ being the gyromagnetic ratio). The pulse causes precession of the spin magnetization of the protons in groundwater around the geomagnetic field, which creates an alternating magnetic field that can be detected by the same loop after the pulse is terminated (the free induction decay method). Oscillating at the Larmor frequency, the SNMR signal has an exponential envelope that depends on the pulse moment q = I 0 τ , with I 0 and τ being the amplitude and duration of the pulse, respectively. Measurements of the magnetic resonance signal are carried out by varying the pulse moment. Distribution of the water content in the subsurface can be derived from the inversion of the SNMR signal. The depth of investigation depends on the loop size and site-specific conditions; it usually ranges from 40 to 120 m [11].
The amplitude of the magnetic resonance signal e 0 can be computed using existing models [10,12,13]. The signal induced in the coincident transmitting/receiving (Tx/Rx) loop is proportional to the sum of the flux of all precessing magnetic moments and thus where M ⊥ is the transverse component of spin magnetization, B 1 is transmitted by the surface loop magnetic field component perpendicular to the geomagnetic field, I 0 is the amplitude of the current in the loop, q is the pulse moment, ω 0 is the Larmor frequency for protons in the geomagnetic field, 0 w(r) 1 is the water content and r = r (x, y, z) is the coordinate vector. The 3D implementation uses overlapping Tx/Rx loops. The 3D imaging consists of measuring the SNMR signal independently in each loop while varying the pulse moment. All individual soundings are incorporated into a data set for 3D inversion. For interpretation, the field measurement equation (1) can be approximated by a matrix equation [14], where A = [ã i, j ] is a rectangular matrix of I × J (i = 1, 2, . . . , I ; j = 1, 2, . . . , J ; I = K × L; J = J x × J y × J z ; K = L l=1 K l ; k l = 1, 2, . . . , K l ; l = 1, 2, . . . , L; j x = 1, 2, . . . , J x , j y = 1, 2, . . . , J y and j z = 1, 2, . . . , J z ), K l is the number of pulse moments for a loop l, L is the number of measuring loops, J x , J y andJ z are the number of cells in the x-, yand z-directions, respectively.
In equation (2), the elements of matrix A represent the amplitude of the magnetic resonance signal generated by water in corresponding cells (assuming a 100% water content in each cell) and computed using equation (1). The set of experimental data is e = (ẽ 1 ,ẽ 2 , . . . ,ẽ i , . . . ,ẽ I ) T , the water content in the corresponding cell is w = (w 1 , w 2 , . . . , w j , . . . , w J ) T and the symbol T denotes the transposition. For simplicity, we assume that the cell size is constant throughout the investigated volume.
Inversion of the 3D-SNMR data is ill-posed. Different methods suitable for the resolution of ill-posed inverse problems can be found in the literature [15]- [19]. For our study, the inversion was conducted according to Tikhonov's regularization method [19]. To find an approximate solution to the matrix equation (2), this method supposes minimization of Tikhonov's functional, where e ε is the vector of the experimental data contaminated by the noise ε = (ε 1 , ε 2 , . . . , ε i , . . . , ε I ) T , w η is the solution vector that minimizes Tikhonov's functional (3), and η > 0 is the smoothing factor. In equation (3), we assume that smoothness is equal in the X-, Yand Z-directions.
To solve this minimization problem, we followed the discrepancy principle introduced by Morozov [20], which assumes that for erroneous data, all solutions that have a residual Aw η −e ε L 2 smaller than the experimental error are acceptable and the best solution is not necessarily the solutions that have the smallest residual. The best solution can be found only by using additional information about the solutions. Tikhonov's regularization method assumes that in addition to having the minimum residual, the solution has to be smooth. Hence, for a given noise ε > 0, we need to find a solution with a residual Aw η −e ε L 2 ε and stabilize it by making (( ∂ ∂ x w η ) 2 + ( ∂ ∂ y w η ) 2 + ( ∂ ∂z w η ) 2 ) small. w η is an approximation of the solution of the matrix equation (2). When ε → 0, η(ε) → 0 and w η → w. For the optimization itself, we used the conjugate gradient method [21].
For estimating the discrepancy, we use the root mean square error (RMSE), where em i and et i are the measured and theoretical amplitudes, respectively, I = L l=1 K l , K l is the number of pulse moments for each loop l and L is the number of loops.
Before measuring the magnetic resonance signal, the magnitude of ambient electromagnetic noise (always non-negative) is measured in the frequency band that is ±100 Hz from the Larmor frequency. The ambient noise ε is characterized by its mean magnitude, where I is the number of readings. Noise has a random phase relative to the signal and hence may increase or diminish the estimate of the signal amplitude. For numerical modeling, we generate synthetic noise with a constant magnitude but random phase and add it to the signal.

Numerical results
It is known that the performance of the SNMR method is site dependent [11] and hence its performance should be studied in conjunction with the site-specific conditions. In this paper, we consider conditions similar to those in our glacier study, which we present in the following section. We studied sensitivity and resolution numerically, using synthetic signals computed for different 3D target volumes located in dry rock (water content 0%). For modeling, we used a square loop with an 80 m side, a Larmor frequency of 2000 Hz, an inclination of the geomagnetic field of 62 • N, a subsurface of 100 .m and a maximum pulse moment of 12 000 A ms.

Sensitivity
It has been reported that the sensitivity of the method varies within the volume investigated with one SNMR loop and that it also depends on the inclination of the geomagnetic field [5,7,13]. Thus, the same volume of water may produce different signals, depending on the water location relative to the loop axis. Figure 1 shows an example of the amplitude of the SNMR signal as a function of target-volume position computed under conditions used for the modeling. As the amplitude of the signal varies within the loop area, the threshold of 3D target-volume detection should depend on the target size and its position relative to the loop. For demonstration purposes, we calculated the maximum amplitude of the magnetic resonance signal generated by a 20 × 20 × 20 m 3 target volume, assuming different depths for the target. We investigated three cases: a target located in the high-sensitivity area corresponding to the maximum signal; a target located in the low-sensitivity area corresponding to the minimum signal; and an intermediate case ( figure 1). The maximum detection depth was estimated by assuming a 10 nV signal detection threshold, which is the detection threshold of the particular SNMR instrumentation (NUMIS) used in this study. Numerical results presented in figure 2 show that the detection threshold may vary significantly. For example, at a depth of 5 m, the minimum volume that SNMR can detect when the target volume is located in the high-sensitivity area is about 560 m 3 . However, when the target is located in the low-sensitivity area, the volume should be four times larger (2300 m 3 ). In practice, a 3D-SNMR setup consists of overlapping loops and the lowsensitivity area of each loop is covered by the higher-sensitivity area of one of the neighboring loops. Thus, the threshold of detection fluctuates between the solid lines shown in figure 2. Electromagnetic noise is one of the major limitations for the method. In most field measurements, noise controls the threshold of signal detection. Figure 3 shows the impact of noise on the sensitivity of SNMR. We calculated the minimum detectable volume assuming different values of the threshold. For this modeling, we set the target volume in the high-sensitivity area of the loop.

Resolution
We performed an inversion of the synthetic signals computed for the one-box and twobox models assuming two experimental setups consisting of 80 × 80 m 2 measuring loops (figure 4). Two basic configurations were investigated: (a) nine half-side overlapped loops;  (b) 12 quarter-side overlapped loops. We use two models: a one-box target volume of 60 × 50 × 20 m 3 and a two-box target volume of 40 × 40 × 20 m 3 with 40% water content in each box. Our target volume is located at a depth of 40 m.
A synthetic signal was computed for the one-box model considering nine half-side overlapped loops. The maximum amplitude of the signal was 89 nV. Results of the inversion of these data are presented in figure 5, where the model is shown as a box. The theoretical signal computed after incorporating the inversion results fits the synthetic signal with RMSE = 4.03 nV. Figure 5 shows that the approximate position of the target was correctly found but the inversion accuracy was relatively low. The estimated maximum water content derived from the 3D-SNMR inversion is 29% instead of 40%, as in the model. When we add noise with a magnitude of 20 nV, the inversion becomes less accurate but stable ( figure 6). For this inversion, the theoretical signal fits the synthetic signal with RMSE = 11.3 nV. The estimated maximum water content derived from 3D-SNMR shows 17% instead of 40%, as in the model. The   influence of noise on the inversion results was studied by adding random noise of different magnitudes to the synthetic signals calculated considering the one-box model and the nine-loop setup. Accuracy of the inversion was estimated using equation (4). The results are presented in figure 7, which shows that increased noise magnitude makes inversion less accurate.
We carried out similar modeling with a two-box model. The maximum amplitude of the signal for this model was 67 nV. For demonstration, we present only inversion results of noiseless data and results with 10 nV random noise added to the synthetic signal. The results of inversion of noiseless data ( figure 8) show that the 3D-SNMR inversion has a higher water content anomaly that corresponds approximately to the position of the model. However, even under favorable conditions, two boxes are not resolved and the water content is underestimated. Inversion of noisy data shows even lower resolution (figure 9). Figures 10 and 11 show the synthetic signal inversion computed for a two-box model but assuming 12 quarter-side overlapped loops. The model is resolved with slightly better accuracy but the increased number of loops requires increased fieldwork time. The improvement obtained in resolution may not always be justified by increased duration of the survey. In both examples, the water content in both boxes is underestimated and the location of the model is defined only approximately even under noiseless conditions.

Experimental results
In this section, we present the experimental verification of the newly developed method, 3D surface magnetic resonance imaging. The study area is in the French Alps, the Mont Blanc Massif ( figure 12). The Tête Rousse glacier is located at an altitude of 3200 m. Water that has accumulated in the glacier represents a hazard for the local population. Thus, authorities estimating the danger need to know the volume of accumulated water. Results of glaciological and ground penetrating radar (GPR) studies carried out by the Laboratoire de Glaciologie et Géophysique de l'Environnement, (LGGE) and Institut des Sciences de la terre (ISTerre)     suggested the possible accumulation of melted water in the glacier [22]. GPR data allow imaging of the internal structure of the glacier (glacier/bedrock border line, pebbles, etc) but cannot reliably identify water in ice; thus additional measurements were required. According to the literature, magnetic resonance measurements can be used in investigations of frozen soils and ice [23]- [26]. GPR results showed that within the glacier the ice body was not homogeneous and thus 1D-SNMR might have provided insufficient resolution. Thus, we decided to use the 3D-SNMR method.
The Laboratoire d'étude des Transferts en Hydrologie et Environnement (LTHE) conducted the 3D-SNMR field study in September 2009. The main goal was to locate the accumulated water and estimate its volume. The 3D-SNMR measurements were carried out using a half-side overlapped loop setup consisting of nine 80 × 80 m 2 square loops ( figure 12). During our study, the ambient noise magnitude ranged from 200 to 800 nV. To improve the signal-to-noise ratio, we measured the magnetic resonance signal using 100 stacks for each pulse moment. The magnetic resonance signal ranged between 18 and 156 nV with an average noise magnitude after stacking of 19.8 nV. During this study, we used the NUMIS plus apparatus available from IRIS Instruments. Field data were inverted using only the signal amplitude.

Internal structure of the Tête Rousse glacier
3D-SNMR inversion requires knowledge of the electrical resistivity of the subsurface that could be estimated knowing the internal structure of the glacier. The ice thickness was known from the previous GPR survey carried out by ISTerre [22]. Small attenuation of the GPR signal outside of the water-filled cavern suggested that the ice body has very high resistivity. According to the literature [27]- [29], the ice resistivity depends on the temperature, water content and other impurities but typically ranges from 0.4 × 10 5 m at −2 • C to 4 × 10 5 m at −58 • C. Direct measurements of water from the glacier showed 96 mkS m −1 and consequently ρ water = 104 m. Numerical modeling showed that the subsurface with resistivity larger than 100 m has a small effect on the 3D-SNMR measurements and the homogeneous subsurface with a resistivity of 100 m can be safely used for the inversion.
The internal structure of the glacier is shown in figure 13. Here, the position of the water caverns is shown schematically. Figures 14 and 15 show water distribution in the Tête Rousse glacier derived from the 3D-SNMR inversion. The theoretical signal fitted the measured amplitudes with RMSE = 7.7 nV. The inversion results show that the water is primarily accumulated in the central part of the glacier in a volume of approximately 40 × 100 × 40 m 3 . A second cavern that contains a smaller amount of water was detected 50 m north of the principal reservoir. The 3D-SNMR results show that outside the central area the water content is close to zero, thus suggesting the absence of water. However, the water content in the main cavern (45%) is lower than the 100% expected for bulk water. The relatively low water content can be explained by insufficient resolution of the inversion. Underestimated water content was also observed with the numerical modeling. The second water-filled cavern north of the principal reservoir is located in an area that lacks sufficient coverage of SNMR loops and may be resolved with lower accuracy. The limited surface area of the glacier prevented us from placing additional loops in the northern direction to improve 3D-SNMR resolution.

Locating the principal reservoir
In June 2010, the 3D-SNMR results were used for the installation of 20 boreholes by LGGE. Figures 14 and 15 show that all of the boreholes that intersected water-filled caverns (pink columns) are located within or very close to the reservoir position shown by the 3D-SNMR inversion. Outside the reservoir, boreholes confirmed the absence of water (black columns).   figure 14). Black columns show boreholes that did not detect water; pink columns show boreholes that intersected water-filled caverns.

Estimating water volume
Numerical modeling shows that 3D-SNMR is not able to resolve small target volumes and provides only results averaged over a larger volume. Moreover, figure 2 shows that the same amount of water may produce different signals depending on the water position relative to the loop axis. For example, if most small water-filled caverns and channels are located in areas of low sensitivity (figure 2), 3D-SNMR may underestimate the water volume. But if water is located only in areas of high sensitivity then the volume may be overestimated. Consequently, accurate measurements of the water volume in small structures are not possible in practice, and 3D-SNMR can provide only estimates of minimum and maximum water amounts. Figure 16. Two three-box models of water distribution in the Tête Rousse. Model 1 is shown by black boxes with the water content in the boxes w 1 = 10%, w 2 = 75% and w 3 = 25%; the total volume of water is V = 50 400 m 3 . Model 2 is shown by red boxes with the water content in the boxes w 1 = 10%,w 2 = 60% and w 3 = 20%; the total volume of water is V = 54 400 m 3 ).
For hazard prevention, knowledge of the minimum volume of accumulated water is essential. 3D inversion provided an approximate location of the water in the glacier. Inversion results allow easy numerical modeling using a box model. The number, location and water content of each box are selected so that the theoretical signal computed by water in the boxes fits the experimental data. Many equivalent models can be created using this approach. To estimate the minimum amount of water, we used water only in areas with high and medium sensitivity of SNMR loops. We have found that the use of three-box models makes it possible to fit the experimental data to an acceptable error. For example, two of these models are shown by black and red lines in figure 16.
Our modeling results suggest that the volume of water accumulated in the glacier cannot be less than 50 000 m 3 . Taking into account a ±10% inaccuracy of the SNMR instrument, we lower our estimate of the minimum volume to 45 000 m 3 . In September 2010, 48 000 m 3 of water have been pumped out of the glacier, and residual water in the cavern is estimated to be a few thousands of cubic meters. In addition, the ice temperature in the eastern part of the glacier is close to 0 • C and this ice may contain 1-2% water that may contribute to the measured SNMR signal but cannot be pumped out. Thus, pumping results are consistent with the lower bound estimate on the subsurface water volume provided by 3D-SNMR.

Conclusions
We have shown that the newly developed 3D-SNMR method allows non-invasive investigation of heterogeneous subsurface formations down to about 80 m. We verified our approach using numerical modeling and experimental measurements. Comparison of 3D-SNMR results with information provided by 20 boreholes showed very good agreement. Thus, this method could be recommended for investigating water-filled karst, hard-rock aquifers and other heterogeneous water-saturated formations.
3D-SNMR is a large-scale method and, similar to most surface geophysical methods, it has limited sensitivity and resolution. The method is limited to investigations of subsurface structures that contain more than a few thousands of cubic meters of water. Performance of the method depends on site-specific conditions and signal-to-noise ratio. The position of the target volume relative to the SNMR loop, the amount of water in the subsurface and the magnitude of ambient electromagnetic noise may have significant effects on 3D-SNMR results. Consequently, the sensitivity and resolution of the method should be investigated for each survey with respect to site-specific conditions. 3D-SNMR is a time-consuming method. To minimize time and labor, we used coincident Tx/Rx loops, which may be important for projects located in extreme conditions, such as high mountains. Under Tête Rousse conditions, the duration of one sounding was about 6 h and about 1 h was required to install the loop. So, only one sounding per day was possible. For our numerical study and inversion, we used a Dell Precision M6500 (Intel ® Core TM i7 CPU X 920 at 2 GHz, 7.92 GB RAM) portable computer. With this computer, one 3D inversion of a nine-loop data set takes about 15 min. Usually five to six inversions are sufficient to optimize Tikhonov's functional. The preparatory work necessary for entering data, making preliminary computations and configuring the inversion also requires a few hours.