1D geological imaging of the subsurface from geophysical data with Bayesian Evidential Learning. Part 2: Applications and software

Journal Pre-proof 1D geological imaging of the subsurface from geophysical data with Bayesian Evidential Learning. Part 2: Applications and software Hadrien Michel, Thomas Hermans, Thomas Kremer, Ann Elen, Frédéric Nguyen PII:

S0098-3004(19)30602-8

DOI:

https://doi.org/10.1016/j.cageo.2020.104456

Reference:

CAGEO 104456

To appear in:

Computers and Geosciences

Received Date: 18 June 2019 Revised Date:

20 December 2019

Accepted Date: 20 February 2020

Please cite this article as: Michel, H., Hermans, T., Kremer, T., Elen, A., Nguyen, Fréé., 1D geological imaging of the subsurface from geophysical data with Bayesian Evidential Learning. Part 2: Applications and software, Computers and Geosciences (2020), doi: https://doi.org/10.1016/j.cageo.2020.104456. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

1D geological imaging of the subsurface from geophysical data with Bayesian Evidential Learning. Part 2: Applications and software Authors1: Hadrien MICHEL1,2,3, Thomas HERMANS2, Thomas KREMER1, Ann ELEN4, Frédéric NGUYEN1 (1) University of Liège, Urban and Environmental Engineering Department, Faculty of Applied Sciences, Liège, Belgium, (2) Ghent University, Department of Geology, Ghent, Belgium, (3) F.R.S.-FNRS (Fonds de la Recherche Scientifique), Brussels, Belgium, (4) KU Leuven, Department of Earth and Environmental Sciences, Leuven, Belgium Corresponding author: Hadrien MICHEL University of Liège Allée de la découverte, 9 B-4000, Liège (Belgium) E-mail: [email protected] Declaration of interest: none Link to the code: https://github.com/hadrienmichel/BEL1D Highlights: -

Application of the Bayesian Evidential Learning 1D imaging framework.

-

A set of MATLAB graphical user interfaces is provided.

-

Consistent results for both SNMR and surface waves.

-

BEL1D can be generalized to any type of 1D geophysical data.

1

Authorship statement: Hadrien Michel developed the BEL1D codes and produced the results presented in this paper. Thomas Kremer has significantly contributed and improved all the SNMR related elements developed in this paper. The Mont Rigi SNMR dataset was acquired and pre-processed by Ann Elen for her Master Thesis. Frédéric Nguyen and Thomas Hermans both co-supervised the whole work. 1

Abstract: Bayesian Evidential Learning 1D imaging (BEL1D) is a framework that enables the rapid interpretation of geophysical data in terms of 1D models of the subsurface with associated uncertainty. In this paper, we present a set of MATLAB toolboxes for the interpretation of SNMR data and dispersion curves from surface waves methods under the BEL1D framework, along with some results on field data. The obtained results are compared to the solutions provided by other approaches for the interpretation of these datasets, to demonstrate the robustness of BEL1D. Even if the method is presented for these two specific applications, it is not constrained to any type of geophysical data, as it only requires the knowledge of the forward model to be used. A generalized toolbox is therefore also proposed, only requiring the user to provide a consistent forward model in order to use BEL1D. Keywords:

1D imaging, uncertainty, Prediction Focused Approach, Bayesian Evidential

Learning,

machine

learning,

inversion,

SNMR,

2

Surface

Waves,

MATLAB.

1

1. Introduction

2

Geophysics offers a large ensemble of methods to investigate subsurface properties.

3

However, building consistent images of the subsurface based on geophysical data remains a

4

challenge for all methods. Often, this process is handled by deterministic inversion algorithms

5

(Aster et al., 2013). Such methods consist in minimizing an objective function (often the sum

6

of the data misfit and a model misfit) using gradient-based approaches. Doing so requires the

7

computation of the (often large) Jacobian matrix of the problem at each iteration, a CPU-

8

costly operation. To establish uncertainty estimates on the models, deterministic inversions

9

typically relies on linear error propagation (Aster et al., 2013; Kemna et al., 2007).

10

More recently, due to the rise of computer performances, stochastic methods have

11

been applied to geophysics and are able to quantify more properly uncertainty. Most

12

algorithms rely on Markov chains Monte Carlos (McMC) methods in order to explore the

13

prior model space (Sambridge, 2002). Those methods offer the advantage to provide a much

14

more complete appraisal of the ensemble of possible solutions. Nevertheless, stochastic

15

methods are often linked to high computation costs, since they require a large number of

16

forward model runs to sample the posterior distribution in the prior model space.

17

For example, deterministic inversion is routinely applied for the interpretation of

18

surface nuclear magnetic resonance data (SNMR), where QT-inversion (an inversion that

19

takes into account the full pulse-moment-time dataset) is the state of the art methodology

20

(Müller-Petke and Yaramanci, 2010). To estimate uncertainty on the obtained model, Müller-

21

Petke et al. (2016) proposed a bootstrap algorithm. On the other end of the spectrum, surface

22

waves dispersion curves are more and more often interpreted using stochastic approaches

23

(Socco et al., 2010), using genetic algorithm (GA), simulated annealing (SA), and

24

neighborhood algorithm (NA) which all sample models randomly under a guidance system.

3

25

In this paper, we present applications of a new stochastic method for imaging the

26

subsurface called Bayesian Evidential Learning 1D imaging (BEL1D), a new framework for

27

