Predictive Modeling of a Paradigm Mechanical Cooling Tower: I. Adjoint Sensitivity Model

Cooling towers discharge waste heat from an industrial process into the atmosphere, and are essential for the functioning of large energy-producing plants, including nuclear reactors. Using a numerical simulation model of the cooling tower together with measurements of outlet air relative humidity, outlet air and water temperatures enables the quantification of the rate of thermal energy dissipation removed from the respective process. The computed quantities depend on many model parameters including correlations, boundary conditions, material properties, etc. Changes in these model parameters will induce changes in the computed quantities of interest (called “model responses”). These changes are quantified by the functional derivatives (called “sensitivities”) 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 needed 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 predictive modeling, including assimilation of experimental measurements and calibration of model parameters to produce optimal predicted quantities (both model parameters and responses) with reduced predicted uncertainties.


Introduction
A mechanical draft cooling tower (MDCT) discharges waste heat from an industrial process into the atmosphere.Cooling towers are essential for the functioning of large energy-producing plants, including nuclear reactors.Using a numerical simulation model of the cooling tower together with measurements of outlet air relative humidity, outlet air and water temperatures enables the quantification of the rate of thermal energy dissipation removed from the respective process.In addition to computing the temperature drop of the cooling water as it passes through the tower, a cooling tower model that derives heat dissipation rates from thermal imagery needs to convert the remotely measured cooling tower throat or area-weighted temperature to a cooling water inlet temperature.Therefore, a cooling tower model comprises two main components, namely: (i) an inner model which computes the amount of cooling undergone by the water as it passes through the tower as a function of inlet cooling water temperature and ambient weather conditions (air temperature and humidity); and (ii) an outer model which uses a remotely measured throat or area-weighted temperature and iterates on the inlet water temperature to match the target temperature of interest.The cooling tower model produces an estimate of the rate at which energy is being discharged to the atmosphere by evaporation and sensible heat transfer.The sensible heat transfer is estimated using the computed change in air or water enthalpy as it passes through the MDCT.If the MDCT fans are on, a prescribed mass flow rate of air and water is used.If the MDCT fans are off, an additional mechanical energy equation is iteratively solved to determine the mass flow rate of air.
The flow regime in the fill section of a cooling tower, which can be cross-flow or counter-flow, determines the type of the respective cooling tower.A model for computing the steady-state thermal performance has been recently presented [1]; this model simulates both cross-flow and counter-flow wet mechanical draft cooling towers.Using as inputs the temperature and mass flow rate of the incoming water together with the temperature and humidity ratio of the incoming ambient air, this model computes the temperature and mass flow rate of the effluent water, as well as the temperature and water vapor content of the exhaust air.The air mass flow rate is specified when the cooling tower operates in the mechanical draft mode.When the fan is turned-off, the cooling tower operates in the natural draft/wind-aided mode, in which case the air mass flow rate is calculated using the numerical model.
In the present work, the counter-flow cooling tower presented in [1] is further developed by applying a considerably more accurate and efficient numerical method for computing the steady state distributions of the following quantities: (i) the water mass flow rates at the exit of each control volume along the height of the fill section of the cooling tower; (ii) the water temperatures at the exit of each control volume along the height of the fill section of the cooling tower; (iii) the air temperatures at the exit of each control volume along the height of the fill section of the cooling tower; (iv) the humidity ratios at the exit of each control volume along the height of the fill section of the cooling tower; and (v) the air relative humidity at the exit of each control volume.In contrast to many cases of non-convergence experienced by the original method [1], the numerical method applied in this work is highly efficient and yields accurate results everywhere, as shown in Section 2.1, below.
More importantly, however, this work presents the development of the adjoint cooling tower sensitivity model for computing efficiently and exactly the sensitivities (i.e., functional derivatives) of the model responses (i.e., quantities of interest) to all 52 model parameters.The adjoint sensitivity model is developed by applying the general Adjoint Sensitivity Analysis Methodology (ASAM) for nonlinear systems, which was originally presented in [2][3][4].Even though the forward cooling tower model is nonlinear in the state functions, the adjoint sensitivity model is linear in the adjoint state functions, which correspond one-to-one to the forward state functions mentioned in the foregoing.Using the adjoint state functions, the sensitivities of each model response to all of the 52 model parameters can be computed exactly using a single adjoint model computation, as opposed to 52 computations, as would be required if the sensitivities would have been computed-approximately as opposed to exactly-using the forward model in conjunction with finite-differences.The sensitivities are needed for (i) ranking the parameters in the order of their importance for contributing to response uncertainties; (ii) propagating the uncertainties (variances and covariances) in the model parameters to quantify the uncertainties (variances and covariances) in the model responses; (iii) performing predictive modeling, which includes assimilation of experimental measurements and calibration of model parameters to produce optimally predicted nominal values for both model parameters and responses, with reduced predicted uncertainties.The development of the adjoint sensitivity model for the counter-flow cooling tower is presented in Section 2 of this work, together with results for the respective adjoint state functions.A discussion of the significance of the results obtained in this work is presented in Section 3. The predictive modeling of the counter-flow cooling tower will be performed by applying the "predictive modeling for coupled multi-physics systems" (PM_CMPS) methodology recently developed in [5]; the results of this predictive modeling will presented in the sequel to this paper [6].

Results
During the period from April, 2004 through August, 2004, a total 8079 measured benchmark data sets for F-area cooling towers (fan-on case) were recorded every fifteen minutes at SRNL for F-Area Cooling Towers [7].Each of these data sets contained measurements of the following (four) Energies 2016, 9, 718 3 of 45 quantities: (i) outlet air temperature measured with the sensor called "Tidbit", which will be denoted as T a,out(Tidbit) ; (ii) outlet air temperature measured with the sensor called "Hobo", which will be denoted as T a,out(Hobo) ; (iii) outlet water temperature, which will be denoted as T meas w,out ; (iv) outlet air relative humidity, which will be denoted as RH meas .Histogram plots of these 7668 measurement sets (each set containing measurements of T a,out(Tidbit) , T a,out(Hobo) , T meas w,out and RH meas ), together with statistical analyses of these measurements, are presented in Appendix A. These measured quantities provide the basis for the choosing the state functions underlying the mathematical modeling of the cooling tower, which is presented in Section 2.1.An accurate and efficient numerical method for solving the equations underlying the counter-flow cooling tower is also presented in this section.Section 2.2 presents the development of the cooling tower adjoint sensitivity model, along with the solution method for computing the adjoint state functions.The numerical accuracy of solving the equations underlying the adjoint sensitivity model will also be verified.

Mathematical Model of the Counter-Flow Cooling Tower
The counter-flow cooling tower considered in this work has been originally developed in [1] and is schematically presented in Figure 1.
Energies 2016, 9, 718 3 of 48 denoted as Ta,out(Hobo); (iii) outlet water temperature, which will be denoted as  ,  ; (iv) outlet air relative humidity, which will be denoted as RH meas .Histogram plots of these 7668 measurement sets (each set containing measurements of Ta,out(Tidbit), Ta,out(Hobo),  ,  and RH meas ), together with statistical analyses of these measurements, are presented in Appendix A. These measured quantities provide the basis for the choosing the state functions underlying the mathematical modeling of the cooling tower, which is presented in Section 2.1.An accurate and efficient numerical method for solving the equations underlying the counter-flow cooling tower is also presented in this section.Section 2.2 presents the development of the cooling tower adjoint sensitivity model, along with the solution method for computing the adjoint state functions.The numerical accuracy of solving the equations underlying the adjoint sensitivity model will also be verified.

