Propagation of Model Uncertainty in the Stochastic Simulations of a Compartment Fire

Model validation and probabilistic simulations are routinely used for quantifying the uncertainties originating from the numerical models and their inputs, respectively. How the two uncertainty types combine in the context of fire risk analyses is not well understood. In this work, we study the propagation of modeling uncertainty to the predicted distributions of probabilistic fire simulations using model validation data representing an uncertain compartment fire scenario. The wall temperatures are predicted in three different ways: one using a coupled model in which the input is the fire heat release rate, and two models using a standalone conduction solver and either experimentally or numerically (CFD) determined heat flux as a boundary condition. Using the predicted wall temperatures, we calculated demonstrative wall failure probabilities assuming different critical threshold temperatures. We propose a simple method for correcting the simulated distributions and probabilities towards the experimentally observed ones. The simulation results with the Fire Dynamics Simulator show that the obtained uncertainties of this particular validation set are similar to the ones reported in the validation guide. In average, the most accurate model over-predicts wall temperature by 5.0% and the prediction uncertainty for both gas phase and solid phase temperature is 10%. The wall temperatures predicted from the measured heat-fluxes show higher modeling uncertainty than the ones predicted by a coupled model of the entire gas-wall system. The proposed correction method is shown to improve the accuracy of the predicted distributions for internal wall temperatures at different times. In practical applications, this would lead to more accurate estimates of the time-dependent failure probabilities.


Introduction
With an aim to enable the performance-based design of fire safety, the development and validation of essential computational tools is underway. In about a halfcentury, in this context, numerous simulation tools emerged. Some of them primarily evolved to solve the fire-related problems, while others are general purpose tools that now include the fire simulation module. Among them, the tools adopting the well-known integration techniques such as the finite difference method that the model uncertainty metrics can be used to statistically compensate for their effect in a probability calculation. Unlike in the previous works, where the uncertainties are usually presented for the peak values, we present the uncertainty pertaining to each time instance of the output.

Parameter Uncertainty
If the inputs of a mathematical model are uncertain then the outputs will be uncertain too. This uncertainty propagation depends upon the characteristics of the model itself. The expression of uncertainty in output, T ¼ f ðXÞ, f being continuous and one time differentiable function, can be derived by Taylor expanding T about its mean and utilizing the definition of standard deviation in T [16]. The first order approximation is, where r 2 T represents variance in T, R X is variance-covariance matrix of the input vector, X, and J ¼ ðJ 1 ; J 2 ; J 3 . . .Þ, J i ¼ @f =@X i . If the input variables, X, are independent of each other then the Eq. 1 would simply reduce to Figure 1 depicts the uncertainty propagation for a simple model, T ¼ pX . A normally distributed output, T $ N p10; p 2 À Á , is obtained for a normally distributed input, X $ N 10; 1 ð Þ. For complex and non-linear problems, such derivation is mathematically challenging, therefore, stochastic methods are adopted. Some examples of stochastic methods are Monte-Carlo (MC), Latin hypercube sampling (LHS) and Fourier amplitude sensitivity test (FAST) [17][18][19].

Combining Model and Parameter Uncertainty
The model uncertainty can be decomposed into two components: systematic bias and random error [7]. The systematic bias is assumed to be a measure of the multiplicative factor by which the observed output is away from the true value. On average, it is the ratio of observed and true output. The random error is assumed to be an additive error that makes the observed output to fluctuate around the true value. We assume that these parameters can be determined for each output parameter, and are constants for a specific type of fire scenario.
The output for a simulation model, T ¼ f ðX Þ, with systematic bias, d, and random error, , iŝ whereT is the simulated quantity and T is the true quantity. Here, the T and are independent and the mean of is zero. For such conditions, the mean and variance of the observed quantity can be written as, Where l T and r 2 T are the mean and variance of the true quantity and r 2 is the variance of the random error. The derivation for these expressions can be found in the Appendix A.
For a normally distributed output, T, Table 1 lists the expressions of distributions in the presence or absence of model uncertainty. Figure 2 shows the histogram plots for specific values of d and r . The left figure compares the effect of only the bias, the middle one compares the effect of only the random error, and the right one compares the effect of both. Figures show that the bias simply shifts the distribution, while the random error widens it. Table 1 The Output Distribution in Presence or Absence of Error S.N.
Distribution ofT Condition Description