modeling the subsurface with geophysical data relying on the constitution of statistics-based

28

relationships between model parameters and data from geophysical experiments. The

29

algorithm consists of following steps (Figure 1Error! Reference source not found.) (see

30

companion paper ‘1D geological modeling of the subsurface from geophysical data with

31

Bayesian Evidential Learning Part 1: Method and demonstration’):

32

•

33 34

(yellow box in Figure 1) and associated data

(blue box in Figure 1) by forward modeling •

Reduction of the dimensionality of the data and model using principal ( )

component analysis (PCA):

35 36

Generation of prior models

•

and

( )

Constitution of statistical relationships between the models parameters and the

37

reduced data (green cloud in Figure 1), using canonical correlation analysis

38

(CCA) :

39

•

(

,

distribution (

41

density estimators •

(

,

)

Generation of posterior distributions of

40

42

)

,

) to field data

by constraining the bivariate (red box in Figure 1) using kernel

Sampling of the constituted distributions and back-transformation of the

43

samples into the original space

44

subsurface

45

(purple box in Figure 1)

, delivering a set of 1D models of the

that are constrained to the knowledge of the geophysical data

46

The method offers the advantage not to require multiple runs of the CPU-costly

47

inversion of the data, as it is based solely on the use of the forward model. BEL1D offers an

48

alternative to some stochastic algorithms, often requiring much larger computation times, as it 4

49

relies on a limited number of independent prior samples to deduce the posterior distribution.

50

One of the key advantage of the approach is that the relationship between the model

51

parameters and the data can be reused for several data sets, as long as they belong to similar

52

prior model, largely reducing the computational cost of BEL1D. On the other hand, the

53

method is also an efficient alternative to deterministic inversion, as the CPU cost remains

54

acceptable and the results are much more informative due to their probabilistic nature. Above

55

these advantages, BEL1D does not require the definition of regularization parameters or of

56

any other form of bias. In the algorithm presented in this paper, we only consider a finite

57

number of layers in the models. This approach,

58

methodologies (e.g. García-Jerez et al., 2016; Günther and Müller-Petke, 2012; Müller-Petke

59

et al., 2016; Wathelet, 2008), reduces significantly the number of parameters to consider.

applied in numerous other inversion

60 61

Figure 1: Illustration of the BEL1D process.

5

62

In this contribution, we validate the effectiveness of BEL1D for two types of

63

geophysical data: a dispersion curve from surface wave analysis and the free induction decay

64

curves from surface nuclear magnetic resonance. Surface waves are seismic waves whose

65

propagation are parallel to the earth’s surface. They propagate according to mechanical and

66

geometric properties of the subsurface. Their amplitude is generally superior to body waves,

67

making their acquisition much easier. The propagation of these waves is governed by

68

geometric dispersion, often modeled by the dispersion curve. In most interpretations, the

69

focus is made on the fundamental mode dispersion curve, hence the strongest one (Socco et

70

al., 2010). The current state-of-the art inversion method for the interpretation of dispersion

71

curves is an improved version of the neighborhood algorithm implemented in Geopsy

72

(Wathelet, 2008). The resulting models are expressed in terms of P-wave velocity (

73

wave velocity ( ) and density ( ). However, surface waves are mostly sensitive to S-wave

74

velocity (Xia et al., 1999). Often, a link between

75

Poisson’s ratio ( ) (Equation 1) as is the case in Geopsy (Wathelet, 2008).

=

and

), S-

is considered through the

&

0.5 − ' × % 1−

(1)

76

Surface Nuclear Magnetic Resonance (SNMR) is a near-surface geophysical method

77

used to retrieve information about the groundwater in the first hundred meters of the

78

subsurface. The method relies on the quantum mechanics behavior of protons composing the

79

water molecules. At rest, the protons are precessing about the Earth’s magnetic field at their

80

Larmor frequency. By applying a magnetic field tuned at this precise frequency, we modify

81

the thermal equilibrium of the proton’s spin. As soon as the injection is stopped, the spins

82

relax back to their thermal equilibrium, emitting a secondary magnetic field that is registered

83

through induced currents in the receiver antenna (free induction decay). The obtained signal

84

can be interpreted in terms of water content and decay time, the latter being correlated with

6

85

hydrogeological parameters of the subsurface such as pore sizes and hydraulic conductivities

86

(Behroozmand et al., 2015). The current state-of-the-art interpretation method is the QT

87

inversion (Müller-Petke and Yaramanci, 2010). This inversion method uses the full QT (pulse

88

moment-time) dataset originating from the field acquisition, in order to estimate both the

89

water content (or partial water contents) and the relaxation time. The method is fully

90

implemented in MRSmatlab (Müller-Petke et al., 2016), a set of toolboxes for the processing,

91

modeling and interpretation of SNMR data.

92

The results for the two tested cases are compared on different aspects. First, the

93

accuracy of the obtained distributions is quantified, either comparing the results to other

94

studies, to benchmark values, or to results originating from state-of-the-art inversion codes.

95

Second, the efficiency is discussed, both in terms of computation time and of memory usage,

96

and compared to the state-of-the-art method for deterministic inversion.

97

2. Implementation of BEL1D

98

2.1. Program workflow

99

The code is built on the Bayesian Evidential Learning 1D imaging (BEL1D) presented in

100

the companion paper “1D geological modeling of the subsurface from geophysical data with

101

Bayesian Evidential Learning Part 1: Method and demonstration“.