Mathematical Model of the Counter-Flow Cooling Tower
The counter-flow cooling tower considered in this work has been originally developed in [1] and is schematically presented in Figure 1.As this figure indicates, forced air flow enters the tower through the "rain section" above the water basin, flows upward through the fill section and the drift eliminator, and exits at the tower's top through an exhaust that encloses a fan.Hot water enters above the fill section and is sprayed onto the top of the fill section to create a uniform, downward falling, film flow through the fill's numerous meandering vertical passages.Film fills are designed to maximize the water free surface area and the residence time inside of the fill section.Heat and mass transfer occurs at the falling film's free surface between the water film and the upward air flow.The drift eliminator above the spray zone removes entrained water droplets from the upward flowing air.Below the fill section, the water droplets fall into a collection basin, placed at the bottom of the cooling tower.The heat and mass transfer processes occur overwhelmingly in the fill section.Modeling the heat and mass transfer processes between falling water film and rising air in the cooling tower's fill section is accomplished by solving the following balance equations: (A) liquid continuity; (B) liquid energy balance; (C) water vapor continuity; (D) air/water vapor energy balance.The assumptions used in deriving these equations are as follows: As this figure indicates, forced air flow enters the tower through the "rain section" above the water basin, flows upward through the fill section and the drift eliminator, and exits at the tower's top through an exhaust that encloses a fan.Hot water enters above the fill section and is sprayed onto the top of the fill section to create a uniform, downward falling, film flow through the fill's numerous meandering vertical passages.Film fills are designed to maximize the water free surface area and the residence time inside of the fill section.Heat and mass transfer occurs at the falling film's free surface between the water film and the upward air flow.The drift eliminator above the spray zone removes entrained water droplets from the upward flowing air.Below the fill section, the water droplets fall into a collection basin, placed at the bottom of the cooling tower.The heat and mass transfer processes occur overwhelmingly in the fill section.Modeling the heat and mass transfer processes between falling water film and rising air in the cooling tower's fill section is accomplished by solving the following balance equations: the air and/or water temperatures are uniform throughout each stream at any cross section; 2.
the cooling tower has uniform cross-sectional area; 3.
the heat and mass transfer occur solely in the direction normal to flows; 4.
the heat and mass transfer through tower walls to the environment is negligible; 5.
the heat transfer from the cooling tower fan and motor assembly to the air is negligible; 6.
the air and water vapor mix as ideal gasses; 7.
the flow between flat plates is unsaturated through the fill section.
This work considers moderately sized towers for which the heat and mass transfer processes in the rain section is negligible.The fill section is modeled by discretizing it in vertically stacked control volumes as depicted in Figure 2. The heat and mass transfer between the falling water film and the rising air in a typical control volume of the cooling tower's fill section is presented in Figure 3. 2. the cooling tower has uniform cross-sectional area; 3. the heat and mass transfer occur solely in the direction normal to flows; 4. the heat and mass transfer through tower walls to the environment is negligible; 5. the heat transfer from the cooling tower fan and motor assembly to the air is negligible; 6. the air and water vapor mix as ideal gasses; 7. the flow between flat plates is unsaturated through the fill section.
This work considers moderately sized towers for which the heat and mass transfer processes in the rain section is negligible.The fill section is modeled by discretizing it in vertically stacked control volumes as depicted in Figure 2. The heat and mass transfer between the falling water film and the rising air in a typical control volume of the cooling tower's fill section is presented in Figure 3.
Energies 2016, 9, 718 4 of 48 2. the cooling tower has uniform cross-sectional area; 3. the heat and mass transfer occur solely in the direction normal to flows; 4. the heat and mass transfer through tower walls to the environment is negligible; 5. the heat transfer from the cooling tower fan and motor assembly to the air is negligible; 6. the air and water vapor mix as ideal gasses; 7. the flow between flat plates is unsaturated through the fill section.This work considers moderately sized towers for which the heat and mass transfer processes in the rain section is negligible.The fill section is modeled by discretizing it in vertically stacked control volumes as depicted in Figure 2. The heat and mass transfer between the falling water film and the rising air in a typical control volume of the cooling tower's fill section is presented in Figure 3.In the mechanical draft mode, the mass flow rate of dry air is specified.With the fan off and hot water flowing through the cooling tower, air will continue to flow through the tower due to buoyancy.Wind pressure at the air inlet to the cooling tower will also enhance air flow through the tower.The air flow rate is determined from the overall mechanical energy equation for the dry air flow.The state functions underlying the cooling tower model (cf., Figures 1-3) are as follows: 1.
the water mass flow rates, denoted as m w (i = 2, . . ., 50), at the exit of each control volume, i, along the height of the fill section of the cooling tower;

2.
the water temperatures, denoted as w (i = 2, . . ., 50), at the exit of each control volume, i, along the height of the fill section of the cooling tower;

3.
the air temperatures, denoted as a (i = 1, . . ., 49), at the exit of each control volume, i, along the height of the fill section of the cooling tower; and 4.
the humidity ratios, denoted as .ω (i) (i = 1, . . ., 49), at the exit of each control volume, i, along the height of the fill section of the cooling tower.
It is convenient to consider the above state functions to be components of the following (column) vectors: In this work, the dagger (+) will be used to denote "transposition", and all vectors will be considered to be column vectors.The governing conservation equations within the total of I = 49 control volumes represented in Figure 2 are as follows [1]: A Liquid continuity equations: 1 (m w , T w , T a , ω; α) m (2) w ,α) T (2) w − ω (1) P atm T (1) a (0.622+ω (1) Control Volumes i = 2,..., I−1: Control Volume i = I:

D
The air/water vapor energy balance equations: a , α) a )H(m a ,α) g,w (T w ,α) Control Volumes i = 2,..., I−1: Control Volume i = I: The components of the vector α, which appears in Equations ( 2)-( 13), comprise the model parameters which are generically denoted as α i , i.e.: where N α denotes the total number of model parameters.These model parameters are experimentally derived quantities, and their complete distributions parameters are not known; however, we have determined the first four moments (means, variance/covariance, skewness, and kurtosis) of each of these parameter distributions, as detailed in Appendix B.
In the original work [1], the above equations were solved using a two stage-iterative method comprising an "inner-iteration" using Newton's method within each control volume, followed by an outer iteration aimed at achieving overall convergence.This procedure, though, did not converge at all of the points of interest.Therefore, after testing several alternatives provided in [8] and [9], we have replaced the original solution method in [1] by Newton's method together with the GMRES linear iterative solver for sparse matrices [10] provided in the NSPCG package [9], which turned out to be the most efficient among those we have tested.This GMRES method [10] approximates the exact solution-vector of a linear system by using the Arnoldi iteration to find the approximate solution-vector by minimizing the norm of the residual vector over a Krylov subspace.The specific computational steps are as follows: (a) Write Equations ( 1)-(13) in vector form as: where the following definitions were used: 4 , . . ., N , u (m w , T w , T a , ω) † (16) (b) Set the initial guess, u 0 , to be the inlet boundary conditions; (c) Steps d through g, below, constitute the outer iteration loop; for n = 0, 1, 2, . .., iterate over the following steps until convergence: (d) Inner iteration loop: for m = 1, 2, . . ., use the iterative GMRES linear solver with the Modified Incomplete Cholesky (MIC) preconditioner, with restarts, to solve, until convergence, the following system to compute the vector δu: where n is the current outer loop iteration number, and the Jacobian matrix of derivatives of Equations ( 2)-(13) with respect to the state functions is the block-matrix: with matrix-components defined in Appendix C. As shown in this Appendix, the Jacobian represented by Equation ( 18) is a non-symmetric sparse matrix of order 196 by 196, with 14 nonzero diagonals.The non-symmetric diagonal storage format is used to store the respective 14 nonzero diagonals, so that the "condensed" Jacobian matrix has dimensions 196 by 14.Since the Jacobian is highly non-symmetric, the cost of the iterations of the GMRES solver grows as O(m 2 ), where m is the iteration number within the GMRES solver.To reduce this computational cost, the GMRES solver is configured to run with the restart feature.The optimized value for the restart frequency is 10 for this specific application.The MIC preconditioner can speed up the convergence of the GMRES solver using the parameters OMEGA and LVFILL [9] in the modified incomplete factorization methods for the MIC preconditioner; for this application the following values were found to be optimal: OMEGA = 0.000000001 and LVFILL = 1.The Jacobian is not updated inside the sparse GMRES solver.The default convergence of GMRES is tested with the following criterion [9]: δu (m) , δu (m)   where z (m) denotes the pseudo-residual at m th -iteration of the GMRES solver, δu (m) is the solution of Equation ( 17) at m th -iteration, and ζ denotes the stopping test value for the GMRES solver.
(e) Set where n is the current outer loop iteration number, and update the Jacobian.(f) Test for convergence of the outer loop until the error in the solution is less than a specified maximum value.For solving Equations ( 2)-(13), the following error criterion has been used: (g) Set n = n + 1 and go to step d.
The above solution strategy for solving Equations ( 2)-( 13) converged successfully for all the 8079 benchmark data sets.As described previously, each of these data sets contained measurements of the following quantities: (i) outlet air temperature measured with the sensor called "Tidbit"; (ii) outlet air temperature measured with the sensor called "Hobo"; (iii) outlet water temperature; (iv) outlet air relative humidity.For each of these benchmark data sets, the outer loop iterations described above (i.e., steps c through g) converge in 4 iterations; for each outer loop iteration, the GMRES solver used for solving Equation (17) converges in 12 iterations.The "zero-to-zero" verification of the solution's accuracy using Equations ( 2) through (13) gives an error of the order of 10 −7 .
In view of the above-mentioned measurements, the responses of interest for this work are as follows: of air temperatures at the exit of each control volume i, (i = 1, . . ., 49); (d) the vector RH RH (1) , . . ., RH (I) † , having components of the air relative humidity at the exit of each control volume i, (i = 1, . . ., 49).
While the water mass flow rates m w , the water temperatures w , and the air temperatures a are obtained directly as the solutions of Equations ( 2)-(13), the air relative humidity, RH (i) , is computed for each control volume using the expression: The bar plots, showing the respective values of the water mass flow rates m w , the water temperatures T (i) w , the air temperatures T (i) a , and the air relative humidity, RH (i) , at the exit of each control volume, are presented in Figures 4-7, below.

