Considerations for parameter optimization and sensitivity in climate models

Multiobjective Approach

A standard optimization approach is to lump variables ϕ_k into a single weighted cost function, such as [1]where f_k is an objective function associated with a particular aspect of the climate system, such as the square error (β = 1) or rms error () in Eq. 1, of a particular climate variable ϕ_k relative to the observed field . Typically, ϕ_k will be time averages (although they could include regression coefficients or other statistics) that depend on space and season, and 〈〉_k is a spatial average (or spatial and seasonal average) over a domain of interest. The domain might be global, but there may also be user-dependent preferences for particular averaging regions, hence the inclusion of subscript k on the spatial average. Optimization of a weighted cost function yields a single optimum (at considerable cost if the procedure involves functional evaluations using the climate model) that is a compromise among these variables. However, the preference vector of weights w_k is arbitrary, because it depends on how a particular user values accuracy in certain regions or in certain variables.

Adopting elements from multiobjective optimization (19–23) provides a natural fit for climate modeling in several respects. Treating f_k as multiple objective functions, optima for each of which are to be considered simultaneously, can identify tradeoffs among variables, the best case for each variable, and dominance relations among potential optima. In addition, much information relevant to decisions on parameter update are qualitative, e.g., indications from a field campaign that a parameterization is less trustworthy in a certain range. Such higher-level information can be subjectively combined in choosing among optima.

Metamodel for Climate Field Parameter Dependence

The climate modeling problem has close parallels to high-dimensional design problems with computationally expensive black-box functions, optimization strategies for which are reviewed in ref. 24. Results below on smoothness in certain measures lead us to adopt a strategy of fitting the parameter dependence, termed metamodeling (25, 26), with a low-order polynomial. An important distinction is that the metamodel here is created for the spatial and seasonal vector field of the climate model output, not merely for the objective functions. Advantages of this approach include analytic insights that aid visualization of the potentially strong regional dependence of the sensitivity and optimization. To the extent that the metamodel is successful, optimization for multiple objective functions is relatively inexpensive. This approach addresses a nagging issue—that a parameter optimization choice based on one set of objective functions is always vulnerable to reevaluation by additional criteria. Here, alternative objective functions are easily added from the same set of climate model runs.

Climate Model Setup

We examine these issues in the International Centre for Theoretical Physics climate model (27, 28), which has been used for a number of climate problems (29, 30). Here a triangular truncation at total wavenumber 30 (T30) version is used, roughly equivalent to a 3.75° × 3.75° spatial resolution. We examine both a set of runs with the atmospheric general circulation model (AGCM) forced by observed sea surface temperature (SST) and the AGCM coupled to a mixed-layer (ML) ocean, using a simple estimate of ocean heat flux transport divergence known as a Q-flux (31). The Q-flux is set such that observed SST is approximated at the standard parameter values and does not vary with parameter. For the imposed-SST case, ensembles of ten 25-y runs (1978–2002), differing in their initial conditions, are carried out for each parameter value considered. For the AGCM-ML runs, 260-y runs are carried out at each parameter value, omitting a 10-y spin up and using ten 25-y averages for analogous statistics.

Parameter Dependence and Fits to Model Output Fields

Here we illustrate the parameter space dependence using four parameters that can have significant impact on the model solution: the subgrid scale wind gustiness that controls the minimum wind speed in the bulk formula for surface fluxes, a viscosity parameter measured as a damping time, a cloud albedo parameter, and the relative humidity parameter from the deep convective parameterization, RH_conv, which represents the observed inhibition of convection unless free troposphere relative humidity is sufficiently high.

AGCM Forced by Observed SST Case.