102

In the proposed codes, the models are described as a succession of finite number of

103

horizontal layers of unknown thickness, with associated parameters that are related to the

104

geophysical experiment (water content, velocities, densities, etc.). This leads to small

105

numbers of uncorrelated parameters. Although the latter is not strictly required by BEL1D, it

106

makes it straightforward to define the prior in a general way. Then, BEL1D requires the

107

generation of geophysical datasets associated with the models, i.e. running the geophysical

108

forward model. 7

109

The method relies on the constitution of statistical relationships between models

110

parameters and the associated datasets. Contrary to classical inversions processes, this method

111

does not require the use of regularization parameters that are known to bias geophysical data

112

interpretation (Aster et al., 2013). Moreover, this method can explore large prior model spaces

113

as efficiently as narrowed ones, due to its construction. It is therefore a very efficient method

114

for the construction of models from geophysical data without the introduction of prior biasing

115

information, contrary to other Bayesian methods that sometimes use unrealistic small prior to

116

converge faster.

117

As BEL1D is ubiquitous, a generalized toolbox has first been developed and can be

118

used with any forward model available. The input for the GUI is the data (A) in a MATLAB

119

‘.mat’ file containing a matrix called data. The user can specify the names and units of the X

120

and Y variables on the graph. Those names are later used in the model display (RMS error

121

scale bar). The user can then describe the prior model space in panel B. There, the names and

122

the ranges of the parameters (up to 3 parameters plus the thickness) can be set for each layer

123

and the sampler selected.

8

124 125 126

Figure 2: Graphical user interface layout in the generalized case.

2.2. Implementing Surface Waves as a toolbox in BEL1D

127

In the MATLAB graphical user interface (GUI) (Figure 3), the models are presented

128

as sets of a finite known number of layers, with unknown thicknesses, seismic P and S

129

velocities and densities. The last layer is considered an infinite half-space. At the user 9

130

requests, it is possible to insert constraints on the velocities through an interval of Poisson’s

131

ratio and to only consider increasing velocities and densities through the layers. Such

132

conditions will relax the independence between the parameters. The used data is the

133

fundamental mode of the dispersion curve (however, it should be stated that BEL1D could be

134

used with multiple modes). This dispersion curve is exploited in the frequency/phase velocity

135

space. The forward modeling is performed via the GPDC command line tool provided by

136

Geopsy (Wathelet, 2008) called through the “jsystem” MATLAB function (Rosenberg, 2018),

137

a function similar to the MATLAB-integrated ‘system’ function.

138

In this case, the user inputs are in the data panel (A), the file describing the dispersion

139

curve is selected. The file formatting is the same as the one from Geopsy (Wathelet, 2008);

140

the prior model space description (B) is similar to the SNMR case, with the choice of the

141

number of layers, the ranges for the different parameters, the type of sampler and the number

142

of models for the prior realizations. However, in this case, the user can also specify if the

143

models should have increasing velocities and densities with depth (checkbox “increasing”)

144

and whether or not to consider Poisson’s coefficient (NaNs to ignore and doubles to consider).

145

Finally, once all the inputs are present, the user can launch BEL1D in panel C.

10

146 147

Figure 3: Graphical user interface layout for the interpretation of dispersion curves from multichannel analysis

148

of surface waves.

149

2.3. Implementing Surface Nuclear Magnetic Resonance as a toolbox in BEL1D

11

150 151

Figure 4: Graphical user interface layout for the interpretation of (multiple-loops) SNMR datasets

152

In this section, we present a new MATLAB graphical user interface (GUI) for the

153

interpretation of SNMR datasets through BEL1D. This GUI is adapted to multiple-loops

154

datasets but is also able to manage more classical single receiver experiments. For SNMR, the

155

models are described through a finite number of layers, with unknown thicknesses (e), water

156

contents (W), and decay times (T2*). The data formatting is the one issued from MRSmatlab

157

(Müller-Petke et al., 2016) as is the kernel function.

158

The input parameters for the GUI are separated in three different panels: data (A), Prior model

159

space (B) and kernels (C). To enter SNMR data, the user must specify first the path to the

160

folder containing the data and then the name(s) of the file(s) containing the experimental data

161

in the table (A1). As the GUI is adapted to multiple loops configurations, the user can enter

162

multiple files corresponding to different acquisitions performed on the same site with 12

163

different configurations. Then, the user must select the configuration that they want to use

164

(A2). To do so, they choose a transmitter size in column 1 and check the chosen receivers

165

used with this transmitter. If multiple files where selected, the configuration can differ for all

166

files, leading to complementary information. The prior model space definition is performed in

167

panel B. The user can choose the number of layers (up to 6), describe the parameter ranges,

168

select the type of sampler and the number of prior realizations. Finally, the kernels for the

169

forward modeling (C) are either computed under the assumption of a resistive Earth or a 100

170

Ohm.m resistivity, either loaded from previously computed files (that are located in the same

171

folder as the data) to handle more complex resistivity models. The software can handle up to

172

