Sensitivity analysis, calibration and validation of a snow indentation model

Available online at www.sciencedirect.com

Journal of Terramechanics Journal of Terramechanics 49 (2012) 315–324 www.elsevier.com/locate/jterra

Sensitivity analysis, calibration and validation of a snow indentation model Jonah H. Lee ⇑, Daisy Huang Department of Mechanical Engineering, University of Alaska Fairbanks, Fairbanks, AK 99775-5905, United States Received 12 September 2012; received in revised form 23 November 2012; accepted 4 December 2012

Abstract Quantification of the mechanical behavior of snow in response to loading is of importance in vehicle-terrain interaction studies. Snow, like other engineering materials, may be studied using indentation tests. However, unlike engineered materials with targeted and repeatable material properties, snow is a naturally-occurring, heterogeneous material whose mechanical properties display a statistical distribution. This study accounts for the statistical nature of snow behavior that is calculated from the pressure-sinkage curves from indentation tests. Recent developments in the field of statistics were used in conjunction with experimental results to calibrate, validate, and study the sensitivity of the plasticity-based snow indentation model. It was found that for material properties, in the semi-infinite zone of indentation, the cohesion has the largest influence on indentation pressure, followed by one of the the hardening coefficients. In the finite depth zone, the friction angle has the largest influence on the indentation pressure. A Bayesian metamodel was developed, and model parameters were calibrated by maximizing a Gaussian likelihood function. The calibrated model was validated using three local and global confidence-interval based metrics with good results. Ó 2012 ISTVS. Published by Elsevier Ltd. All rights reserved. Keywords: Indentation; Snow; Drucker-Prager; Validation; Calibration; Sensitivity; Bayesian; Metrics; Bias; Surrogate

1. Introduction A key component in the modeling of vehicle-terrain interactions for soils and snow [1–3] is the pressure-sinkage relationship of terrain material obtained using indentation testing. Naturally occurring terrain materials have been categorized as random heterogeneous materials [4] such that their properties should be treated statistically. Recent efforts in statistical modeling of vehicle-snow interaction include the interval analysis approach in [5], a metamodeling approach in [6], and a polynomial chaos approach in [7]. The pressure-sinkage relationship used in these efforts were, however, empirically-based. Recently, a snow indentation model in [8] was developed based on plasticity theory such that pressure-sinkage

⇑ Corresponding author. Tel.: +1 907 474 7136.

E-mail address: [email protected] (J.H. Lee).

curves have a physical basis, which gives an improvement over the empirical nature of previous research work. The model uses only a few fundamental material properties such as the cohesion, friction angle and hardening parameters. However, due to limited test data, the material properties of the model were assumed in an ad hoc fashion without consideration of the statistical variations of the pressure-sinkage curves of natural snow. Estimating parameters of a physical model against test data is a statistical process that is usually difficult since it belongs to the class of inverse problems. Indeed, estimating mechanical properties from indentation tests for engineering materials in general poses as an inverse problem. Parameter estimation is also called calibration and is intimately related to the validation of the model. Validation of engineering and scientific models has drawn much attention from academia and industry in recent years resulting in common terminology (e.g., [9]) for validation and related statistical frameworks [10,11]. One accepted

0022-4898/$36.00 Ó 2012 ISTVS. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jterra.2012.12.001

316

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

definition of validation is ‘the process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model’ [9]. Characterization of uncertainties of model and data is an integral part of the validation process. Toward this end, advancement of flexible and reasonably rigorous statistical frameworks such as [10,11] has been made to address the issues of sensitivity, calibration of parameters and validation of models including many applications to road-load interaction in automotive engineering [11–13]. In addition, quantitative validation metrics have been an active research area that provide a more rigorous statistical assessment of the agreement between test results and model predictions [14]. Although headway has been made in the statistical modeling of vehicle-snow interaction, no statistically rigorous efforts have been made to validate the various models developed recently. This paper applies recent results in the field of statistics, to study the sensitivity of the snow indentation model, to calibrate fundamental material properties of the model using newly obtained experimental results, and to validate the model using calibrated material properties and several validation metrics. The paper is organized as follows. Section 2 discusses background in the snow indentation model, statistical methodology in global sensitivity analysis, Bayesian metamodel, calibration, and confidence-interval based validation metrics. Section 3 presents new snow indentation tests. Section 4 discusses results, and Section 5 follows with discussion and conclusions. 2. Background and methodology In the following, we first summarize the snow indentation model where material constants are to be calibrated against test data toward the validation of the model. We then discuss the statistical methods and models used. It should be noted that there are two types of model discussed in this paper, one is the physical snow indentation model, and the other is the statistical model. Unless expressed explicitly, ‘model’ by itself means the snow indentation model. 2.1. Snow indentation model Three deformation zones have been approximately identified in [8]. Zone I is a small, initially linearly elastic region. It is followed by zone II, a strain-hardening region where the pressure bulb developed underneath the indenter has not yet reached the bottom of the snow cover, i.e., it is a zone of semi-infinite depth. Zone III, a finite-depth zone, is a region where the pressure bulb has reached the bottom of the snow cover. The material properties for the plasticity indentation model to be calibrated in this paper are summarized below. Full details of the model can be found in [8] and references therein.