Fig. 1 shows results for slices along the four selected parameter directions for the case of seasonal precipitation. Contrary to our initial expectations, the parameter dependence proves rather smooth in large-scale measures. We underline the caveat that this approximate smoothness will not necessarily hold for regional-scale measures or more exotic statistics. The fact that it tends to hold for measures such as global rms error or spatial correlation to observations for climatological fields can be of great utility in the estimation of sensitivity and optimization. This property motivates use of low-order polynomial fits below. Estimation error is small for the global rms error for the 10-member ensemble used here, but variation among 25-y averages for individual ensemble members is illustrated in Fig. 1D. A single 25-y run would suffice for a first estimate of sensitivity and nonlinearity but would have limitations for statements regarding smoothness.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Root-mean-square error of the ensemble mean AGCM precipitation in June–August (JJA), relative to National Centers for Environmental Prediction reanalysis. The AGCM values are compared to the rms error reconstructed from the quadratic metamodel (Eq. 2) using endpoints for each parameter and to its linear counterpart. Note that the linear metamodel gives quadratic terms (with positive curvature) in the rms error. Parameters are given on the abscissa (nondimensional if no units given). The vertical size of the symbols gives the two standard error estimation range for the ensemble mean. D also shows the rms error for 10 individual ensemble members to illustrate the estimation error that would be associated with a single 25-y run.

Mixed-Layer Case and Sensitivity Under Change of CO₂.

These properties tend to be inherited in the mixed-layer run, as seen in Fig. 2A for the case of RH_conv. The smoothness of the solution for large-scale measures and the reasonable fit of the quadratic metamodel are reproduced, as is the negative curvature of the rms error.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

(A) Root-mean-square error of ensemble mean June–August (JJA) precipitation as in Fig. 1D, but for the mixed-layer coupled case, with rms error reconstructed from linear and quadratic metamodels using endpoints for the relative humidity parameter. (B) Change JJA precipitation for doubled CO₂ minus preindustrial CO₂ runs, shown as rms difference between the spatial field for each value of RH_conv and the case for RH_conv = 0.9 (the value with minimum rms error in the climatology). The vertical size of the symbols gives the two standard error estimation range.

For the case of changes between doubled CO₂ and preindustrial AGCM-ML runs, we cannot of course show rms error because no observations are available. To provide a rough sense of the parameter dependence, Fig. 2B shows rms differences of the model spatial field at each parameter point relative to the spatial field at a reference parameter value, here using the current climate optimum along the RH_conv axis. The metamodel is the quadratic fit as a function of parameter as in Eq. 2 but to the spatial field Δϕ = ϕ_2× - ϕ_pre, where ϕ_2× and ϕ_pre are 250-y averages from the doubled CO₂ and preindustrial cases, respectively, and the two-standard error range is estimated from ten 25-y intervals. Again the rms difference evolves smoothly, exhibiting substantial nonlinearity that is well fit by the quadratic metamodel.

Fig. 3 shows the linear and quadratic contributions of the metamodel (Eq. 2) in the RH_conv direction (fit from endpoints only) for the AGCM-ML preindustrial CO₂ case. The spatial pattern is that of the coefficients a_i, b_ii, and the magnitude is given at parameter endpoint μ_i max (i.e., at RH_conv = 0.9 relative to standard value). The sum of the two fields gives the precipitation change from the standard case at μ_i max. Both the linear sensitivity and the quadratic nonlinearity are sufficiently strong to produce substantial modifications of the precipitation pattern across the parameter range, especially in the tropics where convective precipitation is the dominant contribution. These spatial patterns provide a view of where changes are occurring at the regional scale and of the regional importance of nonlinearity. Analytic results that further exploit these properties are discussed below.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Ensemble-mean measures of June–August precipitation parameter sensitivity in RH_conv for the AGCM-ML, preindustrial CO₂ case (A) linear contribution a_i multiplied by μ_i max to give units of millimeter per day for a value that would occur at the positive end of the feasible range; (B) quadratic contribution b_ii, similarly given as in millimeters per day. The two contributions add to the total difference between μ_i max and the standard case, where subscript i denotes RH_conv.

