MARWAN FAHS § , QIAN SHAO ¶ Thierry A Mara † ‡
email: thierry.mara@ec.europa.eu

Anis Younes

RANDOM SAMPLING FROM JOINT PROBABILITY DISTRIBUTIONS DEFINED IN A BAYESIAN FRAMEWORK *

Keywords: Bayesian framework, model statistical calibration, maximal conditional posterior distribution, posterior uncertainty and sensitivity analyses, numerical drainage experiment AMS subject classifications. 65C60, 62H20

Random variables characterized by a joint probability distribution function (jpdf) defined in a Bayesian framework are generally sampled with Markov chain Monte Carlo (MCMC). The latter can be computationally demanding when the number of variables is high. As an alternative, the maximal conditional probability distribution (MCPD) sampler was recently introduced by some of the authors of the present article to readily and efficiently draw values randomly sampled from the desired jpdf. The MCPD approach provides the probability distribution of a given variable under the condition that the other variables maximized the conditional jpdf. However, contrarily to MCMC, MCPD does not provide enough draws to allow posterior uncertainty and sensitivity analyses of the computer model responses. In the present work, we show how to draw random samples from the MCPD draws under the requirement that the target jpdf possesses a particular dependence structure. Several numerical tests are carried out to prove the efficiency of the new sampling method. The new approach is used to perform the predictive uncertainty and sensitivity analyses of numerical models posterior to their statistical calibration from experimental data.