A simple Drucker-Prager yield function was used to develop the plasticity solutions for the indentation model using two material parameters: the cohesion (pd) and the friction angle (b) which can be related to the absolute tensile and compressive strengths (T and C) as 1 C  C tan b  pd ¼ 0 3 1 T þ T tan b  pd ¼ 0 3

ð1Þ ð2Þ

The hardening of the snow can be expressed by the location of the cap of the Drucker-Prager model, pa, which is parameterized as   3  log10 pa ¼ c1  c2 exp pv  c3 pv ð3Þ p



where pv ¼ 3kk is the volumetric plastic strain, c1, c2 and c3 are constants. 2.2. Statistical methods In this section, we present background of the essential ingredients of the methodology used in this paper: global sensitivity analysis, Bayesian metamodel, calibration, and validation. These include the Gaussian maximum likelihood used in calibration, Gaussian processes used for metamodel and sensitivity analysis, and validation metrics. The methods used are also summarized as flow charts in Figs. 2–4 which will be discussed in detail in Section 4. 2.2.1. Gaussian likelihood Maximum likelihood estimation (MLE) is a common method that estimates the parameters of a statistical model such that there is a maximum probability of the observed data based on the estimated parameters [15]. When the statistical model uses a Gaussian distribution for the random variable x with unknown mean l and standard deviation r, one has ! 2 1 ðx  lÞ f ðx; l; r2 Þ ¼ pffiffiffiffiffiffi exp  ð4Þ 2r2 r 2p To estimate the mean and standard deviation of a sample size of N (xi, i = 1, . . . , N) using MLE, the probability of the data set is first expressed as the product of the probabilities of each point, i.e., the likelihood function L(l, rjx) ! 2 N Y 1 ðxi  lÞ pffiffiffiffiffiffi exp  Lðl; rjxÞ ¼ ð5Þ 2r2 i¼1 r 2p The unknown parameters (l and r) are then obtained by maximizing L which is equivalent in maximizing log(L) such that the normal log-likelihood function is ! 2 N X N ðxi  lÞ 2 logðLÞ ¼  ðlog r þ logð2pÞÞ  ð6Þ 2 2r2 i¼1

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

(a)

317

P

d

H

W

(b)

Fig. 2. Flow chart of global sensitivity analysis.

y, as well as the contribution of the interactions of the parameters to y. In this paper, we adopt the approach of global sensitivities based on variance decomposition as opposed to a local, derivative based approach; see Section 1.2.1 of [16] for more details. The sensitivity measure used in this paper is Tj for the jth parameter (j = 1, . . . , r) defined in [17] as Fig. 1. Indentation test: (a) schematic of the test: diameter of the indenter d = 12.7 mm, diameter of the container W = 46 mm, depth of the snow H = 20 mm, P = force, and (b) experimental set up.

2.2.2. Global sensitivities Analysis of the global sensitivities of a model is such that it does not rely on local derivatives of the function against parameters of interest, but rather uses design of experiments to explore the entire space to efficiently figure out sensitivity [16]. This is especially useful for nonlinear functions such as the snow indentation model of interest in this paper. A statistical model for sensitivity studies can be expressed as y ¼ f ðxÞ þ 

ð7Þ

where y is the scalar output of the physical model, x = [x1, x2, . . . , xr] is a vector of r parameters, and  is an error term due to numerical computation with finite precision.  discussed here should not be confused with strain of snow indentation which uses the same symbol.  will be assumed to be small and will not be subsequently considered. The sensitivity of the parameters x of the model refers to the relative contribution of each of the parameters to the output

Tj ¼

EðVar½f ðxÞjxðjÞ Þ Varðf ðxÞÞ

ð8Þ

