Bayesian calibration of order and diffusivity parameters in a fractional diffusion equation

This work focuses on parameter calibration of a variable-diffusivity fractional diffusion model. A random, spatially-varying diffusivity field with log-normal distribution is considered. The variance and correlation length of the diffusivity field are considered uncertain parameters, and the order of the fractional sub-diffusion operator is also taken uncertain and uniformly distributed in the range (0,1) . A Karhunen-Loève (KL) decomposition of the random diffusivity field is used, leading to a stochastic problem defined in terms of a finite number of canonical random variables. Polynomial chaos (PC) techniques are used to express the dependence of the stochastic solution on these random variables. A non-intrusive methodology is used, and a deterministic finite-difference solver of the fractional diffusion model is utilized for this purpose. The PC surrogates are first used to assess the sensitivity of quantities of interest (QoIs) to uncertain inputs and to examine their statistics. In particular, the analysis indicates that the fractional order has a dominant effect on the variance of the QoIs considered, followed by the leading KL modes. The PC surrogates are further exploited to calibrate the uncertain parameters using a Bayesian methodology. Different setups are considered, including distributed and localized forcing functions and data consisting of either noisy observations of the solution or its first moments. In the broad range of parameters addressed, the analysis shows that the uncertain parameters having a significant impact on the variance of the solution can be reliably inferred, even from limited observations.

Despite recent advances, the application of fractional diffusion models remains challenging. One challenge concerns the numerical solution of the governing equation, which typically involves full matrices that arise from

Governing equation and deterministic solver
We focus on the steady, 1D, two-sided, fractional order diffusion equation with variable diffusivity:  (1), κ denotes the diffusivity, f is the forcing term, whereas ¶ x denotes the first-order derivative, and ¶ a x is the two-sided fractional differential operator of order a Î 0, 1 ( ), ]. We seek the numerical solution of problem(1) by adopting the finite difference method proposed in [17], which combines first-order forward and backward approximations of the Riemann-Liouville derivatives. Therein, the existence and uniqueness of the finite difference solution were established, and a rigorous truncation error analysis was provided. We consider a uniform partition of the domain using a sequence of

Stochastic formulation
A stochastic framework is introduced by considering uncertainty in the order of the fractional exponent α, and in the diffusivity κ. The former is specifically modeled as a random variable, whereas the latter is modeled as a log-normal random field of the form: where m k x ( ) is the median of the log-normal field and q M x, ( )is a real-valued, centered, Gaussian process characterized by the square exponential covariance function: In the covariance expression, L M is the correlation length and s M is the standard deviation of the Gaussian field. In the following L M and s M are either given, or alternatively treated as uniformly-distributed random variables. Without loss of generality, we set throughout the paper m = k x 1 ( ) . The covariance function, C is continuous, symmetric and positive-definite. By Mercer's theorem [78], it can be decomposed spectrally as: where l k are eigenvalues, and f k are the associated L 2 -orthonormal eigenfunctions of the Fredholm integral operator with kernel, C. Thus, the truncated KL expansion of κ is given by [70]: where N KL is the number of the terms retained in the KL expansion, and the x q k ( ) are uncorrelated Gaussian random variables with zero mean and unit variance.

Polynomial chaos surrogates
Surrogate models of the solution, or of QoIs, are used to conduct sensitivity analyses of the solution and to accelerate the Bayesian calibration of κ and/or α [79]. We rely on a PC formalism for this purpose [71]. This section briefly outlines the PC methodologies.
Consider an uncertain model inputs parametrized by an N-dimensional vector of canonical random variables, x, with independent components x 1 , K, x N . The random vector x is sometimes called the stochastic germ. The PC approximation x U ( )  of a generic second-order functional x U ( ) of the model output is expressed as: where N t is the number of terms retained in the PC expansion. Note that Roman subscripts are used to denote discrete spatial locations, whereas Latin subscripts denote PC modes index.
} forms a complete orthogonal polynomial basis of the space of second-order random functionals in x, denoted is the density of i-th random variable, x i . We will rely on non-intrusive approaches [71] to determine the PC coefficients, U β . These methods generally aim to minimize the L 2 error -U U    , based on a discrete set of deterministic model realizations. In this work, both regression and spectral projection approaches are used to perform the minimization. The regression approach [80] calculates the coefficients b   U N 0 t that minimize the least-squares error between the solution x U ( ) and its approximation x U ( )  . Consider a sample set of N s realizations, the corresponding sample set of observations of the QoI x V ( ). The L 2 minimization problem can be written in matrix form as follows: ) is the vector of PC coefficients. In this work, we will also consider the regularized L 2 minimization problem, that is solved using LASSO [81]. In this case, the minimization problem can be expressed as: and g  0 is a parameter that is also optimized by a cross validation procedure. The non-intrusive spectral projection (NISP) approach [82] aims to compute the PC coefficients of the model output through an orthogonal projection onto the PC basis. Hence, the evaluation of the coefficients proceeds as follows: Estimating the expansion coefficients amounts to the computation of N-dimensional integrals in a product space. Classically, the integrals are estimated by employing a numerical quadrature or a sampling approach. Since the Y b are polynomials, we will apply fully-tensorized Gauss quadrature rules when the dimensionality N is small. Otherwise, we rely on an isotropic sparse-grid method (Smolyak formula [83]) based on the nested quadrature rules [71].