Fig. 4 shows corresponding linear and quadratic contributions to the parameter sensitivity of the precipitation change under CO₂ doubling. The strength of the sensitivity is comparable to the magnitudes of seasonal precipitation differences in the tropics seen across a multimodel ensemble in the Coupled Model Intercomparison Project 3 archive (e.g., refs. 14 and 15). The substantial differences in the pattern of change predicted under global warming are here seen in a single model, occurring relatively smoothly under parameter variations. A key feature is the importance of the nonlinear term, which is overall of comparable importance to the linear term and can exceed it in particular regions.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Ensemble-mean measures of parameter sensitivity in RH_conv as in Fig. 3 but for the change in June–August precipitation for doubled CO₂ minus preindustrial case in the AGCM-ML (shown as the value that would occur at μ_i max in millimeters per day) (A) linear contribution a_iμ_i max (millimeters per day) ; (B) quadratic contribution .

Analytic Approximations and Regional Contributions

Climate model output comes in the form of high-dimensional vector fields, but these fields are easily interpreted by a climate modeler as geographic information. It can thus be highly advantageous to have analytic solutions that present aspects of the optimization problem in this form. A common class of question is the difference between an optimization based on a regional spatial average, as opposed to the global average. Analytic approximations to the optimization based on the metamodel permit examination of such differences. Recall that the quadratic fit (Eq. 2) to the AGCM field ϕ gives rms or square-error objective functions f (here omitting subscript k) that have fourth-order terms in μ. Using square error (which has the same extrema as rms), expanding in μ about a reference value, here the standard case, retaining second-order in μ in f, and differentiating yields [3][4][5]where g is the gradient in μ and A the Hessian matrix from the curvature, both evaluated at the standard case, and ϕ_err = ϕ_std - ϕ_obs is the error of the standard case compared to observations. A standard quadratic optimization problem (32) min _μ f, with constraints μ_i min ≤ μ_i ≤ μ_i max would give Aμ = -g for interior solutions that fall within the feasible μ-range, if additional conditions on the curvature are satisfied. For f associated with seasonal precipitation, Figs. 1D and 2A showed that the problem is nonconvex in the RH_conv direction, corresponding to A not being positive semidefinite, when 〈〉 is a global average. When a boundary solution occurs, interior solutions for other variables can be found with the same A_ij and a slight adaptation of g (see SI Appendix). Thus considerable insight into the optimization problem can be gleaned from the approximation (Eq. 3) (noting that this approximation is not used for numerical multiobjective optimization below, which carries all terms).

The off-diagonal terms in A act to change the orientation of the basin in which the optimization occurs. Off-diagonal contributions come from both the linear sensitivity via 〈a_ia_j〉 and from 〈b_ijϕ_err〉. These off-diagonal terms are often small compared to the diagonal (i.e., parameter pairs often do not interact strongly), so focusing on the diagonally dominated problem ∂_μif = g_i + A_iiμ_i yields a sense of basic behavior. The form of g_i is a spatial projection of the sensitivity a_i with the error of the standard case, ϕ_err. Thus regional sensitivity is only effective at yielding a reduction in error to the extent that the spatial pattern of a_i matches the error pattern. For an interior solution, if b_ii is small, the μ_i optimum would simply be . This expression implies that the distance moved in μ_i is a tradeoff between the part of the sensitivity that projects on the error, and thus can reduce it, and the part of the sensitivity that is orthogonal to the error (under the inner product corresponding to the average used in the objective function). When using a global average measure, this tradeoff has substantial implications for regional errors. Suppose a_i has a large value in a small region that does not project on the error. The optimization will yield a solution that reduces global error, at the cost of introducing a large error in that region—a likely explanation for this common modeling experience.