Correction of Output Distribution
If the prior information of d and r is available, one can correct the simulated output towards the true one. The corrected moments, l T and r T , can be derived from Eq. 4. The corrected distribution is then the distribution generated using the corrected moments. The cumulative distribution function (CDF) of the corrected values can be obtained by substituting l T and r T into the general expression of CDF. For Gaussian distribution, the CDF is This method works well for the output distributions that can be represented by the first two moments, e.g., Gaussian and uniform. When the shape of the output distribution cannot be well represented by the first two moments, the output distribution can be corrected using where T is the corrected realization corresponding to the observed realization,T . The derivation for this expression can be found in Appendix A. We illustrate the correction method using two arbitrarily chosen examples [20]. In one of the examples, both the simulated and the true distribution are Gaussian, while in the remaining one, the distribution shape is irregular. First, we calculate the correction parameters, d and r , by comparing the simulated and true values, whereT i and T i are the ith realization of the simulated and the true quantity respectively and N is the sample size. Then, using the correction parameters we estimate the true shape from the simulated one. Figure 3 shows the true, simulated and corrected distributions along with CDF. In the upper plots, the dotted line represents the distribution generated using Eq. 5 and the continuous line represents the distribution generated using Eq. 6. Plots indicate that both methods work well with the normally distributed output. For irregularly distributed outputs, as expected, the estimation is better with Eq. 6. The maximum difference between the CDF of true and the corrected distribution is $ 0.05 and $ 0.01 respectively for Eqs. 5 and 6. Even with Eq. 6, complete trace-backing is not possible because the random error that occurred per realization cannot be known.

Sampling Uncertainty
In the stochastic analysis, the inferred moments and the probabilities depend upon the sample size and sampling method. This is known as sampling uncertainty.  fractiles values, z 95 , are presented for sample sizes N=100, 1000 and 10,000. Higher sample size well represents the distribution and z 95 values increases with the increase in the sample size.
The sampling uncertainty can be presented as AE bounds from the corrected value. For example, if the probability inferred from the corrected distribution is p, then the probability is p AE Dp, where Dp is the sampling uncertainty. The sampling uncertainty for simple MC simulation having sample size N is z a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pð1 À pÞ=N p , where z a is a multiplier number that determines the level of confidence [21]. For 99% level of confidence z a is 2.58. For LHS, such analytical expression is not available, and a separate convergence analysis is needed. Figure 5 shows the result of the convergence analysis carried out for the distributions presented in Fig. 4. The left plot shows z 95 ðN Þ. The right plot shows their difference with the converged value, z 95 (N ¼ 10; 000), and the maximum bound represents the sampling uncertainty. With N = 1000, the corrected z 95 and the sampling uncertainty are 61 and 2 respectively. This means the 95% fractiles value is 61 AE 2.

Validation Experiment
In October of 1998, a series of fire test was carried out at VTT Building Technology with an aim to produce a set of data for validation of fire models [22]. The tests were conducted in a compartment, 10 Â 7 Â 5 m 3 , having one door opening to the large fire testing hall. The walls and ceiling were made of lightweight concrete and the floor was made of normal concrete. Figure 6 depicts one of the test setups with a fire plume and measurement devices. Table 2 lists the material properties and the thickness of the obstructions.
Systematic variations of fire size and locations were made to determine their effect on the fire environment. The selected fire locations are indicated in Fig. 6 and the test series are summarized in Table 3. Test 10 and 14 were for calibrations, and hence not included in the table.
The fire source was n-Heptane circular steel pool placed over a load cell measuring mass loss rate. Water was used under n-Heptane to stabilize the fire. The free height from the water surface to pool edge was 0.13 m. The free height from the fuel surface to pool edge was 0.11 m in most of the tests.
Burning rates, gas temperatures, wall temperatures, and heat fluxes were measured during the tests. There were 30 thermocouples to measure hot gas  layer(HGL) temperature, 46 thermocouples to measure ceiling jet temperature, 25 thermocouples to measure plume temperature, 5 heat-flux gages to measure the heat flux on the wall and 9 thermocouples to measure the inside wall temperature. Figure 7 shows the three locations where the inside wall temperatures were measured. At each location, a light-weight concrete block with three thermocouples was placed in order to measure the temperatures at varying depths from the inner wall surface.