where E() and Var() are the expectation and variance of a function, respectively; x(j) = [x1, . . . , xj1, xj+1, . . . , xr] represents all the parameters except xj. Tj defined here accounts for the contributions of parameter xj and its interaction with other parameters to y. In addition to obtaining the sensitivity measure, it is also desirable to obtain its confidence interval, e.g., via a bootstrapping (resampling) technique [17]. To obtain the sensitivity measure and its confidence interval, many model runs are needed, which can be daunting for even relatively simple models. Consequently, metamodels are oftentimes developed to speed up the calculations. This will be discussed in the next section. In the literature, metamodel is also called surrogate, emulator or simulator. When metamodels are used for sensitivity analysis, the sensitivity measure Tj will be designated as Tb j . 2.2.3. Bayesian metamodel There are two major views of statistics: the frequentist view and the Bayesian view, with the former based on hypothesis testing and the latter based on inference of

318

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

Fig. 4. Flow chart of the calibration process.

is considered as an unknown function modeled statistically as a Gaussian process by its mean function m() where m(x) = E[f(x)], and its covariance function c() where c(x, x0 ) = cov(x, x0 ). Following [20], for a random function f() and a set of points [x1, . . . , xn] in its domain where each xi, i = 1, . . . , n, has r parameters [x1, . . . , xr], the random vector [f(x1), . . . , f(xn)] is assumed to be a multivariate normal distribution with mean T

E½f ðxÞjb ¼ hðxÞ b; Fig. 3. Flow chart of the development of the Time Input (TI) metamodel.

(posterior) probabilities with assumed or known prior knowledge using Bayes’ theorem [18]. See [19] for detailed discussions on these two views where the Bayesian approach was considered to be more flexible and practical. A common approach for a metamodel is to use a Bayesian statistical method where the Gaussian processes are used as the prior probability distribution [10,20]. Gaussian process [21], in contrast to Gaussian distribution which is a distribution over random variables, is a stochastic process and is defined as the probability distribution over functions f(x) such that the set of values of f(x) evaluated at arbitrary points x1, . . . , xn will have a joint Gaussian distribution, i.e., it is a distribution over functions. The Gaussian process has been used widely in many fields of science and engineering including application in the form of Gaussian random field used effectively in characterizing the 3D microstructure of snow [22]. Given data from the simulations of the physical model, used as ‘training’ data, the posterior probability distributions of the model output can be obtained, which is model dependent. In particular, the output from a physical model

ð9Þ

where b is a vector of unknown coefficients, and h() the q unknown regressor functions of the r parameters x = [x1, . . . , xr]; superscript T denotes the transpose of a vector or matrix. The covariance is expressed as cov½f ðxÞ; f ðx0 Þjr2  ¼ r2 cðx; x0 Þ

ð10Þ

0

where c(x, x ), the correlation function between the inputs x and x0 , is assumed to be cðx; x0 Þ ¼ exp½ðx  x0 ÞT Bðx  x0 Þ

ð11Þ

where B is a diagonal matrix of positive ‘roughness’ parameters li, i = 1, . . . , r; in statistical literature, the unknown r2 and roughness parameters li are called hyperparameters. Once the mean function and the covariance function are chosen, one obtains a set of design points [x1, . . . , xn] in parameter space as input to the physical model f(); these points are called the experimental design or the design matrix (with dimensions n  r). In this paper, the design points are generated using the popular Latin Hypercube Sampling (LHS) approach [23] which obtains distributions of parameters from a given multidimensional distribution. The outputs from the design matrix, d = [d1, . . . , dn], are obtained by running the physical model n times. In practice, d are typically appended to the design matrix; for

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

convenience, the design matrix augmented by d will also be referred to collectively as the design matrix. Using Bayes’ theorem, the posterior mean and covariance functions can then be derived and given in Eqs. (3) and (4) in [20] and Eqs. (2) and (3) in [24]. The quality of the metamodel is assessed via cross-validation methods such as the leave-one-out analysis where data from n  1 model runs is used to predict the model data of the remaining run. The metamodel thus obtained can be used to predict new outputs as an emulator of the original physical model f() given a new set of inputs. The Bayesian metamodel approach effectively converts a deterministic physical model into a stochastic model via Bayes’ theorem conditioned with model data. For pressure–displacement relationship, given parameters x, the output of the model is a curve rather than a single pressure. Displacement, or strain, can be considered as a time-like parameter (t). The metamodel discussed previously can then be extended to include an extra input t taking values in 1, . . . , T yðtÞ ¼ f ðx; tÞ þ 