Development of the Cooling Tower Adjoint Sensitivity Model, with Solution Verification
All of the responses of interest in this work, e.g., the experimentally measured and/or computed responses discussed in the previous Sections, can be generally represented in the functional form , where R is a known functional of the model's state functions and parameters.
As generally shown in [2], the sensitivity of such as this response to arbitrary variations in the model's parameters   where the so-called "direct effect" term, direct DR , and the so-called "indirect effect" term, indirect DR , are defined, respectively, as follows: where the components of the vectors are defined as follows: Since the model parameters are related to the model's state functions through Equations ( 2)-(13), it follows that variations in the model parameter will induce variations in the state variables.More precisely, it has been shown in [2][3][4] that to first-order in the parameter variations, the

Development of the Cooling Tower Adjoint Sensitivity Model, with Solution Verification
All of the responses of interest in this work, e.g., the experimentally measured and/or computed responses discussed in the previous Sections, can be generally represented in the functional form R (m w , T w , T a , ω; α), where R is a known functional of the model's state functions and parameters.As generally shown in [2], the sensitivity of such as this response to arbitrary variations in the model's parameters δα ≡ (δα 1 , . . . ,δα N α ) and state functions δm w , δT w , δT a , δω is provided by the response's Gateaux (G-) differential DR m 0 w , T 0 w , T 0 a , ω 0 ; α 0 ; δm w , δT w , δT a , δω; δα , which is defined as follows: DR m 0 w , T 0 w , T 0 a , ω 0 ; α 0 ; δm w , δT w , δT a , δω; δα where the so-called "direct effect" term, DR direct , and the so-called "indirect effect" term, DR indirect , are defined, respectively, as follows: where the components of the vectors R ≡ r (1) , . . ., r (I) , = 1, 2, 3, 4 are defined as follows: Since the model parameters are related to the model's state functions through Equations ( 2)-( 13), it follows that variations in the model parameter will induce variations in the state variables.More precisely, it has been shown in [2][3][4] that to first-order in the parameter variations, the respective variations in the state variables can be computed by solving the G-differentiated model equations, namely: Performing the above differentiation on Equations ( 2) through (13) yields the following forward sensitivity system: where the components of the vectors Q ≡ q (1) , . . ., q (I) , = 1, 2, 3, 4 are defined as follows: The system represented by Equation ( 26) is called the forward sensitivity system, which can be solved, in principle, to compute the variations in the state functions for every variation in the model parameters.In turn, the solution of Equation ( 26) can be used in Equation ( 24) to compute the "indirect effect" term, DR indirect .However, since there are many parameter variations to consider, solving Equation (26) repeatedly to compute DR indirect becomes computationally impracticable.The need for solving Equation (26) repeatedly to compute DR indirect can be circumvented by applying the Adjoint Sensitivity Analysis Procedure (ASAM) formulated in [2][3][4].The ASAM proceeds by forming the inner-product of Equation ( 26) with a yet unspecified vector of the form [µ w , τ w , τ a , o] † , having the same structure as the vector u (m w , T w , T a , ω) † , transposing the resulting scalar equation and using Equation (24).Furthermore, by requiring that the vector [µ w , τ w , τ a , o] † satisfy the following adjoint sensitivity system: it ultimately results that the "indirect effect" term can be expressed in the form The system represented by Equation ( 28) is called the adjoint sensitivity system, which -notablyis independent of parameter variations.Therefore, the adjoint sensitivity system needs to be solved only once, to compute the adjoint functions [µ w , τ w , τ a , o] † .In turn, the adjoint functions are used to compute DR indirect , efficiently and exactly, using Equation (29).As an illustrative example of computing response sensitivities using the adjoint sensitivity system, consider that the model response of interest is the air relative humidity, RH (i) , in a generic control volume i, as given by Equation ( 21).
For this model response, the "direct effect" term, denoted as D RH (i) direct , is readily obtained in the form: where: On the other hand, the "indirect effect" term, denoted as D RH (i)

indirect
, is readily obtained in the form: where: The units of the adjoint functions can be determined from Equation (29) through dimensional analysis.Specifically, the units for the adjoint functions satisfy the following relation: where "[R]" denotes the unit of the response R, while the units for the respective equations are as follows: Table 1 below lists the units of the adjoint functions for four responses: R T a , R T w , R RH (1) and R m w , respectively, in which, T a denotes exit air temperature; T (50) w denotes exit water temperature; RH (1) denotes exit air relative humidity; and m (50) w denotes exit water mass flow rate.

Responses µ
Note that the adjoint sensitivity system represented by Equation ( 28) is linear in the adjoint state functions, so it can be solved by using numerical methods appropriate for large-scale sparse linear systems.In particular, we solved it by using NSPCG, a "Package for Solving Large Sparse Linear w ; (iii) the exit air humidity ratio R RH (1) ; and (iv) the outlet (exit) water mass flow rate R m (50) w , are presented in Figures 8-11. (1) Note that the adjoint sensitivity system represented by Equation ( 28) is linear in the adjoint state functions, so it can be solved by using numerical methods appropriate for large-scale sparse linear systems.In particular, we solved it by using NSPCG, a "Package for Solving Large Sparse Linear Systems by Various Iterative Methods" [9]; 12 to 18 iterations sufficed for solving the adjoint system within convergence criterion of            The numerical accuracy of the computed adjoint functions can be independently verified by first noting that Equations ( 22), ( 23) and (29) imply that: where N α denotes the total number of model parameters, and S j denotes the "absolute sensitivity" of the response R with respect to the parameter α j , and is defined as: All of the derivatives with respect to the model parameter α j on the right side of Equation ( 40) are known quantities.On the other hand, the absolute response sensitivity S j can be computed independently, as follows: 1.
consider a small perturbation δα j in the model parameter α j ; 2. re-compute the perturbed response R α 0 j + δα j , where α 0 j denotes the unperturbed parameter value; 3.
use the finite difference formula: use the approximate equality between Equations ( 41) and ( 40) to obtain independently the respective values of the adjoint function(s) being verified.
The verification procedure described in steps (1)-( 4) above, will be illustrated in the remainder of this section, for verifying the adjoint functions depicted in Figures 8-11

