BAYESIAN SIMULATION-BASED LEARNING: USING NUMERICAL MODELS TO QUANTIFY UNCERTAINTY IN THE SUBSURFACE

Prof. Dr. Ir. Thomas Hermans, Associate professor, Department of Geology

DEPARTMENT GEOLOGY

RESEARCH GROUP HYDROGEOLOGY AND APPLIED GEOLOGY

… & collaborators, postdocs and PhD students

WHO AM I ?

Prof. Dr. Ir. Thomas Hermans

Associate Professor, Department of Geology (since 2017)

Hydrogeology and Applied Geophysics

​

Master in Civil Engineering, mining and geology (ULiege, 2010)

PhD in Applied Geophysics (ULiege, 2014)

​

Research interests: Geophysical inversion, geophysical monitoring, petrophysics, data integration in groundwater models, uncertainty quantification, geothermal energy

​

Current: 8 PhD, 4 postdocs, 2 technicians

​

​

​

​

​

​

http://hydrogeo.ugent.be

UNCERTAINTY IN THE SUBSURFACE: AN UNRESOLVED ISSUE ?

3

LET’S TAKE AN EXAMPLE

4

Calculating the wellhead protection area (WHPA)

THE 1^ST PROBLEM OF THE (HYDRO)GEOLOGIST

Our representation of the subsurface is based on (very) sparse local data

5

What happens between boreholes ?

​

How do I reproduce my observation ?

Well data

Lithology

THE TYPICAL SOLUTION TO CALCULATE THE WHPA ?

6

Data = predictor

Model

Prediction = target

Pumping and Tracing experiments

Distribution of subsurface properties

Inversion or calibration

Prediction

IS THIS APPROACH SATISFACTORY ?

7

Model

Target

To make inversion possible, model(s) are simplified

Parameterization

Boundary conditions

​

The solution is non-unique 🡪 Uncertainty

Models often fail to predict within the correct range

🡪 Uncertainty

Expected pollutant arrival time from model

Observed pollutant arrival time

THE 2^ND PROBLEM OF THE (HYDRO)GEOLOGIST

The inherent model uncertainty propagates in the model prediction 🡪 Decision making based on uncertain outcome

8

N different possible models = N different possible predictions

THE ALTERNATIVE TO THE DETERMINISTIC APPROACH

9

STOCHASTIC APPROACHES

Stochastic inversion…

• Generates an ensemble of possible solutions
• Produces geologically realistic models
• Provides an estimation of uncertainty

​

But stochastic inversion…

• Is extremely computationnally expensive (millions of iterations needed)
• Is tractable only for a limited number of parameters
• Is dependent on the definition of the prior
• Can be difficult to tune to convergence

10

CAN WE LEVERAGE MACHINE LEARNING?

Bayesian Simulation-Based Learning (BaSiL)

• Inversion is complex, but we are interested in the target
• (Machine) learning approach between predictor and target
• Problem: in the subsurface we miss repetitions
• Numerical models as input for learning

11

• Bayesian ?

The learning process should account for uncertainty

• Simulation-based ?

We use numerical models to map how the predictor influences the target.

• Learning ?

Every learning process requires “examples” to derive useful relationships

(Scheidt et al., 2018)

(Hermans et al., 2018, WRR)

(Thibaut et al., 2021, JoH)

Predictor

Models

Target

SOLVING THE WHPA WITH BASIL

12

DIMENSION REDUCTION

Model n samples

Predictor n samples

Target n samples

Statistical learning

Learning Phase

1. Dimension reduction with PCA

Corresponding reconstruction of the predictor (99%) and target (98%)

(Thibaut et al., 2021, JoH)

Solving GW flow and transport with numerical models

Solving GW flow and transport with numerical models

STATISTICAL LEARNING

Model n samples

Predictor n samples

Target n samples

Statistical learning

Learning Phase

2. Canonical Correlation Analysis

Set of independent bivariate linearized relaionships between predictor and target components

(Thibaut et al., 2021, JoH)

PREDICTING UNCERTAINTY

15

Prediction Phase

Observed data

Sampling

Prediction

3. Sampling in the CCA space + backtransformation

Linear regression or Kernel Density Estimation

(Thibaut et al., 2021, JoH)

WHPA PREDICTION FOR 4 DIFFERENT TRUTHS