ð12Þ

to emulate a dynamic simulator resulting in a Time Input (TI) simulator discussed in [25]. The design matrix now has dimensions of nT  r. 2.2.4. Calibration For calibration of model parameters, we used the Gaussian likelihood method. Since the snow indentation model is an approximation of the real phenomenon, there is a bias (or model inadequacy [10]) between the model and test data that needs to be taken into consideration. Let the test data be represented by Yt,i, t = 1, . . . , T, and i = 1, . . . , N where N is the number of replications of test data. The bias between the model f(x, t) and test data is approximated as [13] Bt ¼ Y t  f ðx; tÞ ð13Þ where at time t, Y t is the average of the test data over N replications, and f ðx; tÞ is the average of the surrogate model output over n sets of input parameters x. We then assume the following statistical model for calibration [13] Y t;i ¼ Bt þ f ðx; tÞ þ Eðx; tÞ þ t;i

ð14Þ

where f is the mean function of the surrogate model, E is the surrogate error term having a multivariate Gaussian distribution with a zero mean vector and covariance expressed in Eq. (10),  is an observational random error term with a Gaussian distribution. The likelihood function in Eq. (5) then can be generalized as LðxjYÞ ¼

N Y

1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi detðXÞ i¼1   1  exp  ðY t;i  Bt  fðx; tÞÞT X1 ðY t;i  Bt  fðx; tÞÞ 2 ð15Þ

319

where X is the covariance function of the metamodel in Eq. (12); constant terms are ignored. The calibration parameters x are then obtained by maximizing Eq. (15) or its logarithm. 2.2.5. Validation Once the calibration parameters have been obtained via the Gaussian maximum likelihood approach, they are substituted into the metamodel of Eq. (12) to predict the mean, variance, and confidence interval (called credible interval using Bayesian terminology) of the output of interest. An evaluation of the prediction constitutes the validation of the model. Validation of models can be conducted qualitatively and/or quantitatively. The former often involves expert opinion of, say, the comparison of prediction against test data graphically. To validate models quantitatively, validation metrics are needed which are more stringent than qualitative validation. In this paper, we use three confidence-interval (CI) based validation metrics where the first two were proposed in [14]. The first validation metric is based on the difference in the means between model prediction and test data, and the standard deviation of either test data or model prediction as follows. e at a given time (strain) We define the estimated error E as the difference between the predicted mean of the calibrated model f(xcalibrated, t), and the mean of test data over all replications (Y t in Eq. (13)): e ¼ f ðxcalibrated ; tÞ  Y t EðtÞ

ð16Þ

At a given time t, the confidence interval (CI) of test data is constructed based on the standard deviation of the test data s(t) following standard statistical procedure. A plot e of EðtÞ against the test CI gives a comparison of the estimated error with the error of the test data as a function of t; if the former is smaller than the latter, then one may consider the model to be adequate and one may want to improve the test; otherwise, one may want to improve the model if needed depending upon the intended purpose. A variation of this metric is to construct the CI of model prediction as opposed to that of the test data. Since we are using a Bayesian surrogate model, the CI is constructed e against the model from the surrogate model. A plot of EðtÞ CI gives a comparison of the estimated error with the error of the model prediction as a function of t; if the former is smaller than the latter, then one may consider the model to be adequate for predictive purpose; otherwise, one may want to improve the model. Since the above metric is a function of t, it is possible that the model performs adequately for part of t, but not for others. Consequently, this metric gives a fine-grained (local) evaluation of the quality of the model and can provide guidance on the improvement of the model. The second and third validation metrics are global, i.e., a measure for the entire time interval. They are normalized to provide a high-level summary of the errors in percentage.

320

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

For the second metric, the average relative error metric is first defined as [14]: Z tu EðtÞ f ðxcalibrated ; tÞ  Y ðtÞ 1 e dt ð17Þ ¼ Y t t Y ðtÞ t

u

avg

l

tl

where tu and tl are the upper and lower limits of t. This quantity can be considered as the absolute mean error between test and model normalized by the mean test error. The associated average relative confidence indicator is defined next as: Z xu sðtÞ CI tp;m ¼ p ffiffiffiffi ð18Þ Y ðtÞ dt Y ðt  t Þ N t avg

u

l

xl

where tp,m is the t-statistic of probability p, and degree of freedom m, N is the total number of test data replications. This indicator can be considered as an approximate normalized confidence interval of test data. For this metric, similar to the first metric, if the average relative error is smaller than the average relative confidence indicator, then one may consider the model to be adequate; otherwise, the model may need improvement. The third validation metric is the coverage defined as the percentage of the mean test data points that fall within, say, the 95% CI of the model prediction; this is a measure of the predictability of the model. A larger coverage means a better model.