FDS Model
The fire experiment was modeled using FDS version 6.5.3. Figure 8 shows the 3D representation of the simulation domain with transparent gas region and the gray structural region. The fire driven flows in the gas region were simulated by numerically solving the weakly compressible form of the Navier-Stokes equations. The governing equations are presented in Technical Reference Guide of FDS [23]. The heat-transfer in the structural region was simulated by numerically solving the one-dimensional heat-conduction equation with boundary conditions where x is the wall/ceiling depth from the heat-exposed surface. l represents the wall/ceiling thickness such that the hot-side and cold side surface are at x ¼ 0 and x ¼ l respectively. q 00 is the interface heat-flux, e is the emissivity, T 1 is the ambient temperature and the field variable T(x, t) represents the wall temperature. The density, q, the specific heat capacity, c p , and the thermal conductivity, k, are assumed to be constant. The heat transfer coefficient, h, at the front and the back of the wall is calculated based on a combination of natural and forced convection correlations [23].

Fire Technology 2019
The fire source was modeled as a circular burner with an appropriate heat release rate (HRR), corresponding to a fuel inflow boundary condition. Figure 9 shows the specified HRR for the tests. The text in the figure indicates four different test groups having the same pool diameter. The pool front surface temperature was specified to follow the same trend as HRR, starting from the room temperature at t ¼ 0 and increasing to a peak value of 98.4 C at the time of the peak HRR, and staying in that value until the end of the simulation. To account for the incomplete combustion, soot yield was set to 2%. Figure 10 depicts the discretization of the gas domain. To accurately resolve the cylindrical shape of the heat source, the computational cells near the heat source is refined to 5 cm. The cell size in the rest of the region is 10 cm. Similarly, the right side in Fig. 11 depicts the discretization of the 1D heat-conduction model. To resolve the spacing (0.5 mm) between the thermocouples measuring the wall  temperatures, the compartment wall (0.3 m thick ) is discretized into 1000 grid points with finer spacing close to the surfaces.
The simulation was carried out using a distributed-memory computer. We used two steps to decompose the simulation domain. First, the entire compartment is divided into 6 Â 5 ¼ 30 mesh regions. The right side in Fig. 10 shows the bottom six of them, the view from the top. Then, for the bottom layer, one of the mesh regions is further divided into nine regions. This results all together 29 þ 9 ¼ 38 individual mesh regions. It applies to all tests except the ones in which the pool is located in the middle of the compartment. For these cases, the gas domain is divided into 9 Â 5 ¼ 45 mesh regions with one of them further divided into nine, resulting all together 44 þ 9 ¼ 53 individual mesh regions. In both cases, there are two additional mesh region covering the outside geometry, see Fig. 8. Each mesh region used one core and 500 MB memory from a CPU. The models representing all 17 tests were computed at the same time and the total computation time was $ 6 h.