16

(Thibaut et al., 2021, JoH)

HOW MANY MODELS ARE NEEDED FOR LEARNING

17

(Thibaut et al., 2021, JoH)

CONCLUSIONS

What are the advantages ?

• Only forward numerical models are used to generate the target and the predictor
• Once the learning phase is done, generating the posterior is cheap

🡪 No approximation needed for the posterior (e.g. ensemble smoother)

• No need to invert many data sets
• No approximation needed for the forward model (e.g., proxy-model)
• Number of prior samples is limited

🡪 Computationnally affordable

​

What are the limitations ?

• Consistency of the prior assumptions with the data
• Back to the 1^st problem of the (hydro)geologist
• Dimension reduction and learning with uncertainty

🡪 No guarantee of success

18

EXPERIMENTAL DESIGN WITH BASIL

19

STATING THE EXPERIMENTAL DESIGN PROBLEM

1. The most informative data set = the one reducing the most the prediction uncertainty (= data utility function)
2. The most informative data set depends on the subsurface uncertainty and the uncertain data sets
3. It is necessary to identify the most informative data sets for the whole range of possible data
4. This requires to stochastically inverse many data sets (not tractable)
5. This is commonly solved by approximating averaging techniques or proxy for the forward solver

20

BASIL FOR EXPERIMENTAL DESIGN

BaSiL is a perfect candidate to solve the ED problem

​

1. The statistical model is learned from the prior models
2. Once known this model can be used for any consistent (= within prior) data set
3. Prediction for a new dat set is straightforward and not expensive (sampling + back-transform)
4. It can be repeated as many time as needed (with data sets unseen during the prior) to estimate the posterior
5. The average data utility function is computed

​

!!! A new statistical model is required for any newly proposed data set !!!

21

WHICH WELL IS THE MOST INFORMATIVE

1. Prediction of WHPA from 400 training models
2. Data utility function = modified Hausdorff distance with the true prediction and Structural similarity index (SSIM)
3. Average over 250 unseen data sets (based on K-fold cross-validation)
4. Most informative well is deduced

22

WHICH WELL IS THE MOST INFORMATIVE

23

1. Prediction of WHPA from 400 training models
2. Data utility function = modified Hausdorff distance with the true prediction and Structural similarity index (SSIM)
3. Average over 250 unseen data sets (based on K-fold cross-validation)
4. Most informative well is deduced

​

🡪 Wells 4 and 6 are the most informative

TESTING COMBINATIONS OF WELLS

Which (combination of tests) is the most informative ?

We just need to predict for an independent (not used for training) set of models and assess uncertainty reduction

​

24

CONCLUSIONS

What are the advantages ?

• Once the learning phase is done, generating posterior is cheap
• Additional effort for ED is limited
• Possibility to test many data combinations
• No simplification needed to solve the ED problem
• ED problem is prior data acquisition
• No problem of consistency

​

What are the limitations ?

• Sequential optimization (when real data is acquired) requires model parameter estimation
• Validity of conclusion in other contexts ?

25

SOLVING FOR MODEL PARAMETERS AND GEOPHYSICAL INVERSION

26

IMAGING – MODEL PARAMETERS

27

(Michel et al., 2020, Comp. & Geos)

The target h = the model m

We replace inversion by statistical learning

(Michel et al., 2022, GJI)

GEOPHYSICAL INVERSION

28

BaSiL for 1D Geophysics

(Michel et al., 2020, Comp. & Geos)

Learning phase

(Michel et al., 2022, GJI)

GEOPHYSICAL INVERSION

29

BaSiL for 1D Geophysics

(Michel et al., 2020, Comp. & Geos)

Prediction phase

(Michel et al., 2022, GJI)

THE SNMR EXAMPLE

30

sNMR = surface Nuclear Magnetic Resonance

• Directly sensitive to the water content of the soil
• Relies on the response of water to electromagnetic perturbations (quantic level)

​

REDUCTION OF DIMENSIONALITY

31

PCA

STATISTICAL RELATIONSHIP

32

SAMPLING THE POSTERIOR

33

BACK-TRANSFORMATION

34

COMPARISON WITH MCMC

• Validated against DREAM (Vrugt, 2016)
• Uncertainty slightly oversestimated

​