in a clear, insulative enclosure maintained at a temperature of 20 °C using a low-temperature forced air stream. During an indentation test, snow was placed into a small sample holder underneath the indenter pin, and the indenter pin was manually moved down until it was close to, but not touching, the snow surface. Then the pin was moved down at the selected rate while resisting force was measured by the load cell, which was coupled to the top of the pin. Resisting force versus indentation displacement was output by the tribometer. Tests were carried out with an indenter of 12.7 mm diameter and an indentation speed of 5 mm/s for a snow 20 mm deep. A 20-kg load cell with 0.01 N resolution was used. Two fine-grained snows were used in this paper, one was collected on December 29, 2009 and tested on September 14, 2010; the other was collected on November 16, 2010 and tested on March 18, 2011. 4. Results The implementation of the methodology discussed in Section 2.2 uses computer programs coded in the R language [27]. Global sensitivity analysis uses the R package in [28] which also employs metamodels to represent the physical model; the Bayesian metamodel uses the R package discussed in [24]. 4.1. Global sensitivity analysis

3. Experiments To validate the model, indentation tests were performed on virgin, natural snows. Snow was collected in the Fairbanks, Alaska area where air temperature, snow surface temperature, and snow pit temperature under the surface were recorded at each snow collection. The experimental procedures used in this paper follow those in [26] where much more details can be found. Snow was characterized in situ using a high-resolution snow micropenetrometer, which gave information about penetration resistance and average grain size. Then samples were removed and measured, weighed, and sieved on-site to obtain density and grain size distribution. Specimens were also collected and brought to the laboratory and stored in freezers which were maintained at 30 °C. There, they were additionally characterized via inspection with an optical microscope, and scanning using 3D X-ray MicroTomography. From these measurements, the snow densities were determined to range from about 170 to 240 kg/ m3, and the average grain sizes were estimated to be 1 mm. Indentation tests were performed using a CETR UMT tribometer shown in Fig. 1. The testing region was enclosed

For global sensitivity analysis, the ranges of parameters, selected by expanding the ranges used in [8], are shown in Table 1. The design points were generated using Latin Hypercube Sampling with a uniform distribution where the parameters are normalized in [0, 1]. Note that the strain  refers to the logarithmic strain. Out of the many metamodels available in [28], we used the Gaussian process metamodel. To gain physical insight of the indentation process, we examine the sensitivities for zone II (hardening or semi-infinite zone) and zone III (finite-depth zone) separately. The dimensions of the design matrix (n  r) for zone II are 149  7, and those for zone III are 199  7. Fig. 2 presents the procedure of sensitivity analysis as a flow chart. The sensitivity measure Tb j and its 95% confidence interval for the parameters are shown in Tables 2 and 3 for zone II and zone III, respectively. In interpreting Tb j , it should be noted that Tb j accounts for the contribution of the j-th parameter and its interaction with all other parameters; consequently, the sum of Tb j is larger than 1. For zone II, Table 2 shows that pd has the largest influence on indentation pressure and explains between 58.5%

Table 1 Range of parameters for global sensitivity analysis. pd (kPa)

b (deg)

c1

c2

c3

Depth (m)



(1, 50)

(10, 50)

(0.4, 0.8)

(3, 5)

(0.1, 0.5)

(0.2, 0.6)

(0, 0.7)

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324 Table 2 Sensitivity measure and associated 95% confidence interval using the Gaussian process metamodel in [28] for zone II.

321

Table 4 Ranges of parameters to generate the design matrix; c1 = 0.66, c3 = 0.18; depth = 20 mm.

Tb

95% CI of Tb

Snow sample

pd (kPa)

b (deg)

c2

pd b c1 c2 c3 Depth 

0.676 0.068 0.010 0.106 0.000 0.028 0.220

(0.585, 0.757) (0.049, 0.110) (0.000, 0.043) (0.055, 0.130) (0.000, 0.028) (0.019, 0.056) (0.164, 0.286)

December 2009 November 2010

(0.1, 6.3) (0.1, 7.3)

(10, 20) (10, 20)

(4, 30) (4, 30)

0.000 0.217 0.045 0.012 0.032 0.234 0.732

(0.000, 0.039) (0.156, 0.261) (0.013, 0.090) (0.000, 0.050) (0.000, 0.069) (0.214, 0.300) (0.670, 0.774)