∂R/∂T
(1) a = 1.Thus, the adjoint functions corresponding to the outlet air temperature response T a are computed by solving the adjoint sensitivity system given in Equation (28) using r (1) 3

∂R/∂T
(1) a = 1 as the only non-zero source term; for this case, the solution of Equation (28) has been depicted in Figure 8.
(a) Verification of the adjoint function o (49)   Note that the value of the adjoint function o (49) obtained by solving the adjoint sensitivity system given in Equation ( 28) is o (49) = −4.687× 10 −4 [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 (40) yields the following expression for the sensitivity of the response R = T a to T a,in : /∂T a,in = 1.0309 × 10 3 [J/ (kg • K)] is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation (41) can be used to compute the approximate sensitivity S FD 46 ; subsequently, this value can be used in conjunction with Equation (43) 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: Numerically, the inlet air temperature T a,in (= T db ) has the nominal ("base-case") value of T 0 a,in = 299.11[K].2)-( 13) with the value of T pert a,in yields the "perturbed response" value T = 978.20 [K], as indicated in Figure 8. Now select a variation δω in in the inlet air humidity ratio ω in , and note that Equation (40) yields the following expression for the sensitivity of the response R = T Re-writing Equation (45) in the form: Energies 2016, 9, 718 18 of 45 indicates that the value of the adjoint function τ (49) a could be computed independently if the sensitivity S 48 were available, since the o (49) has been verified in (the previous) Section 2.2.1 (a) and the quantity h (50) g,a (T a,in , α) is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation (41) can be used to compute the approximate sensitivity S FD 48 ; subsequently, this value can be used in conjunction with Equation (46) to compute a "finite-difference sensitivity" value, denoted as τ (49) a SFD , for the respective adjoint, which would be accurate up to second-order in the respective parameter perturbation: Numerically, the inlet air humidity ratio ω in has the nominal ("base-case") value of a,nom of the response T a is T = 978.20 [K] obtained by solving the adjoint sensitivity system given in Equation (28), cf. Figure 8.When solving this adjoint sensitivity system, the computation of τ Note that the values of the adjoint functions τ (1) w and µ (1) w obtained by solving the adjoint sensitivity system given in Equation (28) are as follows: τ (1) w = −1.49067.× 10 −6 [K/(J/s) ] and µ (1) w = 3.5735 [K/ (kg/s)], respectively, as indicated in Figure 8. Now select a variation δT w,in in the inlet water temperature T w,in , and note that Equation (40) yields the following expression for the sensitivity of the response R = T (1) a to T w,in : a to m w,in : Numerically, the inlet water temperature, T w,in , has the nominal ("base-case") value of T 0 w,in = 298.79[K], while the nominal ("base-case") value of the inlet water mass flow rate is m 0 w,in = 44.021229[kg/s].As before, the corresponding nominal value T a,nom of the response T a is T

∂R/∂T
(50) w = 1.. as the only non-zero source term; for this case, the solution of Equation (28) has been depicted in Figure 9.
(a) Verification of the adjoint function o (49)   Note that the value of the adjoint function o (49) obtained by solving the adjoint sensitivity system given in Equation ( 28) is o (49) = −2.214× 10 −4 [K/ (J/kg)], as indicated in Figure 9. Now select Energies 2016, 9, 718 20 of 45 a variation δT a,in in the inlet air temperature T a,in , and note that Equation (40) yields the following expression for the sensitivity of the response R = T (50) w to T a,in : Re-writing Equation (52) in the form: indicates that the value of the adjoint function o (49) could be computed independently if the sensitivity S 46 were available, since the quantity ∂N (49) 4 /∂T a,in = 1.0309 × 10 3 [J/ (kg • K)] is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation ( 41) can be used to compute the approximate sensitivity S FD 46 ; subsequently, this value can be used in conjunction with Equation (53) 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: w,pert − T to ω in : g,a (T a,in , α) . (55) Re-writing Equation (55) in the form indicates that the value of the adjoint function τ (49) a could be computed independently if the sensitivity S 48 were available, since the o (49) has been verified in (the previous) Section 2.2.2(a) and the quantity h (50) g,a (T a,in , α) is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation ( 41) can be used to compute the approximate sensitivity S FD 48 ; subsequently, this value can be used in conjunction with Equation (56) to compute a "finite-difference sensitivity" value, denoted as τ (49) a SFD , for the respective adjoint, which would be accurate up to second-order in the respective parameter perturbation: Numerically, the inlet air humidity ratio ω in has the nominal ("base-case") value of The corresponding nominal value T = −76.12[K] obtained by solving the adjoint sensitivity system given in Equation (28), cf. Figure 9.When solving this adjoint sensitivity system, the computation of τ to T w,in : [K/ (J/kg)] can also be considered as being accurate, since they constitute the starting point for solving the adjoint sensitivity system in Equation (28).Hence, the unknowns in Equation ( 58) are the adjoint functions µ Numerically, the inlet water temperature, T w,in , has the nominal ("base-case") value of T 0 w,in = 298.79[K], while the nominal ("base-case") value of the inlet water mass flow rate is m 0 w,in = 44.021229[kg/s].As before, the corresponding nominal value T    w = −5.730× 10 −7 [K/(J/s)], respectively, which are obtained by solving the adjoint sensitivity system given in Equation (28), cf. Figure 9.
(a) Verification of the adjoint function o (49)   Note that the value of the adjoint function o (49) obtained by solving the adjoint sensitivity system given in Equation ( 28) is o (49) = 1.860 × 10 −5 (J/kg) −1 , as indicated in Figure 10.Now select a variation δT a,in in the inlet air temperature T a,in , and note that Equation (40) yields the following expression for the sensitivity of the response R = RH (1) to T a,in : /∂T a,in = 1.0309 × 10 3 [J/ (kg • K)] is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation (41) can be used to compute the approximate sensitivity S FD 46 ; subsequently, this value can be used in conjunction with Equation (65) 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 parameter perturbation: Numerically, the inlet air temperature T a,in (= T db ) has the nominal ("base-case") value of T 0 a,in = 299.11[K].The corresponding nominal value RH (1) nom of the response RH (1) is RH = 1.815 × 10 −5 (J/kg) −1 .This result compares well with the value o (49) = 1.860 × 10 −5 (J/kg) −1 obtained by solving the adjoint sensitivity system given in Equation (28), cf., Figure 10.When solving this adjoint sensitivity system, the computation of o (49) depends on the previously computed adjoint functions o (i) , i = 1, . . ., I − 1; hence, the forgoing verification of the computational accuracy of o (49) also provides an indirect verification that the functions o (i) , i = 1, . . ., I − 1, were also computed accurately.