• But McMC …
• difficult to tune to convergence (data representation)
• Computation time is large !!!

35

OVERESTIMATION OF THE UNCERTAINTY

36

BEL1D

DREAM

• No forward model run during prediction
• Prediction based on approximation in a lower dimensional space
• Some models do not fit the data
• Possiblity to filter models based on their RMS error

Noise level is 35 nV

OVERESTIMATION OF THE UNCERTAINTY

37

• Possibility to track back the “bad” models
• Origin related to models out of the trend
• Happens when prior uncertainty is large
• Is it possible to improve the approximation ?

Noise level is 35 nV

ITERATIVE PRIOR RESAMPLING

38

Concept

​

• Inspired by iterative spatial resampling (Mariethoz et al., 2010, WRR) and Sampling importance resampling (Dosne et al., 2016, J. Pharmacokinet phar)

​

• 3 phase process
• Adding the posterior samples to the prior
• Re-running learning operations
• Repeat iteratively until convergence

​

(Park and Caers, 2020, Comp. & Geos.;

Michel et al., 2022, GJI)

HOW MANY PRIOR SAMPLES DO WE NEED

• No improvement observed when more than 1000 samples are used

​

​

​

​

​

​

​

​

​

​

• Total CPU time ~1000 forward model runs + PCA + CCA

39

In this case, BaSiL is actually faster than deterministic inversion !

FIELD APPLICATION

Mont Rigi natural reserve

• Metric peat layer
• Cambrian bedrock

40

GPR estimated peat thickness map (Modified from Wastiaux and Schumacker, 2003)

• Transmitter = Receiver (20 m diameter)
• Noise estimated: 18 nV

FIELD APPLICATION: POSTERIOR MODELS

41

Results compared to QT-inversion (gray line)

​

Non-unicity of the solution illustrated by the RMS error range

​

Equivalence between thickness and water content of layer 1

​

Uncertainty on the relaxation time

​

ITERATIVE PRIOR RESAMPLING

42

Field case Mont Rigi

• Large uncertainty for the first layer

​

• Strongly reduced after IPR (6 iterations)

​

• All samples within noise level

CONCLUSION

• Efficient method to estimate the posterior distribution
• Combined with IPR: similar results as McMC (DREAM)
• Only relies on forward model runs
• Computationnally affordable
• No regularization
• Any noise model can be added (including correlated noise)

43

SEQUENTIAL OPTIMIZATION

44

STATING THE SEQUENTIAL ED PROBLEM

1. The most informative data set = the one reducing the most the prediction uncertainty (= data utility function)
2. Once the first data is known, a posterior distribution can be computed, including model parameters. Using the BaSiL framework, it requires to predict the target and the model parameters simultaneously
3. Update the learning samples to derive the next data to be collected and learn a new relationship
4. The process is repeated sequentially until uncertainty cannot be reduced further or budget is not available

45

OPTIMIZING TEMPERATURE MEASUREMENT

Assuming 1D vertical flow and advective-dispersive heat transport, it is possible to estimate vertical groundwater flux based on temperature measurement

​

It is using fixed boundary conditions and fixed values for thermal parameters

46

OPTIMIZING TEMPERATURE MEASUREMENT

If the set of parameters are uncertain, a stochastic approach must be used and a posterior distribution for the flux is obtained

​

The uncertainty can be reduced by choosing the optimized depth to record the temperature

​

Bayesian Optimal Experimental Design

1. Joint estimation of flux and subsurface parameters
2. Sequential optimization of sample depth

47

PROBABILISTIC BAYESIAN NEURAL NETWORK

Predictor = set of temperatures (The first 2 points are always measured)

​

Target = flux + subsurface parameters

​

Target = pdf (Gaussian mixture model)

48

RESULTS

Low fluxes are more difficult to predict

​

The model is not sensitive to the thermal conductivity

​

The bottom temperature can be accurately predicted (deep sampling)

49

CONCLUSION

• The optimum sequence is different for each test case
• Sequential design is needed

​

• On average, a shallow point is needed
• confirm the gradient close to the surface (flux)

​

• On average, a deep point is needed
• Estimate the bottom temperature

​

• 4 points are sufficient, measuring a 5^th points does not lead to significant uncertainty reduction

50

CONCLUSION AND PERSPECTIVES