10000 8000 6000

pd b c1 c2 c3 Depth 

Pressure (Pa)

95% CI of Tb

4000

Tb

2000

Parameter

0

Table 3 Sensitivity measure and associated 95% confidence interval using the Gaussian process metamodel in [28] for zone III.

12000

Parameter

0.0

and 75.7% of the indentation pressure with a 95% confidence; other parameters can be interpreted similarly. The strain  has the second largest influence on indentation pressure which is not surprising since indentation pressure increases with strain. Next in importance is c2 and b with the rest of the parameters of little importance. That the depth has little influence over the indentation pressure for zone II is also expected since it is a semi-infinite zone. In summary, the significance of material properties for zone II, in decreasing order of significance is pd, c2, b. For zone III, Table 3 shows that the significance of parameters is ranked by  > depth > b. For a finite-depth zone, indentation pressure can increase dramatically with deformation which explains the significance of ; similarly for depth. It is interesting to note that b becomes the most important material property for zone III with pd playing almost no role in indentation pressure, and with c2 playing an insignificant role–almost the opposite of the results for zone II. 4.2. Bayesian metamodel The Bayesian metamodel is used for the calibration, prediction and validation of the physical model; the steps in building the metamodel are given in Fig. 3. Based on the results from Section 4.1, only a reduced set of the parameters, pd, b, and c2 will be allowed to vary with the rest fixed. A prerequisite of the values of parameters for building the metamodel is that the output from the physical model due to the parameters need to cover the domain of interest, i.e., the test data in our case. Using the Time Input model from [25] in Eq. (12), the design matrix is again created using Latin Hypercube Sampling with initial dimensions of n  r where n = 50, and r = 3. For each model run, one

0.1

0.2

0.3

0.4

Strain Fig. 5. Indentation pressure from model data (dashed lines) and test data (solid lines) for November 2010 snow.

uses T values of the model output such that the final design matrix has dimensions of nT  r where T = 5, i.e., 250  3. The parameters are adjusted until the test data are covered by the model output. During the step in creating the design matrix, it was found that a slight modification of the original snow indentation model for stage I was needed to improve the coverage of test data by model data. In particular, the strength at the end of zone I is quite small and was set to a small threshold value of 100 Pa in all subsequent steps. The ranges of parameters thus arrived at are given in Table 4 for both snow samples. Model output and the test data are shown in Fig. 5 for the November 2010 snow with N = 6 data replications, as an example. After the design matrix is established, the hyperparameters of the metamodel need to be estimated. Initial values are assumed and an MLE procedure was used in [24] to find optimized values. A leave-one-out cross-validation was then conducted to evaluate the quality of the metamodel shown in Fig. 6. The closer the data points to the main diagonal, the better and the performance of the metamodel seems adequate. 4.3. Calibration, prediction and validation 4.3.1. Calibration After the Bayesian metamodel was developed, calibration of the parameters was conducted by maximizing the Gaussian likelihood function in Eq. (15). The steps for calibration are summarized in Fig. 4. The calibrated con-

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

6000 2000

4000

Pressure (Pa)

8000 6000 4000

0

2000

Predicted pressure (Pa)

10000

8000

12000

322

0

0.0

0.1

0.2

0.3

0.4

Strain 4000

6000

8000

10000

12000

Actual model output (Pa) Fig. 6. Cross-validation results of the surrogate model for December 2009 snow.

Fig. 7. Predicted indentation pressure (solid line) with 95% confidence interval (dashed lines) using calibrated material constants compared with experimental data for December 2009 snow.

pd (kPa)

b (deg)

c2

T (kPa)

C (kPa)

December 2009 November 2010

3.4 3.9

11.9 11.7

13.6 25.6

3.1 3.6

3.6 4.7

0

stants for two snow samples are given in Table 5 where the compressive (C) and tensile strength (T) have been converted from pd and b using Eq. (1) for reference. The calibrated material properties of the two snows are similar even though they were collected about a year apart and tested at different times.

2000

Snow sample

Pressure (Pa)

8000

Table 5 Calibrated material parameters.

6000

2000

4000

0