The linear contribution has positive curvature, consistent with the possibility of having an interior minimum. However the quadratic term 2〈b_iiϕ_err〉 can reverse this curvature if sufficiently negative, as in Fig. 2a. Again only the projection on ϕ_err matters. Recall that the quantity acting like an inner product, 〈〉, is typically a spatial average, but can be user dependent, e.g., a regional measure instead of a global average. For a quick view of whether objective functions defined for spatial averages on different regions will tend to have properties similar to the global average or strongly differing, one can examine maps of the spatial structure of these contributions such as those shown in Fig. 5. The quantity (a_iϕ_err) normalized by the diagonal contribution (Fig. 5A) is proportional to contributions to the objective function gradient g_i in RH_conv (at μ = 0). The fact that most points on the map are negative indicates that optimization using various regional averages in the objective function will most commonly yield optima that move in the same, positive direction in RH_conv as occurs when the global average was used. The quantity (Fig. 5B) is the contribution of nonlinearity in the model field to the diagonal term. The global average of this term is negative, with magnitude larger than one, so it reverses the curvature. The fact that this term has large negative contributions almost everywhere indicates that this influence of nonlinear terms in the RH_conv direction tends to hold consistently for regional averages. For other parameters, regional properties need not be so consistent with those of the global optimization, for instance containing regions of opposite sign. Examination of the spatial patterns in the analytic solution can provide a low-cost way of identifying such regional issues and can be easier to interpret than a large number of regional objective functions.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Contributions to the approximate analytic solution for the ensemble-mean June–August precipitation preindustrial CO₂ case of Fig. 3 (A) The spatial pattern (a_iϕ_err) (nondimensionalized by ), the area-weighted spatial average of which yields g_i in Eq. 4. (B) The spatial pattern 2(b_iiϕ_err) (nondimensionalized by ) whose area-weighted spatial average appears in Eq. 5.

Multiobjective Optimization

The success of the metamodel over the feasible parameter range, permitting the construction of various objective functions, facilitates multiobjective optimization, with separate objective functions for each major climate variable. Repeated optimizations using standard packages for constrained optimization (here the Nonlinear Interior point Trust Region Optimization package, ref. 32) prove quite practical, although advanced multiobjective techniques (20, 21) could be applied.

Fig. 6 quantifies the contradiction among objective functions, i.e., the tradeoffs in optimizing for one variable versus another, as the location of the optima in parameter space for different climate variables, illustrating with the June–August season. The spread of the optima in parameter space is substantial in two of the parameter directions shown. In the RH_conv direction, the optima all occur at or near the same boundary of the feasible range. This agreement favors an update of the parameter value, but the fact that it occurs at the end of the feasible range suggests careful scrutiny of the physical processes involved in the parameterization in this range. The sensitivity to high values of water vapor in the free troposphere required for convective onset is closely related to other work pointing to the importance of a moist environment in the lower troposphere for development of deep convective plumes (33–38). In the other parameters in Fig. 6, boundary optima are also common. Boundary optima are likewise found in other measures we have examined. In selecting parameter updates, these boundary optima would tend to be distrusted to a degree that depends on the climate modeler’s confidence in the parameterization at the end of the range.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Global optima calculated using an order-N² (large sphere) or order-N (small sphere) procedure for global-average rms error objective functions for several model variables. The vertical pressure level of the variable, e.g., T200, is given in hectopascal. The projection of each 3D location on the lower plane is shown as a visual aid. Precip, precipitation; T, temperature; U and V, wind components; Ω, vertical velocity; MSLP, mean sea level pressure; LSTA, land surface temperature.

The optima in Fig. 6 are a nondominated set, i.e., comparing the objective functions f_k between any two optima, one cannot improve one f_k without making others worse. Additional information about the tradeoffs among optima are provided by constructing from the proxy the set of values of f_k when each is at its own optimum (the ideal objective vector or utopia point) (19, 22). These values provide a useful normalization for the objective function f_k. Examining normalized f_k at each of the other optima indicates that some tradeoffs are only a few percent worse than the optimum value and highlights certain others. For instance, a significant tradeoff in this model involves land surface temperature, which is roughly 30% off its optimum value at the optima for precipitation, vertical velocity, and low-level wind.

Also shown in Fig. 6 is a comparison of a metamodel based on order-N climate model evaluations versus a procedure that uses order N². The cheaper order-N procedure, using the approximation that parameters do not interact strongly in the quadratic terms, works well enough for most variables that examination of the properties would be useful prior to investing in the additional runs, and could guide the choice of these.