Bayesian Simulation-Based Learning is an efficient framework for prediction, inversion and experimental design

• Accurate uncertainty estimation
• Basic statistical learning or advanced NN depending on the complexity
• Computationally more efficient than McMC
• Learning is based on numerical forward modeling of involved physics

​

Limitations

• Dependence on the prior model
• Consistency with field data (a model is always a simplification)
• Complex targets require more training samples and more advanced learning techniques

51

CONCLUSION AND PERSPECTIVES

Dimension reduction is the key to success

• Layered models (independent parameters)
• Smooth distribution (high spatio-temporal correlation)
• Parameterization of complex 1D-2D-3D model distributions ?

Deep Neural Network (Lopez et al., 2021, Comp. & Geos., 2022, JGR:SE)

​

52

Representation of complex structures (8325 pixels) with only 20 dimensions

​

Problem: highly non-linear mapping

True model

DNN representaion

CONCLUSION AND PERSPECTIVES

Dimension reduction is the key to success

• Layered models (independent parameters)
• Smooth distribution (high spatio-temporal correlation)
• Parameterization of complex 1D-2D-3D model distributions ?

Parameterization of geometries (grid independent) (e.g., Goebel et al., 2019, AGU Fall Meeting, H13A-04)

​

53

CONCLUSION AND PERSPECTIVES

Dimension reduction is the key to success

• Layered models (independent parameters)
• Smooth distribution (high spatio-temporal correlation)
• Parameterization of complex 1D-2D-3D model distributions ?

Global and local variables (Park and Caers, 2020, Comp. & Geos.; Oware et al., 2019, Geophysics)

​

54

Samples

Eigenimages

Hermans Thomas�Associate Professor��DEPARTMENT OF GEOLOGY��E thomas.hermans@ugent.be�T +32 9 264 46 60�M +32 499 13 88 53��www.ugent.be�

Ghent University�@ugent�Ghent University

​

PREDICTION RELATED TO GROUNDWATER SYSTEMS: ATES

56

Summer

Water table

Warm well 17°C

Cold well 7°C

Heat pump

25°C

12°C

The prediction of interest is the capacity of the aquifer to store heat

Approximation by the temperature in the hot well during seasonal cycles

storage

storage

recovery

recovery

HOW DO WE GET THIS PREDICTION ?

57

Data

Models

Prediction

“Standard Method”

Thermal response test, Push-pull test, Tracing experiment, pumping test etc.

storage

storage

recovery

recovery

Inversion/Calibration

Simulation/Forecasting

IS THIS APPROACH SATISFACTORY ?

58

Models

Prediction

Parameterization : zonation, layered model, simplification to reduce the number of unknowns, etc.

​

Choice of the boundary conditions, type of parameters (flow, heat transport, etc.)

​

Uncertainty ?

Calibration

Prediction

Real

Estimated

Models often fail to predict within the correct range

​

Uncertainty ?

EXAMPLE: HEAT STORAGE IN A SHALLOW AQUIFER

59

0

3

7

Coarse gravel

Sandy gravel

Loam

Injection well

Depth (m)

Prediction

Simulation of an ATES system

​

Temperature in the hot well during the recovery phase

(Hermans et al., WRR, 2018)

GENERATING PRIOR MODELS

60

Models

What do we know, what do we ignore ?

500 realizations = prior models

SENSITIVITY ANALYSIS OF THE PREDICTION

61

Prediction

Identification of the most sensitive parameters

IDENTIFICATION OF INFORMATIVE DATA SET(S)

62

Data

Designing an informative experiment

A heat tracing experiment monitored with ERT

Identification of transport processes and heterogeneity

SENSITIVITY ANALYSIS OF THE DATA

63

Identification of the most sensitive parameters

Overlapping with predictions ?

Data

CAN HEAT TRACING PREDICT HEAT STORAGE?

64

models from the prior field observations

Data d

Prediction h

1. Dimension reduction (PCA)

2. Linearization

(CCA)

Learning phase: Finding a direct relationship between data and prediction

CAN HEAT TRACING PREDICT HEAT STORAGE?

65

Prediction phase: Sampling the posterior in the low-dimension space

1. Modeling the relationship

​

• Sampling at the location of the observed data

Linear regression, kernel density, transport maps, etc.

ESTIMATING THE PREDICTION + UNCERTAINTY