Standalone Analysis
Excluding the gas phase computation, we carry out a separate analysis to predict the compartment wall thermal response. In this case, the heat-diffusion in the wall is simulated in response to a pre-defined boundary heat-flux that represent the possible fire scenario. Figure 11 shows a schematic diagram of the standalone model. The left figure shows the Gauge heat-flux, q 00 , measured during the experiment and a ¼ k=ðqc p Þ indicates the wall material property. i ¼ 1; 2; . . . N and j ¼ 1; 2; . . . represent the spatial and temporal discretization with N nodes and t time steps respectively. The cold-side boundary condition is both convective and radiative heat flux with a heat transfer coefficient h and emissivity e respectively. The wall temperatures, T i;j , were predicted in response to the heat-flux obtained in two different ways; (i) measured during the experiment, q 00 Exp , (ii) predicted from the (CFD+FDM) coupled analysis, q 00 FDS . We compare the solutions of the coupled model and the standalone model to find out how the different error types propagate to the wall temperature predictions. If, for instance, q 00 Exp was error-free and there was no error in interpretation of q 00 Exp as a boundary condition, then the standalone model with q 00 Exp boundary conditions should be more accurate than the coupled model. This is due to the fact that, in the coupled model, both input uncertainties (most importantly fuel mass loss rate) and the gas phase model uncertainty propagate to the wall temperature prediction. In general, the measurement uncertainty of q 00 Exp is higher than the measurement uncertainty of the fuel mass loss rate [22], and the relative performance of the different modeling techniques is not obvious. In addition, comparing the coupled model uncertainties against the standalone model with q 00 FDS boundary condition will indicate how much error is generated by the process of interpreting specified (measured or predicted) heat fluxes as a boundary condition for the numerical model. For this test, a noticeable discrepancy can be seen in the beginning and around 8 min. The experimental uncertainty in this test may also be higher than average as the temperatures were significantly higher, and as there was only one repetition of this particular scenario. Figure 13 show the measured and predicted inside wall temperatures for the three versions of the heat flux boundary condition. Figure 14 show the times at which the inside wall temperature exceeds a given threshold, T cr C. The maximum value of the temperature measured during the experiment was 500 C, and the horizontal axis of the lower plots is normalized by this maximum value, i.e., T cr=500. Most of the curves showing the threshold times do not reach the end because the peak temperature value for those test is below 500 C. The results indi- cate that the coupled model predictions overlap more to the measured values than the predictions from the standalone model with q 00 Exp boundary condition. Although the discrepancy in temperatures is small, the discrepancy in times to reach a threshold temperature are sometimes very high, especially when the threshold is close to a semi-steady temperature of the particular experiment. Table 4 lists the bias, d, and the second central moments of random errors, r , calculated following the methods explained in [7]. The d represents the average deviation of model prediction from the measured value. The random errors are presented as a relative term, i.e., e r ¼ r =lT . e r fEg represent the random experimental errors and e r fMg represent the random model errors. For the calculation, we used all the measurement points mentioned in Sect. 3.1. Appendix B shows the scatter plots. The uncertainty values obtained from the current experiment are close to the ones reported in the FDS Validation Guide, except for the Gauge Heat Flux output quantity. The model uncertainty and the systematic bias for the Gauge Heat flux are higher than those reported in the Validation Guide.

Modeling Uncertainty
In addition to the discretization scheme explained in Sect. 3.2, we studied how the mesh configurations affected the uncertainties ( Table 5). The corresponding

Fire Technology 2019
scatters plots are shown in Appendix C. The first configuration is the one that we explained in Sect. 3.2. In the second configuration, the 10 cm mesh region was made coarse, down to 20 cm. In the remaining two configurations, uniform cell size was used in the entire gas domain and the cylindrically shaped burner was simplified to a rectangle shape. For cylindrical shaped burner and multi-mesh configuration, the bias remains unchanged. For the rectangular burner, the bias increases despite the mesh refinement. This is due to the imperfect modeling of the burner vent area. The vent area is poorly represented when the burner surface is not perfectly aligned with the mesh face. For the coarse mesh, the effective vent area can be lower than the specified value. This results in a lower HRR and ultimately the lower bias. Figure 15 visualizes the model uncertainty as a function of time. d and e r were calculated, using Eq. 7 at each time, by comparing the measured and predicted temperatures presented in Fig. 12. The plots show that on average, Gauge heatfluxes are underestimated, while Ceiling Jet and Plume temperatures are overestimated. The Gauge heat-fluxes have higher random components than the Ceiling jet temperatures and Plume temperatures. Most importantly, we see that the   Figures 16 and 17 show the d and e r calculated based on the wall temperatures and threshold times presented in Figs. 13 and 14 respectively. The plot indicates that the wall temperatures on average are overestimated. Due to this, the predicted times are underestimated. In the models based on the FDS gas phase, the early transient temperatures are underestimated, and threshold times hence overestimated. The effect is much less in the right-most plot (standalone model with q 00 Exp boundary condition), indicating that either the HRR boundary condition or CFD solution introduces a temporal delay in the early phase. Comparison of the coupled and standalone analysis predictions over the entire time period, however, indicates that the modeling uncertainty is higher for the latter one. Furthermore, the modeling uncertainty is higher for boundary flux q 00 Exp . We therefore conclude that, for the wall temperature prediction, the propagation of gas-phase modeling uncertainty is less harmful than the propagation of heat-flux measurement uncertainty. Better predictions can be achieved with the coupled analysis.

