Sensitivity Analysis, Bayesian Calibration and Forward Propagation of Uncertainties Using quoFEM
Pedro Arduino  University of Washington
SangRi Yi  SimCenter, UC Berkeley and
Aakash Bangalore Satish  SimCenter, UC Berkeley
key Words: quoFEM, OpenSees, Tapis, Python
Resources
This example makes use of the following DesignSafe resources:
Step  Notebook 

Sensitivity analysis  
Bayessian calibration  
Forward propagation 
The notebooks, and required scripts, are available in the Community Data folder and can be executed without any modification. Users are invited to try this notebook and use any parts of it.
Description
This illustrative use_case demonstrates several UQ techniques using the parameters of the PM4Sand constitutive model, a liquefactioncapable soil model in OpenSees. This complex material model is often calibrated using a small number of experimental results which yields imperfect information about its parameters. This leads to uncertain model predictions. Quantifying such uncertainties and inspecting the uncertainty bounds of model predictions can provide more information about the importance of each model parameter. Recognizing these uncertainties can incentivize more sophisticated modeling and calibration techniques that can better utilize the available data from experiments to reduce these bounds and provide more robust and higher fidelity simulations.
In this use case, the amount of reduction in the uncertainty in PM4Sand parameters calibrated to Cyclic Direct Simple Shear (CyDSS) test data is inspected, and the resulting uncertainty is propagated in an earthquake excitation simulation of a soil column. Three steps of UQ analyses, schematically shown in Fig. 1, are presented:

Global sensitivity analysis to get an insight into which parameters are critical in triggering liquefaction. This is an important first step to decide which parameters need to be included in the calibration process.

Bayesian calibration to obtain the posterior probability distribution of the PM4Sand parameters based on the CyDSS test dataset.

Forward propagation to investigate how the uncertainty that remains after the Bayesian calibration (characterized by the posterior probability distribution) affects the prediction of an earthquake response.
Fig.1  Probabilistic calibration of soil model (step 2) with sensitivity analysis (step 1) and prediction of uncertainty in estimation of lateral sperading (step 3)
Implementation
Uncertainty Quantification Using quoFEM
Accurate quantification of uncertainty requires wellestablished workflows that incorporate sophisticated UQ techniques with advanced simulation models and frameworks. The SimCenter quoFEM tool streamlines this process by offering comprehensive workflows in a single tool, which can be accessed locally or remotely through a web browser using the DCV client on DesignSafe. Furthermore, users can utilize the Jupyter Hub environment on DesignSafe to manage the same, or additional, runs via Python scripts, defining job variables and submitting jobs through the Tapis system. This allows for seamless collaboration and efficient job management, resulting in faster and more effective UQ analysis.
In this context, the notebooks included in this usecase complement input generated by quoFEM and therefore must be considered together. Details on how to run this example using the quoFEM desktop can be found here[]. In this document three complementary notebooks are discussed that correspond to each of the steps mentioned above. A link to each notebook is included at the beginning of each section.
To connect SimCenter applications and Jupyter notebooks in DesignSafe, it's essential to ensure that all required tools are accessible from both frameworks. The SimCenter's vision is wellaligned with this concept and offers all of the necessary functionality through backend applications installed in DesignSafe that can be accessed via Tapis apps. This is schematically shown in Fig 2. Additionally, all SimCenter workflows are stored in JSON files that represent all steps in a workflow. This file is readily accesible using a JSON parser.
Fig.2  Running quoFEM analysis using remote computing resources at DesignSafe
In quoFEM, the workflow data is stored in a tmp.SimCenter folder that can be accessed from the quoFEM desktop, sent to an HPC, or accessed from a notebook. For the examples presented in this document, the information included in this folder is sufficient to run all cases.
In order to facilitate the discussion of each notebook, it is helpful to first identify common aspects that are present in all workflows run from Jupyter. These include:
 Setup Tapis App job
 Display job workflow
 Run Tapis job
 Postprocess results