(b) Verification of the adjoint function τ (49) a
Note that the value of the adjoint function τ (49) a obtained by solving the adjoint sensitivity system given in Equation ( 28) is τ (49) a = −67.047,as indicated in Figure 10.Now select a variation δω in in the inlet air humidity ratio ω in , and note that Equation (40) yields the following expression for the sensitivity of the response R = RH (1) to ω in : Re-writing Equation (67) in the form indicates that the value of the adjoint function τ (49) a could be computed independently if the sensitivity S 48 were available, since the o (49) has been verified in (the previous) Section 2.2.3(a) and the quantity h (50) g,a (T a,in , α) is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation (41) can be used to compute the approximate sensitivity S FD 48 ; subsequently, this value can be used in conjunction with Equation (68) to compute a "finite-difference sensitivity" value, denoted as τ (49) a SFD , for the respective adjoint, which would be accurate up to second-order in the parameter perturbation: Numerically, the inlet air humidity ratio ω in has the nominal ("base-case") value of nom of the response RH (1) is RH w = −0.2926(kg/s) −1 , respectively, as indicated in Figure 10.Now select a variation δT w,in in the inlet water temperature T w,in , and note that Equation (40) yields the following expression for the sensitivity of the response R = RH (1) to T w,in : = −172.515and o (1) = 4.816 × 10 −5 (J/kg) −1 can also be considered as being accurate, since they constitute the starting point for solving the adjoint sensitivity system in Equation (28).Hence, the unknowns in Equation ( 70) are the adjoint functions µ w .A second equation involving solely these adjoint functions can be derived by selecting a perturbation, δm w,in , in the inlet water mass flow rate, m w,in , for which Equation (40) yields the following expression for the sensitivity of the response R = RH (1)  to m w,in : Numerically, the inlet water temperature, T w,in , has the nominal ("base-case") value of T 0 w,in = 298.79[K], while the nominal ("base-case") value of the inlet water mass flow rate is m 0 w,in = 44.021229[kg/s].As before, the corresponding nominal value RH (1) nom of the response RH (1) is RH

∂R/∂m
(50) w = 1.as the only non-zero source term; for this case, the solution of Equation ( 28) has been depicted in Figure 11.
(a) Verification of the adjoint function o (49)   Note that the value of the adjoint function o (49) obtained by solving the adjoint sensitivity system given in Equation ( 28) is o (49) = 1.603 × 10 −5 [(kg/s) / (J/kg)], as indicated in Figure 11.Now select a variation δT a,in in the inlet air temperature T a,in , and note that Equation (40) yields the following expression for the sensitivity of the response R = m (50) w to T a,in : /∂T a,in = 1.0309 × 10 3 [J/ (kg • K)] is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation ( 41) can be used to compute the approximate sensitivity S FD 46 ; subsequently, this value can be used in conjunction with Equation (75) 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 parameter perturbation: to ω in : Re-writing Equation (77) in the form: indicates that the value of the adjoint function τ (49) a could be computed independently if the sensitivity S 48 were available, since the o (49) has been verified in (the previous) Section 2.2.4(a) and the quantity h (50) g,a (T a,in , α) is known.To first-order in the parameter perturbation, the finite-difference formula given in Equation (41) can be used to compute the approximate sensitivity S FD 48 ; subsequently, this value can be used in conjunction with Equation (78) to compute a "finite-difference sensitivity" value, denoted as τ (49) a SFD , for the respective adjoint, which would be accurate up to second-order in the parameter perturbation: Numerically, the inlet air humidity ratio ω in has the nominal ("base-case") value of w obtained by solving the adjoint sensitivity system given in Equation (28) are as follows: τ (1) w = 2.6710 × 10 −7 (J/kg) −1 and µ (1) w = 0.3377, respectively, as indicated in Figure 11.Now select a variation δT w,in in the inlet water temperature T w,in , and note that Equation (40) yields the following expression for the sensitivity of the response Since the adjoint functions τ (49) a and o (49) have been already verified as described in Section 2.2.4(a) and (b), it follows that the computed values of adjoint functions τ (1) a = −3.1344[kg/s] o (1) = 8.1328 × 10 −7 [(kg/s) / (J/kg)] can also be considered as being accurate, since they constitute the starting point for solving the adjoint sensitivity system in Equation (28).Hence, the unknowns in Equation ( 80) are the adjoint functions µ Numerically, the inlet water temperature, T w,in , has the nominal ("base-case") value of T 0 w,in = 298.79[K], while the nominal ("base-case") value of the inlet water mass flow rate is m 0 w,in = 44.021229[kg/s].As before, the corresponding nominal value m  These values compare well with the values µ (1) w = 0.3377 and τ (1) w = 2.6710 × 10 −7 (J/kg) −1 , respectively, which are obtained by solving the adjoint sensitivity system given in Equation (28), cf. Figure 11.

Discussion
Based on the counter-flow cooling tower presented in [1], this work introduced a considerably more accurate and efficient numerical method for computing the steady state distributions of the following quantities: (i) the water mass flow rates at the exit of each control volume along the height of the fill section of the cooling tower; (ii) the water temperatures at the exit of each control volume along the height of the fill section of the cooling tower; (iii) the air temperatures at the exit of each control volume along the height of the fill section of the cooling tower; (iv) the humidity ratios at the exit of each control volume along the height of the fill section of the cooling tower; and (v) the air relative humidity at the exit of each control volume.Subsequently the adjoint cooling tower sensitivity model was conceived in this work for computing efficiently and exactly the sensitivities (i.e., functional derivatives) of the model responses (i.e., quantities of interest) to all 52 model parameters.The adjoint cooling tower model was developed by applying the general adjoint sensitivity analysis methodology (ASAM) for nonlinear systems, which was originally presented in [2][3][4].The adjoint sensitivity model is linear in the adjoint state functions, which correspond one-to-one to the forward state functions.Using the adjoint state variables, the sensitivities of each model response to all of the 52 model parameters can be computed exactly using a single adjoint model computation, as opposed to 52 computations, as would be required if the sensitivities had been computed using the forward model in conjunction with approximate finite-differences.The various adjoint state functions were computed and their accuracy was verified, thus setting the stage for the specific numerical computation of the respective response sensitivities.By applying the "predictive modeling for coupled multi-physics systems" (PM_CMPS) methodology recently developed in [5], the sequel to this paper [6] will compute numerically the sensitivities, which will be used for the following purposes (i) ranking the parameters in the order of their importance for contributing to response uncertainties; (ii) propagating the uncertainties (variances and covariances) in the model parameters to quantify the uncertainties (variances and covariances) in the model responses; (iii) performing predictive modeling, which includes assimilation of experimental measurements and calibration of model parameters to produce optimally predicted nominal values for both model parameters and responses, with reduced predicted uncertainties. is considered to be saturated.Applying this criterion to all of the 8079 benchmark data sets ("fan-on", with the average exhaust air velocity at the shroud equal to 10 m/s) provided in [7] for F-area cooling towers leads to the identification of 7668 measured data sets that fall into the "unsaturated" case analyzed in this work.Histogram plots of these 7668 measurement sets (each set containing measurements of T a,out(Tidbit) , T a,out(Hobo) , T meas w,out and RH meas ), together with statistical analyses thereof are presented in the remainder of this Appendix.
The measured outlet (exit) air relative humidity, RH meas , 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% to ±5.0% for relative humidity from 95% to 100%.The 7668 measured values of the outlet (exit) air relative humidity, RH meas , considered according to the results produced by CTTool code [1] to be "unsaturated," are presented in the histogram plot shown in Figure A2.As shown in this figure, although the computed relative humidity for each of the 7668 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, 6975 data sets have their respective RH meas less than 100% while the other 693 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 7668 benchmark data sets plotted in Figure A2 were considered as "unsaturated", since their respective RH meas was less than 105%.This plot, as well as all of the other histogram plots in this work, have their total respective areas normalized to unity.The 7668 measured values of the outlet (exit) air relative humidity, RH meas , considered according to the results produced by CTTool code [1] to be "unsaturated," are presented in the histogram plot shown in Figure A2.As shown in this figure, although the computed relative humidity for each of the 7668 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, 6975 data sets have their respective RH meas less than 100% while the other 693 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 7668 benchmark data sets plotted in Figure A2 were considered as "unsaturated", since their respective RH meas was less than 105%.This plot, as well as all of the other histogram plots in this work, have their total respective areas normalized to unity.
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 7668 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 other 693 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 7668 benchmark data sets plotted in Figure A2 were considered as "unsaturated", since their respective RH meas was less than 105%.This plot, as well as all of the other histogram plots in this work, have their total respective areas normalized to unity.The statistical properties of the (measured air outlet relative humidity) distribution shown in Figures 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 7668 data sets for the other measurements, namely for: the outlet air temperature [Ta,out(Tidbit)] measured using the "Tidbit" sensors; the outlet air temperature [Ta,out(Hobo)] measured using the "Hobo" sensors; and the outlet water temperature [ , 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  Comparing the results in Equations (A1) and (A2) shows that eliminating the second column and row in Equation (A1) yields a 3-by-3 matrix which has entries essentially equivalent to the covariance matrix in Equation (A2).In turn, this result indicates that the temperature distributions measured by the "Tidbit" and "Hobo" sensors, respectively, need not be treated as separate data sets for the purposes of uncertainty quantification and predictive modeling.
The sensors' standard deviations (namely: σ sensor = 0.2K for each of the responses T w , and σ sensor = 2.8% for the response RH (1) ) have been taken into account for the data at the 100%-saturation point, by including the 693 data sets that have their respective measured relative humidity, RH meas , between 100% and 104.1%.In addition, the respective sensors' uncertainties (standard deviations) must also be taken into account for the 6975 data sets that have their respective RH meas less than 100%.Since the various measuring methods and devices are independent of each other, the standard deviation, σ statistic , stemming from the statistical analysis of the 7668 benchmark data sets and the standard deviation, σ sensor , stemming from the instrument's uncertainty are to be combined according to the well-known formula "addition of the variances of uncorrelated variates", namely: )

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 α = 52) , 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.