Measured and Predicted Moments
The model uncertainty metrics presented in Figs. 15 and 16 are relative quantities and do not visualize well the quality of parameter uncertainties. Figure 18 com-

Fire Technology 2019
pares the measured and the predicted outputs in terms of their first two moments. The line or dot represents the first moment and the half-length of the error bar represent the second moment. Plots show that the measured and the predicted values are close to each other. In average, the first moments of the ceiling jet temperatures and plume temperatures are slightly over-predicted and the first moments of Gauge Heat-fluxes are slightly under-predicted. In most of the plots, the error bar after 8 minutes extends towards the negative axis. This is because after 8 minutes, the mean is close to zero and the standard deviation is high, see Figs. 12 and 13. Figures 19 and 20 respectively show the first two moments of wall temperatures and the times at which the wall crosses the threshold temperature, T cr C. In average, the wall temperatures are slightly overpredicted and because of this, the times are underpredicted. Overall, the simulated first two moments are close to the observed one.

Temperature and Probability Correction
Assuming that the wall fails when it crosses a given temperature threshold, the failure probability would be the fraction of the number of the test cases in which the wall temperature rises above this threshold. We now try to understand how the modeling uncertainty in temperatures propagates to such a probability and  how it can be corrected. As the temperature predictions of the previous section were indistinctly close to the measurements, it became very difficult to demonstrate the corrections of probabilities. We, therefore, used the results corresponding to the mesh configuration with the highest modeling uncertainty.
In Figure 21, the upper plots show the predicted, measured and the corrected wall temperatures at three different times for each of the tests. The corrected temperature were obtained using Eq. 6 and single values of model unceratinty parameters (d ¼ 1:15, e r fMg ¼ 0:16, see Table 5). We see that where the prediction and measurement are apart, the corrected value is usually closer to the measurement.
The lower plots of Figure 21 show the failure probabilities at different times for three different threshold temperatures. Initially, the walls are at ambient temperature and probabilities are zero. The probabilities increase as the number of tests exceeding the given threshold increases. The upper middle plot shows the temperatures at 4 min. At this time, the number of temperatures above 100 C is 13 for the predicted, corrected as well as the measured quantities. Therefore the probability is, 13/17 $ 0.8 (lower left plot). The number of points crossing 300 C, however, is 7 for the predicted, 3 for the corrected and 2 for the measured quantities, The predicted probabilities are higher than the measured ones. The corrected probability values are closer to the measured ones. Even though the real modeling uncertainty varies with respect to time, the probability correction carried out using the generalized (constant) value was effective at each time instance. This indicates that the model uncertainty values can be generalized for the failure probability correction. Here we used the uncertainty parameters obtained from the same campaign that we used for testing the method. In the validation guide, however, the uncertainty parameters are calculated from the result of numerous fire experiments, hence representing more generalized values.

Stochastic Analysis
In this study, the variation of fire size, the pool area, pool location and the width of the opening door represent the input parameter uncertainty. For stochastic inputs listed in Table 6, we carry out MC simulation using the model having uniform cell of size of 10 cm. The sampling size, N, is 100 and the sampling method is LHS. The selected fire type is t-square fire. For such fire, HRR is calculated using fire growth time, t g , and peak HRR as where t is time in second.   Table 5). The plot shows that the correction assuming Gaussian shape is not appropriate for wall temperatures, i.e., the distribution reaches negative axis. The correction using Eq. 6, however, narrows the width of the distribution without deviating significantly from the observed shape. Figure 23 shows the contour plot for the CDF, U, of wall temperatures. The vertical axis shows the temperature range, the horizontal axis shows the time and the embedded text show the U values. The left plot shows the predicted values.  The middle and right plots show the corrected values calculated according to Eqs. 5 and 6 respectively. Assuming that the wall fails when it crosses a given temperature threshold, the failure probability would be the fraction of the number of the test cases in which the wall temperature rises above this threshold. From the above CDF plots one can infer the failure probability. For example, the predicted probability that the wall temperature rises above 100 C before 6 min is 1-0.1$ 0.9, where as the corrected probability is 1-0.2 $ 0.8. Similarly the predicted probability for wall to rise above 200 C before 6 min is 1-0.6 $ 0.4 and the corrected probability is 1-0.7 $ 0.3. The predicted probabilities are higher than the measured ones. This is due to bias in the temperature prediction.

Discussion
The study propose correction, Eqs. 5 and 6, for the stochastically simulated output,T , based on the requirement that the corrected quantity, T, and the random error, , are independent of each other and the mean of is zero. In general, the output and the total error are dependent and the mean of the total error may not be zero. The significance of the uncertainty model presented in this study is that the total error is decomposed into a dependent constant, i.e., the ratio of simulated and corrected mean, d ¼ lT =l T , and a random component, ¼T À d Á T , which implies that the mean of must be zero.
In this study, the correction method is illustrated using the true and observed data that are perfectly aligned in time. Figure 24 demonstrate the applicability of the method when data are shifted in time. It is found that the corrected probability density for wall temperatures shifted up to 5 s backward or forward perfectly overlaps with the corrected probability density for wall temperatures that are not shifted in time. Upper and lower plots compare the probability density of the wall temperatures shifted backward by 10 s and 1 min respectively. Plots show that there is an acceptable difference in probability density when the wall temperatures are shifted by 10 s. For 1 min shift, however, the difference is significant. The method, therefore, may not work when the data are significantly shifted in time. In Sect. 4.2 we conclude that the coupled analysis approach for wall temperatures prediction is more accurate than the standalone approach. This could be confusing as the random model errors, e r fMg, presented in Table 4 is lowest for standalone model with boundary q 00 Exp . The basis for this reasoning is the comparatively high bias values for standalone model with boundary q 00 Exp presented in Table 4 and the first two plots in Figure 16. Similarly, the average model uncertainty values presented in Table 4 may not fully comply with the values presented in Figs. 15 and 16. This is because the values in in Table 4  Finally, the study presents predicted and corrected CDF of wall temperatures calculated for the input stochastics listed in Table 6. The proposed correction method handles only one type of different uncertainties appearing in a probabilistic simulation with deterministic models. Other uncertainty types, input uncertainty and sampling uncertainty deserve their own studies when aiming at accurate fire risk analyses. Figure 25 presents an overall procedure for uncertainty management in the stochastic simulation. Estimation of input uncertainty distribution is crucially important for the simulation outcome and can require significant effort if the number of uncertain parameters is high. Luckily, in a nonlinear system, such as fire, the number of dominating input parameters is usually small [24]. For sampling uncertainty, the convergence of the distribution moments can be studied, as explained in Sect. 2.4. This would be very expensive if a complex numerical method such as CFD is being used. Means to quantify the sampling convergence in LHS could possibly be developed using surrogate models, such as the response surface method.

Conclusion
In this work, we show that the model uncertainties based on the peak outputs and the current experimental data are similar to the ones estimated from the FDS validation database. The coupled analysis (FDS alone) had the smallest model uncertainties in wall temperatures. The higher uncertainties in the standalone analyses were caused by the high uncertainties of the heat flux, i.e. additional uncertainty propagation. The model uncertainties were found to vary over time, however, the probability correction using the generalized uncertainty parameters was effective at each time instance. The model uncertainties reported in the context of a model validation can be, therefore, used for correcting the output distributions resulting from parameter (input) uncertainty. Nevertheless, the proposed method for the model uncertainty compensation may not be effective when the model uncertainty cannot be generalized. Further work is needed to study the effect of Latin hypercube sampling uncertainty in failure probability calculation. Also, validation using larger experimental datasets and a wider range of output quantities would be valuable.

Open Access
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.