Instructions for (1) setting up and (2) running Tapis jobs can be found here. These steps are generally applicable for launching any Tapis app from a Jupyter notebook in DesignSafe. The most significant step in this process is determining the appropriate Tapis app to utilize. To perform uncertainty quantification within SimCenter backend applications, the following Tapis app is used:
#Select tapisapp
app_name = 'simcenteruqfrontera'
app_id = 'simcenteruqfrontera4.0.0u4'
storage_id = 'designsafe.storage.default'
# Get Tapis app
app = ag.apps.get(appId=app_id)
Postprocessing of results is specific to the problem being solved and can be achieved using Python to access output data stored in archived files. For this purpose it is important to identify the location of data files. This is explained here.
Displaying a quoFEM job workflow is useful for understanding the data and steps followed in the workflow. These include: simulation tools, input variables, UQ methods used, and remote directories/folders. An schematic of a typical quoFEM workflow is shown in Fig. 3.
Fig.3  Elements of quoFEM workflow (only relevant elements for launching notebook from JupyterHub
To display the JSON file the IPython.display module can be used:
python
# Display Workflow in JSON file
import IPython.display
import json
jsonPath = os.getcwd()+"/tmp.SimCenter/templatedir/"+parameters["inputFile"]
with open(jsonPath) as f:
jsonInfo = json.load(f)
IPython.display.JSON(jsonInfo)
To modify the workflow, the user can either manually change the workflow files within the tmp.SimCenter folder or regenerate the workflow using quoFEM. Regenerating the workflow using quoFEM is the preferred approach, as the quoFEM desktop is specifically designed to facilitate workflow creation. On the other hand, Jupyter notebooks offer more flexibility in terms of postprocessing, generating plots, and manipulating data.
Step 1 – Global Sensitivity Analysis
The PM4Sand constitutive model has 24 parameters. Among them, apparent relative density Dr, shear modulus coefficient Go, and contraction rate parameter hpo, are known to be important for predicting liquefaction responses [2]. Therefore, these three parameters theta = {Dr, Go, hpo} are considered in the UQ analyses and their prior distributions are assumed to be uniform distributions with the ranges shown in Table 1. These prior distributions shall capture a plausible wide range that includes all possible parameter values for the target soils. The experimental data will be used to constrain this wide range to the domain that best describes the behavior exhibited by the specimen during the experiments.
Table 1.  Prior distributions of PM4Sand parameter
Parameter  Distribution  Range 

Dr  Uniform  0.10.6 
Go  Uniform  200  2000 
hpo  Uniform  0.01  5 
The sensitivity analysis is performed for a simulation model that reproduces the CyDSS test shown in Figs. 4 and 5. The output quantity of interest is the number of cycles until the onset of liquefaction (denoted as Y). The onset of liquefaction is defined as the time step when the shear strain shown in Fig. 4 exceeds 3.5%. Liquefaction capacity is affected by the initial shear stress typically characterized by the cyclic shear stress ratio (CSR; i.e., ratio of horizontal cyclic shear stress to vertical consolidation stress). In this sensitivity analysis, a CSR of 0.175 is considered. Two variancebased global sensitivity indices are evaluated:
where theta_i is the parameter of interest ( i.e., one of the {Dr, Go, hpo} ) , theta~i denotes the other two parameters, E_X[.] and Var_X[.] denote mean and variance of function over X, respectively, and the vertical bar denotes ‘conditional on’. The former index, called the maineffect index, quantifies how much of the variance of Y is attributed to the parameter theta_i, while the latter index, called the totaleffect index, also considers the joint contributions of theta_i and other parameters [3].
Fig.4  Single element FE model used in sensitivity analysis and Bayesian calibration
Fig.5  (a) simulated cyclic stressstrain curve; (b)stress path during the simulated cyclic direct simple shear test; (c) evolution of pore water pressure ratio during the simulated CyDSS test
The sensitivity analysis is performed using the algorithm in Weirs et al. (2012) through the Dakota engine that interfaces with quoFEM [3]. 2500 simulations were performed using the prior distributions in Table 1. The resulting sensitivity is shown in Fig. 6(a) which indicates that Dr is the dominating parameter for the response Y. This is also confirmed by inspecting the scatter plot of Fig. 6(b): Dr (horizontal axis) demonstrates a stronger influence on the output (vertical axis) compared to the influence of the other parameters shown in (c) and (d). Based on this, we can expect that the CyDSS observations will help constrain the uncertainty in Dr, while the reduction of uncertainty in hpo and Go will be relatively limited. Additionally, different types of experiments would be needed to better characterize those other parameters.
Fig.6  (a) Sensitivity analysis results for the critical number of cycles given CSR = 0.172; (b)– (d) Individual inputoutput scatter plots
Step 2 – Bayesian Parameter Calibration
Consider now the observations of the CyDSS experiment in Table 2, that are publicly available on the DesignSafe data depot [4]. We assume that the observed count of cycles at different CSR values, denoted as Y_i^m (i = 1,…,6), is given by the simulation model predictions and an added Gaussian noise. The latter captures various inaccuracies such as inherent uncertainty in the phenomenon, the imperfection of our simulation model, and measurement error. Given the above assumptions, we can denote the relationship between the data and model prediction, Y_i (theta), as
where noise epsilon_i is assumed to have zeromean and unknown variance sigma^2_{epsilon,i}. Given the six measurement values, we can use a Bayesian approach to evaluate the posterior distribution of the parameters of PM4Sand and the unknown noise variances:
where p(∙) denotes the (joint) probability distribution, and c is the normalization constant that ensures the area under the posterior distribution is one. From Eq. (3),
p(Y^m_i  theta, sigma^2) is a Gaussian distribution with mean Y_i(theta) and variance of sigma^2. The prior distribution of theta is in Table 1. Following best practices, inverse Gamma priors with the shape parameter alpha = 3 and scale parameter beta = 2 are introduced for the sigma^2 measurement variances [5]. The posterior sample of theta in this example is obtained using the transitional Markov chain Monte Carlo (TMCMC) sampling technique [6] that is available in quoFEM through the UCSDUQ engine. This is an expensive calculation that greatly benefits from the available HPC resources at DesignSafe.
Table 2. Cyclic direct simple shear (CyDSS) test experimental data
Cyclic Shear stress ratio (CSR)  Number of cycles to onset of liquefaction 

0.105  26 
0.105  21 
0.130  13 
0.151  5 
0.172  4 
0.200  3 
Fig.7  Comparison of calibrated model predictions and experimental data
Figure 7 compares the experimental data with the calibrated model predictions of the loadcycle counts, while Fig. 8 shows the calibrated parameter sample from the joint posterior distribution. Figure 8 shows that uncertainty in all variables is reduced by calibrating to the observed data, but the reduction was most apparent in Dr. This is in line with our expectations from the earlier sensitivity analysis. The results also highlight a strong dependency between Dr and hpo, indicating that multiple combinations of Dr and hpo produce nearoptimal solutions. None of these features are captured by a deterministic estimator that results from a conventional errorminimizing optimization approach (e.g., red diamond marker shown in the same figure). It is also important to recognize that a nonnegligible amount of uncertainty remains in the parameter estimates, and this produces substantial uncertainty in the model predictions. The dark blue bounds in Fig. 7 show the level of uncertainty in the estimated number of cycles to liquefaction, but this simulation model was prepared to reproduce the experimental setup. When the calibrated constitutive model is applied in another simulation, the responses can exhibit different scales of uncertainties.
Fig.8  PM4Sand model parameters sampled from the joint posterior distribution
Step 3 – Forward Propagation
A forward propagation analysis is helpful to characterize uncertainties in a simulation model. For this purpose it is good practice to run such an analysis and characterize the effect of uncertainties on applicationspecific quantities of interest before practically applying these parameter values in a simulation for decision making.
The obtained samples of the soil parameters in Fig. 8 are used to predict the uncertainty in the lateral spreading response of a site subjected to an earthquake (Loma Prieta Gilroy Array #2) with peak ground acceleration of 0.37 g. The soil column model shown in Fig. 9 is introduced in which the liquefiable layer in the middle is modeled using PM4Sand and the other parts are assumed to remain elastic throughout the shaking.
Fig.9  Schematic of 1D soil layer with liquefiable soil used in the forward propagation analysis.
The results of 500 simulations are shown in Fig. 10. The mean and standard deviation of the residual displacement at the surface level (6 m) are 0.24 m and 0.02 m, respectively. Depending on the application, the uncertainty in these results can be considered reasonably low. The sample of the predictive distribution shown on the top of the vertical profile can further be utilized in reliability and risk assessment workflows
Fig.10  Predicted earthquake response of soil column
References

McKenna, F., et al.: NHERISimCenter/quoFEM: Version 3.0.0. Zenodo (2022).

Boulanger, R.W., Ziotopoulou, K.: PM4Sand (Version 3.1): A sand plasticity model for earthquake engineering applications. Department of Civil and Environmental Engineering, University of California, Davis, Davis, CA, Report UCD/CGM17/01 (2017).

Weirs, V.G., et al.: Sensitivity analysis techniques applied to a system of hyperbolic conservation laws. Reliab. Eng. Syst. Saf. 107, 157–170 (2012).

Morales, B., Humire, F., Ziotopoulou, K.: Data from: Cyclic Direct Simple Shear Testing of Ottawa F50 and F65 Sands, 1 February 2021. Distributed by Design SafeCI Data Depot. https://doi.org/10.17603/ds2eahz9466. Accessed 28 June 2021.

Ching, J., Chen, Y.C.: Transitional Markov chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging. J. Eng. Mech. 133(7), 816–832 (2007).

Carlin, J.B., Vehtari, A., Stern, H.S., Rubin, D.B., Gelman, A., Dunson, D.B.: Bayesian Data Analysis, 3rd edn. Taylor & Francis, United Kingdom (2014).
Citation and Licensing

Please cite Aakash B. Satish et al. (2022) to acknowledge the use of resources from this use case.

Please cite SangRi Yi et al. (2022) to acknowledge the use of resources from this use case.

Please cite Rathje et al. (2017) to acknowledge the use of DesignSafe resources.

This software is distributed under the GNU General Public License.