4.3.2. Prediction Predictions of the mean and variance of the metamodel are obtained by substituting the calibrated parameters to the metamodel with results shown in Figs. 7 and 8 where the solid red line indicates the predicted mean and the dashed blue lines indicate the 95% confidence interval. The predictions appear to be good qualitatively especially in view of the large variations of the test data which is the rule rather than the exception. A more rigorous evaluation of the predictions is given below using validation metrics introduced previously. 4.3.3. Validation Validation results using the first (time-dependent) metric are shown in Fig. 9a and b for the December 2009 snow. Fig. 9a shows that the estimated error lies within the CI of the test data indicating that the calibrated model has good performance over the entire range of strain. The estimated error increases with strain, following the trend of the CI of the test data, as a consequence of the calibration procedure.

0.0

0.1

0.2

0.3

0.4

Strain Fig. 8. Predicted indentation pressure (solid line) with 95% confidence interval (dashed lines) using calibrated material constants compared with experimental data for November 2010 snow.

Fig. 9b shows that the CI of the model prediction, constant over the entire range of strain, is almost half of that of the CI of the test data. Except at large strains, the estimated error still lies within the CI of the model. However, improvement of the model at large strains could be considered if one is concerned about the prediction of the model at that level of strain. Fig. 10a and b show the results in applying the first validation metric for the November 2010 snow. For both figures, the estimated errors lie within the CI of the test data and model prediction. Consequently, the results are slightly better than those for the December 2009 snow. The results in applying the two global metrics to the two snows are shown in Table 6.

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

0.10

0.15

0.20

0.25

0.30

0.35

1000 500 0 −500 −1500 −1000

Estimated error and experiment confidence interval (Pa)

2000 1000 0 −1000 −2000

Estimated error and experiment confidence interval (Pa)

0.05

1500

(a)

(a)

0.05

0.40

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.30

0.35

0.40

Strain

Strain

0.10

0.15

0.20

0.25

0.30

0.35

200 0 −200 −400

Estimated error and model confidence interval (Pa)

1000 500 0 −500 −1000 0.05

400

(b)

(b) Estimated error and model confidence interval (Pa)

323

0.05

0.40

0.10

0.15

0.20

0.25

Strain

Strain Fig. 9. (a) Estimated error and confidence interval for test, and (b) estimated error and confidence interval for model; December 2009 snow; mean error (solid lines), confidence interval (dashed lines).

Fig. 10. (a) Estimated error and confidence interval for test, and (b) estimated error and confidence interval for model; November 2010 snow; mean error (solid lines), confidence interval (dashed lines).

Since the relative error is much smaller than the relative confidence indicator for both snows, the performance of the calibrated model is good. These summary metrics show clearly that the December 2009 snow has a larger confidence indicator than the November 2010 snow. This is due to the top curve in Fig. 7 which is farther apart from other curves. The coverage for the December 2009 snow is slightly less than that for the November 2010 snow indicating the model for the latter is slightly better than the former. In addition, the 94–100% coverage of both snows shows that the calibrated model performs quite well for the indentation tests. To the best of our knowledge, this is the first time that fundamental mechanical properties of snow have been estimated from indentation tests with rigorous statistical prediction and validation of the indentation model that performs well.

Table 6 Results of global validation metrics. Snow type

Coverage (%)

Relative error

Relative confidence indicator

December 2009 November 2010

94 100

0.12 0.12

0.73 0.44

5. Discussion and conclusions In this paper, we used a few interrelated statistical methods to calibrate parameters of the snow indentation model, and infer statistically the pressure-sinkage relationship based on new test data with replications. The Bayesian metamodel is a key component of the methodology serving the following roles: (1) an efficient emulator of the physical

324

J.H. Lee, D. Huang / Journal of Terramechanics 49 (2012) 315–324

model for the calibration, prediction and validation of the model, and (2) as a non-intrusive stochastic model of the physical model such that statistical inference can be made. Global sensitivity analysis was conducted rigorously and efficiently for several purposes: (1) to gain physical insight of the influence of the multidimensional parameters to model output, (2) to select only influential parameters for calibration which utilizes the Bayesian metamodel, (3) to avoid the possibility of getting physically unreasonable results for the calibration process [13]. The results presented here are from the slightly-modified indentation model, as a consequence of applying the statistical methodology during the design of experiment stage, and before the validation metrics were applied, demonstrating another aspect of the usefulness of the methodology employed. In summary, our results indicated that, in terms of material properties, pd, c2 influence significantly the pressuresinkage relationship for zone II (semi-infinite hardening), whereas b has the most important contribution for zone III (finite-depth). The three interval-based validation metrics used indicated that the calibrated model is quite adequate for the two snows used in the paper. Acknowledgments The authors gratefully acknowledge the support of this work by the US Army TACOM Life Cycle Command under Contract No. W56HZV-08-C-0236, through a subcontract with Mississippi State University. This work was performed in part for the Simulation Based Reliability and Safety (SimBRS) research program. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the US Army TACOM Life Cycle Command. References [1] Wong JY. Theory of ground vehicles. 3rd ed. Wiley Interscience; 2001. [2] Muro T, O’Brien J. Terramechanics. Land locomotion mechanics. A.A. Balkema Publishers; 2004. [3] Lee JH. An improved slip-based model for tire-snow interaction. SAE Int J Mater Manuf 2011;4:278–88. http://dx.doi.org/10.4271/2011-010188. [4] Torquato S. Random heterogeneous materials: microstructure and macroscopic properties. 1st ed. Springer; 2001. [5] Lee JH, Liu Q, Mourelatos ZP. Simulation of tire-snow interfacial forces for a range of snow densities with uncertainty. SAE Trans J Mater Manuf 2006:408–18. [6] Li J, Mourelatos ZP, Lee JH, Liu Q. Prediction of tire-snow interaction forces using metamodeling. SAE congress. Paper 2007-011511; 2007.

[7] Li L, Sandu C, Lee J, Liu B. Stochastic modeling of tire-snow interaction using a polynomial chaos approach. J Terramech 2009;46:165–88. [8] Lee JH. A new indentation model for snow. J Terramech 2009;46:1–13. [9] ANSI/ASME V& V 10-2006: guide for verification and validation in computational solid mechanics. American Society of Mechanical Engineers; 2006. ISBN: 079183042X. [10] Kennedy MC, O’Hagan A. Bayesian calibration of computer models. J R Stat Soc: Ser B (Stat Methodol) 2001;63(3):425–64. [11] Bayarri MJ, Paulo R, Cafeo JA, Cavendish J, Lin C, Tu J. A framework for validation of computer models. Technometrics 2007;49(2):138–54. [12] Drignei D. Functional ANOVA in computer models with time series output. Technometrics 2010;52:430–7. [13] Drignei D, Mourelatos ZP. Parameter screening in dynamic computer model calibration using global sensitivities. ASME J Mech Des 134:081001. doi:10.1115/1.4006874. [14] Oberkampf D, Barone MF. Measures of agreement between computation and experiment: validation metrics. J Comput Phys 2006;217:5–36. [15] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the art of scientific computing. 3rd ed. Cambridge University Press; 2007. [16] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global sensitivity analysis: the primer. Wiley-Interscience; 2008. ISBN: 0470059974. [17] Storlie CB, Swiler LP, Helton JC, Sallaberry CJ. Implementation and evaluation of nonparametric regression procedures for sensitivity analysis of computationally demanding models. Reliab Eng Syst Safety 2009;94(11):1735–63. [18] Bayes T, Price R. An essay towards solving a problem in the doctrine of chance. By the late Rev. Mr. Bayes, communicated by Mr. Price, in a letter to John Canton, M.A. and F.R.S.. Philos Trans R Soc Lond 1763;1:370–418. [19] Bayarri MJ, Berger JO. The interplay of bayesian and frequentist analysis. Stat Sci 2004;19(1):58–80. [20] Oakley J, O’Hagan A. Bayesian inference for the uncertainty distribution of computer model outputs. Biometrika 2002;89(4):769–84. [21] Rasmussen CE, Williams CKI. Gaussian processes for machine learning. MIT Press; 2006. ISBN: 026218253X. [22] Yuan H, Lee JH, Guilkey J. Stochastic reconstruction of the microstructure of equilibrium form snow and computation of effective elastic properties. J Glaciol 2010;56:405–14. http://dx.doi.org/ 10.3189/002214310792447770. [23] McKay MD, Beckman RJ, Conover WJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 1979;21:239–45. [24] Hankin RKS. Introducing BACCO, an R bundle for bayesian analysis of computer code output. J Stat Softw 2005;14(16):1–21. [25] Conti S, O’Hagan A. Bayesian emulation of complex multi-output and dynamic computer models. J Stat Plann Infer 2010;140:640–51. [26] Huang D, Lee JH. Mechanical properties of snow using indentation tests: size effects. J Glaciol 59. doi:10.3189/2013JoG12J064. [27] R Development Core Team. R: A language and environment for statistical computing. R foundation for statistical computing, Vienna, Austria; 2005. . ISBN: 3-900051-07-0. [28] Storlie C. CompModSA: an R package for the sensitivity analysis for complex computer models. R package version 1.3.1; 2011. .