Thermal Properties (Functions of State Variables) Math. Notation Defining Equation or Correlation
h f (T w ) = saturated liquid enthalpy (J/kg) h f (T w , α) a 0 f + a 1 f T w H g (T w ) = saturated vapor enthalpy (J/kg) h g,w (T w , α) 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 parameters α 1 through α 4 (i.e., the dry bulb air temperature, dew point temperature, inlet water temperature, and atmospheric pressure) were measured at the SRNL site at which the F-area cooling towers are located.Among the 8079 measured benchmark data sets [7], 7688 data sets are considered to represent "unsaturated conditions", which have been used to derive the statistical properties (means, variance and covariance, skewness and kurtosis) for these model parameters, as shown below in Figures B1-B4 and Tables B4-B7 .    (B1) 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.ω in = 0.622P vs (T dp , α) P atm − P vs (T dp , α) = 0.622e Note 4: The Reynold's number is defined (equivalently) as follows: Note 5: The Schmidt number is defined as follows: Note 6: The Sherwood number is defined as follows: (B5) Note 7: The Nusselt number is defined as follows: Re d > 10000 (B6) Note 8: The overall uncertainty of the Nusselt number is estimated to be 25% in the laminar flow region, 40% in the turbulent flow region, and the uncertainty is assumed to increase linearly in the transition region.This is a more conservative estimate of the overall uncertainty than assuming 25% for all three regimes.
For subsequent use, the above quantities are considered to be the components of the I × I matrix A 1 defined as follows: The derivatives of the "liquid continuity equations" [cf.Equations ( 2)-( 4)] with respect to T The derivatives of the "liquid continuity equations" [cf.Equations ( 2)-( 4)] with respect to T ; i = 1, . . ., I; j = i.

(C13)
For subsequent use, the above quantities are considered to be the components of the I × I diagonal 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 diagonal matrix D 1 defined as follows: The derivatives of the liquid energy balance equations [cf.Equations ( 5)-( 7)] with respect to m For subsequent use, the above quantities are considered to be the components of the I × I matrix A 2 defined as follows:  The derivatives of the liquid energy balance equations [cf.Equations ( 5)- (7)] with respect to T The derivatives of the liquid energy balance equations [cf.Equations ( 5)- (7)] with respect to T The derivatives of the liquid energy balance equations [cf.Equations ( 5)-( 7)] 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 A 3 defined as follows: The derivatives of the water vapor continuity equations [cf.Equations ( 8)- (10)] with respect to w are as follows: ∂N For subsequent use, the above quantities are considered to be the components of the I × I matrix The derivatives of the water vapor continuity equations [cf.Equations ( 8)- (10)] with respect to a are as follows: ∂N For subsequent use, the above quantities are considered to be the components of the I × I matrix The derivatives of the water vapor continuity equations [cf.Equations ( 8)- (10)] with respect to ω (j) are as follows: ∂N

Figure 1 .
Figure 1.Flow through a counter-flow cooling tower.

Figure 1 .
Figure 1.Flow through a counter-flow cooling tower.
(A) liquid continuity; (B) liquid energy balance; (C) water vapor continuity; (D) air/water vapor energy balance.The assumptions used in deriving these equations are as follows:

Figure 2 .
Figure 2. Control volumes   1,..,  i I comprising the counter-flow cooling tower, together with the symbols denoting the forward state functions   ( ) ( ) ( ) ( ) , , , , 1,..,  i i i i w w a T T i I m 

Figure 3 .
Figure 3. Heat and mass transfer between falling water film and rising air in a typical control volume of the cooling tower's fill section.

Figure 2 .
Figure 2. Control volumes (i = 1, .., I) comprising the counter-flow cooling tower, together with the symbols denoting the forward state functions m (i) w , T

Figure 2 .
Figure 2. Control volumes   1,..,  i I comprising the counter-flow cooling tower, together with the symbols denoting the forward state functions   ( ) ( ) ( ) ( ) , , , , 1,..,  i i i i w w a T T i I m 

Figure 3 .
Figure 3. Heat and mass transfer between falling water film and rising air in a typical control volume of the cooling tower's fill section.

Figure 3 .
Figure 3. Heat and mass transfer between falling water film and rising air in a typical control volume of the cooling tower's fill section.
flow rates at the exit of each control volume i, (i = 1, . . ., 49); (b) the vector T w T at the exit of each control volume i, (i = 1, . . ., 49); (c) the vector T a T

Figure 4 .
Figure 4. Bar plot of the water mass flow rates () i w m , ( 2,...,50)  i , at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 6 .
Figure 6.Bar plot of the air temperatures () i a T , ( 1,...,49) i  , at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 4 .
Figure 4. Bar plot of the water mass flow rates m (i) w , (i = 2, . . ., 50), at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 4 .
Figure 4. Bar plot of the water mass flow rates () i w m , ( 2,...,50)  i , at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 6 .
Figure 6.Bar plot of the air temperatures () i a T , ( 1,...,49) i  , at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 5 .
Figure 5. Bar plot of the water temperatures T (i) w , (i = 2, . . ., 50), at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 4 .
Figure 4. Bar plot of the water mass flow rates () i w m , ( 2,...,50)  i , at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 6 .
Figure 6.Bar plot of the air temperatures () i a T , ( 1,...,49) i  , at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 6 .
Figure 6.Bar plot of the air temperatures T (i) a , (i = 1, . . ., 49), at the exit of each control volume along the height of the fill section of the cooling tower.

Figure 7 .
Figure 7. Bar plot of the air relative humidity

Figure 7 .
Figure 7. Bar plot of the air relative humidity RH (i) , (i = 1, . . ., 49), at the exit of each control volume along the height of the fill section of the cooling tower.
Iterative Methods"[9]; 12 to 18 iterations sufficed for solving the adjoint system within convergence criterion of ζ = 10 −12 .The bar plots of the adjoint functions corresponding to the four measured responses of interest, namely: (i) the exit air temperature R T (1) a ; (ii) the outlet (exit) water temperature R T (50)

.
The bar plots of the adjoint functions corresponding to the four measured responses of interest, namely: (i) the exit air temperature (1) a RT; (ii) the outlet (exit) water temperature (50) w RT ; (iii) the exit air humidity ratio (1) R RH ; and (iv) the outlet (exit) water mass flow rate (50) w Rm , are presented in Figures 8-11.

Figure 8 .
Figure 8. Bar plots of adjoint functions for the response

Figure 9 .Figure 9 .
Figure 9. Bar plots of adjoint functions for the response R T (50) w as functions of the height

Figure 10 .
Figure 10.Bar plots of adjoint functions for the response

Figure 11 .
Figure 11.Bar plots of adjoint functions for the response

Figure 11 .
Figure 11.Bar plots of adjoint functions for the response

2. 2 . 1 .
Verification of the Adjoint Functions for the Outlet Air Temperature Response T (1) a When R = T (1) a , the quantities r (i) defined in Equation (25) all vanish except for a single component, namely: r (1) 3

( 1 )
a,pert = 297.6073[K].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD 46 T (1) a,pert −T (1) a,nom δT a,in = 0.4802.Using this value together with the nominal values of the other quantities appearing in the expression on the right side of Equation (44) yields o (49) SFD = −4.658× 10 −4 [K/ (J/kg)].This result compares well with the value o (49) = −4.687× 10 −4 [K/ (J/kg)] obtained by solving the adjoint sensitivity system given in Equation (28), cf., Figure 8.When solving this adjoint sensitivity system, the computation of o (49) depends on the previously computed adjoint functions o (i) , i = 1, . . ., I − 1; hence, the forgoing verification of the computational accuracy of o (49) also provides an indirect verification that the functions o (i) , i = 1, . . ., I − 1, were also computed accurately.(b) Verification of the adjoint function τ (49) a Note that the value of the adjoint function τ (49) a obtained by solving the adjoint sensitivity system given in Equation (28) is τ (49) a

( 1 )
a,nom = 297.4637[K].Consider next a perturbation δω in = (0.0063) ω 0 in , for which the perturbed value of the inlet air humidity ratio becomes ω pert in = ω 0 in + δω in = 0.0138612.Re-computing the perturbed response by solving Equations (2)-(13) with the value of ω pert in yields the "perturbed response" value T (1) a,pert = 297.4824[K].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD 48 T (1) a,pert −T (1) a,nom δω in = 216.32[K].Using this value together with the nominal values of the other quantities appearing in the expression on the right side of Equation (47) yields τ (49) a SFD = 978.25 [K].This result compares well with the value τ (49) a a , i = 1, . . ., I − 1; hence, the forgoing verification of the computational accuracy of τ (49) a also provides an indirect verification that the functions τ (i) a , i = 1, . . ., I − 1 were also computed accurately.(c)Verification of the adjoint functions τ

4 ∂T
a and o(49) have been already verified as described in Section 2.2.1 (a) and (b), it follows that the computed values of adjoint functions τ (1)a = 2410.83[K] o (1) = −9.5142× 10 −4 [K/ (J/kg)]can also be considered as being accurate, since they constitute the starting point for solving the adjoint sensitivity system in Equation (28).Hence, the unknowns in Equation (48) are the adjoint functions µ (1) w and τ (1) w .A second equation involving solely these adjoint functions Energies 2016, 9, 718 19 of 45 can be derived by selecting a perturbation, δm w,in , in the inlet water mass flow rate, m w,in , for which Equation (40) yields the following expression for the sensitivity of the response R = T (1)

T
= 294.579029[K].Consider next a perturbation δω in = (0.0063) ω 0 in , for which the perturbed value of the inlet air humidity ratio becomes ω pert in = ω 0 in + δω in = 0.0138612.Re-computing the perturbed response by solving Equations (2)-(13) with the value of ω pert in yields the "perturbed response" value T (50) w,pert = 294.634438[K].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD 48 δω in = 639.98[K].Using this value together with the nominal values of the other quantities appearing in the expression on the right side of Equation (57) yields τ (49) a SFD = −75.64[K].This result compares well with the value τ (49) a

aw
, i = 1, . . ., I − 1; hence, the forgoing verification of the computational accuracy of τ (49) a also provides an indirect verification that the functions τ (i) a , i = 1, . . ., I − 1 were also computed accurately.(c)Verification of the adjoint functions τ obtained by solving the adjoint sensitivity system given in Equation (28) are as follows: τ (1) w = −5.730× 10 −7 [K/(J/s) ] and µ (1) w = 1.3996 [K/ (kg/s) ], respectively, as indicated in Figure 9. Now select a variation δT w,in Energies 2016, 9, 718 22 of 45 in the inlet water temperature T w,in , and note that Equation (40) yields the following expression for the sensitivity of the response R = T (50) w

4 ∂T
a and o(49) have been already verified as described in Section 2.2.2(a) and (b), it follows that the computed values of adjoint functions τ (1) a = −0.88745[K] o (1) = −1.38335× 10 −6 second equation involving solely these adjoint functions can be derived by selecting a perturbation, δm w,in , in the inlet water mass flow rate, m w,in , for which Equation (40) yields the following expression for the sensitivity of the response R = T (50) w to m w,in :

3 T
= 294.579[K].Consider now a perturbation δT w,in = (0.00033) T 0 w,in , for which the perturbed value of the inlet air temperature becomes T pert w,in = T 0 w,in + δT w,in = 298.89[K].Re-computing the perturbed response by solving Equations (2)-(13) with the value of T pert w,in yields the "perturbed response" value T (50) w,pert = 294.5895[K].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD a perturbation δm w,in = (0.001) m 0 w,in , for which the perturbed value of the inlet air temperature becomes m pert w,in = m 0 w,in + δm w,in = 44.065252[kg/s].Re-computing the perturbed response by solving Equations (2)-(13) with the value of m pert w,in yields the "perturbed response" value T (50)w,pert = 294.5804[K].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" now all of the numerical values of the known quantities in Equations (58) and (59) yields the following system of coupled equations for obtaining µ

( 1 )
nom = 86.11678%.Consider next a perturbation δT a,in = (0.001) T 0 a,in , for which the perturbed value of the inlet air temperature becomes T pert a,in = T 0 a,in + δT a,in = 299.40911[K].Re-computing the perturbed response by solving Equations (2)-(13) with the value of T pert a,in yields the "perturbed response" value RH (1) pert = 85.55717%.Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" this value together with the nominal values of the other quantities appearing in the expression on the right side of Equation (66) yields o (49) SFD

( 1 )=w
nom = 86.11678%.Consider next a perturbation δω in = (0.0063) ω 0 in , for which the perturbed value of the inlet air humidity ratio becomes ω pert in = ω 0 in + δω in = 0.0138612.Re-computing the perturbed response by solving Equations (2)-(13) with the value of ω pert in yields the "perturbed response" value RH (1) pert = 86.28967%.Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" 19.6252.Using this value together with the nominal values of the other quantities appearing in the expression on the right side of Equation (69) yields τ (49) a SFD = −67.034.This result compares well with the value τ (49) a = −67.047obtained by solving the adjoint sensitivity system given in Equation (28), cf.Figure 10.When solving this adjoint sensitivity system, the computation of τ (49) a depends on the previously computed adjoint functions τ (i) a , i = 1, . . ., I − 1; hence, the forgoing verification of the computational Energies 2016, 9, 718 25 of 45 accuracy of τ (49) a also provides an indirect verification that the functions τ (i) a , i = 1, . . ., I − 1 were also computed accurately.(c) Verification of the adjoint functions τ Note that the values of the adjoint functions τ obtained by solving the adjoint sensitivity system given in Equation (28) are as follows: τ (1) w = −1.168× 10 −8 (J/s) −1 and µ (1)

4 ∂T
the adjoint functions τ (49) a and o (49) have been already verified as described in Section 2.2.3(a) and (b), it follows that the computed values of adjoint functions τ (1) a

( 1 ) 3 RHRH
nom = 86.11678%.Consider now a perturbation δT w,in = (0.00033) T 0 w,in , for which the perturbed value of the inlet air temperature becomes T pert w,in = T 0 w,in + δT w,in = 298.89[K].Re-computing the perturbed response by solving Equations (2)-(13) with the value of T pert w,in yields the "perturbed response" value RH (1) pert = 86.13838%.Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD a perturbation δm w,in = (0.001) m 0 w,in , for which the perturbed value of the inlet air temperature becomes m pert w,in = m 0 w,in + δm w,in = 44.065252[kg/s].Re-computing the perturbed response by solving Equations (2)-(13) with the value of m pert w,in yields the "perturbed response" value RH (1) pert = 86.11725%.Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD 45 nom δm w,in = 0.000107.Inserting now all of the numerical values of the known quantities in Equations (70) and (71) yields the following system of coupled equations for obtaining µ (1) w and τ (1) w : 0.00216 = − (0.0161) µ (1) w + (223597) τ

