Predictive Modeling of a Buoyancy-Operated Cooling Tower under Unsaturated Conditions: Adjoint Sensitivity Model and Optimal Best-Estimate Results with Reduced Predicted Uncertainties

Nuclear and other large-scale energy-producing plants must include systems that guarantee the safe discharge of residual heat from the industrial process into the atmosphere. This function is usually performed by one or several cooling towers. The amount of heat released by a cooling tower into the external environment can be quantified by using a numerical simulation model of the physical processes occurring in the respective tower, augmented by experimentally measured data that accounts for external conditions such as outlet air temperature, outlet water temperature, and outlet air relative humidity. The model’s responses of interest depend on many model parameters including correlations, boundary conditions, and material properties. Changes in these model parameters induce changes in the computed quantities of interest (called “model responses”), which are quantified by the sensitivities (i.e., functional derivatives) of the model responses with respect to the model parameters. These sensitivities are computed in this work by applying the general adjoint sensitivity analysis methodology (ASAM) for nonlinear systems. These sensitivities are subsequently used for: (i) Ranking the parameters in their importance to contributing to response uncertainties; (ii) Propagating the uncertainties (covariances) in these model parameters to quantify the uncertainties (covariances) in the model responses; (iii) Performing model validation and predictive modeling. The comprehensive predictive modeling methodology used in this work, which includes assimilation of experimental measurements and calibration of model parameters, is applied to the cooling tower model under unsaturated conditions. The predicted response uncertainties (standard deviations) thus obtained are smaller than both the computed and the measured standards deviations for the respective responses, even for responses where no experimental data were available.


Introduction
Nuclear and other energy-producing large-scale industrial installations must include systems that guarantee the safe discharge of residual heat from the industrial process into the atmosphere. In the case of nuclear power plants, the residual heat is discharged into the atmosphere by using cooling towers. The amount of heat released into the external environment can be quantified by using a numerical simulation model of the physical processes that take place within the respective cooling tower, complemented by experimentally measured data that characterize external conditions such following (four) measured quantities: (i) Exhaust air temperature according to the "Tidbit" sensor, which will be denoted in the following as T a,out(Tidbit) ; (ii) Exhaust air temperature according to the "Hobo" sensor, which will be denoted in the following as T a,out(Hobo) ; (iii) Outlet water temperature, which will be denoted in the following as T meas w,out ; (iv) Exhaust air relative humidity, which will be referred to as RH meas . For a data set to be intended as unsaturated is required that the value of the exhaust air relative humidity (RH), computed by making use of the SNRL simulation code CTTool [1], is smaller than 100%, while the threshold is set on correspondence of the saturation line (RH = 100%). With this standard set, 6717 data sets among the total 8079 measured have matched the requirement above and have been therefore identified as "unsaturated", leading to include them in the present study. Histogram plots and analyses of the statistical moments of the selected 6717 data sets have been gathered in Appendix A. The state functions chosen for the analysis of the cooling tower mathematical model hereby considered have been selected according to the quantities comprised in the experimentally measured data sets. A numerical method, which has been presented in detail in [2], has been implemented to improve both the accuracy and the efficiency in the solution of the constitutive equations system governing the counter-flow cooling tower model. The cooling tower model analyzed in this work has been originally presented in [1] and has been thoroughly described in [2], from which Figures 1-3 have been taken. Natural draft air enters the tower just below the fill section, through which it passes flowing upwards along the height of the tower finally exiting at the top, through an exhaust comprising a fan. Inlet water enters the tower below the drift eliminator, and is sprayed downwards over the fill section, leading to the creation of a uniform film flow through the fill.
Energies 2016, 9,1028 3 of 57 "Tidbit" sensor, which will be denoted in the following as , ( ) a out Tidbit T ; (ii) Exhaust air temperature according to the "Hobo" sensor, which will be denoted in the following as , ( ) a out Hobo T ; (iii) Outlet water temperature, which will be denoted in the following as , meas w out T ; (iv) Exhaust air relative humidity, which will be referred to as meas RH . For a data set to be intended as unsaturated is required that the value of the exhaust air relative humidity ( ) RH , computed by making use of the SNRL simulation code CTTool [1], is smaller than 100%, while the threshold is set on correspondence of the saturation line . With this standard set, 6717 data sets among the total 8079 measured have matched the requirement above and have been therefore identified as "unsaturated", leading to include them in the present study. Histogram plots and analyses of the statistical moments of the selected 6717 data sets have been gathered in Appendix A. The state functions chosen for the analysis of the cooling tower mathematical model hereby considered have been selected according to the quantities comprised in the experimentally measured data sets. A numerical method, which has been presented in detail in [2], has been implemented to improve both the accuracy and the efficiency in the solution of the constitutive equations system governing the counter-flow cooling tower model. The cooling tower model analyzed in this work has been originally presented in [1] and has been thoroughly described in [2], from which Figures 1-3 have been taken. Natural draft air enters the tower just below the fill section, through which it passes flowing upwards along the height of the tower finally exiting at the top, through an exhaust comprising a fan. Inlet water enters the tower below the drift eliminator, and is sprayed downwards over the fill section, leading to the creation of a uniform film flow through the fill. The magnitude of the heat and mass transfer processes occurring in the counter-flow cooling tower of interest can be mathematically determined by obtaining the solution of the nonlinear system comprising the following balance equations: (A) Liquid continuity; (B) Liquid energy balance; (C) Water vapor continuity; (D) Air/water vapor energy balance; (E) Mechanical energy balance. In the derivation process of these equations, the cooling tower model has undergone several assumptions, such as: 1. air and water stream temperatures are uniform at any cross section; 2. the cross-sectional area of the cooling tower is assumed to be uniform; 3. the heat and mass transfer only occur in the direction normal to flows;    Here and in the following, the above state functions are assumed to be components of the following (column) vectors: † (2) ,..., ,..., In this paper's notation, the dagger ( ) † is used to denote "transposition", and all vectors in this work are column vectors. Because of the aforementioned similarities with the partially saturated case in [2], the constitutive equation system governing the cooling tower model of interest is closely comparable with the one provided in Equations (2)-(15) in [2]. For this unsaturated analysis, the governing equations within the total of I = 49 control volumes represented in Figure 2 are as follows [1]: A. Liquid Continuity Equations: (i) Control Volume i = 1:     Here and in the following, the above state functions are assumed to be components of the following (column) vectors: , ω ≡ ω (1) , ..., ω (I) † , m a (1) In this paper's notation, the dagger ( †) is used to denote "transposition", and all vectors in this work are column vectors. Because of the aforementioned similarities with the partially saturated case in [2], the constitutive equation system governing the cooling tower model of interest is closely comparable with the one provided in Equations (2)-(15) in [2]. For this unsaturated analysis, the governing equations within the total of I = 49 control volumes represented in Figure 2 are as follows [1]: Liquid Continuity Equations: (i) Control Volume i = 1: 1 (m w , T w , T a , ω, m a ; α) m (2) w − m w,in + M(m a ,α) R P (2) vs (T (2) w ,α) T (2) w − ω (1) P atm T (1) a (0.622+ω (1) Control Volumes i = 2,..., I − 1: Control Volume i = I: Liquid Energy Balance Equations: w )h (2) g,w (T (2) w , α) = 0; Control Volumes i = 2,..., I − 1: (iii) Control Volume i = I: C. Water Vapor Continuity Equations: (ii) Control Volumes i = 2,..., I − 1: (iii) Control Volume i = I: D.
The Air/Water Vapor Energy Balance Equations: (i) Control Volume i = 1: a , α) Control Volumes i = 2,..., I − 1: The vector α, which appears in Equations (2)- (14), comprises as its components the 47 model parameters, which will be denoted in the following as α i , i.e., where N α = 47 represents the total number of model parameters. These model parameters value and standard deviations have been experimentally derived, causing their statistical distributions not to be completely known; the first four statistical moments (means, variance/covariance, skewness, and kurtosis) of each of these parameter distributions have nevertheless been quantified, as presented in Appendix B. As discussed in [2], the numerical method used in the original work [1] to solve Equations (2)-(14) failed to achieve convergence for all the considered data sets, and for this reason it was substituted with a more accurate and efficient one, based on Newton's method together with the GMRES linear iterative solver for sparse matrices [9] comprised in the NSPCG package [10]. This method has been thoroughly described in [2], Section 2.1.
The Jacobian matrix of derivatives of Equations (2)- (14) with respect to the state functions is denoted as: The majority of the components of this block matrix remain unaltered from the Jacobian matrix in [2], Equation (20), whose components are detailed in [2], Appendix C; the components which are caused by the unsaturated operating conditions to be different from the study in [2] are separately detailed in Appendix C of this paper. The above matrix J (u n ) is a non-symmetric sparse matrix of order 197 by 197; the diagonal storage method used in NSPCG to relevantly reduce the matrix size, together with the approximation introduced by setting the column vectors (E 1 , ..., E 4 ) and the row vectors (A 5 , ..., D 5 ) to zero, are discussed in [2], Section 2.1.
As detailed in Appendix A, each of these data sets comprises measurements of the following quantities: (i) Outlet air temperature measured with the "Tidbit" sensor; (ii) Outlet air temperature measured with the "Hobo" sensor; (iii) Outlet water temperature; (iv) Outlet air relative humidity.
Accordingly to the aforementioned quantities contained in the benchmark data sets, the following responses of interest have been chosen for this work:  It is important to note that the water mass flow rates m (i) w , the water temperatures T (i) w , the air temperatures T (i) a and the air mass flow rate m a are computed by solving Equations (2)- (14), while the air relative humidity value, RH (i) , is obtained through the following expression: For the unperturbed base case, all model parameters (α i ) are used in solving Equations (2)-(14) with their nominal values, which are listed in Appendix B. It is worth specifying that the nominal values for parameters α 1 through α 5 (i.e., the dry bulb air temperature, dew point temperature, inlet water temperature, atmospheric pressure and wind speed) are statistically computed by averaging the values of the respective quantities in the 6717 data sets considered for this study.
The bar plots displayed below in Figures

Development of the Cooling Tower Adjoint Sensitivity Model
The development of the cooling tower adjoint sensitivity model follows the same path as in [2], Section 3.1, where the topic is discussed more in detail. The total sensitivity of a model response R (m w , T w , T a , ω, m a ; α), with respect to arbitrary variations in the model's parameters δα ≡ (δα 1 , ..., δα N α ) and state functions δm w , δT w , δT a , δω, δm a , around the nominal values m 0 w , T 0 w , T 0 a , ω 0 , m 0 a ; α 0 of the parameters and state functions, is obtained by means of the G-differential of the model's response to these changes. This G-differential is referred to as DR m 0 w , T 0 w , T 0 a , ω 0 , m 0 a ; α 0 ; δm w , δT w , δT a , δω, δm a ; δα , and introducing the adjoint sensitivity functions it becomes: where the "indirect effect" term, DR indirect , is obtained as: and where the vector [µ w , τ w , τ a , o, µ a ] † is required to be the solution of the following adjoint sensitivity system:  The vectors R ≡ r and where R 5 is defined as follows: while the vectors Q ≡ q (1) , ..., q (I) , = 1, 2, 3, 4 in Equation (19) comprise the derivatives of the model's equations with respect to model parameters, i.e.,: and where Q 5 is defined as follows: It is worth reminding that the adjoint sensitivity system in Equation (20) is independent of parameter variations. This feature allows the selected adjoint functions [µ w , τ w , τ a , o, µ a ] † to be computed by solving the adjoint sensitivity system just once.
Dimensional analysis allows determining the units of the adjoint functions from Equation (19). Namely, the units for the adjoint functions must satisfy the following relations: where [R] denotes the unit of the response R, while the units for the respective equations are as follows:  w ; (iii) The exit air humidity ratio R RH (1) ; (iv) The outlet (exit) water mass flow rate R m (50) w ; and (v) the air mass flow rate R m a . Let S j denote the "absolute sensitivity" of the response R with respect to the parameter α j , and is defined as:

Responses µ
An independent method to compute absolute response sensitivities S j is by making use of perturbative finite difference methods, such as: considering an arbitrarily small perturbation δα j to the model parameter α j ; 2. re-computing the perturbed response R α 0 j + δα j , where α 0 j denotes the unperturbed parameter value; 3.
using the finite difference formula 4. using the approximate equality between Equations (27) and (28) to obtain independently the respective values of the adjoint function(s) being verified.
The independent verification methodology discussed in steps (1)-(4) above will be clearly illustrated in Appendix D, where the adjoint functions depicted in Figures 8-12 will be verified.
4. using the approximate equality between Equations (27) and (28) to obtain independently the respective values of the adjoint function(s) being verified. ,..., ,...,  a R T , the value of the adjoint function μ a is −0.12651.      R RH , the value of the adjoint function μ a is −0.00743.   w R m , the value of the adjoint function μ a is −0.0306.
w , the value of the adjoint function µ a is −0.0306.
For the response  a R m , the value of the adjoint function μ a is 5.805.
The independent verification methodology discussed in steps (1)-(4) above will be clearly illustrated in Appendix D, where the adjoint functions depicted in Figures 8-12 will be verified.

Sensitivity Analysis Results and Rankings
As has been discussed above, there are a total of 8079 measured benchmark data sets for the cooling tower model operated in "fan-off" regime. As it has also been mentioned, 6717 benchmark data sets (out of the total of 8079 data sets) are considered to correspond to the "unsaturated conditions" which are analyzed in this work. The nominal values for boundary and atmospheric conditions used in this work were obtained from the statistics of these 6717 benchmark data sets corresponding to "unsaturated conditions". In turn, these "unsaturated" boundary and atmospheric conditions were used to obtain the sensitivity results reported in this Subsection. Sections 3.1.1-3.1.5, below, provide the numerical values and rankings, in descending order, of the relative sensitivities . For the response R m a , the value of the adjoint function µ a is 5.805.

Sensitivity Analysis Results and Rankings
As has been discussed above, there are a total of 8079 measured benchmark data sets for the cooling tower model operated in "fan-off" regime. As it has also been mentioned, 6717 benchmark data sets (out of the total of 8079 data sets) are considered to correspond to the "unsaturated conditions" which are analyzed in this work. The nominal values for boundary and atmospheric conditions used in this work were obtained from the statistics of these 6717 benchmark data sets corresponding to "unsaturated conditions". In turn, these "unsaturated" boundary and atmospheric conditions were used to obtain the sensitivity results reported in this Subsection. Sections 3.1.1-3.1.5, below, provide the numerical values and rankings, in descending order, of the relative sensitivities computed using the adjoint sensitivity analysis methodology for the five model responses T w , RH (1) , and m a . Note that the relative sensitivity RS (α i ) of a response R (α i ) to a parameter α i is defined as . Thus, the relative sensitivities are unit-less and are very useful in ranking the parameters to highlight their relative importance for the respective response. Thus, a relative sensitivity of 1.00 indicates that a change of 1% in the respective parameter will induce a 1% change in a response that is linear in the respective sensitivity. The higher the relative sensitivity, the more important the respective parameter to the respective response.
3.1.1. Sensitivity Analysis Results and Rankings for the Outlet Air Temperature, T (1) a Table 2 lists the sensitivities, computed using Equations (18) and (19), of the air outlet temperature with respect to all of the model's parameters. The parameters have been ranked according to the descending order of their relative sensitivities. As the results in Table 2 indicate, the first parameter (i.e., T w,in ) has a relative sensitivity around 90%, and is therefore the most important for the air outlet temperature response, T a , since that means that a 1% change in T w,in would induce a 0.91% change in T (1) a . The next four parameters (i.e., T db , T a,in , a 0 , T dp ) have relative sensitivities between 1% and 6%, and are therefore somewhat important. Parameters #6 through #9 (i.e., a 1 , V w , D h , D f an ) have relative sensitivities between 0.1% and 0.8%. The remaining 38 parameters are relatively unimportant for this response, having relative sensitivities smaller than 1% of the largest relative sensitivity (with respect to T a,in ) for this response. Positive sensitivities imply that a positive change in the respective parameter would cause an increase in the response, while negative sensitivities imply that a positive change in the respective parameter would cause a decrease in the response. The results and ranking of the relative sensitivities of the outlet water temperature with respect to the most important nine parameters for this response are listed in Table 3. The largest sensitivity of T (50) w is to the parameter T w,in , and has the value of 0.5055; this means that a 1% increase in T w,in would induce a 0.5055% increase in T (50) w . The sensitivities to the remaining 38 model parameters have not been listed since they are smaller than 1% of the largest sensitivity (with respect to T w,in ) for this response. The results and ranking of the relative sensitivities of the outlet water mass flow rate with respect to the most important 12 parameters for this response are listed in Table 4. This response is most sensitive to m w,in (a 1% increase in this parameter would cause a 1.01% increase in the response) and the second largest sensitivity is to the parameter T w,in (a 1% increase in this parameter would cause a 0.214% decrease in the response). The sensitivities to the remaining 35 model parameters have not been listed since they are smaller than 1% of the largest sensitivity (with respect to m w,in ) for this response. The results and ranking of the relative sensitivities of the outlet air relative humidity with respect to the most important 29 parameters for this response are listed in Table 5. The first three sensitivities of this response are the most relevant; in particular, an increase of 1% in T db or T a,in would cause an increase in the response of 0.27% or 0.25%, respectively. On the other hand, an increase of 1% in T w,in would cause a decrease of 0.32% in the response. The sensitivities to the remaining 18 model parameters have not been listed since they are smaller than 1% of the largest sensitivity (with respect to T w,in ) for this response. The results and ranking of the relative sensitivities of the air mass flow rate with respect to the most important 14 parameters for this response are listed in Table 6. The first three sensitivities of this response are very large (relative sensitivities larger than unity are customarily considered to be very significant). In particular, an increase of 1% in T a,in or T db would cause a decrease in the response of 38.51% or 38.49%, respectively. On the other hand, an increase of 1% in T w,in would cause an increase of 36% in the response. The sensitivities to the remaining 33 model parameters have not been listed since they are smaller than 1% of the largest sensitivity (with respect to T a,in ) for this response.
Overall, the air mass flow rate, m a , displays the largest sensitivities, so this response is the most sensitive to parameter variations. The other responses, namely the outlet air temperature, the outlet water temperature, the outlet water mass flow rate and the outlet air relative humidity display sensitivities of comparable magnitude.

Cross-Comparison of Sensitivity Results
In Tables 7-11, the ranked relative sensitivities for each response are compared side-by-side between three operating conditions, i.e., the two subcases discussed in [2] (partially saturated subcase I, completely saturated subcase II) and the unsaturated case analyzed in this paper. Among the three operating conditions, the "unsaturated case" is defined as a working condition in which air is unsaturated from the inlet to outlet of the cooling tower; while in the saturated subcase II, on the contrary, air is saturated from inlet to outlet of the cooling tower; the saturated subcase I is the combination of the these two cases, i.e., air in the lower portion of the fill section of the cooling tower is in unsaturated conditions, reaching saturation at some point along the height of the tower and remaining saturated in the upper part of the cooling tower. Cross-comparison of sensitivity results reveals the sensitivity variations between the three operating conditions.
The relative sensitivities and corresponding parameters listed in Table 7 are extracted from  Tables 1 and 6 in [3], and Table 2 in this paper. As shown in Table 7, for all three operating conditions, the first most sensitive parameters of the response of air outlet temperature, T a , is the same (i.e., T w,in ). The 2nd and 3rd most sensitive parameter are inverted in the unsaturated case with respect to the two subcases of the saturated case, but with values very close between the two parameters. The parameters that ranks in 4th place for this response is the same for all cases (i.e., a 0 ). The 5th parameter is a 0 for Subcases I and II and T dp for the unsaturated case.
For the first parameter (i.e., T w,in ), the unsaturated case displays the largest sensitivity for this response; subcase II has the smallest sensitivity; while subcase I has an intermediate value of sensitivity between the two. This is expected since subcase I is a mixed case between the unsaturated case and the saturated subcase II, as explained above. For all the remaining parameters in the table the situation is reversed, with Subcase II showing the largest sensitivity values and the unsaturated case presenting the smallest ones, with Subcase I still in the middle. Generally, the sensitivity magnitude of subcase I is slightly closer to that of subcase II. This can be explained by the fact that air remains unsaturated less than half of the height of the fill section, and flows in saturated conditions for more than half of the height of the fill section, as analyzed in [2].
The relative sensitivities and corresponding parameters listed in Table 8 are extracted from  Tables 2 and 7 in [3], and Table 3 in this paper. As shown in Table 8, for the response of water outlet temperature, T (50) w , both the unsaturated case and subcase I are most sensitive to the parameter T w,in , whereas subcase II is most sensitive to the parameter T a,in . As a comparison, the response of water outlet temperature to the parameter T w,in ranks in 3rd place, with a value comparable to the other two cases. The next two most sensitive parameters that rank from 2nd to 3rd places of this response are also different between the operating conditions: for both the unsaturated case and subcase I, parameters T a,in and T db rank in 2nd and 3rd places, respectively; however, for subcase II, parameters that take the 2nd and 3rd places are T db and T w,in , respectively. The parameters that take the 4th and 5th places are also different between the operating conditions, as shown in the table. Overall, for the response w , the sensitivity behavior of subcase I is more similar to that of the unsaturated case.
The relative sensitivities and corresponding parameters listed in Table 9 are extracted from Tables 3  and 8 in [3], and Table 4 in this paper. As shown in Table 9, for all three operating conditions, the first two most sensitive parameters of the response of water outlet mass flow rate, m (50) w , are the same (i.e., m w,in and T w,in , respectively). In addition, for each of the first two parameters, all three operating conditions have comparable sensitivity magnitudes. This indicates that the sensitivities of the first two parameters are insensitive to the operating condition change. The third most sensitive parameter of this response is different between the operating conditions: for both the unsaturated case and subcase I, this parameter is T dp ; whereas for subcase 2, this parameter is T a,in . Similarly, the parameters that take the 4th and 5th places are also different between the operating conditions, as shown in the table.
The relative sensitivities and corresponding parameters listed in Table 10 are extracted from Tables 4 and 9 in [3], and Table 5 in this paper. As shown in Table 10, for Subcases I and II, the first three most sensitive parameters of the response of air outlet relative humidity, RH (1) , are the same (i.e., T a,in , T db and T dp , respectively); the order is different for the unsaturated case. The next two most sensitive parameters that rank the 4th and 5th places of this response are different between the operating conditions.
For each of the first three parameters, all three operating conditions are sensitive to the parameter changes. In which, subcase II is the most sensitive case; and the unsaturated case is the least sensitive case comparatively. For instance, 1% change in T a,in , T db or T dp will cause around 0.2% change in RH (1) for the unsaturated case, around 2% change in RH (1) for subcase I; and nearly 15% change in RH (1) for subcase II, respectively. Overall, for the response of air outlet relative humidity, RH (1) , the sensitivity behavior of subcase I is also more similar to that of subcase II, as also the signs of most of the sensitivity values, inverted in the unsaturated case with respect to Subcase I and II, show in Table 10.
The relative sensitivities and corresponding parameters listed in Table 11 are extracted from  Tables 5 and 10 in [3], and Table 6 in this paper. As shown in Table 11, for all the operating conditions, the first three most sensitive parameters of the response of air mass flow rate, m a , are the same (i.e., T db , T a,in and T w,in , respectively) with the order of the first two being swapped for the unsaturated case. P atm is the 4th more sensitive parameter in all operating conditions, and with values comparable between the three cases; the parameters ranking in 5th place are different for the three operating conditions.
For each of the first three parameters, all three operating conditions are sensitive to the parameter changes. Differently from the response RH (1) , subcase II is this time the least sensitive case, while the unsaturated case is the most sensitive case comparatively. For instance, 1% change in T a,in , T db or T w,in will cause around 38% change in m a for the unsaturated case, around 24% change in m a for subcase I; and nearly 22% change in m a for subcase II, respectively. Overall, for the response of air mass flow rate, m a , the sensitivity behavior of subcase I is also more similar to that of subcase II.     where the measured correlated parameters are: α 1 T db , α 2 T dp , α 3 T w,in , α 4 P atm , and α 5 V w .
The a priori parameter covariance matrix, C αα , is: Energies 2016, 9,1028 24 of 52 The a priori covariance matrix of the computed responses, C comp rr , is given below: The best-estimate nominal parameter values have been computed as follows: in conjunction with the a priori matrices given in Equations (29)-(32) and the sensitivities presented in Tables 1-5. The resulting best-estimate nominal values are listed in Table 12, below. The corresponding best-estimate absolute standard deviations for these parameters are also presented in this table. These values are the square-roots of the diagonal elements of the matrix C pred αα . For comparison, the original nominal parameter values and original absolute standard deviations are also listed. As the results in Table 12 indicate, the predicted best-estimate standard deviations are all smaller or at most equal to (i.e., left unaffected) the original standard deviations. The parameters are affected proportionally to the magnitudes of their corresponding sensitivities: the parameters experiencing the largest reductions in their predicted standard deviations are those having the largest sensitivities.

Predicted Best-Estimated Response Values with Reduced Predicted Standard Deviations
The predicted response covariance matrix, C pred rr , is as follows: The non-zero elements with the largest magnitudes of best-estimate response-parameter correlation matrix, C The resulting best-estimate predicted nominal values are summarized in Table 13. To facilitate comparison, the corresponding measured and computed nominal values are also presented in this w . For this response, therefore, the predicted best-estimate nominal value has been obtained by a forward re-computation using the best-estimate nominal parameter values listed in Table 12, while the predicted best estimate standard deviation for this response has been obtained as follows: The results presented in Table 13 indicate that the predicted standard deviations are smaller than either the computed or the experimentally measured ones. This is indeed the consequence of using the PM_CMPS methodology in conjunction with consistent (as opposed to discrepant) computational and experimental information. Often, however, the information is inconsistent, usually due to the presence of unrecognized errors. Solutions for addressing such situations have been proposed in [10]. It is also important to note that the PM_CMPS methodology has improved (i.e., reduced, albeit not by a significant amount) the predicted standard deviation for the outlet water flow rate response, for which no measurements were available. This improvement stems from the global characteristics of the PM_CMPS methodology, which combines all of the available simultaneously on phase-space, as opposed to combining it sequentially, as is the case with the current state-of-the-art data assimilation procedures [11,12]. Table 13. Computed, measured, and optimal best-estimate nominal values and standard deviations for the outlet air temperature, outlet water temperature, outlet air relative humidity, outlet water mass flow rate and air mass flow rate responses.

Discussion
The original numerical method presented in [1] for the model solution has been replaced in this work with a considerably more accurate and efficient one which guarantees convergence of the computations for all of the available data sets. The adjoint model of the cooling tower has been implemented to compute exactly and efficiently the sensitivities of the model responses to all the 47 model parameters. The adjoint sensitivity model yields the adjoint state functions which are used to compute the sensitivities of each model response to all of the 47 model parameters by means of just one adjoint model computation. These adjoint state functions have been computed and their numerical accuracy has been independently verified. The response sensitivities to all model parameters have been computed for the following responses: (i) the outlet air temperature; (ii) the outlet water temperature; (iii) the outlet water mass flow rate; (iv) the air outlet relative humidity; and (v) air mass flow rate. Thes sensitivities have been subsequently used within the "predictive modeling for coupled multi-physics systems" (PM_CMPS) methodology [4] to obtain: (a) optimal best-estimate prediction for the model parameter values; (b) optimal best-estimate nominal values of the model responses; (c) reduced predicted standard deviations for the best-estimate calibrated model parameter values; and (d) reduced predicted standard deviations for the best-estimate predicted response values.
The results presented in this work show that the PM_CMPS methodology reduces the predicted standard deviation to values that are smaller than both the standard deviations of the measured and the computed response, respectively, even for responses, such as for the air mass flow rate, for which no experimentally measured values are available. This reduction stems from the fact that the PM_CMPS methodology simultaneously combines all the available data in the phase-space; customary data assimilation methodologies [11,12] only allow a sequential combination of the available information. All in all, the application of the PM_CMPS methodology has produced and improved, calibrated and validated model for simulating the functioning of a buoyancy-operated cooling tower under unsaturated conditions. Ongoing work aims at using second-order sensitivities, to be computed by applying the 2nd-ASAM presented in [13,14]. The availability of second-order response sensitivities will enable the computation of non-Gaussian features, such as skewness and kurtosis, of the response distributions of interest.

Appendix A. Statistical Analysis of Experimentally Measured Responses for SRNL F-Area Cooling Towers
Histogram plots of the 6717 measurement sets considered in this work (each set containing measurements of α , ( ) The measured outlet (exit) air relative humidity, meas RH , was obtained using Hobo humidity sensors. The accuracy of these sensors is depicted in Figure A1, which indicates the following tolerances (standard deviations): ±2.5% for relative humidity from 10% to 90%; between ±2.5% and ±3.5% for relative humidity from 90% to 95%; and ±3.5%-±4.0% from 95% to 100%. However, when exposed to relative humidity above 95%, the maximum sensor error may temporally increase by an additional 1%, so that the error can reach values between ±4.5% and ±5.0% for relative humidity from 95% to 100%. Figure A1. Humidity sensor accuracy plot (adopted from the specification of HOBO Pro v2).
As shown in this Figure A2, although the computed relative humidity for each of the 6717 data sets is less than 100%, the measured relative humidity Consequently, all the 6717 benchmark data sets plotted in Figure A1, were considered as Figure A1. Humidity sensor accuracy plot (adopted from the specification of HOBO Pro v2).
As shown in this Figure A2, although the computed relative humidity for each of the 6717 data sets is less than 100%, the measured relative humidity RH meas actually spans the range from 33.0% to 104.1%; in this range, 4925 data sets have their respective RH meas less than 100% while the other Energies 2016, 9,1028 28 of 52 1792 data sets have their respective RH meas over 100%. This situation is nevertheless consistent with the range of the sensors when their tolerances (standard deviations) are taken into account, which would make it possible for a measurement with RH meas = 105% to be nevertheless "unsaturated". Consequently, all the 6717 benchmark data sets plotted in Figure A1, were considered as "unsaturated", since their respective RH meas was less than 105%. The statistical properties of the (measured air outlet relative humidity) distribution shown in Figure A2 have been computed using standard packages, and are presented in Table A1. These statistical properties will be needed for the uncertainty quantification and predictive modeling computations presented in the main body of this work.   The statistical properties of the (measured air outlet relative humidity) distribution shown in Figure A2 have been computed using standard packages, and are presented in Table A1. These statistical properties will be needed for the uncertainty quantification and predictive modeling computations presented in the main body of this work. The histogram plots and their corresponding statistical characteristics of the 6717 data sets for the other measurements, namely for: the outlet air temperature [T a,out(Tidbit) ] measured using the "Tidbit" sensors; the outlet air temperature [T a,out(Hobo) ] measured using the "Hobo" sensors; and the outlet water temperature [T meas w,out ] are reported below in Figures A3-A6, and Tables A2-A5, respectively. The statistical properties of the (measured air outlet relative humidity) distribution shown in Figure A2 have been computed using standard packages, and are presented in Table A1. These statistical properties will be needed for the uncertainty quantification and predictive modeling computations presented in the main body of this work.  Figure A3. Histogram plot of the air outlet temperature measured using "Tidbit" sensors, within the 6717 data sets collected by SRNL from F-Area cooling towers. Figure A3. Histogram plot of the air outlet temperature measured using "Tidbit" sensors, within the 6717 data sets collected by SRNL from F-Area cooling towers.              Ordering the above-mentioned four measured responses as follows: (i) outlet air temperature T a,out(Tidbit) ; (ii) outlet air temperature T a,out(Hobo) ; (iii) outlet water temperature T meas w,out ; and (iv) outlet air relative humidity RH meas out , yields the following "measured response covariance matrix", denoted as Cov T a,out(Tidbit) , T a,out(Hobo) , T meas w,out , RH meas out : Cov (A1) For the purposes of uncertainty quantification, data assimilation, model calibration and predictive modeling, the temperatures measurements provided by the "Tidbit" and "Hobo" sensors can be combined into an "averaged" data set of measured air outlet temperatures, which will be denoted as T meas a,out . The histogram plot and corresponding statistical characteristics of this averaged air outlet temperature are presented in Figure A6 and A comparison between the results in Equations (A1) and (A2) makes clear that the elimination of the second column and row in Equation (A1) yields a 3-by-3 matrix which has entries basically equivalent to the covariance matrix shown in Equation (A2). Therefore, this means that the temperature distributions measured by the "Tidbit" and "Hobo" sensors do not need to be dealt with as separate data sets for the purposes of uncertainty quantification and predictive modeling.
The standard deviation of the humidity sensor utilized for the measurements (σ sensor = 5.0% for the response RH (1) ) have been already considered by including in the category of the "unsaturated" data sets those that have their respective measured relative humidity, RH meas , up to 105.0%. In addition to that, the respective uncertainties of the temperature sensors (standard deviations, σ sensor = 0.2K for both responses T (1) a and T (50) w ) must also be taken into consideration for the 6717 data sets. The measuring methods and devices are not dependent with respect to each other, therefore the data standard deviation σ statistic , stemming from the statistical analysis of the 6717 benchmark data sets, and the sensor standard deviation, σ sensor , stemming from the instrument's uncertainty, must stack according to the well-known formula of "addition of the variances of uncorrelated variates", i.e.,: (A3) Coupling the above relation with the result presented in Equation (A2) will lead to incremented values of the variances on the diagonal of the respective "measured covariance matrix"; this new form of the covariance matrix which will be denoted as Cov T meas a,out , T meas w,out , RH meas out . The obtained result is: Cov T meas a,out , T meas w,out , RH meas out =    In the predictive modeling formalism (which includes uncertainty quantification, data assimilation, and model calibration) the covariance matrix between the measured parameters and responses is required as an input. In the case of interest, all the parameters and responses can be considered as uncorrelated, except for the measured responses considered in this Appendix and the measured parameters listed in Appendix B. The "parameter-response" covariance matrix in Equation (A5), indicated as Cov T meas a,out , T meas w,out , RH meas , α 1 , ..., α 47 , refers to the above mentioned parameters (namely: dry-bulb air temperature, T db ; dew-point air temperature, T dp , inlet water temperature, T w,in , atmospheric pressure, P atm , and wind speed V w ) and responses (i.e., average outlet air temperature, outlet water temperature, and outlet air relative humidity):

Appendix B. Model Parameters for the SRNL F-Area Cooling Towers
The mean values and standard deviations for the independent model parameters α i , ( i = 1, ..., N α = 47) ,presented in Table B1, below, have been derived in collaboration with Dr. Sebastian Aleman of SRNL (private communications, 2016).  The above independent model parameters are used for computing various dependent model parameters and thermal material properties, as shown in Tables B2 and B3, below.

Dependent Scalar Parameters Math. Notation Defining Equation or Correlation
Mass diffusivity of water vapor in air (m 2 /s) D av (T a , α) a 0,dav T 1.5 a 1,dav +(a 2,dav +a 3,dav T)T Air velocity in the fill section (m/s) v a (m a , α) A out π r f an 2 Table B3. Thermal Properties (Dependent Scalar Model Parameters).

Thermal Properties (Functions of State Variables) Math. Notation Defining Equation or Correlation
a 0g + a 1g T w H g (T a ) = saturated vapor enthalpy (J/kg) h g,a (T a , α) a 0g + a 1g T a C p (T) = specific heat of dry air (J/kg K) C p (T, α) a 0,cpa + (a 1,cpa + a 2,cpa T)T P vs (T w ) = saturation pressure (Pa) P vs (T w , α) P c · e a 0 + a 1 Tw , in which P c = 1.0 Pa P vs (T a ) = saturation pressure (Pa) P vs (T a , α) P c · e a 0 + a 1 Ta , in which P c = 1.0 Pa Note 1: The measurements of parameters α 1 -α 5 (i.e., the dry bulb air temperature, dew point temperature, inlet water temperature, atmospheric pressure and wind speed) were taken at the SRNL site, where the F-area cooling towers are located. Out of the 8079 total benchmark data sets [8], 6717 data sets have been considered in this study, since "unsaturated"; through these data sets the statistical properties (means, variance and covariance, skewness and kurtosis) for these model parameters have been derived, as shown in Figures B1-B5 and Tables B4-B8. αα 5 (i.e., the dry bulb air temperature, dew point temperature, inlet water temperature, atmospheric pressure and wind speed) were taken at the SRNL site, where the F-area cooling towers are located. Out of the 8079 total benchmark data sets [8], 6717 data sets have been considered in this study, since "unsaturated"; through these data sets the statistical properties (means, variance and covariance, skewness and kurtosis) for these model parameters have been derived, as shown in Figures B1-B5 and Tables B4-B8. Figure B1. Histogram plot of dry-bulb air temperature data collected by SRNL from F-Area cooling towers. Figure B1. Histogram plot of dry-bulb air temperature data collected by SRNL from F-Area cooling towers. Figure B1. Histogram plot of dry-bulb air temperature data collected by SRNL from F-Area cooling towers. Figure B2. Histogram plot of dew-point air temperature data collected by SRNL from F-Area cooling towers.         Figure B5. Histogram plot of wind speed data collected by SRNL from F-Area cooling towers. The 5-by-5 covariance matrix for the above experimental data has also been computed and is provided below, with the four model parameters ordered as follows: dry-bulb air temperature T db , dew-point air temperature T dp , inlet water temperature T w,in , atmospheric air pressure P atm , and wind speed V w .

Minimum Maximum Range Mean Std. Dev. Variance Skewness Kurtosis
Cov T db ; T dp ; T w,in ; P atm ; V w = The covariance matrix (above) neglects the uncertainty associated with sensor readings throughout the data collection period. When combining uncertainties by adding variances, the contribution from the sensors is 0.04 K for each of the first three parameters, which accounts for a maximum of ca. 1% of the total variance (for the inlet water temperature, specifically). The uncertainty in the atmospheric pressure sensor is at this time unknown. For these reasons, their contribution to overall uncertainty is considered insignificant at this time.

Appendix C. Derivative Matrix (Jacobian) of the Model Equations with Respect to the State Functions
As mentioned in Section 2, the Jacobian matrix presents similarities with the Jacobian matrix detailed in [2], Equation (C1). More precisely, the sub-matrices (A i , B i , C i , D i , E i ; i = 2, 3, 4, 5) whose elements represents the derivatives of Equations (6)- (14) with respect to the vector valued state function u (m w , T w , T a , ω, m a ) † remain the same as in [2]; for reasons of brevity, they have not been reported in this paper. On the other side, the sub-matrices (A i , B i , C i , D i , E i ; i = 1) whose elements represents the derivatives of Equations (2)-(4) with respect to the vector valued state function u (m w , T w , T a , ω, m a ) † , are different from their respective formulations in [2], and therefore they will be hereby detailed with the following notation: For subsequent use, the above quantities are considered to be the components of the I × I matrix A 1 defined as follows: For subsequent use, the above quantities are considered to be the components of the I × I matrix B 1 defined as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to T (j) a are as follows: ∂N For subsequent use, the above quantities are considered to be the components of the I × I matrix C 1 defined as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to ω (j) are as follows: ∂N For subsequent use, the above quantities are considered to be the components of the I × I matrix D 1 defined as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to m a are as follow follows: (1) For Re d < 2300 (2) For 2300 ≤ Re d ≤ 10, 000 (3) For Re d > 10, 000 For subsequent use, the above quantities are considered to be the components of the I column vector E 1 defined as follows: (C22)

Appendix D. Derivatives of Cooling Tower Model Equations With Respect to Model Parameters
The differences between the governing equations for this study and for subcase I in [2] concern only the "liquid continuity equations". Other governing Equations (i.e., liquid energy balance equations; water vapor continuity equations; and the air/water vapor energy balance equations) are the same for both cases.
For this reason, the derivatives of Equations (6)- (14) with respect to the model parameters remain the same as in [3], Equation (A3); for reasons of brevity, they have not been reported in this paper. On the other side, the derivatives of Equations (2)-(4) with respect to the model parameters are different from their respective formulations in [3], and therefore they will be hereby detailed with the following notation: For the sake of brevity, only the nonzero derivatives have been reported in this appendix. The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (1) : T db are as follows: where: ∂M(m a , α) ∂D av (T db , α) The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (4) : P atm are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (7) : µ are as follows: (D6) where: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (8) : υ are as follows: where: ∂M(m a , α) ∂υ The derivatives of the "liquid continuity equations" (cf. Equations (A1)-(A4)) with respect to the parameter α (11) : f mt are as follows: where: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (13) : a 0 are as follows: where: ∂P The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (14) : a 1 are as follows: where: ∂P The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (18) : a 0,dav are as follows: ∂D av (T db ,α) was defined previously in Equation (D3), and The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (19) : a 1,dav are as follows: where ∂M(m a ,α) ∂D av (T db ,α) was defined previously in Equation (D3), and The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (20) : a 2,dav are as follows: ∂D av (T db ,α) was defined previously in Equation (D3), and The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (21) : a 3,dav are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (26) : a 0,Nu are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (27) : a 1,Nu are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (28) : a 2,Nu are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (29) : a 3,Nu are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (39) : D h are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (40) : A f ill are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (41) : A sur f are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (42) : Pr are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (43) : w tsa are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (44) : m w,in are as follows: The derivatives of the "liquid continuity equations" (cf. Equations (2)-(4)) with respect to the parameter α (47) : Sc are as follows: where ∂M(m a , α) ∂Sc Re-writing Equation (E1) in the form indicates that the value of the adjoint function µ a could be computed independently if the sensitivity S 5 were available, since the quantity ∂N 5 /∂V w = −2.1795 J/ m 4 /s is known. To first-order in the parameter perturbation, the finite-difference formula given in Equation (28) can be used to compute the approximate sensitivity S FD 5 ; subsequently, this value can be used in conjunction with Equation (E2) to compute a "finite-difference sensitivity" value, denoted as [µ a ] SFD , for the respective adjoint, which would be accurate up to second-order in the respective parameter perturbation: Numerically, the wind speed V w has the nominal ("base-case") value of V 0 w = 1.859 [m/s]. The corresponding nominal value T The same parameter perturbation was utilized to perform the same verification procedure for the adjoint function µ a with respect to the other four responses; Table E1 displays the obtained results, which compare well with the values in the bar plots in Figures 9-12. Note that the value of the adjoint function o (49) obtained by solving the adjoint sensitivity system given in Equation (20) is o (49) = −1.313 × 10 −5 [K/ (J/kg)], as indicated in Figure 8. Now select a variation δT a,in in the inlet air temperature T a,in , and note that Equation (27) yields the following expression for the sensitivity of the response R = T (1) a to T a,in : /∂T a,in = 1.0309 × 10 3 [J/ (kg · K)] and ∂N 5 /∂T a,in = 0.40491 J/ m 3 · K are known. To first-order in the parameter perturbation, the finite-difference formula given in Equation (28) can be used to compute the approximate sensitivity S FD 45 ; subsequently, this value can be used in conjunction with Equation (E5) to compute a "finite-difference sensitivity" value, denoted as o (49) SFD , for the respective adjoint, which would be accurate up to second-order in the respective parameter perturbation:      g,a (T a,in , α) is known. To first-order in the parameter perturbation, the finite-difference formula given in Equation (28) can be used to compute the approximate sensitivity S FD 46 ; subsequently, this value can be used in conjunction with Equation (E8) to compute a "finite-difference sensitivity" value, denoted as τ