1. Introduction. Sampling random variables from a given joint probability distribution function (jpdf) is a challenging issue if the latter is not defined in a closedform. This is the case when modellers wish to statistically calibrate the input variables of their computer model in a Bayesian framework (e.g. [START_REF] Campbell | Statistical calibration of computer simulations[END_REF][START_REF] Bayarri | A framework for validation of computer models[END_REF]). Sampling from the jpdf is necessary to perform, for instance, the predictive uncertainty and sensitivity analyses of the computer model responses (e.g. [START_REF] Saltelli | Sensitivity analysis[END_REF][START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF][START_REF] Mara | Non-parametric methods for global sensitivity analysis of model output with dependent inputs[END_REF]. Many authors privilege sampling the input values from a given jpdf with Markov chain Monte Carlo (MCMC) [START_REF] Vrugt | Validity of first-order approximations to describe parameter uncertainty in soil hydrologic models[END_REF][START_REF] Qian | Understanding predictive uncertainty in hydrology modeling: The challenge of identifying input and structural errors[END_REF][START_REF] Marzouk | A stochastic collocation approach to Bayesian inference in inverse problems[END_REF][START_REF] Schoups | A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic and non-Gaussian errors[END_REF]34,[START_REF] Laloy | Efficient posterior exploration of a high-dimensional groundwater model from two-stage markov chain Monte Carlo simulation and polynomial chaos expansion[END_REF].

MCMC relies on a rejection/acceptance sampling to generate random draws sampled from the jpdf [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. It has been subject to several developments and improvements during the last two decades aiming at accelerating its convergence (e.g. [START_REF] Grenander | Representations of knowledge in complex system[END_REF][START_REF] Haario | An adaptive Metropolis algorithm[END_REF][START_REF] Green | Delayed rejection in reversible jump Metropolis-Hastings[END_REF][START_REF] Haario | DRAM: Efficient adaptive MCMC[END_REF][START_REF] Braak | Differential evolution markov chain with snooker updater and fewer chains[END_REF][START_REF] Dostert | Efficient uncertainty quantification techniques in inverse problems for Richards' equation using coarse-scale simulation models[END_REF] among others). However, MCMC sampling is challenging in high-dimensions and remains computationally demanding.

Recently, a new approach has been introduced for the statistical calibration of computer models [START_REF] Mara | Inversion and uncertainty of highly parameterized models in a Bayesian framework by sampling the maximal conditional posterior distribution of parameters[END_REF]. The authors named it the maximal conditional probability distribution sampler. The first step of this approach is to seek all the probable local optima of the jpdf (assuming the existence of a finite number of local optima). Then, several maximizations of the conditional jpdf are performed for different prescribed values of one selected variable. The values assigned to the selected variable are successively drawn around each local optimum while the values of the other variables are investigated by maximizing the conditional jpdf. This provides what is called the maximal conditional probability distribution (MCPD) of the selected variable.

The evaluation of the MCPDs is independent of each other. Consequently, they can be evaluated simultaneously by distributing the calculations over several computers (or a multi-core computer). This feature drastically decreases the computation time and makes the inversion of highly parameterized problems feasible. For instance, in [START_REF] Mara | Inversion and uncertainty of highly parameterized models in a Bayesian framework by sampling the maximal conditional posterior distribution of parameters[END_REF], a flow model with d = 104 parameters has been calibrated with the MCPD approach by distributing the calculations over eight cores (each core assessing 13 MCPDs). An overall of 54 560 model calls was necessary to assess the MCPDs, but thanks to the parallelization the actual waiting time (also called computational time unit) was about 6 500 model calls. Several comparisons between MCPD and a MCMC sampler have been carried out in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF] which demonstrated the efficiency of the MCPD approach.

However, the MCPD draws represent a few probabilistic set of values sampled from the jpdf. If stochastic samples distributed over the target jpdf are desired (for the purpose of uncertainty and sensitivity analyses of computer models for instance), the MCPD draws are not sufficient. In the present work, we extend the approach in order to generate Monte Carlo samples from the MCPD draws. This extension is based on the orthogonalization procedure introduced in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] which assumes a particular correlation structure between the variables. The new sampling method is tested on two numerical target distributions of increasing complexity before application to the identification of soil hydraulic parameters from a synthetic multi-step outflow experiment and to the calibration of a reactive transport model from field observations. In these last two applications, the use of Monte Carlo samples for uncertainty and sensitivity analyses is illustrated.

The paper is organized as follows: First, we briefly discuss the different methods to generate random samples from a given jpdf in Section 2. Then, we recall the concept of MCPD sampling and describe succinctly its assessment in Section 3. The approach for sampling random variables from the desired distribution is explained in Section 4. The new approach is evaluated through several numerical exercises of different dimensions and complexities in Section 5. In Section 6 we apply the new sampling approach to the calibration of a drainage model from synthetic multi-step outflow experiment. Finally, the statistical calibration of a reactive transport model from field data is undertaken in Section 7 before concluding in Section 8.

2. Sampling from joint distribution: Existing approaches. Let x = (x 1 , . . . , x d ) T be a vector of random variables distributed over the joint probability density function p(x), where the superscript T denotes the transpose operator. We are asked to generate x, a sample of size N , of the random vector. The most straightforward approach to perform this task is to use the inverse Rosenblatt transformation [START_REF] Rosenblatt | Remarks on the multivariate transformation[END_REF], that is,

                   x 1 = F -1 1 (u 1 ) x 2 = F -1 2|1 (u 2 |x 1 ) x 3 = F -1 3|1,2 (u 3 |x 1 , x 2 ) . . . x d = F -1 d|1,...,d-1 (u d |x 1 , . . . , x d-1 ) (2.1)
where u is a random vector uniformly and independently distributed within the unit hypercube [0, 1] d and F -1 i|1,...,i-1 is the inverse cumulative distribution of x i conditioned onto (x 1 , . . . , x i-1 ). By sampling u ∈ [0, 1] d and knowing the set of inverse conditional cumulative distribution functions (CDF) {F -1 1 , F -1 2|1 , . . . , F -1 d|1,...,d-1 }, the desired sample of x is obtained by applying Eq. (2.1). Unfortunately, in many situations and especially in model statistical calibration, the conditional cumulative distribution functions are unknown.

More often, random variables are defined by their unconditional marginal CDFs {F 1 (x 1 ), . . . , F d (x d )} and a copula density c : [0, 1] d → [0, +∞[. The copula density contains the dependence/correlation structure of the jpdf. It is convenient to define random variables in this way because many algorithms are proposed in the literature to generate random samples [START_REF] Nelsen | An introduction to copulas[END_REF]. For instance, the inverse Nataf transform is designed to cope with Gaussian copulas [START_REF] Nataf | Détermination des distributions dont les marges sont données[END_REF]. We note that defining random variables from copula theory is tantamount to define the jpdf in a closed-form because of the following equality:

p(x) = c(F 1 (x 1 ), . . . , F d (x d ))p 1 (x 1 ) . . . p d (x d ), where p i (x i ) = dF i (x i )/dx i
is the marginal pdf of x i . However, there are many problems for which such an explicit definition is not possible as when, for model calibration purposes, the target jpdf is derived in a Bayesian framework. In such a situation, Markov chain Monte Carlo samplers are preferred.

MCMC samplers are based on acceptance-rejection algorithms. They employed the Metropolis-Hastings algorithm (or a variant of it) to sample x from p(x) thanks to a given proposal distribution q(•|•). There are several MCMC variants proposed in the literature. Amongst them, we can cite the Gibbs sampler [START_REF] Gelfand | Sampling-based approaches to calculating marginal densities[END_REF] which requires the knowledge of the conditional distributions p i|∼i (x i |x ∼i ), ∀i = 1, . . . , d (x ∼i stands for all the x-variables except x i ). When the conditional distributions are unknown, the Langevin [START_REF] Grenander | Representations of knowledge in complex system[END_REF] or Hamiltonian [START_REF] Hoffman | The No-U-Turn sampler: Adaptatively setting path lengths in Hamiltonian Monte Carlo[END_REF] MCMC sampler can be a good choice to generate x. Notably, the efficiency of these samplers is enhanced if the Jacobian of log (p(x)) is provided. When the latter is not available or not differentiable, the DiffeRential Evolution Adaptive Metropolis sampler is a good alternative, in particular the version exploiting the archive of past states called DREAM ZS [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (zs) and high-performance computing[END_REF][START_REF] Vrugt | Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation[END_REF]. DREAM ZS is particularly efficient, since instead of drawing new candidates from a prescribed proposal distribution q, new candidates are drawn thanks to an archive of past states.

Despite of the obvious improvements in MCMC sampling, the use of this method to generate a sample of x ∼ p(x) can still be computationally demanding. This is particularly embarrassing when the computational time required to evaluate the target jpdf is high. In this case, the two stages surrogate-based MCMC approach of [START_REF] Cui | Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis-Hastings algorithm[END_REF] seems a promising alternative.

The aim of the present work is to introduce a new fast approach to generate random draws from the desired jpdf based on the MCPD sample. This work capitalizes upon recent works of some of the authors of the present article [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF][START_REF] Mara | Inversion and uncertainty of highly parameterized models in a Bayesian framework by sampling the maximal conditional posterior distribution of parameters[END_REF]. The approach, described in the next sections, is particularly efficient if the Jacobian of log (p(x)) is provided. It also requires a particular correlation structure between the random variables. Although these constraints narrow its field of applications, it is believed that the proposed approach remains a valuable tool in many model calibration problems.

3. The maximal conditional posterior distribution. The maximal conditional posterior distribution of x i writes

P i (x i ) = max x-i (p(x -i |x i )) × p(x i ).
(3.1)

P i (x i
) is interpreted as the posterior probability function that maximizes the conditional posterior distribution p(x -i |x i ). The MCPD of x i is assessed in a discrete form by evaluating Eq. (3.1) at different values of x i . In practice, the sampled values of x i (denoted x i * ) are picked around each probable local optimum (estimated beforehand) within its prior uncertainty range. This gives,

x * -i = argmax x-i {p(x -i |x i = x i * )} (3.
2)

P i (x i * ) = p(x * -i |x i * ) × p(x i * ) = p(x * ) (3.3)
The MCPD approach makes sense provided that the target jpdf p(x) admits a finite number of modes. Notably, the success of MCPD sampling relies on the ability to retrieve all the probable local optima of p(x). In the present work, all the local optima are searched by use of gradient-based methods with multiple starting points. Multiple tries ensure that the algorithm starts searching the local optima from different regions of the input space. Of course, the selection of the optimization technique is a matter of choice.

The algorithm used to compute the MCPDs is thoroughly explained in [START_REF] Mara | Inversion and uncertainty of highly parameterized models in a Bayesian framework by sampling the maximal conditional posterior distribution of parameters[END_REF] and [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF]. For convenience, it is succinctly recalled here. The algorithm is divided into three parts: in part 1, all the probable optima of p(x) are investigated. In part 2, the posterior range of variation of each parameter around each probable optimum is roughly estimated by preliminarily evaluating its MCPD. In part 3, the evaluation of the MCPDs is refined. It is worth reminding that fast computation is possible because these three parts can take advantage of parallel computing.

Let us denote the MCPD draws in the vicinity of the local optimum x opt,m by x ki,m d i=1 , ∀k i = 1, . . . , N i . We set p(x opt,1 ) ≥ p(x opt,2 ) • • • ≥ p(x opt,m ) which means that x opt,1 is the global optimum. The subset x ki,m , P ki,m i = p(x ki,m )/p(x opt,1 ) Ni ki=1 ∀m = 1, . . . , M represents the discretized MCPD of variable x i . As defined, the discretized MCPD representation is scaled within [0, 1]. As explained above, the MCPD draws of x i are estimated by setting x i * = x i ki,m and solving Eq. (3.2). The value x i ki,m is successively drawn in the vicinity of the current optimal value of x i . Figure 3.1 depicts an example of the MCPD draws for a unimodal target distribution. On the diagonal, the discretized MCPDs are plotted in blue crosses. On the lower off-diagonal, at row #j and column #i, the pairs (x i ki,m , x j ki,m ) and (x i kj ,m , x j kj ,m ) are plotted. The first pair is obtained while P i is assessed whereas the second one stems from the assessment of P j . This scatterplot shows possible correlations between the x i and x j draws. For instance, x 1 and x 2 show no correlation (row #2, column #1) because virtually two orthogonal curves are observed while the MCPD draws of x 2 and x 3 are correlated (row #3, column #2).

In [START_REF] Mara | Non-parametric methods for global sensitivity analysis of model output with dependent inputs[END_REF], it was demonstrated that the MCPD draws were (by far) faster to sample than the MCMC draws. Nevertheless, the two samples are not comparable because the probabilistic MCPD draws are sampled with constraints (maximization of the conditional jpdf, see Eq. (3.2)) while the stochastic MCMC draws are obtained from the same jpdf but unconditionally. Therefore, posterior uncertainty and sensitivity analyses of computer models cannot be performed with the MCPD draws. For this purpose, a Monte Carlo (MC) sample that covers the posterior input space would be preferable. Hence, our aim in this work is to propose an approach to generate MC samples from the MCPD draws. To this end, an approach that can impose a desired correlation structure amongst independent MC samples is needed. This is described in the next section. (4.2). Subplots on the diagonal display the MCPD draws of the original variable x i and its transformed one xi (blue crosses and red squares resp.). The lower off-diagonal plots show the pairwise correlations between the original variables. The draws decorrelated with the transformation Eq. ( 4.2) are depicted on the upper off-diagonal. [START_REF] Cui | Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis-Hastings algorithm[END_REF]. The new sampling approach. In this section we describe the algorithm to generate stochastic Monte Carlo samples from a set of probabilistic MCPD sample. The success of the method heavily relies on the decorrelation procedure proposed in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] . Here we recall the procedure in § 4.1. In the case of multimodal target jpdf, a balanced sampling scheme is necessary which is described in § 4.2. We discuss the validity of the approach and a possible alternative in case of failure in § 4.3. The algorithm of the proposed approach, called MCPD-MC, is given in § 4.4.

Decorrelation of the MCPD draws.

Assume that the vector of dependent random variables x = (x 1 , . . . , x d ) T is function of a set of independent latent random variables x = (x 1 , . . . , xd ) T . Let us further assume that the dependence structure of the x-variables is strictly due to their relationships with the latent variables which are of the form,

                         x i1 = x1 x i2 = x2 + f i21 (x 1 ) x i3 = x3 + f i31 (x 1 ) + f i32 (x 2 ) . . . x i d = xd + d-1 j=1 f i d j (x j ) (4.1)
where, in the j-th equation, R f ij k (x k )p(x k )dx k = 0, ∀k = 1, . . . , j -1 which ensures that the set of functions {f ij 1 (x 1 ), . . . , f ij j-1 (x j-1 )} is orthogonal (according to [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]). We note that each equation in Eq. (4.1) is related to the generalized additive model representation [START_REF] Hastie | Generalized additive models[END_REF].

It is worth noticing that Eq. (4.1) may not always be possible for any set {i 1 , . . . , i d } = {1, . . . , d}, for instance, when the relationship between two variables is non-monotonic (see an example in § 5.2). However, we assume that it is at least possible for one of the sets (as for the non-monotonic case).

Proposition 4.1. Under assumption (4.1), the transformation procedure proposed in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] that writes

             x1 = x i1 x2 = x i2 - E [(x i2 -E [x i2 ]) |x 1 ] x3 = x i3 - 2 j=1 E [(x i3 -E [x i3 ]) |x j ] . . . xd = x i d - d-1 j=1 E [(x i d -E [x i d ]) |x j ] (4.2)
where E [•] is the mathematical expectation and E [•|•] is the conditional expectation, provides a set of independent variables x.

Eq. (4.2) slightly differs from the original set of equations in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] because in their paper it was assumed that x was a vector of standardized random variables (i.e. centered and reduced). The modified version here considers the conditional expectations, x ij -E x ij j = 2, . . . , d, which are centered. Indeed, the reduction of the x-variables (V [x i ] = 1) is not mandatory for the independence of the x-variables. The proposed approach to generate random samples w.r.t. the target jpdf, first, requires to apply transformation (4.2) to the MCPD draws. Numerically this is achieved by approximating

E [(x i -E [x i ]) |x j ] = f ij (x j
) with a regression method. In the present work, polynomial regressions are used in conjunction with Schwartz's criterion for the automatic selection of polynomial degree [START_REF] Schwartz | Estimating the dimensions of a model[END_REF].

Applying the previous transformation to the set x ki,m Ni ki=1 and ∀i = 1, . . . , d, which is the MCPD sample drawn in the vicinity of the local optimum x opt,m , yields the independent sample xki,m Ni ki=1 , the pairs xki,m i , Pi (x ki,m i ) as well as the overall set of f ij functions around the optimum. An example of such a transformation is depicted in Figure 3.1 (red squares). One can note that the distribution of x1 is the same as the one of x 1 (row #1, column #1). The one of x2 is similar to the one of x 2 except that it is centered (row #2, column #2). P3 (x 3 ) is centered and narrower as compared to P 3 (x 3 ) (row #3, column #3). The shrinkage of x3 's distribution is due to the correlation between x 3 and x 2 . Pairwise analysis shows no correlation amongst the transformed draws (upper off-diagonal plots) contrarily to the MCPD draws depicted on the lower off-diagonal which shows a strong correlation between x 2 and x 3 (row #3, column #2). 4.2. Sampling from the target distribution. Given the estimated set of f ij functions and the set of transformed MCPD Pi , we are now ready to propose a new approach to generate x from the joint probability distribution without further evaluating p(x). First, given that x is a vector of independent random variables in Eq. ( 4.1), we start by drawing the x-variables from the Pi (x i )'s. For this purpose, we use the latin hypercube sampling [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] that requires the individual cumulative distribution functions (CDF) around each optimum. The CDF of xi around the m-th optimum is defined as follows,

A m = upi loi Pi (x i )dx i F m i ( Xi ) = 1 A m Xi loi Pi (x i )dx i , ∀i = 1, . . . , d (4.3)
with

lo i = min{x 1,m i , . . . , xNi,m i }, up i = max{x 1,m i , . . . , xNi,m i } and Xi ∈ [lo i , up i ]. Numerically, F m i ( Xi )
is estimated with the Simpson's quadrature rule for numerical integration. Let us denote {x k,m } nm k=1 the independent latin hypercube sample of size n m . Then, the desired sample x of the original random variables are simply obtained by transforming the previous draws with Eq. (4.1). The f ij functions are replaced by their polynomial approximations already computed in the decorrelation procedure. By repeating this procedure around all optima, one gets N = M m=1 n m draws of the x-variables sampled from the jpdf that was evaluated to get the MCPD draws. In the sequel, MCPD-MC refers to the MCPD-based Monte Carlo approach described in this section.

It has to be noted that, if a final sample of size N is sought, the sample size n m for each optimum must respect the balance area of the modes, that is,

n m N = A m M k=1 A k . ( 4.4)
This ensures that each optimum is sampled in good proportion. For the sake of completeness, we show in Figure 4.1 the Monte Carlo draws obtained with the proposed approach applied to the MCPD draws depicted in Figure 3.1. Note that the MCPD-MC draws are randomly sampled with the desired correlation structure.

4.3. Discussion. Contrarily to MCMC, the MCPD sampler relies on the use of an optimizer. The computational effort of the algorithm resides in the optimization steps. In the numerical exercises below, the optimizations are performed with a gradient-based algorithm. The convergence of these algorithms is accelerated if the Jacobian matrix/vector of the target distribution is also provided. In our studies, the latter is systematically computed, unless the contrary is mentioned. As we underline above, the optimization processes can be distributed over several computers (or parallel sessions). This advantage is also systematically exploited in the current studies.

In practice, it is not known beforehand whether the correlation structure of the input variables satisfies Eq. (4.1) (the strong assumption of the method). Thereby, it is recommended to evaluate the target distribution with the generated sample. Our experience suggests that, if for several candidates {x k } N k=1 of the MCPD-MC sample, the target distribution produces the value zero, then it is likely that the assumption is not valid.

Alternatively, if evidence yields to reject assumption Eq. ( 4.1), one can adopt a less restricted assumption by considering the following orthogonalization procedure,

             x1 = x i1 x2 = x i2 - E [(x i2 -E [x i2 ]) |x 1 ] x3 = x i3 -E [(x i3 -E [x i3 ]) |x 1 , x2 ] . . . xd = x i d -E [(x i d -E [x i d ]) | x∼d ] (4.5)
which relies on a dependence structure of the form,

                         x i1 = x1 x i2 = x2 + f i21 (x 1 ) x i3 = x3 + f i31 (x 1 ) + f i32 (x 2 ) + f i312 (x 1 , x2 ) . . . x i d = xd + d-1 j1=1 f i d j1 (x j1 ) + d-1 j2>j1 f i d j1j2 (x j1 , xj2 ) + • • • + f i d 12...d-1 (x 1 , . . . , xd -1 )
(4.6) with still the orthogonality constraints on f i k j1j2...js (x j1 , . . . , xjs ) with {j 1 , . . . , j s } ⊆ {1, . . . , k-1} in each equation to ensure its uniqueness. For the interested readers, the expansions in Eq. (4.6) are related to the analysis of variance decomposition [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The difficulty in this case is to approximate the multidimensional functions in the different equations. This can be achieved, for instance, with the non-parametric polynomial chaos expansion method of [START_REF] Shao | Bayesian sparse polynomial chaos expansion for global sensitivity analysis[END_REF]. Adopting (4.6) has the drawback to slow down the MCPD-MC method as one has to build (d -2) expansions (needless to perform such a complex expansion for x i1 and x i2 ). This alternative approach is not considered in the numerical examples treated in § 5, § 6 and § 7 as assumption (4.1) is fulfilled.

The upper off-diagonal plots show no correlation amongst the draws (red dots). The lower diagonal plots show the desired draws sampled from the target distributions (gray dots).

Note that step 4 requires to evaluate the jpdf, and thus the likelihood function. This step might be computationally demanding. It is recommended to start with a small sample size N (say N = 50) and if the test is successful (i.e. n/N > 1%), then increase the Monte Carlo sample size (which is computationally cheap). In all the following exercises, we find n = 0 and therefore this issue is not discussed much.

Numerical exercises.

5.1. Eleven-dimensional multimodal target distribution. We start with the following multimodal target distribution,

p(x|µ 1 , µ 2 , µ 3 , C) = 1 6 N (µ 1 , 5C) + 2 6 N (µ 2 , 5I d ) + 3 6 N (µ 3 , 5I d ) (5.1)
where The case d = 11 was studied in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF] in which the authors demonstrated that the MCPD draws were faster to generate than the MCMC draws (the case d = 25 was considered in [START_REF] Mara | Inversion and uncertainty of highly parameterized models in a Bayesian framework by sampling the maximal conditional posterior distribution of parameters[END_REF]). Indeed, the MCPD sampler only required around one thousand evaluations of p(x|µ 1 , µ 2 , µ 3 , C) to obtain an accurate estimate of the variables' distributions. Thanks to the parallel computation, the computational time units (CTU) were about 190. This means that actually the waiting time corresponded to 190 model calls. This is very short compared to the CTU of the DREAM ZS sampler which was 10 000 (with 11 chains in parallel). The obtained MCPD sample is considered hereafter for generating the MCPD-MC sample.

N (µ i , 5I
Eq. (5.1) has been assessed with the MCPD-MC draws. The results are depicted in Figure 5.1 for a sample of size N = 4 096. On the diagonal, the MCPD estimate of x i (row #i, column # i) and the Monte Carlo draws (x i k,m , p(x k,m )) , ∀k = 1, . . . , N and m = 1, 2, 3 are reported. Note that the MCPD-MC draws are located below the MCPD curves. This is because the MCPD-MC approach samples probable solutions while the MCPD sampler draws the most probable variable set for prescribed values of one of them.

The lower off-diagonal scatterplots show the pairwise correlations of MCPD-MC draws. The estimated Pearson correlation coefficients for the first mode are also reported. They are close to the analytical correlation coefficients, r

12 ≈ C 1,2 = -0.5 and r

(1) 13 ≈ C 1,3 = 0.8. The upper off-diagonal plots compare the empirical densities with the analytical marginal densities. It can be inferred that each mode has been sampled in good proportion.

Ten-dimensional twisted Gaussian target distribution.

Let us now consider another challenging problem which was also analyzed in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF]. We target the twisted Gaussian density proposed in [START_REF] Haario | Adaptive proposal distribution for random walk Metropolis algorithm[END_REF] and defined as follows, p(x) = p(x 1 )p(x 2 |x 1 ) 10 i=3 p(x i ) with p(x 1 ) = N (0, 100), p(x 2 |x 1 ) = N (-0.1x 2 1 + 10, 1) and p(x i ) = N (0, 1), ∀i = 3, . . . , 10. This jpdf is very challenging because it exhibits a non-monotonic relationship between x 2 and x 1 . Indeed, the conditional expectation of x 2 on x 1 is E [x 2 |x 1 ] = -0.1x 2 1 + 10. Generating the MCPD draws required approximately 1 900 evaluations of p(x) while DREAM ZS required at least 20 000 evaluations to converge (see the aforementioned paper for more details). We can notice that the relationship between E [x 1 |x 2 ] cannot be written in the form of f 12 (x 2 ) since for one value of x 2 two values of x 1 are possible. Hence, Eq. (4.1) is not valid for the variable set (x 2 , x 1 , . . . , x 10 ) so ordered, but it is for the set (x 1 , x 2 , x 3 , . . . , x 10 ). Indeed, it is possible to write E [x 2 |x 1 ] = f 21 (x 1 ) with f 21 the non-monotonic function mentioned above. Consequently, the transformation in Eq. (4.2) starts with

x1 = x 1 , then x2 = x 2 -E [(x 2 -E [x 2 ]) |x 1 ], x3 = x 3 -E [(x 3 -E [x 3 ]) |x 1 ] -E [(x 3 -E [x 3 ]) |x 2 ]
and so on.

The results are depicted in Figure 5.2 for the first three variables. Note that, on the first diagonal plot (row #1, column #1) the MCPD and MCPD-MC draws of x 1 are depicted, while on row #2 and column #2 the draws of x 2 are represented. The MCPD-MC draws are located beneath the MCPD curves. We note that the uncertainty ranges of x 1 and x 2 are very large which is not an issue for the MCPD sampler. The upper off-diagonal plots prove that the MCPD draws have been successfully orthogonalized (they are actually independent). The lower off-diagonal scatterplot of

The on-diagonal plots represent the MCPDs estimate (blue line) as well as the MCPD-MC draws. The lower off-diagonal scatterplots represent the pairwise correlation of the MCPD draws (blue crosses) and the MCPD-MC draws (gray dots). r

(1) ij is the estimated pairwise correlation coefficient of x i and x j in the vicinity of the first mode. The upper off-diagonal confirms that the MCPD-MC sample of each variable fits its analytical marginal density (red curve).

x 2 versus x 1 (row #2, column # 1) shows that, because of their strong dependence, the MCPD-MC draws (gray dots) are located close to the MCPD draws (blue dots).

6. Calibration of a drainage model. [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF].1. The model and the dataset. The MCPD-MC approach proposed in the present work has been specifically developed to address the issue of model calibration in a Bayesian framework. In this section, we consider the calibration of a soil drainage model at the laboratory scale. Once again, this problem was also studied in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF] to compare the performance of the MCPD approach with the one of the MCMC sampler. Here, we repeat the same numerical experiment in order to illustrate the MCPD-MC approach.

We model a laboratory multistep outflow drainage experiment in which a column of length L = 6 cm and diameter D = 8.5 cm is filled with sand and initially saturated with water. The column is drained by imposing at the lower boundary of the column multistep prescribed negative pressure heads. Flow in the column is modelled by the non-linear one-dimensional Richard's equation,

∂θ ∂t = ∂ ∂z K(h) ∂h ∂z -1 , (6.1)
in conjunction with the Mualem-van Genuchten (MvG) retention curve [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF][START_REF] Van Genuchten | A closed form equation for predicting the hydraulic properties of unsaturated soils[END_REF],

K(S e ) = k s • S λ e 1 -1 -S 1/m e m 2 , (6.2)
where t (min) is time, z (cm) is the vertical coordinate (positive downward) and m = 1 -1/n. The water content θ (cm 3 .cm -3 ) and the pressure head h (cm -1 ) are the state variables. The effective saturation S e is defined as follows,

S e = θ -θ r θ s -θ r =    1 |αh| n h < -1/α 1 h ≥ -1/α . (6.3)
and K (cm 3 .min -1 ) is the unsaturated hydraulic conductivity. The soil hydraulic parameters to be estimated are: k s (cm 3 .min -1 ) the saturated hydraulic conductivity, θ s (cm 3 .cm -3 ) the saturated water content, θ r (cm 3 .cm -3 ) the residual water content and the MvG fitting coefficients α (cm -1 ), n (-) and λ (-).

Eqs.(6.1-6.3) are solved with a standard Galerkin finite element method in conjunction with the Newton linearization method. Synthetic data were obtained by running the numerical model for a given input parameter set and noising the model responses with independent Gaussian random noises. The responses of interest are the pressure head h and the soil water content θ at the center of the column. The inverse modelling is performed from the observed pressure head and water content denoted by y h and y θ respectively.

Setting the calibration problem in a Bayesian framework yields the following posterior jpdf,

p(x, σ h , σ θ |y h , y θ ) ∝ 1 σ N h h σ N θ θ exp - 1 2 SS h (x) σ 2 h + SS θ (x) σ 2 θ , (6.4)
where N h and N θ are the number of observed data of pressure and water content respectively and SS h and SS θ are the sum of square errors of pressure head and water content respectively. This posterior jpdf was derived by assuming normal errors for h and θ and independent uniform priors for the hydraulic parameters within large plausible ranges (see [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF] for more details). Eight unknowns were sought in the inverse problem, including the vector of hydraulic parameters x = (k s , θ r , θ s , α, n, λ) as well as the error variances σ 2 h (cm 2 ), σ 2 θ (cm 6 .cm -6 ). The different maximization processes were performed with the Levenberg-Marquardt algorithm [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] as the partial derivatives of the model responses w.r.t. the hydraulic parameters were computed by the numerical model. The performance of the MCPD sampler was discussed in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF] and a comparison with the DREAM ZS MCMC sampler of [START_REF] Laloy | Efficient posterior exploration of a high-dimensional groundwater model from two-stage markov chain Monte Carlo simulation and polynomial chaos expansion[END_REF] was also carried out. We recall that the MCPD sampler required an overall of 7 500 model calls. But thanks to the parallelization, the computational time unit (related to the real waiting time) was about 2 000 model calls which corresponded to the assessment of θ s 's MCPD. An MCPD sample of size 185 was obtained. DREAM ZS required 64 000 model calls but with eight chains in parallel the computational time unit was 8 000. In the next subsections, we discuss the MCPD-MC sampling of the calibrated model parameters and the posterior predictive uncertainty of the model responses. Of particular interest is the posterior uncertainty of the predicted cumulative outflow at the end of the simulated experiment. In effect, the cumulative outflow was not considered as a measurement in the calibration problem.

Parameter uncertainty quantification.

The MCPDs of the hydraulic parameters as well as those of the likelihood hyperparameters are depicted in Figure 6.1 (the on-diagonal plots). They show bell-shaped parameters' posterior marginal pdf with a clear optimal value. The support of the MCPDs are quite narrow except for the MvG parameter λ. This might be explained by a lack of sensitivity of h(t) and θ(t) to this parameter. From the 185 MCPD draws, a MCPD-MC sample of size 512 was generated. They are plotted in Figure 6.1 (the scatterplots).

The MCPD-MC draws were then propagated through the model (Eq. (6.1)) in order to estimate the probability assigned to each draw (Eq. (6.4)). They are plotted with the MCPD curves (the on-diagonal scatterplots in Figure 6.1). Note that after evaluating Eq. (6.4), the probability value of each draw was scaled between [0,1] by dividing by the probability value of the MAP. The MCPD-MC draws are assigned high probability values, proving that they stem from the region of high probability of the posterior parameter space. Actually, none of the MCPD-MC draw was assigned a probability zero which is an indication of the reliability of the generated draws as explained in Step 4 of the algorithm defined in § 4.4.

Pairwise scatterplots of the MCPD-MC draws are depicted on the lower offdiagonal. They show either strong linear or strong non-linear correlations between the hydraulic parameters, which is a sign of the overparameterization of the calibration problem. This was already noted in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF]. Indeed, k s is positively correlated to (θ s , α, λ) and negatively correlated to (θ r , n). This means that, for instance, an increase of k s -value would cause a deviation of the predicted state variables from the observations that could be compensated by increasing (θ s , α, λ) and simultaneously decreasing (θ r , n). The error variances are not correlated with the hydraulic parameters. Given that for this problem, Eq. ( 6.4) admits one single optimum, in the sequel the reference to the superscript m in the MCPD draws is dropped (see Section 3). 6.3. Predictive uncertainty of the observed variables. Except the fact that the MCPD-MC draws are assigned high probabilities, so far, there is no clue that they are really sampled from the target posterior jpdf. One way to check this assumption is to assess the 95% credible intervals assigned to the observed state variables with the MCPD sample. In this way, one can compare the predictive uncertainty obtained with the 95% credible intervals evaluated with the few MCPD draws. Indeed, as explained in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF], it is possible to assess the posterior uncertainty of the observed model responses because, for the latter, the likelihood functions have been defined before calibration. Let us denote by y * out a new observation of either the pressure head (out = h) or the water content (out = θ), then the posterior pdf of y * out reads,

p(y * out |y h , y θ ) = d i=1 Ni ki=1 P i (x ki )p(y * out |y h , y θ , x ki , σ ki h , σ ki θ ) d i=1 Ni ki=1 P i (x ki ) (6.5)
where one recognizes P i (x ki ) the probability assigned to the MCPD draw x ki (i.e. the posterior jpdf Eq. (6.4) evaluated at x ki ). The likelihood function at y * out evaluated at the k i -th MCPD draw is denoted by p(y * out |y h , y θ , x ki , σ ki h , σ ki θ ). These two quantities are obtained by evaluating the model responses of interest at each MCPD draw x ki ,

k i = 1, . . . , N i , i = 1, . . . , d.
The two predicted 95% credible intervals are represented on Figure 6.2. We can note that they match very well. Furthermore, they encompass most of the observations. In our viewpoint, this is a clear indication that the MCPD-MC sample was generated with respect to the posterior jpdf. We recall that in [START_REF] Mara | A comparison of two Bayesian approaches for uncertainty quantification[END_REF], MCPD and MCMC yielded similar predicted uncertainties for this calibration problem.

6.4. Posterior uncertainty analysis of the cumulative outflow. In subsurface hydrology, the cumulative outflow prediction is of high importance to predict water recharge of aquifers. While it is an issue to measure in-situ, this quantity is easily measured in laboratory experiment. In the present study, it corresponds to the amount of water exiting the column per unit surface in the end of the drainage experiment. The cumulative outflow is computed as follows,

C s = L 0 (θ(z, 0) -θ(z, t f )) dz,
where t f = 240 min is the duration of the drainage experiment. Predicting the uncertainty of this model response directly with the 185 MCPD draws is not possible because it requires the knowledge of the likelihood of this variable. The uncertainty of cumulative outflow a posteriori (i.e. after calibration) can be assessed by propagating the MCPD-MC sample through this model response variable.

The predicted posterior density of the final cumulative outflow is depicted in Figure 6.3. It shows a very narrow support with C s ∈ [0.69, 0.75] cm with a mode around 0.71 cm. This indicates that, as far as the prediction of the cumulative outflow is concerned, the hydraulic parameters have been satisfactorily calibrated. On the opposite, a larger support of the predicted posterior density would indicate a poorly calibrated model. In that case, a sensitivity analysis would help identifying which input variables were responsible for the lack of accuracy [START_REF] Saltelli | On the relative importance of input factors in mathematical models: Safety assessment for nuclear waste disposal[END_REF]. It is worth mentioning that sensitivity analysis can be easily carried out with the MCPD-MC draws at hand. For this purpose, one can use the method developed in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] which is valid under the assumption represented by Eq. (4.1). Although in the present case such an analysis is not necessary as the cumulative outflow is predicted accurately, we have undertaken a sensitivity analysis exercise for the sake of completeness.

The sensitivity analysis setting addressed in this exercise is the following: What is the smallest subset of hydraulic parameters x 1 ⊂ x = (k s , θ r , θ s , α, n, λ) that mostly explained the predicted variance of C s ? To answer this question, we consider the following Sobol' index [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]:

S closed x1 = V x1 E x2|x1 [C s |x 1 ] V x [C s ] (6.6)
where V is the variance operator and x = (x 1 , x 2 ). We investigate the smallest subset x 1 with the highest Sobol' index and such that S closed x1 > 0.99 which means that x 1 explains more than 99% of the variance of C s . From the MCPD-MC sample of size 512 already generated, the Sobol' indices of different subsets of various cardinality were computed with the method introduced in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF]. Obviously several subsets of equal cardinality that satisfied S closed x1 > 0.99 were found. But the one with the highest Sobol' index was attributed to x 1 = (k s , θ r , n). This result indicate if narrower predictive uncertainty is required for the cumulative outflow, one should devote further effort to reduce the uncertainty in the input set x 1 = (k s , θ r , n). 7. Calibration of a reactive transport model.

7.1. Problem setting. We consider the field observations of Valocchi et al. [START_REF] Valocchi | Transport of ion exchanging solutes in ground water: chromatographic theory field simulation[END_REF] that was used as a simulation benchmark by several authors ([48, 49, 6, 7]). The problem is a field experiment involving nonlinear ion exchange reactions in Palo Alto Baylands, California, where a multi-ion solution was injected into a well. For aqueous components, namely Mg 2+ , Ca 2+ , Na + and Cl -, were monitored at an observation well (named S23) located 16 m from the injection well. In this work, we only focus on the reactive species Mg 2+ and Ca 2+ . The experiment is modelled by a set of one-dimensional advective, dispersive and reactive transport equations,

∂C sp ∂t = -U ∂C sp ∂x + D ∂ 2 C sp ∂x 2 + R sp (x, t) (7.1)
where C sp (mmol.l -1 ) is the concentration of cation sp ∈ (Mg 2+ , Ca 2+ , Na + ), U (m.h -1 ) is the groundwater velocity, D (m 2 .h -1 ) is the dispersion coefficient of the homogeneous medium and R sp (mmol.l -1 .h -1 ) is a sink/source rate of cation due to chemical reactions. The different chemical reactions and associated log-constants of reactions are,

Na + + S - ⇋ NaS, log(K N a ) Ca 2+ + 2S -⇋ CaS 2 , log(K Ca ) Mg 2+ + 2S -⇋ MgS 2 , log(K Mg )
The chemical reactions are parameterized by two selectivity coefficients as follows, log(K CN ) = log(K Ca ) -2 log(K N a ) log(K MN ) = log(K Mg ) -2 log(K N a ) with log(K N a ) = 4. and S -the surface charge of the soil. This set of equations is completed by initial conditions and is solved numerically with the sequential noniterative approach [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF]. The domain of length 16 m is discretized with a uniform mesh of 400 elements. The simulated time T = 2000 h is discretized with a uniform time step ∆t = 0.2 h.

The model is parameterized by the vector of input parameters x = (log(K CN ), log(K MN ), T f , U, D), where T f (-) is the total fixed concentration. The observed data are the concentrations of Ca 2+ and Mg 2+ (mmol/l). We assume log-normal likelihood for each concentration data and independent uniform prior uncertainty ranges for x. This leads to the following target jpdf,

p(x, σ Ca , σ Mg |y Ca , y Mg ) ∝ 1 σ NCa Ca σ NMg CMg exp - 1 2 SS Ca (x) σ 2 Ca + SS Mg (x) σ 2 Mg , (7.2)
which has the same form as Eq. ( 6.4) except that y Ca and y Mg are the vectors of observed log-concentrations while SS Ca and SS Mg are the sum of squares of the differences between predicted and observed log-concentration of cations. Mg ≃ 0.14. Indeed, the MCPD curves are unimodal (Figure 7.1, on the diagonal), rather Gaussian for the parameters in vector (log(K CN ), log(K MN ), T f ) and skewed for those in (U, D, σ 2 Ca , σ 2 Mg ). The parameters in the former subset are well identified despite of the strong collinearity between (log(K CN ), log(K MN )) (row #2, column #1). Those in the second subset have large posterior uncertainty ranges with a negative strong correlation between (U, D) (row #5, column #4). The hyperparameters are not correlated to the parameters (see the last two rows).

From the MCPD sample, a MCPD-MC sample of size N = 512 was generated. To check the validity of Eq. (4.1), for each MCPD-MC draw the predicted logconcentrations of Ca 2+ and Mg 2+ were collected after running the transport model. Then, with the collection of model responses the jpdf in Eq. (7.2) was evaluated. Therefore, each MCPD-MC draw is assigned a probability value which has been scaled properly within [0, 1] and plotted on Figure 7.1 (dot plots). On the diagonal plots, we notice that the MCPD-MC draws are assigned high probability values. The dot plots filled the space between the MCPD curves and the x-axis. None of the MCPD-MC draw was assigned a probability zero. This is an indication that the MCPD-MC is successfully sampled from the desired jpdf. This is confirmed by the comparison of the predictive uncertainties obtained with the MCPD draws (by using Eq. (6.5)) and the Monte Carlo sample respectively (see Figure 7.2). Notably, the predictive uncertainties are narrow and encompass the observations. This indicates that the choice of the log-normal likelihood is reasonable.

In our study, cations N a + were not involved in the calibration process, only Ca 2+ and Mg 2+ were. Hence, it can be questioned whether the calibration of the model from the concentration measurements of Ca 2+ and Mg 2+ are sufficient to estimate the total mass of Na + (denoted M N a ) accurately. A Monte Carlo sample of size N = 512

Conclusion.

In this paper, we have extended the maximal conditional posterior distribution sampling to generate Monte Carlo draws from a joint probability distribution function. The proposed approach requires that the correlation structure amongst the random variables be such that Eq. (4.1) is valid at least for one possible ordering of the random variables in the set (x 1 , . . . , x d ). Generating Monte Carlo samples is essential, for instance, to perform uncertainty and sensitivity analyses of model responses. This is illustrated in our work.

The numerical exercises have shown that the proposed approach, named MCPD-MC sampler, was able to sample from distribution functions with possibly multimodality and complex correlation structure. The MCPD-MC approach relies on the orthogonalization procedure introduced in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] that produces independent samples under assumption (4.1). Fast generation of random samples from the MCPD draws is performed without evaluating the jpdf further. However, evaluation of the jpdf with each generated draw is necessary to verify that the latter is located in the region of high probability. This verification step is required to check the validity of Eq. (4.1).

Applications to model calibration have demonstrated the interest of the proposed approach. With the MCPD-MC sample, uncertainty and sensitivity analyses of the model responses that were not used in the calibration step could be performed. Furthermore, sensitivity analysis can be easily carried out with the method developed in [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] since the latter relies on the same assumption as the MCPD-MC sampler. However, it has to be underlined that Eq. (4.1) may not be satisfied in some situations. Therefore, great care must be taken when applying the MCPD-MC approach. Otherwise, another alternative (more expensive) relying on Eq. (4.6) can be employed.

It has to be underlined that the current version of the MCPD-MC sampler for statistical model calibration only handles Gaussian likelihood functions. Current development of the sampler aims at including a less restricted likelihood function by accounting for the Generalized Error Distribution (also called Exponential Power distribution [START_REF] Box | Bayesian inference in statistical analysis[END_REF]).

Figure 3 . 1 .

 31 Figure 3.1. Example of MCPD draws prior and posterior to the transformation procedure Eq. (4.2). Subplots on the diagonal display the MCPD draws of the original variable x i and its transformed one xi (blue crosses and red squares resp.). The lower off-diagonal plots show the pairwise correlations between the original variables. The draws decorrelated with the transformationEq. (4.2) are depicted on the upper off-diagonal.

4. 4 .Figure 4 . 1 .

 441 Figure 4.1. Same as Figure 3.1 with the Monte Carlo draws obtained with the new approach.The upper off-diagonal plots show no correlation amongst the draws (red dots). The lower diagonal plots show the desired draws sampled from the target distributions (gray dots).

  d ) is the multiGaussian distribution of mean vector µ i and covariance 5I d . I d is the d-dimensional identity matrix which indicates that the parameters (x 1 , . . . , x d ) are independent in the second and third Gaussian distributions in Eq. (5.1). In the present work, we consider the case d = 11 with the correlation matrix C having null off-diagonal elements except for C 1,2 = C 2,1 = -0.5 and C 1,3 = C 3,1 = 0.8. These non-null terms impose, for the first Gaussian distribution in Eq. (5.1), a negative correlation between x 1 and x 2 and a strong positive correlation between x 1 and x 3 . The three modes of each variable are grouped in the vectors of means µ 1 = (-5, -4, . . . , 4, 5) T , µ 2 = (1, 2, . . . , 11) T and µ 3 = (11, 10, . . . , 1) T .

Figure 5 . 1 .

 51 Figure 5.1. Sampling from the multimodal target distribution. Results of first three variables for a sample of size N = 4 096 are depicted.The on-diagonal plots represent the MCPDs estimate (blue line) as well as the MCPD-MC draws. The lower off-diagonal scatterplots represent the pairwise correlation of the MCPD draws (blue crosses) and the MCPD-MC draws (gray dots). r

2 Figure 5 . 2 .

 252 Figure 5.2. Sampling from the banana-shaped distribution. Results of first three variables for a sample of size N = 4 096 are depicted.

Figure 6 . 1 .

 61 Figure 6.1. The generated MCPD-MC draws of the numerical drainage experiment. On the diagonal, the red curves are the MCPDs, the black dots represent the MCPD-MC draws.

3 ]Figure 6 . 2 .

 362 Figure 6.2. The dots represent the data used for soil hydraulic parameter identification. The predictive uncertainty ranges (95% credible interval) estimated with the MCPD and MCPD-MC samples match satisfactorily well, confirming that the MCPD-MC draws are reliable.

Figure 6 . 3 .

 63 Figure 6.3. Predictive uncertainty of the cumulative outflow in the end of the experiment.

7. 2 .

 2 Uncertainty quantification. The statistical calibration of the transport model involves seven variables, namely the model parameters x = (log(K CN ), log(K MN ), T f , U, D), and the hyperparameters (σ 2 Ca , σ 2 Mg ) of the likelihood function. Their MCPD estimates are depicted in Figure 7.1. The optimization algorithm found a single set of optimum values which is, log(K opt,1 CN ) ≃ 8.69, log(K opt,1 MN ) ≃ 8.43, T opt,1 f ≃ 735.64, U opt,1 ≃ 0.74, D opt,1 ≃ 4.35, σ opt,1 Ca ≃ 0.09 and σ opt,1

Figure 7 . 1 .

 71 Figure 7.1. The generated MCPD-MC draws of the reactive transport experiment. On the diagonal, the red curves are the MCPDs, the black dots represent the MCPD-MC draws.

Figure 7 . 2 .

 72 Figure 7.2. Breakthrough curves of cations in the reactive transport experiment. The predictive uncertainty ranges (95% credible interval) estimated with the MCPD and MCPD-MC samples are comparable, confirming that the MCPD-MC draws are reliable.

Figure 7 . 3 .

 73 Figure 7.3. Predictive uncertainty of the total mass of Na + involved in the chemical reactions during the flow transport experiment.

Acknowledgements. The MCPD-MC sampler developed in MATLAB is available upon request from the corresponding author (mara@univ-reunion.fr). The authors are grateful to the anonymous referees for their insightful comments that helped improving the manuscript. Qian SHAO acknowledges the support from the National Science Foundation of China (Grant No. 11702199).

* This work was supported by the French National Research Agency through the research project RESAIN n • ANR-12-BS06-0010-02.

was obtained by propagating the MCPD-MC draws through the model response M N a . The estimated pdf is depicted in Figure 7.3. The support of the latter is large, ranging from 10 to 19 mol.m -2 with a mode around 16 mol.m -2 .

Identifying those uncertain input variables responsible for the large uncertainty in the model response is the role of sensitivity analysis. In the present analysis, we quantify the amount of the variance of M N a explained by each individual input parameter. This statistic is called the Pearson correlation ratio or the first-order Sobol' index and is defined as follows [START_REF] Pearson | On the general theory of skew correlation and non-linear regression[END_REF][START_REF] Sobol | Estimating the approximation error when fixing unessential factors in global sensitivity analysis[END_REF],

which is the same as Eq. (6.6) but for individual variable x i such that (x i , x ∼i ) = x.

Using the sensitivity analysis method of Mara & Tarantola [START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF] as previously, we find: S KCN ≃ 0.74, S KMN ≃ 0.73, S T f ≃ 0.18, S U ≃ 1 and S D ≃ 0.99. These results simply claim that the large uncertainty in M N a prediction is due to the inaccurate assessment of the groundwater velocity U . The high sensitivity index of D is merely due to its strong correlation with U (see row #5, column #4 in Figure 7.1) which was poorly estimated from the measurements of Ca 2+ and Mg 2+ concentrations. Including measurements of chloride concentrations in the calibration dataset did not improve much the flow velocity estimate even though Cl -is a tracer (not shown).