16 different configurations at the same time (up to four transmitters and four receivers. Both

173

the data and the kernels are in the MRSmatlab format (Müller-Petke et al., 2016). Finally,

174

once all the inputs are present, the user can launch BEL1D in panel D.

175

3. Results

176

The different toolboxes that are presented are automatically generating most of the

177

different figures that are presented in the first part of those companion papers: the CCA space,

178

the load of each parameter in the different CCA space dimensions, the posterior parameters

179

and models distributions. However, in this section, we will mainly focus on the posterior

180

(parameters and models distributions) analysis and compare it with our knowledge of the

181

sites. Nonetheless, a brief analysis of the CCA space will be proposed as it may give an

182

insight on the sensitivity of the experiment towards the parameters of interest.

183

3.1. Case 1: InterPacific – Surface Waves

184

The InterPacific project (Intercomparison of methods for site parameter and velocity

185

profile characterization) aimed at assessing the reliability of in-hole and surface wave

186

methods for the estimation of shear wave velocity (Garofalo et al., 2016b). In this section, we

13

187

will use the mean dispersion curve that arose from these experts’ analyses of the raw data to

188

demonstrate BEL1D applicability to retrieve the shear wave velocity profile. We will focus on

189

the case proposed by Garofalo et al. (2016b) in the Mirandola test-site located in the Po plain.

190

The site is characterized by soft sediments (sand, clays and gravels) overlaying a Pliocene

191

bedrock (sand and marine pelite), found at a depth between 50 and 150 meters (Garofalo et

192

al., 2016b, 2016a). We will here assume that the soil is constituted of three layers: two soft

193

layers, representing the alluvial plain and its variations and a third layer representing the

194

bedrock. A full description of the prior model space is given in Table 1, with the additional

195

information that the parameters

196

coefficient is always within 0.2 and 0.5 for each layer. Thickness (e) [m]

,

and

P-wave (

are increasing with depth and that the Poisson’s

velocity S-wave

) [m/s]

velocity Density

( ) [m/s]

( )

[kg/m³]

Min

Max

Min

Max

Min

Max

Min

Max

Layer 1

5

50

200

4000

100

500

1500

3500

Layer 2

45

145

200

4000

100

800

1500

3500

200

4000

300

2500

1500

3500

Layer 3 197

Table 1: Mirandola – Prior model space

198

We generated 5000 prior realizations and simulated their fundamental mode dispersion

199

curve (the associated data). To sample the models, we used uniformly distributed variables

200

with a rejection sampler in order to satisfy the increase of the parameters values with depth

201

and the Poisson’s coefficient interval. From the original 64 dimensions of the data (the 64

202

couples frequency-phase velocity of the sampled dispersion curve), PCA reduced the

203

dimensionality to 11 (the minimum to be able to perform further steps of BEL1D, as the

204

dimension of

205

– the (PCA reduced) model space). In the CCA model space, the first two dimensions (where

(

– the PCA reduced dataset – must be larger or equal to the one of

14

(

(=

)

206

shear-wave velocities and thicknesses are gathered) show high correlations but all the 9

207

remaining dimensions presents scattered values (Figure 5). This implies that the field

208

experiment may enable a significant reduction of uncertainty for the shear-wave velocities

209

and thicknesses, whereas no substantial reduction is expected for the other parameters. Due to

210

the lack of existing error models for dispersion curves (to our knowledge), BEL1D is applied

211

without specific noise analysis. This is further discussed later on. The process results in the

212

production of 1000 models in the posterior.

213 214

Figure 5: Mirandola – Data versus model space in the CCA space (first 5 dimensions out of 11). The blue dots

215

represent each prior realizations and the orange line the position of the field data.

15

216 217

Figure 6: Mirandola - Comparison of the prior (blue) and posterior (orange) parameters distributions.

218

As expected, the main reductions of uncertainty are obtained for shear-wave velocities

219

(Figure 6), which are known to be the most sensitive to surface waves. Then, the reduction of

220

uncertainty is substantial for compression waves velocities (especially for the third layer), due

221

to the relationship constrained by the Poisson’s coefficient interval. Finally, the model is

222

nearly insensitive to the density in agreement with previous studies (e.g.: Xia et al., 1999). No

223

particular correlation is observed between the parameters as all but one Pearson’s correlation

224

coefficient are below 0.5 in absolute value, the only exception being observed between

16

,)

,)

225

and

226

coefficient.

(0.59) due to the limitation of possible values for

, constrained by the Poisson’s

227

The depth of the interface with the bedrock still remains very uncertain, but the most

228

probable value matches borehole (Garofalo et al., 2016a) measurements (Figure 7) whereas

229

the prior space did not imply such values (borehole information were not taken into account to

230

build the prior model space). Moreover, our uncertainty covers the range proposed by the

231

different experts (Figure 8). The comparison of the results obtained through BEL1D with the

232

results from the InterPacific project (Figure 8) shows that BEL1D provides a posterior

233

distribution that encompass the InterPacific results for the alluvial plain. However, the

234

bedrock seems to present larger

235

This may originate from the use of the average dispersion curve from the aggregated data

236

issued by the whole set of analysts, instead of the individual dispersion curves of each analyst.

values compared to the results of 5 out of the 14 analysts.

237 238

Figure 7: Mirandola – Estimated depth to the bedrock. The borehole value originates from Garofola (2016a).

17

239 240

Figure 8 : Mirandola – Posterior models distributions and their associated RMS error. The black lines in the

241

profile represent the models obtained from surface wave methods by the experts that participated in the

242

InterPacific project (Garofalo et al., 2016b).

243

Figure 9 shows that BEL1D was able to reject unlikely areas in the data space. From

244

this figure, we can see that the posterior distribution of dispersion curves all exhibits a

245

stepwise decrease, also observed in the prior distribution, instead of the steady decrease of the

246

experimental phase velocity observed from 2 to 15 Hz. We attribute this to the limited number

247

of layers used in our prior, inherent to a blocky imaging framework. From Figure 9, it appears

248

that most of the dispersion curves picked by the experts teams are fitting inside the most

249

probable area (lowest RMS). If we assume that the range of dispersion curves picked by the

250

expert group is representative of the noise level, BEL1D samples posterior models within this

18

251

range. Our result is therefore valid, even when accounting for the intrinsic uncertainty on the

252

dispersion curve.

253 254

Figure 9 : Mirandola – Comparison of the prior (grey) and posterior (blue to red) data spaces. The colors of the

255

posterior dispersion curves represent the RMS error with the same scale as the one used in Figure 8 (blue = low

256

RMSE and red = high RMSE). The white dispersion curve represents the average dispersion curve as defined in

257

Garofalo (2016b). The black dots represent picked dispersion curves from the different teams of experts.

258

To further validate our results, we used the same dataset in Geopsy (Wathelet, 2008)

259

in order to produce models that explain the dataset (Figure 10). We used the same prior model

260

space to run the inversion. The results show that the models that explain best the dataset are

261

similar to the one resulting from BEL1D.

19

262 263 264

Figure 10: Mirandola – Results from Geopsy on the dataset.

3.2. Case 2: Mont Rigi (Belgium) – SNMR

265

The Mont Rigi is located in the Belgian Fagnes region, in the Eastern part of the

266

country. This site presents an ideal case study for SNMR as it is remote and far from any

267

electromagnetic noise sources, due to its natural reserve classification.

20

268 269

Figure 11: Peat thickness at the Mont Rigi test site according to GPR prospection. The coordinates are in

270

Belgian Lambert 72 (EPSG:31370). (Modified from Wastiaux and Schumacker, 2003)

271

Geologically, the site is characterized by a metric peat layer on top of a Cambrian

272

bedrock (La Venne formation) known as an aquiclude with a very low water content (Gilson

273

et al., 2017). In contrast, peat is known to present very high water contents with observed total

274

porosities around 90% (Wastiaux, 2008). According to previous GPR exploration (Wastiaux

275

and Schumacker, 2003), the peat layer should have a thickness between 2.5 and 4.5 meters

276

(Figure 11). The SNMR response should therefore allow to easily distinguish the two layers

277

with a properly designed experiment. The designed on-site experiment consisted of a classical

278

coincident loops transmitter/receiver couple with a diameter of 20 meters.

21

279 280

Figure 12 : Extract of the obtained SNMR signal from Mont Rigi after processing in MRSmatlab (Müller-Petke

281

et al., 2016).

282

Figure 12 shows the measured data set, after signal processing operations (despiking,

283

harmonic modelling, reference noise cancellation and despiking again). We see that the noise

284

amplitude is quite low (18 nV), and a clear decay is observed for the smallest pulse moments

285

(those most sensitive to the first meters where the peat layer is located).

286

First, the signals have been truncated to 0.2 seconds to decrease the impact of noise,

287

mostly present at large timeframes, where the NMR signal is absent. Then, the signals have

288

been inverted using the QT Inversion (Müller-Petke and Yaramanci, 2010), to constitute a

289

benchmark to compare to the results of BEL1D. We used a smooth-mono model description

290

(smooth models with a single relaxation time for a given depth) and a regularization

291

parameter of 6000 for both water content and relaxation time using the L-curve criteria. Then, 22

292

BEL1D was applied with the uniformly distributed prior model space described in Table

293

2Error! Reference source not found.. We sampled 1000 prior realizations and produced

294

1000 models for the posterior.

295

Thickness [m]

Water content [%]

Relaxation time [ms]

Minimum

Maximum

Minimum

Maximum

Minimum

Maximum

Layer 1

0

7.5

30

80

0

200

Layer 2

/

/

0

15

100

400

Table 2: Prior model space for the Mont Rigi SNMR experiment

296

The CCA space analysis (Figure 13) shows that we could expect a significant

297

reduction of uncertainty for some parameters, as the first three dimensions show high

298

correlations. The first dimension, dominated by the water content of the half-space shows an

299

especially narrow relationship between the reduced data and model space. On the other hand,

300

the last dimension, mostly represented by the water content of the first layer, shows no

301

specific correlation. We should therefore expect a much better reduction for the parameters of

302

the half-space than for the first layer.

23

303 304

Figure 13: Mont Rigi – Data versus model space in the CCA space. The blue dots represent each prior

305

realizations and the orange line the position of the field data.

306

The noise level (after denoising) impact on BEL1D was assessed by analyzing noise

307

propagation in the PCA space (Hermans et al., 2016) for a random sample of 50 models from

308

the prior. The values of the covariance matrix

309

of the order of 10*&+ ). However, once transformed into CCA space, it resulted in significant