=
= 43.598097[kg/s].Consider next a perturbation δω in = (0.0063) ω 0 in , for which the perturbed value of the inlet air humidity ratio becomes ω pert in = ω 0 in + δω in = 0.0138612.Re-computing the perturbed response by solving Equations (2)-(13) with the value of ω pert in yields the "perturbed response" value m (50) w,pert = 43.603424[kg/s].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" 61.53 [kg/s].Using this value together with the nominal values of the other quantities appearing in the expression on the right side of Equation (79) yields τ (49) a SFD = −102.39[kg/s].This result compares well with the value τ (49) a = −102.42[kg/s] obtained by solving the adjoint sensitivity system given in Equation (28), cf. Figure 11.When solving this adjoint sensitivity system, the computation of τ (49) a depends on the previously computed adjoint functions τ (i) a , i = 1, . . ., I − 1; hence, the forgoing Energies 2016, 9, 718 28 of 45 verification of the computational accuracy of τ (49) a also provides an indirect verification that the functions τ (i) a , i = 1, . . ., I − 1 were also computed accurately.(c) Verification of the adjoint functions τ Note that the values of the adjoint functions τ second equation involving solely these adjoint functions can be derived by selecting a perturbation, δm w,in , in the inlet water mass flow rate, m w,in , for which Equation (40) yields the following expression for the sensitivity of the response R = m (50) w to m w,in :

3 m
w,nom = 43.598097[kg/s].Consider now a perturbation δT w,in = (0.00033) T 0 w,in , for which the perturbed value of the inlet air temperature becomes T pert w,in = T 0 w,in + δT w,in = 298.89[K].Re-computing the perturbed response by solving Equations (2)-(13) with the value of T pert w,in yields the "perturbed response" value m (50) w,pert = 43.591565[kg/s].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD 06531 [(kg/s) /K].Next, consider a perturbation δm w,in = (0.001) m 0 w,in , for which the perturbed value of the inlet air temperature becomes m pert w,in = m 0 w,in + δm w,in = 44.065252[kg/s].Re-computing the perturbed response by solving Equations (2)-(13) with the value of m pert w,in yields the "perturbed response" value m (50) w,pert = 43.641962[kg/s].Using now the nominal and perturbed response values together with the parameter perturbation in the finite-difference expression given in Equation (41) yields the corresponding "finite-difference-computed sensitivity" S FD 45 m (50) w,pert −m (50) w,nom δm w,in = 0.99646.Inserting now all of the numerical values of the known quantities in Equations (80) and (81) yields the following system of coupled equations for obtaining µ (1) w and τ (1) w : −0.06531 = − (0.0161) µ (1) w + (223597) τ

Figure A1 .
Figure A1.Humidity sensor accuracy plot (adopted from the specification of HOBO Pro v2).

Figure A2 .
Figure A2.Histogram plot of the measured air outlet relative humidity, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A1 .
Figure A1.Humidity sensor accuracy plot (adopted from the specification of HOBO Pro v2).

Figure A2 .
Figure A2.Histogram plot of the measured air outlet relative humidity, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A2 .
Figure A2.Histogram plot of the measured air outlet relative humidity, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).
below in Figures A3 through A6, and Tables A2 through A5, respectively.

Figure A3 .
Figure A3.Histogram plot of the air outlet temperature measured using "Tidbit" sensors, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A3 .
Figure A3.Histogram plot of the air outlet temperature measured using "Tidbit" sensors, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A4 .
Figure A4.Histogram plot of the air outlet temperature measured using "Hobo" sensors, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A4 .
Figure A4.Histogram plot of the air outlet temperature measured using "Hobo" sensors, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A5 .
Figure A5.Histogram plot of water outlet temperature measurements, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A6 .
Figure A6.Histogram plot of air outlet temperatures averaged from Figures A3 and A4.

Figure A5 .
Figure A5.Histogram plot of water outlet temperature measurements, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure A6 .
Figure A6.Histogram plot of air outlet temperatures averaged from Figures A3 and A4.

Figure B1 .
Figure B1.Histogram plot of dry-bulb air temperature data collected by SRNL from F-Area cooling towers (unsaturated conditions).

Figure B1 .
Figure B1.Histogram plot of dry-bulb air temperature data collected by SRNL from F-Area cooling towers (unsaturated conditions).

Note 2 :
Temperature and pressure values are initially input in units of [C] and [mb], respectively, but are internally converted to [K] and [Pa] for computational purposes.Note 3: Inlet air humidity ratio is defined as follows:

Table 1 .
Units of the adjoint functions for different responses.
The corresponding nominal value T a,nom = 297.4637[K].Consider next a perturbation δT a,in = (0.001) T 0 a,in , for which the perturbed value of the inlet air temperature becomes T pert a,in = T 0 a,in + δT a,in = 299.40911[K].Re-computing the perturbed response by solving Equations ( 49) 49)

Table A1 .
Statistics of the air outlet relative humidity distribution [%].

Table A1 .
Statistics of the air outlet relative humidity distribution [%].

Table A2 .
Statistics of the air outlet temperature distribution [K], measured using "Tidbit" sensors.

Table A2 .
Statistics of the air outlet temperature distribution [K], measured using "Tidbit" sensors.Histogram plot of the air outlet temperature measured using "Tidbit" sensors, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Table A2 .
Statistics of the air outlet temperature distribution [K], measured using "Tidbit" sensors.

Table A4 .
Water outlet temperature distribution statistics [K].

Table A5 .
Statistics of the averaged air outlet temperature distribution [K].Histogram plot of water outlet temperature measurements, within the 7688 data sets collected by SRNL from F-Area cooling towers (unsaturated conditions).

Table A4 .
Water outlet temperature distribution statistics [K].

Table A5 .
Statistics of the averaged air outlet temperature distribution [K].Histogram plot of air outlet temperatures averaged from Figures A3 and A4.

Table A4 .
Water outlet temperature distribution statistics [K].

Table A5 .
Statistics of the averaged air outlet temperature distribution [K].
[6]ng the relation in the above Equation (A3) in conjunction with the result presented in Equation (A2) will lead to an increase of the variances on the diagonal of the respective "measured covariance matrix", which will be denoted as Cov T meas a,out , T meas w,out , RH meas out .The final result obtained is:As indicated in the accompanying PART II[6], the predictive modeling formalism (which includes uncertainty quantification, data assimilation, and model calibration) also requires as "input" the covariance matrix between the measured parameters and responses.All of the parameters and responses are uncorrelated, except possibly for the measured responses considered in this Appendix and the measured parameters considered in Appendix B. The following "parameter-response" covariance matrix, denoted as Cov T meas a,out , T meas w,out , RH meas , α 1 , . . ., α 52 , is obtained for the respective parameters (namely: dry-bulb air temperature, T db ; dew-point air temperature, T dp , inlet water temperature, T w,in , and atmospheric pressure, P atm ) and responses (i.e., average outlet air temperature, outlet water temperature, and outlet air relative humidity): RH meas , α 1 , . . ., α 52 =

Table B1 .
Parameters for SRNL F-area Cooling Towers.

Table B4 .
Statistics of the dry-bulb temperature (set to air inlet temperature) distribution [K].

Table B4 .
Statistics of the dry-bulb temperature (set to air inlet temperature) distribution [K].

Table B7 .
Statistics of the atmospheric pressure distribution [Pa].The 4-by-4 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 , and atmospheric air pressure P atm .Cov T db ; T dp ; T w,in ; P atm =