Setup
This section outlines the scope of the forward uncertainty quantification, sensitivity analysis and inference studies. As outlined in table 1, four cases are considered in section 6. Case1 concerns only one uncertain parameter, namely the fractional order α. In Case2, we use fixed L M , s M , and α, and first study the solution sensitivity on the diffusivity field's KL coordinates. We then tackle the problem of calibrating κ. Case3 also considers the impact of the leading diffusivity field's KL coordinates but also includes variability in α. Case4 explores the sensitivity of the first two solution's moments, namely the mean, m x ( ), and standard deviation, s x ( ), to the correlation length, L M , and standard deviation, s M , defining the log-normal diffusivity field, and the fractional order, α. When considered uncertain, the parameters are assumed to be independent and follow uniform priors: Each case also considers two types of forcing functions, distributed and localized. In the first type, the forcing function is In the second type, the source has a sharp peak at x 0 in the interior of the domain and follows As outlined in table 1, different PC methodologies are used depending on the case considered. In Case1, we use the NISP approach to approximate the dependence of the solution u on α. In Case2, we rely on a LASSO regression approach to build the PC surrogate of u as a function of the KL coordinates. As for Case3, the surrogate of u is built via a fully-tensorized NISP machinery. In Case4, we use a sparse-grid, pseudo-spectral projection (PSP) technique to construct functional approximations of the solution moments, in terms of α, L M and s M .
In all cases, we consider a discretization of the physical domain, 0, 1 [ ], with P=300 nodes, and we apply the finite-difference scheme introduced in section 2. This discretization follows a careful grid refinement study (not shown), which revealed that this resolution is sufficient to obtain accurate predictions for the entire range of physical parameters considered. Also, unless otherwise stated, the KL is truncated by retaining the first five dominant modes (i.e.
, which enables us to capture 99.9% of the underlying Gaussian field variability at the lowest value of the correlation length considered, L M =0.5.

Numerical experiments
This section presents for each case the surrogate construction, the results of the sensitivity analysis, and, finally, the inference results. For Cases 2-4, we perform a global sensitivity analysis [84] by calculating the Sobol indices [85]. These indices provide a variance-based importance measure of the effects of individual parameters on the variability of the solution or the QoI. After that, we focus on the Bayesian calibration problem based on synthetic data. The quality of inferred parameters is assessed based on a limited amount of noisy observations for each case scenario. Specifically, particular values of the parameters, called hereafter the true values, are picked within their respective prior ranges, and noisy observations are generated from the corresponding surrogates using a multiplicative Gaussian noise. The Maximum Likelihood Estimation (MLE) and Maximum A Posteriori (MAP) values of the parameters are retrieved, and the posterior densities are explored using an adaptive Markov Chain Monte Carlo (MCMC) algorithm [86]. Sufficiently long MCMC chains (not shown) ensure that stationary densities are adequately sampled. Kernel density estimation (KDE) [87,88] was subsequently applied to the MCMC samples to estimate the posterior densities.

CASE1
We build a surrogate for the solution u considering one uncertain parameter, α. For this purpose, we used fixed and m = k 1. Figure 1 shows the corresponding diffusivity field.

Surrogate error
For the PC construction, we utilize a Gauss-Legendre quadrature rule with N t =8 nodes, enabling us to determine a seventh-order expansion. In figure 2 we compare the solution u and its surrogate u  for 16 randomly sampled values of α in the prior range. The figure demonstrates the surrogate accuracy for both distributed and localized forcings. For a more quantitative accuracy assessment, the following relative error estimate, is the surrogate prediction and  i ʼs are i.i.d. Gaussian random variables with zero mean and standard deviation r = 0.01 t . The log-likelihood of the observations is consequently expressed as is the observation vector and ρ is a hyper-parameter that we also infer from the data. A Jeffrey's prior for ρ is assumed, with density r r = p 1 2 ( ) . Noisy observations are generated at m d =5 observation points, located at x=0.1, 0.3, 0.5, 0.7 and 0.9. For the distributed forcing case, the noisy observations are generated from a single experiment, observing the solution at the specified stations. In the localized forcing case, the observations come from m d independent experiments involving a localized forcing centered at the corresponding observation station and using the same 'system' (medium). The 'true' solution profiles (u ) and the noisy observations are depicted in figure 3.
The marginal posterior densities of α for the distributed and localized forcing scenarios are plotted in figure 4. In both instances, the posterior of α is a narrow distribution with a well-defined peak. The MAP estimates are a = 0.4177 for the distributed forcing experiment, and a = 0.4788 for the localized forcing   figure 3. A very close agreement between the true and inferred profiles is observed in the distributed case, whereas we notice some discrepancies for the localized forcing experiments. These higher discrepancies are consistent with the inferred MAP and MLE of α, which are closer to the truth in the distributed forcing instance. Additional insight into the behavior of the two instances can be gained from figure 5 and 6 which depict the solution, u, together with the 95% bounds of the pushed-forward prior and the 50% bounds of the posterior of α, considering distributed and localized forcing, respectively. The plots indicate that the response of the stochastic solution and the inferred posterior can depend significantly on the imposed forcing and provide insight into the computed posterior distributions. Specifically, one notes that though the posterior bounds in the distributed forcing case are tight, they clearly reflect the bias between the MAP and true parameter (see figure 4). for localized forcing setups, as indicated. Also plotted are curves showing the 95% bounds of the solution prior (blue dash) and the 50% bounds of the solution posterior (red dash). The noisy measurements used to perform the inference are also depicted. Note that a single inference step is performed using all five data points simultaneously.  On the other hand, the 50% bounds in the distributed forcing case are slightly broader but clearly capture the true solution. Of course, a more detailed analysis based on repeated realizations would be needed to provide a full characterization of the expected information gain associated with different strategies and its dependence on  the number of observations. Such analysis is outside the scope of the present exploratory study and is left for future work. ,..., s then, for each sample point x i we evaluate the corresponding solution u i using the deterministic model, to constitute the training set of solutions used in the LASSO method. Moreover, surrogates using NISP were also established and produce surrogates (not shown) of similar quality. Figure 7 shows a comparison between the deterministic solution u and its surrogate for different forcing at randomly selected values of x. A good agreement is observed between the model and its surrogate. To provide a quantitative estimate of surrogate errors, we plot in figure 8, profiles of local L 2 errors between the surrogate and deterministic realizations generated using an independent coarse sample. Shown are results obtained for distributed forcing and forcing functions localized at x 0 =0.1 and 0.5. As expected, the localized profiles peak at the same locations as those of the forcing function. However, in all cases, the errors are small and fall below 2% of the peak values.

Global sensitivity analysis
We first explore the variability of u due to the KL coordinates x i of the diffusivity field, for the case of distributed forcing. Figure 9 displays the first-order and total-order partial variances of u due to x i . The plots show clearly that the first KL coordinate, x 1 , is responsible for most of the variance. The second coordinate, x 2 , exhibits a small contribution to the variability of u, between the boundaries and midpoint of the domain. One also concludes from figure 9 that other coordinates have an insignificant effect on the variability of the solution. For the sensitivity of u to the KL coordinates, figure 10 reports first-order and total-order partial variances of u due to the KL coordinates x i when f is localized at x 0 =0.1, and at x 0 =0.5. The figure indicates that most variability in u is due to the first mode x 1 for the two forcing locations, while x 2 has an appreciable effect only when x 0 =0.1. Other coordinates are seen to be unimportant regardless of the forcing location.

( )
The noisy observations are generated using where the  i ʼs are i.i.d. Gaussian with zero mean and standard deviation r = 0.01 t and x I II , denotes one of the two instances. The inference relies on the following expression of the log-likelihood: is the observation vector, and ρ is the noise hyper-parameter that is assumed to follow a Jeffrey's prior.
We first comment on the inference results for InstanceI. Table 2 presents the MAP and MLE values of the inferred x i and noise level ρ, for the two forcing types. The table shows that the MLE values of the first three coordinates are in better agreement with the true values compared to the MAP. The MAP value of ρ is significantly higher than the actual value used to corrupt the synthetic data. Note that, in light of the sensitivity analysis presented in the previous section, one cannot expect to recover the coordinates of the higher modes from limited noisy observations because these coordinates have an insignificant impact on the solution. The MLE and MAP estimators are further studied in figure 11 in terms of solution profiles and diffusivity fields. Consistently with the results of table 2, the plots reveal a closer agreement between the true solution and the corresponding MLE estimators than is obtained for the MAP. This is seen for both the distributed and localized forcing setups. The discrepancies between the MLE estimator of κ and the truth are much more significant than the discrepancies between the MLE estimator of u and the true solution. Further, we remark that the localized forcing experiment gives a better MLE estimate of the diffusivity field, while the MAP estimator remains far from its true value in all forcing scenarios. Figure 12 shows the posterior marginals of the field coordinates. The first coordinate, x 1 , exhibits a marked peak, while the other coordinates essentially retain their prior shapes. One can interpret the present results in light of the dominant sensitivity of x 1 , and the insignificant contribution of the remaining coordinates. Further, one notes that the second and third coordinates of the true field have large magnitudes, with quite unlikely a priori value; it is well-known [89] that, for an informative prior such as the Gaussian one, the Bayesian calibration needs a substantial amount of data to localize the posterior in a region of low prior probability. These two remarks explain why, in the present case, the Bayesian inference leads to broad predictive distributions, as further discussed below.
The results outlined above for InstanceI are in contrast with those obtained for InstanceII. Specifically, the MLE and MAP values of the coordinate x i reported in table 3 are now in much closer agreement with their true values, with differences not exceeding a fraction of the a priori standard deviation. Consistently, the corresponding MLE and MAP solution profiles, reported in figure 13, are also closer to the true profile, for both  the distributed and localized forcings. Also, the MAP and MLE estimates of the diffusivity fields are also closer to the true κ, compared to the case of InstanceI shown in figure 11, although differences are still significant.
Finally, the posterior marginals of the coordinates presented in figure 14 reveal a more appreciable information gain for x 2 , compared to InstanceI in figure 12, with a clear peak value close to, but distinct from, the origin. The improvement in information gain for x 2 is more pronounced for the localized forcing than for distributed f; this suggests that the asymmetrical structure of the localized forcing can better sense the antisymmetry of the mode associated with this coordinate.
To provide additional insight into the inferred solutions, we plot in figure 15 and 16 the solution, u, the 95% bounds of the pushed-forward prior, and the 50% bounds of the pushed-forward posterior, respectively for distributed and localized forcing. Contrasted in each of these figures are results obtained for Instances I and II. Unlike Instance II, one can observe that the true solution in Instance I can lie appreciably outside the bounds depicted. This is consistent with the previous discussion concerning the fact that a parameter with low prior probability was drawn in Instance I. A larger amount of data would be needed to better estimate the field in these situations.   Figure 18. Case3: first-order (left) and total-order (right) partial variances of u due to α, x 1 and x 2 , for the distributed forcing case. Figure 19. Case3: first-order and total-order partial variances of u due to α, x 1 and x 2 . Plotted are results obtained using localized forcing with x 0 =0.1 (top row) and x 0 =0.5 (bottom row).

CASE3
In Case3, we aim to infer α and coordinates of the diffusivity field κ. Based on the results of section 6.2, we restrict the parametrization of κ to the first two KL modes, setting x = 0 i for > i 2. The prior correlation length and standard deviation of the log-normal distribution of κ are held fixed, L M =0.75 and s = 0.5 M , respectively. We first construct an approximation of u in terms of α, x 1 , and x 2 , using a NISP approach to estimate the PC coefficients. Here, the multidimensional quadrature points resulting from the tensorization of 8 Gauss-Legendre (α uniformly distributed) and Gauss-Hermite (x x , 1 2 i.i.d. Gaussian) nodes. Similarly, the PC basis consists of fully tensorized, Legendre and Hermite polynomials with maximum (partial) order 7. In figure 17, the PC surrogate, û, of u is compared with the solution for a selected value of the random input vector, considering both distributed and localized forcing. A close agreement between u and û is seen for all instances depicted in figure 17.

Global sensitivity analysis
The Sobol decomposition of the variance of u is presented in figure 18 for distributed forcing, and in figure 19 for localized forcing. First-order and total-order variances are plotted to illustrate the contributions of the three uncertain inputs, α, x 1 , and x 2 . In all cases, the results indicate that α has the dominant role on the variability of the solution, followed by x 1 , whereas the impact of x 2 is comparatively negligible. First and total-order partial variances due to x 1 peak in the middle of the domain for the distributed forcing, and at x 0 for the localized one. One can also deduce from the curves that mixed interactions between α and x 1 are appreciable at locations around the maximum of the forcing.

Inference of α and κ
To estimate α, x 1 and x 2 , we consider noisy observations of the solution, is the observation vector, and ρ is a hyper-parameter equipped with Jeffrey's prior. Table 4 summarizes the MAP and MLE estimates obtained from the Bayesian calibration for the distributed and local forcing setups. The table reveals that the MLE values of α, x 1 , and x 2 are in better agreement with their respective true values than the MAP estimates. In particular, the MAP estimators of α overestimate the true values for both setups. Further, the noise level is correctly estimated by the MLE value but overestimated by the MAP value. Similar to Case2, more observations would undoubtedly help to improve the MAP estimators. Figure 20 compares the true solution profiles and the diffusivity field with their MAP and MLE counterparts in the distributed and localized forcing experiments. We see that the MLE profiles are in excellent agreement with the corresponding true profiles. In contrast, substantial discrepancies are observed between the MAP profiles, especially for the distributed forcing experiment. The differences between the MLE and MAP profiles reflect the calibration results in table 4. Concerning the diffusivity fields, figure 20 depicts the closer agreement of κ using MLE, especially in the distributed forcing scenario. The present results are also consistent with previous ones, which indicated that, at least for the present framework, sampling data at several locations in a single distributed forcing experiment provides a better means to recover the parameters than collecting a single response from the same number of localized forcing experiments.
The posterior marginals of the inferred parameters are plotted in figure 21. The curves indicate that for both distributed and localized experiments, the posterior marginals for α have a sharp peak. The marginals of x 1 and x 2 exhibit discernible peaks, although only the marginal of x 1 appears to differ significantly from its (Gaussian) prior. Again, this behavior was expected from the results of the global sensitivity analysis.

CASE4
In the last case, we focus on the inference of the statistical characteristics of the stochastic diffusivity field jointly with the fractional order α. Specifically, besides α, we try to infer the correlation length, L M , and standard deviation, s M , of the Gaussian field underlying κ. The inference uses the mean and standard deviation of the solution to calibrate the parameters. At given values of α, L M , and s M , the mean, μ, and standard deviation, σ, of the solution are estimated, sampling the first five KL coordinates of the corresponding κ. In this work, we rely on a large LHS strategy to obtain well-converged estimates of μ and σ. Figure 22 compares the estimates of μ and σ obtained on LHS sets with sizes = N 15, 000 s and 25,000. The differences are not significant and, consequently, we used = N 25, 000 s for the surrogate construction. To build the surrogates for the dependencies of μ and σ on α, L M , and s M , we considered a sparse-grid quadrature with Féjer rules [71]. A detailed convergence study of the surrogates was conducted to assess the suitability of the PC approximation; a level 5 quadrature was found sufficient in this case. Figure 22 illustrates this study, contrasting the surrogates prediction for level 4 and level 5 in the sparse grid construction, considering a distributed forcing with the same values of L M , s M and α as before. A similar exercise for the localized forcing functions, with peaks located at x 0 =0.1, 0.3, 0.5, 0.7, and 0.9, showed an excellent agreement between the direct estimates of μ and σ and their sparse grid surrogate predictions, see figure 23 for an illustration.
The sparse grid projection methods are susceptible to evaluation noise. Therefore, a high number of LHS samples are necessary to ensure a sufficiently converged estimation of μ and σ at each of the sparse grid nodes. To alleviate this computational burden, one can instead rely on an alternative approach presenting a more robust behavior against evaluation noise. As a final check of our surrogates for μ and ν, we reduced the LHS set size to = N 1, 000 s for the estimation of the μ and σ, and relied on a regression approach using a LHS set of 300 triplets s a L , , M M ( ) . Figure 24 compares a particular realization of μ and σ and their regression surrogate approximations m  and s  for distributed and localized forcing at x 0 =0.1 and x 0 =0.5. A detailed analysis (not shown) proved that the regression surrogates achieve an accuracy comparable to the sparse grid ones, but for just a fraction of the computational cost. These experiments also demonstrated that the sparse grid surrogates are not affected by the estimation noise, and they will be used in the remainder of the section.

Global sensitivity analysis
Results of the variance decomposition are provided in figure 25, considering distributed and localized forcings. Concerning the solution mean μ (top row), we only show the total-order partial variances due to α, L M and s M , because the contribution of α is highly dominant, L M and s M have weak impacts, and interaction effects are negligible for the three cases considered. In contrast, the variance decomposition of the solution's standard deviations shows non-negligible contributions and interactions between α and s M . At the same time, the correlation length has an insignificant impact on σ, with a total-order variance due to L M substantially smaller than the others.

Inference of L M , s M and α
We consider the problem of inferring L M , s M , and α from observations of the solution statistics, namely the solution's mean μ and standard deviation σ. Noisy observations of the mean, m ĩ , and standard deviation, s ĩ , are generated at selected points x i , from their respective PC surrogates, according to: respectively denote the vector of observations of the solution mean and standard deviation, and r r r = , ) is the vector of hyper-parameters, whose components are assumed to follow independent Jeffrey's prior distributions.
The calibration is performed considering distributed and localized forcing. Table 5 provides the computed MLE and MAP values of L M , s M , α, and noise parameters; the last row of table 5 also shows the true values used for the generation of the noisy observations. The MAP and MLE values of s M and α are in excellent agreement and close to the true values. In contrast, the estimates of the correlation length differ greatly depending on the experiment. This finding is consistent with the negligible impact of L M on the mean and standard deviation of the solution, reported in the sensitivity analysis. The inferred hyper-parameters, r 1 and r 2 , differ significantly from their true values, the MAP estimates being closer than the MLE values, which consistently underestimate the noise level. Figure 26 shows the true profiles of μ and σ and their corresponding MAP and MLE estimates using the values reported in table 5. Also shown are the m d =16 noisy observations used for the calibration. The inferred MAP and MLE profiles are in excellent agreement with the corresponding truth profiles for both the mean and standard deviation and all three forcing functions considered.
Finally, figure 27 presents the marginal posteriors of L M , s M and α. As expected, the marginal posteriors of s M and α are concentrated and exhibit unique, well-defined maxima. On the other hand, we note that the marginal pdf of L M remains flat and exhibits several equivalent local maxima. The flatness of the marginal posterior of L M is consistent with the findings of the sensitivity analysis, discussed above, which indicates that it has an insignificant impact on μ and σ. Consequently, one expects the present observations to be essentially uninformative about this parameter, such that it essentially retains its flat uniform prior.

Conclusions
In this paper, we have investigated a steady, 1D, fractional diffusion equation involving uncertain fractional order and uncertain, spatially-varying, diffusivity field. Four different cases were analyzed, corresponding to different combinations of random inputs, QoIs, and localized or distributed forcing terms For all cases, we have constructed surrogate models for the QoIs, using non-intrusive PC approaches. With these PC surrogates, we have conducted sensitivity analyses of the QoIs. In particular, these analyses revealed the dominant role of the fractional order, α, followed by the first KL coordinates involved in the diffusion field's parametrization (largescale fluctuations). The sensitivity analyses also showed that the correlation length of the Gaussian process underlying the diffusion field has a weak impact on the solution mean and variance, at least for the range of values considered.
We also explored the calibration of the uncertain inputs, using direct observations of the solution or some of its statistics. A Bayesian formalism was specifically adopted with synthetic observation derived from the surrogates. In particular, the experiments showed that the dominant parameters, such as the fractional order, α, can be robustly inferred using limited information. The investigations also showed that using the same number of noisy measurements in different experimental setups may lead to appreciably different inverse solutions. This suggests the application of optimal experimental design methodologies [90,91] to quantify the utility of different setups. This will be addressed in future work, which will also consider the extension of the existing approaches to higher dimensions, as well as the incorporation of fast solvers [47,53,54] to enhance the efficiency of the surrogate construction.