310

changes in bandwidths for the kernel density estimation, with between 0.04 and 0.13 instead

311

of the default 0.01 value.

(

in PCA space are relatively low (maximum

24

312 313

Figure 14 : SNMR at Mont Rigi - Prior (blue) and posterior (orange) distributions of the parameters.

314

The analysis of the posterior parameters model space (Figure 14) shows that the

315

reduction of uncertainty from the prior model space to the posterior is significant and

316

produces consistent results, when compared to the geological context and the QT inversion

317

solution (dashed gray lines in Figure 15). The thickness of the peat layer tends to be smaller

318

than the maximum estimate (4.5 m) from to the GPR interpretation (Wastiaux and

319

Schumacker, 2003). This is explained by the equivalence of the produced signal from a high

320

water content in a thin layer and a lower water content in a thick layer. The use of a narrower

321

prior taking into account the GPR information would prevent such a behavior.

25

322

No specific correlation between the parameters is arising from BEL1D. However, we

323

observe that some relations are created. For example, the maximum value of the thickness for

324

the first layer tend to decrease when the water content of the aforementioned layer increases.

325

This behavior arises from the non-uniqueness of the SNMR response stated above.

326

Analyzing the set of models produced by BEL1D (Figure 15), it is observed that a

327

large majority of the models show a root mean square error lower than the noise level in the

328

dataset, implying that the set of probable models provides an efficient estimation of the sites

329

parameters. The RMS errors of the posterior samples show that they all explain the data to a

330

similar level; the uncertainty is thus intrinsic to the data set and the non-unicity of the

331

solution. The non-uniqueness of the solution is easily observable, with two distinct behaviors:

332

one with a median water content and a thick first layer and another with a very high water

333

content but a thinner first layer. Those two cases show similar RMS errors, hence are

334

explaining the data at the same level. The probability of each posterior model is shown in the

335

RMS error colorbar for which each color represents the same number of models. Therefore, it

336

is observed that the low RMS error models also correspond to the most probable ones.

26

337 338

Figure 15: SNMR at Mont Rigi - Posterior models distributions with their associated misfit and results of the QT

339

inversion (dashed gray line). The black line in the RMS error graph represents the level of noise in the original

340

dataset.

341

In the data space, it is observed that the simulated data (exponential decays) from the

342

posterior space are corresponding to the experimental data (Figure 16). We show that, even if

343

the raw data are noisy, BEL1D is able to dramatically reduce the uncertainty on the posterior

344

distributions, centering the posterior data on the experimental ones. However, the fit is better

345

for the first pulse moments, hence the most sensitive to this particular soil configuration.

27

346 347 348

Figure 16: SNMR at Mont Rigi - Comparison of the prior (gray), posterior (blue) and experimental (red) data.

3.3. CPU timing and memory usage

349

In this section, we monitor computation time and RAM usage using the MATLAB®

350

integrated profiler (Table 3). The reported computation time only accounts for BEL1D run in

351

the GUI (‘pushbutton_Run_Callback’). The RAM usage corresponds to the peak memory

352

observed when using all the possibilities of the GUI (run BEL1D, show parameters, models

353

and save results). All the tests have been performed on a workstation equipped with an Intel® 28

354

Core™ i7-6800K hexacore CPU with 64 Go of RAM at 2400MHz and a GTX 1050 GPU in

355

MATLAB R2018a. Even if the GUIs enable the optional use of the parallel computing

356

toolbox of MATLAB, this toolbox was not used for these tests. Computation time [sec] Peak Memory [Gb] (forward run time/noise impact estimation) Case 1 (5000 models in the prior)

139.640

(90.273/Not 0.178

Available) Neighborhood algorithm on case 1

<5

Not available

Case 2 (1000 models in the prior)

46.951 (0.575/16.917)

0.118

QT inversion on case 2

7.717

0.012

L-curve on case 2 (QT inversion)

71.905

0.043

(2977 models visited)

357

Table 3: Computation time and RAM usage

358

When comparing the results of BEL modeling on the second case to QT inversion in

359

MRSMatlab (Müller-Petke et al., 2016) performances, it is observed that QT inversion is

360

faster and uses way less memory. However, this timing assumes that we know the optimal

361

regularization parameter in advance. If we want to account for the computation of the L-curve

362

in the CPU time, we should add 71.905 seconds to the computation time. Moreover, the

363

solution obtained through QT-inversion is deterministic, whereas BEL1D provides a

364

stochastic solution in a timeframe that is very similar.

365

In the case of surface waves, it is observed that the neighborhood algorithm (Wathelet,

366

2008) performs way better in terms of CPU time than BEL1D. This is mainly explained by

367

the highly optimized forward model used in Geopsy, that is strongly hindered by multiple

29

368

system calls from MATLAB (jsystem function (Rosenberg, 2018)) followed by read and write

369

operations from a hard drive in our case.

370

From Table 3, it is observed that the performances of BEL1D highly depends on the

371

efficiency of the forward model run. The other operations mainly depend on the complexity of

372

the dataset (PCA is more demanding for SNMR than for surface waves as the datasets are

373

much larger) and the number of models in the prior (mainly hindering the kernel density

374

estimation). At equal performance of the forward model, the efficiency of BEL1D is thus

375

directly proportional to the number of models sampled from the prior.

376

Overall, the use of BEL1D appears as an efficient method to model the earth based on

377

geophysical data, both in terms of accuracy, CPU usage and memory cost. In the case of large

378

campaigns, a set of prior realizations and the associated CCA spaces could be created once

379

and reused for the interpretation of all the data, removing the costly forward model runs from

380

the computation times, as well as the dimensionality reduction steps thus leading to very rapid

381

stochastic interpretations (less than a minute). For example, a SNMR survey at Mont-Rigi to

382

determine the thickness of the peat layer could use the same prior for all the interpretation of

383

the sounding curves.

384

4. Conclusion

385

We successfully applied BEL1D to two different field cases using completely different

386

geophysical methods: surface waves analysis and surface nuclear magnetic resonance. The

387

first case showed the applicability of BEL1D to dispersion curves obtained by analysis of

388

surface waves. We showed that the resulting posterior models were coherent with the range of

389

results independently obtained by experts (Garofalo et al., 2016b) and is inherent to the

390

uncertainty around the dispersion curve. Moreover, we compared with the results from a

391

direct search inversion technique, the Neighborhood Algorithm, and obtained similar results. 30

392

In addition, we observed that the most likely depth to the bedrock as originated from the

393

process corresponds to the depth measured in the borehole.

394

For the second case, BEL1D has been applied to SNMR data acquired in the Belgian

395

Fagnes. We showed that the process outputs reasonable ranges of uncertainty on the

396

parameters and includes the QT Inversion results. The careful analysis of the RMS error of

397

the proposed models indicates that the range of uncertainty delivered by BEL1D is intrinsic to

398

the non-unicity of the solution for geophysical inverse problems: many models explain the

399

data to the same level.

400

Those two applications demonstrate the versatility of BEL1D and its applicability to

401

multiple geophysical methods with very few adaptations. Provided that a forward model

402

already exists, the adaptation of the codes to any other 1D method would be straightforward.

403

For example, the forward operator from HV-Inv (García-Jerez et al., 2016) could be used for

404

the interpretation of H/V data. System calls may also help further extend the applicability of

405

BEL1D. This may be the case for the use of all the 1D forward operators implemented in

406

pyGIMLi (Rücker et al., 2017) – vertical electrical sounding, frequency- and time-domain

407

electromagnetic and magnetotelluric.

408

BEL1D also gives a total freedom in the definition of the prior model as it is efficient

409

with both large and narrow prior. As demonstrated with the surface wave, constraints between

410

parameters can be easily implemented.

411

In term of performance, the computation cost of BEL1D is directly proportional to the

412

efficiency of the forward model and the number of samples in the prior model space. It must

413

be noted that the prior samples being independent, the simulation of synthetic data can be

414

fully parallelized. The performance of BEL1D can thus be easily estimated from the available

415

computing facilities. The cost is therefore much smaller than other stochastic methods

31

416

requiring ten to hundreds thousands of runs to converge towards a solution. In the SNMR

417

case, we show that BEL1D provides the posterior distribution in a computation time similar to

418

the state-of-the-art deterministic QT inversion. Moreover, computation timing could be

419

dramatically reduced by pre-field forward modelling operations for large geophysical

420

campaigns where the prior would be similar. Prior realizations and the associated datasets

421

would then be reused (PCA and CCA operations also), letting only the kernel density

422

estimation, sampling and back-transformation operations to be applied after field acquisition,

423

leading to extremely rapid imaging.

424

In its present version, BEL1D is dedicated to the 1D inversion of geophysical data in

425

layered Earth model where the contrast between layer is supposed to be sharp. In the future,

426

we plan to extend the method to smoothly varying systems accounting for a large number of

427

thin, correlated layers.

428

5. Acknowledgment

429

The authors would like to thank the contributors to the development of MRSmatlab,

430

the cornerstone of SNMR data interpretation in BEL1D. Marc Wathlet directed us towards the

431

InterPacific project for the surface wave application. Professor S. Foti and Federico Passeri

432

shared with us the raw dispersion curves from the different teams of experts for the Mirandola

433

case. Hadrien Michel is a fellow funded by the F.R.S.-FNRS.

434

6. Computer codes availability

435

The different MATLAB codes that have been used to produce the results presented in

436

this paper are available at github.com/hadrienmichel/BEL1D. They are released with this

437

paper under the BSD-2-Clause license. The codes consist in a series of graphical user

438

interfaces for the interpretation of SNMR data (BEL1DSNMR), dispersion curves from

439

surface waves analyses (BEL1DSW) and for global cases, assuming the user provides a 32

440

forward model (BEL1DGENERAL). All those functions are standalone (no specific toolbox

441

required) in MATLAB but can run faster when the Parallel Computing toolbox is available.

442

The main codes have been developed by H. Michel (ULiège, Belgium):

443

Hadrien MICHEL

444

Urban and Environmental Engineering Department

445

Allée de la découverte, 9

446

B-4000 Liège

447

Belgium

448

The forward models that are used are Open Source MATLAB toolboxes for SNMR

449

(MRSMatlab: Müller-Petke et al., 2016) and command line tools for surface waves (Geopsy:

450

Wathelet, 2008).

451

In order to run the codes, a MATLAB license is required. Minimum requirements for

452

the computations are the ones for your MATLAB license, at the exception of RAM, for which

453

a minimum of 8 Go is recommended and more could be required for more complex

454

computations. If using the parallel capabilities of MATLAB (parallel computing toolbox), a

455

large number of cores is preferred and, globally, a 4-core (or more) processor is advised.

456

7. References

457

Aster, R., Borchers, B., Thurber, C., 2013. Parameters Estimation and Inverse Problems, 2nd

458

ed. Academic Press, New York. 376pp.

459

Behroozmand, A.A., Keating, K., Auken, E., 2015. A Review of the Principles and

460

Applications of the NMR Technique for Near-Surface Characterization. Surveys in

461

Geophysics 36, 27–85. https://doi.org/10.1007/s10712-014-9304-0

33

462

García-Jerez, A., Piña-Flores, J., Sánchez-Sesma, F.J., Luzón, F., Perton, M., 2016. A

463

computer code for forward calculation and inversion of the H/V spectral ratio under

464

the

465

https://doi.org/10.1016/j.cageo.2016.06.016

diffuse

field

assumption.

Computers

&

Geosciences

97,

67–78.

466

Garofalo, F., Foti, S., Hollender, F., Bard, P.Y., Cornou, C., Cox, B.R., Dechamp, A.,

467

Ohrnberger, M., Perron, V., Sicilia, D., Teague, D., Vergniault, C., 2016a.

468

InterPACIFIC project: Comparison of invasive and non-invasive methods for seismic

469

site characterization. Part II: Inter-comparison between surface-wave and borehole

470

methods.

471

https://doi.org/10.1016/j.soildyn.2015.12.009

Soil

Dynamics

and

Earthquake

Engineering

82,

241–254.

472

Garofalo, F., Foti, S., Hollender, F., Bard, P.Y., Cornou, C., Cox, B.R., Ohrnberger, M.,

473

Sicilia, D., Asten, M., Di Giulio, G., Forbriger, T., Guillier, B., Hayashi, K., Martin,

474

A., Matsushima, S., Mercerat, D., Poggi, V., Yamanaka, H., 2016b. InterPACIFIC

475

project: Comparison of invasive and non-invasive methods for seismic site

476

characterization. Part I: Intra-comparison of surface wave methods. Soil Dynamics

477

and

478

https://doi.org/10.1016/j.soildyn.2015.12.010

Earthquake

Engineering

82,

222–240.

479

Gilson, M., Briers, P., Ruthy, I., Dassargues, A., 2017. Carte hydrogéologique de Wallonie,

480

Planchettes Elsenborn - Langert - Dreiherrenwald n°50/3-4, 50A/1. Service public de

481

Wallonie, DGO 3, Belgique. [In french]

482

Günther, T., Müller-Petke, M., 2012. Hydraulic properties at the North Sea island of Borkum

483

derived from joint inversion of magnetic resonance and electrical resistivity

484

soundings.

485

https://doi.org/10.5194/hess-16-3279-2012

Hydrology

and

Earth

System

34

Sciences

16,

3279–3291.

486

Hermans, T., Oware, E., Caers, J., 2016. Direct prediction of spatially and temporally varying

487

physical properties from time-lapse electrical resistance data. Water Resources

488

Research 52, 7262–7283. https://doi.org/10.1002/2016WR019126

489

Kemna, A., Nguyen, F., Gossen, S., 2007. On linear model uncertainty computation in

490

electrical imaging. Presented at the SIAM Conferance on Mathematical and

491

Computational Issues in Geosciences, Santa Fe.

492

Müller-Petke, M., Braun, M., Hertrich, M., Costabel, S., Walbrecker, J., 2016. MRSmatlab —

493

A software tool for processing, modeling, and inversion of magnetic resonance

494

sounding data. GEOPHYSICS 81, WB9–WB21. https://doi.org/10.1190/geo2015-

495

0461.1

496

Müller-Petke, M., Yaramanci, U., 2010. QT inversion — Comprehensive use of the complete

497

surface

NMR

data

set.

498

https://doi.org/10.1190/1.3471523

GEOPHYSICS

75,

WA199–WA209.

499

Rosenberg, A.A., 2018. jsystem: Fast drop-in replacement for matlab’s slow “system”

500

command: avivrosenberg/matlab-jsystem. https://github.com/avivrosenberg/matlab-

501

jsystem (Accessed on 16th of June 2019).

502

Rücker, C., Günther, T., Wagner, F.M., 2017. pyGIMLi: An open-source library for

503

modelling and inversion in geophysics. Computers & Geosciences 109, 106–123.

504

https://doi.org/10.1016/j.cageo.2017.07.011

505 506

Sambridge, M., 2002. Monte Carlo methods in geophysical inverse problems. Reviews of Geophysics 40 (3), 1009. https://doi.org/10.1029/2000RG000089

507

Socco, L.V., Foti, S., Boiero, D., 2010. Surface-wave analysis for building near-surface

508

velocity models — Established approaches and new perspectives. GEOPHYSICS 75,

509

A83-A102. https://doi.org/10.1190/1.3479491

35

510 511

Wastiaux, C., 2008. Les tourbières sont-elles des éponges régularisant l’écoulement ? Bulletin de la Société géographique de Liège 50, 57-66. [In french]

512

Wastiaux, C., Schumacker, R., 2003. Topographie de surface et de subsurface des zones

513

tourbeuses des réserves naturelles domaniales des Hautes-Fagnes (No. C60/1-2-3).

514

Université de Liège, Service de phytosociologie, flore et végétation des Hautes-

515

Fagnes. [In french]

516

Wathelet, M., 2008. An improved neighborhood algorithm: Parameter conditions and

517

dynamic

518

https://doi.org/10.1029/2008GL033256

519

scaling.

Geophysical

Research

Letters

35,

L09301.

Xia, J., Miller, R.D., Parl, C.B., 1999. Estimation of near-surface shear-wave velocity by

520

inversion

of

Rayleigh

521

https://doi.org/10.1190/1.1444578

waves.

522

36

GEOPHYSICS

64,

691–700.

-

A new approach for 1D imaging from geophysical data. A set of MATLAB graphical user interfaces is provided. Estimation of the uncertainty on the model parameters. Consistent results for SNMR. BEL1D can be generalized to any type of 1D geophysical data.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: