On the Need to Determine Accurately the Impact of Higher-Order Sensitivities on Model Sensitivity Analysis, Uncertainty Quantification and Best-Estimate Predictions

: This work aims at underscoring the need for the accurate quantification of the sensitivities (i.e., functional derivatives) of the results (a.k.a. “responses”) produced by large-scale computational models with respect to the models’ parameters, which are seldom known perfectly in practice. The large impact that can arise from sensitivities of order higher than first has been highlighted by the results of a third-order sensitivity and uncertainty analysis of an OECD/NEA reactor physics benchmark, which will be briefly reviewed in this work to underscore that neglecting the higher-order sensitivities causes substantial errors in predicting the expectation and variance of model responses. The importance of accurately computing the higher-order sensitivities is further highlighted in this work by presenting a text-book analytical example from the field of neutron transport, which impresses the need for the accurate quantification of higher-order response sensitivities by demonstrating that their neglect would lead to substantial errors in predicting the moments (expec-tation, variance, skewness, kurtosis) of the model response’s distribution in the phase space of model parameters. The incorporation of response sensitivities in methodologies for uncertainty quantification, data adjustment and predictive modeling currently available for nuclear engineering systems is also reviewed. The fundamental conclusion highlighted by this work is that confidence intervals and tolerance limits on results predicted by models that only employ first-order sensitivities are likely to provide a false sense of confidence, unless such models also demonstrate quantitatively that the second- and higher-order sensitivities provide negligibly small contributions to the respective tolerance limits and confidence intervals. The high-order response sensitivities to parameters underlying large-scale models can be computed most accurately and most efficiently by employing the high-order comprehensive adjoint sensitivity analysis methodology, which overcomes the curse of dimensionality that hampers other methods when applied to large-scale models involving many parameters.


Introduction
The state-of-the-art computational tool for sensitivity analysis, uncertainty quantification and "data adjustment" of model responses (effective multiplication factors, reaction rates) in reactor physics and reactor criticality safety is embodied in the module TSURFER [1] of the code system SCALE, developed at the Oak Ridge National Laboratory. The methodology underlying TSURFER is the so-called "generalized linear least squares adjustment (GLLSA)" procedure [2], which is limited to incorporating just the first-order sensitivities of the response with respect to the model parameters. The firstorder response sensitivities to the many thousands of model parameters (cross-sections, The distribution of neutrons in the PERP system is modeled using the standard form of the time-independent Boltzmann neutron transport equation, namely: where ( ) , Q E r is a spontaneous fission source, having the following expression: The neutron angular flux ( ) where n denotes the unit outward normal vector at any point s V ∈ ∂ r on the body's outer surface V ∂ , and f E denotes the highest neutron energy considered for the system under consideration. The terms on the left side of Equation (1) represent the neutron streaming and total collision loss processes, while the terms on the right side of Equation (1) represent the neutron scattering and fission production processes within the PERP benchmark.
The neutron flux is computed by solving numerically Equations (1) -(3) using the multigroup discrete ordinates particle and radiation transport code PARTISN [23], which solves the following multi-group approximation of Equation (1) where 0 μ ′ Ω Ω   , where d r is the radius of the PERP sphere and where the boundary condition shown in Equation (5) imposes no incoming particle flux. The vector α , which appears in the arguments of the various quantities in Equation (4), represents the "vector of imprecisely known model parameters," including atomic number densities, microscopic group cross-sections, etc. The detailed definition of the components of the vector α for the PERP benchmark is provided in [6][7][8][9][10][11], but this detailed definition is not used in this work. The convention used in PARTISN [23] for obtaining the multigroup transport equation shown in Equation (4) The group fluxes ( ) , g r ϕ Ω are computed numerically using the MENDF71X [24] 618-group cross-sections collapsed to 30 energy groups, in conjunction with a P3-Legendre expansion of the scattering cross-section, an angular quadrature of S256 and a fine-mesh spacing of 0.005 cm (comprising 759 meshes for the plutonium sphere of radius of 3.794 cm and 762 meshes for the polyethylene shell of thickness of 3.81 cm). The other quantities that appear in Equations (4) - (11) are described in detail in [6][7][8][9][10][11].
The response of interest for the PERP benchmark is the leakage of neutrons out of the sphere's outer surface. The leakage response, denoted as ( ) L α , is defined mathematically as the following integral over the outer surface of the spherical benchmark, comprising neutrons of all energies leaking through the surface: As shown in [6][7][8][9][10][11], the distribution of the PERP leakage in each energy interval is depicted in Figure 1. The value of the total leakage computed using Equation (1) for the PERP benchmark is 1.7648 × 10 6 neutrons/sec. The vector α , which appears in the arguments of the various quantities in Equations (4) - (12), represents the "vector of imprecisely known model parameters." For the PERP benchmark, the vector α comprises 21,976 imprecisely known (uncertain) model parameters, as follows: 180 group-averaged total microscopic cross-sections, 21,600 group-averaged scattering microscopic cross-sections, 120 parameters describing the fission process, 60 parameters describing the fission spectrum, 10 parameters describing the system's sources and six isotopic number densities. The exact and efficient computation of the 21,976 first-order and 482,944,576 second-order sensitivities of the leakage response to the PERP model parameters was reported in [6][7][8][9][10][11].
Thus, the third-order relative sensitivity,  The results presented in Table 2, reproduced (with permission) from [6], illustrate the impact of uncertainties in the total group cross-sections on the response's expected value and variance. The group total cross-sections were considered to be uncorrelated and normally distributed, having uniform relative standard deviations of 5% (moderate) and 10% (large), respectively. The impact of the other uncertain cross-sections and nuclear parameter underlying the PERP benchmark are presented in [6][7][8][9][10][11], where they are shown to accentuate the impacts presented in Table 2 and Figure 2, below. In Table 2, the expected value of the leakage response is denoted as has the following mathematical expression:  denotes the second-order contributions to the expected value of the leakage response and is defined as follows: In Equation (13), the quantity ( ) 0 L α represents the leakage response computed using the expected cross-section values; the superscript "U" denotes "uncorrelated" and the superscript "N" denotes "normally distributed." The quantity , g t i s , which appears in Equation (14), denotes the standard deviation associated with the total microscopic crosssection for energy group , 1,..., , g g G = and isotope , 1,..., The variance of the leakage response of the PERP benchmark is denoted as Table 2 and has the following mathematical expression: In Equation (15), the quantity [ ] ( ) ( ) ( )   3  2  2  3, ,   1 , , , var ( ) denotes the contributions stemming from the third-order sensitivities, for 30 G = energy groups and 6 I = isotopes. The results presented in the second column of Table 2 are illustrated graphically in Figure 2, which is reproduced (with permission) from [6].
, , As illustrated in Figure 2, the relationship indicates that the contribution from the second-order sensitivities exceed the contributions from the firstorder ones, and even more importantly, the contributions from the third-order sensitivities are dominant, being much larger than the contributions from either the first-order or the second-order sensitivities. Hence, neglecting the third-order sensitivities would cause a significant error in quantifying the standard deviation of the leakage response. The results in Table 2 and Figure 2 also indicate that the contributions from the second-order sensitivities cause a significant shift, by about 40%, of the expected value,  Table 2 also underscore the fact that the larger the parameters' standard deviations, the larger the impact of the second-order and third-order sensitivities. It was shown in [6][7][8][9][10][11] that correlations among the cross-sections enhance the effects of the high-order sensitivities.
The results obtained in [6][7][8][9][10][11] and summarized in this section, which indicated that the PERP leakage response was most sensitive to the hydrogen cross-sections, have motivated the high-order sensitivity and uncertainty analysis to be presented next, in Section 3.

High-Order Sensitivity and Uncertainty Analysis of Neutron Slowing Down in Hydrogen
The multigroup fluxes and cross-sections defined in Equations (6) -(10) are used not only in PARTISN [23] but are used in all of the deterministic solvers (e.g., SCALE [2]) of the Boltzmann multi-group neutron transport equation. Although the multigroup Boltzmann equation shown in Equation (4) is formally exact (i.e., no approximations have been made in deriving it), this equation cannot be solved as it stands since the energy-dependent flux ( ) , ,E ϕ r Ω , which enters into the definitions of the multigroup fluxes and crosssections defined in Equations (6) -(10), is not known at that stage. The first step in the usual procedure for determining the energy-dependent flux ( ) , ,E ϕ r Ω is to partition the energy interval of interest into very many ("fine" or "ultra-fine") energy groups and solve the "infinite-medium" energy-dependent Boltzmann quasi-analytically to obtain the en-ergy-dependent "infinite-medium" flux [18][19][20]. In view of the results presented in Section 2, it is of interest to investigate the sensitivities of the energy-dependent neutron flux with respect to the total microscopic cross-section of hydrogen in a hydrogen-moderated system, since the energy-dependent neutron flux is directly related, via Equation (12), to the leakage response.
The main features of the energy-dependent neutron flux are determined by considering a medium of large ("infinite") extent, in which neutrons emitted by a source (which may be an external source or an internal fission source) will lose energy ("slow down") following scattering collisions and will disappear following absorption collisions with the nuclei in the medium. A widely used model [17][18][19][20], which includes the scattering and absorption processes of neutrons emitted by an energy-dependent source ( ) Q E in a hydrogen-moderated system (which includes absorbing nuclide(s) having atomic mass(es) much larger than unity), has the following form: Recalling that The Taylor series in Equation (22) The results obtained in Equations (22) and (24) are essential for using the probability distribution of ( ) s E σ , which is seldom known completely in practice, to construct statistical indicators (i.e., mean value, standard deviation, etc.) for the distribution of ( ) The uniform distribution defined in Equation (26)

Sensitivity Analysis
It follows from the expression given in Equation (20) that the absolute sensitivities of ( ) It follows from Equation (31) that: In turn, it follows from Equation (32) that the smallest order, n, which is required to ensure that the Taylor-series expansion shown in Equation (22) represents the exact flux given in Equation (19) within an error ( ) n ε Φ , must satisfy the following inequality: In particular, Equations (33) and (24) As an example of the use of the Taylor-series expansion in sensitivity analysis, consider that it is desired to predict the value of ( ) is "one standard deviation away from the mean," i.e., when: Replacing the result obtained in Equation (35) into the definition provided in Equation (23) and taking 2. If only the first-order sensitivities of ( ) then the first-order expansion ( 1 n = ) in Equation (23) yields the following result: It follows from Equations (37) and (38) that the error of the first-order approximation is: Thus, the result obtained by keeping only the first-order sensitivities in the Taylor expansion underestimates the exact flux by 33.34%.
It follows from Equations (40) and (37) that the error of the second-order approximation is as follows: ; 19.24% ; Thus, the result obtained by keeping the first-and second-order sensitivities in the Taylor expansion overestimates the exact flux by 19.24%.
It follows from Equations (42) and (37) that the error of the third-order approximation is as follows: ; 11.11% ; Thus, the result obtained by keeping the first-, second-and third-order sensitivities in the Taylor expansion underestimates the exact flux by 11.11%.
It follows from Equations (44) and (37) that the error of the fourth-order approximation is as follows: ; 6.42% ; Thus, the result obtained by keeping the first-, second-, third-and fourth-order sensitivities in the Taylor expansion overestimates the exact flux by 6.
The n th -order approximation, denoted as provided by the following truncated Taylor series: Using Equations (19) and (46) in Equation (49) yields the following result: As indicated by the result obtained in Equation (50) The skewness of the distribution of the response ( ) and is defined as follows: 3   2  2  2  2   3  3  3  2  2 2   2  2   1  2  3  1  1  1  ln  ln  2  1  1  1 2 1 Skewness indicates the direction and relative magnitude of a distribution's deviation from the normal distribution. Note that, for this illustrative model, The n th -order approximation of the third-order moment, denoted as ( It follows from Equations (56), (23) and (48) that the first-order approximation of the third-order moment, denoted as ( ) [ ] 1 3 H μ Φ , has the following expression: It also follows from Equations (56), (23) and (48)  μ Φ , has the following expression: It is evident that high-order approximations for the third-order moment are very laborious to obtain. On the other hand, since the definition of skewness involves not only the response's third moment but also the response's variance, it is possible to compute approximate values of the skewness, denoted as , by using approximations of different orders for the response's variance and third moment, respectively, as follows: The n th -order approximation of the fourth-order moment, denoted as ( ) It follows from Equations (62), (23) and (48) that the first-order approximation of the fourth-order moment, denoted as ( ) [ ] 1 4 H μ Φ , has the following expression: It also follows from Equations (62), (23) and (48)  μ Φ , has the following expression: The kurtosis of a probability distribution indicates the propensity of the respective distribution to produce outliers or heavy tails.
The approximate kurtosis, denoted as , can be computed by using approximations of different orders for the response's variance and fourth moment, respectively, as follows: The flexible definition for the approximate skewness Equation (66) aims at obtaining higher accuracy for the kurtosis with less computational effort than would be required if the fourth moment were computed at the same high order = m n as the variance. Thus, the variance would be computed relatively inexpensively at an order of approximation > n m, while the fourth moment would be computed only up to the order of approximation < m n, expecting that Kurt for a sufficiently high order m . Section 3.1 has presented the investigation of the order of sensitivities that would be needed to achieve a target accuracy for predicting changes in the flux due to changes in the cross-section ( ) s E σ by using the Taylor-series expansion given in Equation (23).
Analogously, it is of interest to investigate the accuracy provided by the approximations The first-order approximate expectation has the following expression: The expectation is accurate only to the zeroth-order term, since the contribution from the first-order sensitivities vanishes after the integration of the respective first-order (odd) term over the symmetric-interval around the expectation of the parameter Replacing the expression given in Equation (67) into Equation (52) yields: is accurate up to and including the contributions stemming from the first-order sensitivities.
2. If first-order and second-order sensitivities of ( ) available, then the second-order expansion takes on the following particular form of Equation (23): The expectation has the following expression: The expectation is accurate up to and including the contributions stemming from the second-order sensitivities.
Replacing the expression given in Equation (70) into Equation (52) and using Equation (71) yields: is accurate up to and including the contributions stemming from the first-order sensitivities but the contributions from the second-order sensitivities are incomplete because the term β will arise from the third-order sensitivities, as will be shown below.
3. If all sensitivities up to and including the third-order sensitivities of ( ) takes on the following particular form of Equation (23): The expectation has, therefore, the following expression: The expectation is only as accurate as the lower-order approximation , it is accurate up to and including the contributions stemming from the second-order sensitivities, because the contribution from the third-order sensitivities vanishes after the integration of the respective third-order (odd) term over the symmetric interval around the expectation of the parameter Replacing the expression given in Equation (73) into Equation (52) and using Equation (74) yields: The variance ( ) [ ]  (23): The expectation The expectation is accurate up to and including the contributions stemming from the fourth-order sensitivities.
Replacing the expression given in Equation (77) into Equation (52) and using Equation (77) yields: The normalized exact standard deviation is defined as follows: The results obtained from using the exact expression provided in Equations (79)      ; H E σ Φ also increases. With pronounced skewness, standard statistical inference procedures such as a confidence interval for a mean will be not only incorrect, in the sense that the true coverage level will differ from the nominal (e.g., 95%) level, but they will also result in unequal error probabilities on each side. Very importantly, the results presented in bold in Table 6 indicate that the first-order approximation for the third moment ( ) [ ] 1 3 H μ Φ of the predicted response distribution is always zero. Thus, using just the first-order response sensitivities will always be 100% incorrect unless the distribution of the response under consideration happens to be symmetric. Conversely, if only first-order sensitivities are considered, the third-order moment of the response is always zero. Thus, the response distribution is always predicted to be symmetrical if only first-order sensitivities are considered. Hence, a "first-order sensitivity and uncertainty quantification" will always produce an erroneous third moment (and hence skewness) of the predicted response distribution, unless the distribution happens to be symmetrical. At least second-order sensitivities must be used in order to estimate the third-order moment (and, hence, the skewness) of the response distribution. The overall conclusion that follows from the results presented in Table 6 is that loworder (i.e., first-and second-order) approximations for the third-order moment, of the response distribution can provide the correct sign (i.e., positive or negative) of the response skewness, but cannot provide reliably correct values for the skewness. Significantly higher-order approximations, which would require the computation of high-order response sensitivities to model parameters, are needed in order to obtain reliable values for the third-order response moment

[ ]
3 H μ Φ and, hence, for the skewness of the response distribution. Comparing the results in Table 6 for increasing (positive) values for the parameter β , the low-order (first-and second-order) approximations become increasingly unreliable as β increases toward unity. This trend indicates that, as the response distribution becomes wider, sensitivities of an increasingly higher-order would be needed in order to improve the predicted values for the response skewness. Using a loworder (in Table 6: second-order) approximation for the third moment higher-order (in Table 6: fourth-order) approximation for the variance (so that a comparable amount of computational effort would be needed for computing the respective approximations) does not significantly improve the final results. If accurate values for the response skewness are required, it is not possible to circumvent the need for computing high-order sensitivities, which in turn enable the computation of high-order approximations to the third-order moment and, hence, skewness. The accuracies of the first-order and second-order approximations for the fourth-order moment are illustrated by the results presented in Table 7.
The results presented in Table 7 also indicate (similar to that indicated by the results presented in Table 6) that for small values of the parameter β , the first-order approximation,  Table 7 are similar to those that were drawn from the results presented in Table 6 for the response's third-order moment and skewness, namely: low-order (i.e., first-and second-order) approximations for the fourth-order moment, [ ] and, hence, for the kurtosis of the response distribution. The results in Table 7 indicate that the low-order (first-and second-order) approximations become increasingly unreliable as β increases (with positive values) toward unity, which means that, as the response distribution becomes wider, sensitivities of an increasingly higherorder would be needed in order to improve the predicted values for the response kurtosis.
If accurate values for the response kurtosis are required, it is not possible to circumvent the need for computing high-order sensitivities. It is known [18][19][20] that the form showed in Equation (18) of the slowing down neutron flux in a hydrogenous medium also holds for the asymptotic (below the cut-off energy of the fission source) slowing down flux in an infinite non-hydrogenous medium, where it takes on the following form: In Equation (82), the quantity ( ) where It follows from Equations (83) and (84) which has the same functional form as the series for in Equation (20). Therefore, the analysis that has been performed in the foregoing for is also applicable to perform sensitivity and uncertainty analyses aimed at quantifying the impact of the uncertainties afflicting the other parameters (cross-sections, number densities, source parameters, etc.) on the distribution of the slowing down flux ( ) ; E Φ α in a scattering and absorbing medium, but these additional analyses are beyond the scope of this work.

Use of Model Response Sensitivities in First-Order Data Adjustment and High-Order Predictive Modeling Methodologies
This section reviews the incorporation of sensitivities (of various orders) in data adjustment and predictive modeling methodologies in order to indicate the level of accuracy that could be expected when using these methodologies. Section 4.1 reviews the first-order "generalized linear least squares adjustment (GSSLA)" methodology [2] incorporated in the data adjustment software system TSURFER [1], developed at Oak Ridge National Laboratory (ORNL), which represents the current state of the art in data adjustment used in nuclear engineering for nuclear criticality safety and reactor physics/shielding applications. Section 4.2 compares the TSURFER-GSSLA data adjustment methodology to the high-order predictive modeling (HO-BERRU-PM) methodology recently developed by Cacuci [13,21], highlighting the current capabilities and limitations of these two methodologies.

First-Order Generalized Least Squares Data Adjustment (GLLSA) Methodology in TSURFER
The current state of the art in "data adjustment" is represented by the so-called "generalized linear least squares adjustment (GLLSA)" procedure [2] implemented in the soft-ware module TSURFER [1] of the ORNL-SCALE software system used for reactor criticality safety applications. The GSSLA procedure in TSURFER is based on concepts proposed in the 1970s [25][26][27] and is limited to incorporating just the first-order sensitivities of the response with respect to the model parameters. The end results for the "adjusted" values produced by the TSURFER-GLLSA procedure for the nuclear cross-sections (considered as "model parameters") and effective multiplication factors (considered as "model responses") are obtained by minimizing a (user-chosen) quadratic cost functional, which represents the discrepancies ("least-squares") between measurements and computations. The GSSLA procedure in TSURFER considers the following a priori information, which is described in the paragraphs labeled (A) through (C), below: The relation in Equation (97) indicates that the expected value of the calculated response in TSURFER is the same as the TSURFER-calculated response itself. As has been shown in Section 3 and will also be shown in Subsection 4.2, however, the assumption made in TSURFER that Equation (97) is correct holds only if the second-and higher-order terms in the Taylor series of ( ) i k α are discarded. Of course, Equation (97) is invalid when second (and/or higher) order sensitivities are not discarded.
The linear first-order expression provided in Equation (96) is also used in TSURFER to obtain expressions for the variances/covariances, where the symbol [ ] ( ) The matrix kk C  is called in TSURFER the "I by I absolute covariance matrix for prior calculated responses," and Equation (102) is called the "sandwich rule." The GLLSA procedure [2] implemented in TSURFER "consolidates" the a priori information provided above in items (A) through (C) in order to obtain "adjusted" (i.e., "improved estimates of") mean values and covariances for the parameters (cross-sections), as well as "adjusted" values for the means and covariances of the responses and the (target) application under consideration. The GLLSA procedure computes the adjusted values for the means and covariances of the parameters and responses by performing a constrained minimization of the user-chosen "chi-square functional" defined below Furthermore, the quantities ′ α and ′ m in Equation (103) represent the "adjusted values" for the model parameters and responses, which are to be determined as the end results of performing the constrained minimization procedure. This minimization procedure is performed by constructing an augmented Lagrangian functional, which is obtained by using Lagrange multipliers to append the "hard constraints" represented by Equation (96) to the user-defined "chi-square functional" shown in Equation (103) and, subsequently, determining the values/expressions for ′ α and ′ m at which the first-order derivatives of the resulting augmented Lagrangian vanish.
The TSURFER formulas for the adjusted quantities that result from the GLLSA procedure are reproduced below, written in what TSURFER calls "absolute format": 1. The formula for computing the M -dimensional vector of adjusted cross-sections ("model parameters"), denoted in TSURFER as ′ α , is as follows: where: † + = + dd kk mm kα αα kα mm Additionally, 2. The formula for computing the M M × -dimensional covariance matrix for the adjusted cross-sections, denoted in TSURFER as ′ ′ α α C , is as follows: 3. The M -dimensional vector of adjusted measured responses produced by TSURFER's GLLSA procedure is denoted as ′ m , and the formula for computing it is as follows: 5. The GLLSA procedure in TSURFER forces the following relationships to hold: where ′ m denotes the I -dimensional vector of adjusted measured responses produced by the GLLSA procedure and ( ) ′ ′ k α denotes the I -dimensional vector of adjusted calculated responses obtained using the adjusted nuclear parameters ′ α . Consequently, the following relationships are also forced to hold in TSURFER: 6. TSURFER also computes the "consistency indicator between the calculations and measurements," which is given in TSURFER by the following "chi-square formula": 2 ? .
The application response bias in TSURFER is denoted as α β and defined as the expected deviation of the original calculated application response, denoted as ( ) a k α , from the best estimate of the measured response, denoted as a m′ , which is unknown but is assumed to obey some probability distribution. If the application response did actually have a prior measured value a m , then the best estimate for the experiment value would be the adjusted value a m′ provided by the GLLSA procedure. Mathematically, TSURFER defines the application bias as the expectation of the difference between ( ) a k α and a m′ , i.e., The first-order linear relation shown in Equation (113) is used in TSURFER as the basis for constructing "confidence intervals" and "tolerance limits" for the predicted multiplication factors (predicted model responses) for which there are no measurements. TSURFER assumes a first-order approximation of the dependence of the effective multiplication factors on the underlying cross-sections and, therefore, assumes a linear relationship between the initial experimental biases Of all of the above assumptions made in TSURFER for constructing "confidence intervals" and "tolerance limits" for the predicted multiplication factors, the most questionable (and severe) is the assumption of linearity shown in Equation (113), which holds only if all sensitivities higher than the first order are discarded. However, TSURFER cannot compute second-(or higher) order sensitivities, so the validity of Equation (113) is unproven. Evidently, if second-and/or higher-order sensitivities were not negligible, Equation (113) would be severely in error, and consequently, the "confidence intervals" and "tolerance limits" computed in TSURFER would be severely misleading. In conclusion, the "confidence intervals" and "tolerance limits" computed in TSURFER are not trustworthy until TSURFER implements capabilities for computing higher-order sensitivities and proves quantitatively that these are negligibly small. Otherwise, all of the uncertainty analysis and data adjustment results, including tolerance limits and confidence intervals, produced by TSURFER must be interpreted with due caution.

High-Order Sensitivities in the HO-BERRU-PM Methodology
The HO-BERRU-PM methodology (high-order best-estimated results with reduced uncertainties predictive modeling) conceived by Cacuci [13,21] combines experimental and computational information in the joint phase space of responses and model parameters, including second-order and third-order response sensitivities to model parameters. The HO-BERRU-PM methodology uses the maximum entropy principle [22] to eliminate the need for introducing and "minimizing" a user-chosen "cost functional quantifying the discrepancies between measurements and computations," thus yielding results that are free of subjective user-interferences, while generalizing and significantly extending the current dynamic data assimilation procedures [28,29]. Incorporating correlations, including those between the imprecisely known model parameters and computed model responses, the HO-BERRU-PM also provides a quantitative metric, constructed from sensitivity and covariance matrices, for determining the degree of agreement among the various computational and experimental data while eliminating discrepant information.
In contradistinction to TSURFER, which uses the linear (first-order) model shown in Equation (96), Cacuci's HO-BERRU-PM framework [13,21] uses the third-order multivariate Taylor In Equations (114) Terms involving response sensitivities higher than the fourth order in the parameters' standard deviations have been discarded in Equation (122). The expression in Equation (122) involves implicitly the fourth-order moment, 4 ijkl μ , of the multivariate parameter distribution function ( ) p α α and, explicitly, the fourth-order parameter correlation, ijkl q ; these fourth-order moments are defined as follows: In particular, the variance of a response ( ) 1 c i r α is obtained by setting 1 2 i i = in Equation (122). When the parameters follow a multivariate Gaussian distribution (as it is often assumed in practical applications), the following relations hold: 0; ; , , , 1,..., .
The third-order cumulant for three responses, The covariance of a response, The covariances between the computed responses and the model parameters defined in Equation (126) are considered to be the components of an ( ) r N N α × -dimensional matrix denoted as rα C and defined as follows: ( ) ( ) 1. The best-estimate posterior expectation values for the vector of predicted model parameters; this vector is denoted as be α and is given by the following formula: where the vector d denotes the vector of deviations between the expected value of the computed response and the expected value of the measured response, namely: 2. The best-estimate posterior parameter covariance matrix, denoted as be α C , for the bestestimate parameters be α : In Equation (138) 3. The best-estimate posterior expectation values for the vector of predicted responses, which is denoted as be r and which has the following expression: 4. The best-estimate posterior covariance matrix for the best-estimate responses be r , which is denoted as be r C and which has the following expression: C , which means that the variances contained on the diagonal of the best-estimate matrix be r C will be smaller than the experimentally measured variances contained in m C . Hence, the addition of new experimental information will reduce the predicted best-estimate response variances in be r C by comparison to the measured variances contained a priori in m C .
5. The posterior covariance matrix comprising the best-estimate correlations between the best-estimate parameters be α and the best-estimate responses be r , which is denoted as be r α C , has the following expression: 6. The HO-BERRU-PM methodology also yields a "chi-square" indicator, which measures the agreement/disagreement between the experimental responses and the computed responses. This indicator is denoted as 2 HO χ (since it includes higher-order sensitivities and correlations) and has the following expression: As the expression obtained in Equation (142)

Comparison: HO-BERRU-PM Methodology versus TSURFER-GLLSA Methodology
The correspondences between the notations used in TSURFER and HO-BERRU-PM, respectively, are as follows:  At the most fundamental conceptual level, the TSURFER-GLLSA methodology relies fundamentally on minimizing a user-chosen "chi-square" functional subject to the requirement that the following hard constraint be satisfied: "the first-order Taylor-expansion of the calculated response as a linear function of the model parameters." In contradistinction, the HO-BERRU-PM methodology eliminates the need for a "user-chosen functional to be minimized" by employing the maximum entropy principle to combine, without imposing any "hard constraints," the available information regarding the measured and computed responses, including a nonlinear third-order Taylor-series representation of the underlying model, connecting the model parameters to the model's calculated responses. The distinctions between the best-estimate parameter values, be α , predicted by the HO-BERRU-PM methodology, and the adjusted parameter values, ′ α , produced by TSURFER-GLLSA, are obtained by replacing Equations (129) The difference between the best-estimate covariances be α C for the best-estimate parameter values, be α , predicted by the HO-BERRU-PM methodology and the adjusted covariances The difference between the best-estimate response values, be r , predicted by the HO-BERRU-PM methodology and the adjusted response values, ′ m , produced by TSURFER, is obtained by replacing Equations (129) The difference between the high-order indicator, 2 HO χ , produced by the HO-BERRU-PM methodology and the indicator 2 TSURFER χ produced by TSURFER is obtained by replacing Equations (129) and (143) in Equation (142) and subtracting Equation (112) from the resulting equation. These operations yield the following relation: The distinctions between the posterior adjusted quantities produced by TSURFER-GLLSA versus the HO-BERRU-PM framework are summarized in Table 9. It is evident from Equations (144) − (148) that the expressions of the quantities predicted by the HO-BERRU-PM methodology comprise, as particular cases, the quantities adjusted by TSURFER-GLLSA. The quantities predicted by the HO-BERRU-PM methodology will become identical to the quantities adjusted by TSURFER-GLLSA if (and only if) all of the sensitivities higher than the first order are ignored. Notably, the computation within the HO-BERRU-PM methodology of the best-estimate parameter and response values, together with their corresponding best-estimate covariance matrices, only requires the computation of ( ) as the corresponding matrix, dd C  , which needs to be inverted in TSURFER, even though the matrices in the HO-BERRU-PM methodology comprise not only the first-order response sensitivities (as in TSURFER-GLLSA) but also comprise all of the second-and third-order response sensitivities to the model parameters. This important HO-BERRU-PM characteristic is essential for applications to large-scale practical problems, since in the overwhelming majority of practical situations, the number of model parameters far exceeds the number of responses of interest.

Discussion and Conclusions
This work has underscored the need for computing higher-order (i.e., higher than first-order) sensitivities (functional derivatives) of model responses with respect to the model parameters. The significant potential impact of the higher-order sensitivities on the expected values and variances/covariances for the calculated and predicted model responses have also been highlighted. In particular, if only first-order sensitivities are considered, the third-order moment of the response is always zero. Hence, a "first-order sensitivity and uncertainty quantification" will always produce an erroneous third moment (and, hence, skewness) of the predicted response distribution, unless the unknown response distribution happens to be symmetrical. At least second-order sensitivities must be used in order to estimate the third-order moment (and, hence, the skewness) of the response distribution. Skewness indicates the direction and relative magnitude of a distribution's deviation from the normal distribution, while kurtosis indicates the propensity of the predicted response distribution to have heavy tails and/or outliers. With pronounced skewness, standard statistical inference procedures such as constructing a confidence interval for the mean (expectation) of a computed/predicted model response will be not only incorrect, in the sense that the true coverage level will differ from the nominal (e.g., 95%) level, but the error probabilities will be unequal on each side of the predicted mean.
The evident conclusion that can be drawn from this work is that the consideration of only the first-order sensitivities is insufficient for making credible predictions regarding the expected values and uncertainties (variances, covariances, skewness) of calculated and predicted/adjusted responses. At the very least, the second-order sensitivities must also be computed in order to enable the quantitative assessment of their impact on the predicted quantities. Since the second-order sensitivities impact decisively the expected values and the skewness of the calculated/predicted responses, they will also impact the computation of confidence intervals and tolerance limits for the predicted expectation of these responses. For large-scale nonlinear systems, the second-order response sensitivities can be computed exactly and efficiently by applying the adjoint methodology developed by Cacuci [5]. Extending the methodology presented in [5] to enable the exact and efficient computation of third-order (and higher-order) sensitivities for large-scale nonlinear systems is currently in progress. For large-scale (response-coupled) linear forward/adjoint systems, all of the response sensitivities up to and including the fourth-order sensitivities can be computed exactly and most efficiently by applying the methodology (fourth CASAM), recently developed by Cacuci [30].