66

Prediction phase: Back-transforming in the original space

IS IT WORKING ? VALIDATING WITH INDEPENDENT DATA

67

Data = 1 cycle push-pull test

Prediction = 2 cycles push-pull test

(Hermans et al., 2019, Hydrogeol. J.)

VALIDATION

We verify if we can predict the second experiment using the first one

68

TESTING POTENTIAL DATA SETS: 1-DAY VS 5-DAY

1-day = increasing part of the breakthrough curve

5-day = whole breakthrough curve

​

​

69

1-day

5-day

Same results, similar uncertainty

​

Is the 1-day experiment « sufficient »?

TESTING POTENTIAL DATA SETS: SINGLE OR MULTIPLE PUSH-PULL TESTS

70

Data 1

Data 2

Test 1

Test 2

CONCLUSIONS

What are the advantages ?

• No need for inverting model parameters
• Only forward modelling is needed for learning
• Computing cost limited to the forward model
• Full parallelization possible
• Limited number of samples is sufficient
• Decoupled learning and prediction phases
• Allows testing different data sets

​

What are the limitations ?

• Statistical learning requires dimension reduction and assumptions
• The prior must be consistent with the data (always the case)
• Does not directly provide insights on model parameters (at this stage)

71

A 2D + TIME CASE: TIME-LAPSE ERT

72

A heat tracing experiment

• Alluvial aquifer
• Heterogeneous with sand (top 4.5 m) and coarse gravel (bottom 2.5 m)
• ERT cross-section perpendicular to flow

(Hermans et al., 2016, 2018, WRR)

A 2D CASE: TIME-LAPSE ERT

73

(Hermans et al., 2016, 2018, WRR)

A 2D CASE: TIME-LAPSE ERT

74

Dimension reduction

• 99.5 % of the data in 50 first PCA dimensions, 95% of the variance in 14 dimensions for the prediction
• CCA
• Linear correlation is high
• Possibility to apply linear

Regression (Tarantola, 2005)

​

(Hermans et al., 2016, 2018, WRR)

A 2D CASE: TIME-LAPSE ERT

75

Prediction

​

(Hermans et al., 2016, 2018, WRR)

mean

StD

A 2D CASE: TIME-LAPSE ERT

76

Field application

• Same experimental lay-out
• Data sets filtered to reduce noise: 410 resistance, 6 time-steps (30 h)
• Prior based on fixed sand layer + sequential Gaussian simulations in gravel

(Hermans et al., 2016, 2018, WRR)

A 2D CASE: TIME-LAPSE ERT

77

Posterior

(Hermans et al., 2016, 2018, WRR)

A 2D CASE: TIME-LAPSE ERT

78

Validation

​

Observation from piezometers

• Flow dominated by the gravel
• Heat plume split in 2
• Higher temperature on the left side

​

• Validation of T at two locations (T loggers)
• Large uncertainty in the middle of the section

(Hermans et al., 2016, 2018, WRR)

A 3D CASE: TIME-LAPSE ERT

79

A 3D CASE: TIME-LAPSE ERT

80

A 3D CASE: TIME-LAPSE ERT

81

Dimension reduction

​

• We keep 30 PCA dimensions (95% of the data, 90% of the prediction)
• CCA
• Linear correlation is high
• But distributions are not

Gaussian

• Use of KDE

​

A 3D CASE: TIME-LAPSE ERT

82

Prediction

​

• The median sample is very close to the true distribution

​

• BEL captures both the spatial and temporal behavior of the heat

​

• The amplitude is much closer than with deterministic inversion

​

A 3D CASE: TIME-LAPSE ERT

83

Validation

• BEL perfectly captures the evolution of the average temperature
• The relative error of estimation of T is only significant where the change of T is low, i.e. where the absolute error is low

​

PRIOR UNCERTAINTY

84

Description of the prior model space

2 layer model

… AND SYNTHETIC DATA MODELLING

85

Single transmitter/receiver loop

50 m of diameter

2 turns

HOW DO WE DEAL WITH NOISE ?

• Learning is done on synthetic noise free data
• But field data are noisy !
• The impact on PCA scores can be estimated and translated into a covariance matrix in CCA

86

Noise propagation

Uncertainty on d

Increase of uncertainty

Controlled via KDE