- Email: [email protected]
Optimization of a middle atmosphere diagnostic scheme
Solar-Terresrriol Physics. Vol. 59, No. 9. pp 98.&992. 1997 @ 1997 Elsevier Science Ltd. All rights reserved. Rimed in Great Britain 1364.6826/97 $17.00+0.04 S1364-6826(96)00034-7
Journal
Pergamon PII:
of Amosphrric
and
Optimization of a middle atmosphere diagnostic scheme Rashid A. Akmaev * High Altitude Observatory, National Center for Atmospheric Research, Boulder, Colorado, U.S.A. (Received in final form 2 January 1996; accepted 20 January 1996)
Abstract-A new assimilative diagnostic scheme based on the use of a spectral model was recently tested on the CIRA-86 empirical model. It reproduced the observed climatology with an annual global rms temperature deviation of 3.2 K in the 15-l IO km layer. The most important new component of the scheme is that the zonal forcing necessary to maintain the observed climatology is diagnosed from empirical data and subsequently substituted into the simulation model at the prognostic stage of the calculation in an annual cycle mode. The simulation results are then quantitatively compared with the empirical model, and the above mentioned rms temperature deviation provides an objective measure of the ‘distance’ between the two climatologies. This quantitative criterion makes it possible to apply standard optimization procedures to the whole diagnostic scheme and/or the model itself. The estimates of the zonal drag have been improved in this study by introducing a nudging (Newtonian-cooling) term into the thermodynamic equation at the diagnostic stage. A proper optimal adjustment of the strength of this term makes it possible to further reduce the rms temperature deviation of simulations down to approximately 2.7 K. These results suggest that direct optimization can successfully be applied to atmospheric model parameter identification problems of moderate dimensionality. @ 1997 Elsevier Science Ltd
1.
ing data assimilation technique (e.g. Anthes, 1974; Ghil and Malanotte-Rizzoli, 1991) is iteratively integrated towards a steady state. The final result of these iterations is the total zonal mean zonal drag obtained in such a way as to keep the model climatology in the vicinity of empirical data. The close connection of the diagnostic scheme to the simulation model facilitates an objective assessment of the retrieved zonal drag. At the second stage, the zonal forcing is substituted back into the unmodified model, and the model climatology is then quantitatively compared with the initial empirical climatology. Akmaev (1994) tested this scheme on the CIRA-86 empirical model (Fleming et al., 1990). The simulations using the retrieved zonal forcing were shown to reproduce the observed climatology with an annual global rms temperature deviation of 3.2 K or about 2% in the 15-I IO km layer. The closeness of the simulations to the empirical model suggested for the first time that the diagnosed forcing yielded reliable quantitative estimates of the intensity of real physical processes governing the global-scale circulation in the middle atmosphere.
INTRODUCTION
Since the work of Murgatroyd and Singleton (1961), diagnostics of observed climatologies have been used to study the mechanisms maintaining the gross features of the global-scale circulation and tracer transport in the middle atmosphere. The zonal drag or forcing imposed by dissipating and breaking waves of all scales from global to mesoscale, is of particular importance and interest (e.g. Shine, 1989; Medvedev and Fomichev, 1994). Akmaev (1994) recently proposed a new assimilative middle atmospheric diagnostic scheme and compared it with more traditional diagnostic procedures (see also, Shine, 1989; Medvedev and Fomichev, 1994, and references therein). The scheme employed in this study and detailed by Akmaev (1994) is based on the use of a spectral model developed by Akmaev et al. (1992). At the diagnostic stage the model with the zonal momentum equation modified in the spirit of the nudg-
* Present address: Department of Aerospace Engineering Sciences, Campus Box 429, University of Colorado, Boulder, CO 80309-0429, U.S.A.
983
R. A. Akmaev
984
Both the spectral model and the diagnostic scheme as a whole have some adjustable parameters. The diagnostic scheme is also quite flexible and it may be subjected to further modification via the introduction of additional parameters. This flexibility and the existence of the above mentioned quantitative criterion make it possible to apply thorough optimization procedures to the simulation model, the diagnostic module, or to the whole diagnostic-simulation complex. The first results of such optimization experiments are reported in this paper. It has been proposed, for example, that estimates of the zonal drag could be improved by introducing a nudging, or Newtonian-cooling term into the thermodynamic equation at the diagnostic stage (Akmaev, 1994). As is demonstrated below, a proper adjustment of the strength of this Newtonian-cooling term makes it possible to further reduce the rms temperature deviation down to approximately 2.7 K. These results strongly support the objective quantitative approach to middle atmospheric modeling and diagnostics. The results also clearly demonstrate that direct optimization can successfully be applied to atmospheric model parameter identification problems of moderate dimensionality as opposed, for example, to the more demanding adjoint model technique (e.g. Daley, 1991; Ghil and Malanotte-Rizzoli, 1991) (see also the bibliography collected by Courtier et al., 1993), the method of choice for large-dimension optimization problems. Section 2 briefly describes the diagnosticsimulation scheme and its modifications introduced in this study. Section 3 presents direct optimization routines suitable for the problem. In Section 4 the results of a series of optimization experiments are presented and discussed.
2. DIAGNOSTICSIMUL.ATlON
SCHEME
The diagnostic-simulation scheme employed in this study is based on the middle atmospheric spectral model developed and described in more detail by Akmaev et a/. (1992) and Akmaev (1994). A global three-dimensional spectral model with a triangular truncation T13 corresponding to a moderate latitude-longitude resolution of about 8.6” x 9” is used. The model covers the altitude range between 100 mb (-15 km) and
17 pressure scale heights (- 120 km) with 20 layers of equal thickness in log-pressure coordinates with a vertical resolution of approximately 5 km. Zonally-symmetric diagnostics are performed as predetermined by the zonal symmetry of empirical data, so the model has been run in a two-dimensional mode with the longitudinal dependency of model variables ‘turned off.’ As was suggested by Akmaev (1994) the radiation transfer module has been updated by replacing the parameterization of the ozone ultraviolet absorption of Strobe1 (1978) and Apruzese et al. (1982) with that of Shine and Rickaby (1989), providing better results in the stratosphere due to enhanced heating primarily in the Huggins bands. The diagnostic scheme and the approximations involved are presented and extensively discussed by Akmaev (1994). At the diagnostic stage the model is initialized with the CIRA-86 empirical temperatures and zonal winds (Fleming et af., 1990) for each month. It is then integrated towards a quasi-equilibrium state and the total zonal drag, aZ, necessary to maintain the empirical climatology is determined from the following iterative procedure. The zonal momentum equation can be written in the form:
au = R, -t-a=, at
-
where u is the zonal mean zonal wind and R, represents all other terms in the equation that can be explicitly computed in a zonally symmetric model. At every nth time step aZ is estimated as
z = 0.5 (ZT-” -R:). The averaging of two consecutive iterations in equation (2) is used to prevent the model from spurious oscillations. The initial estimate 2 may be set to zero, or taken from previous diagnostic results. This procedure generally yields reasonable estimates of aZ if the integration continues for IO-15 model days from the initial zero conditions and if the model remains in close vicinity to the empirical climatology during this integration. Akmaev (1994) suggested that these estimates could potentially be improved if an additional nudging term is inserted into the thermodynamic equation to more strongly constrain the model climatology towards observations at the diagnostic stage. Specifically,
Optimization
aT
,t=RT-a(T-T“),
of a middle
atmosphere
(3)
where RT is the right-hand side of the standard thermodynamic equation, 01 is the nudging (Newton-cooling) coefficient, and To is the ‘observed’ temperature. In this study, o( = o((z) is assumed to depend on the vertical coordinate only. A conceptually straightforward and computationally efficient nudging or relaxation technique that may be considered a degraded form of the Kalman filter (Daley, 1991; Ghil and MalanotteRizzoli, 1991) was first introduced for atmospheric data assimilation by Anthes (1974) and continues to be extensively investigated for use in meteorology and oceanography. It is usually applied to prognostic model equations, in a fashion similar to equation (3), to force the model to ‘account’ for available observations, or, in other words, to blend observations and simulations in an ‘optimal’ manner. Generally, o( is a matrix which may depend on the time differential between times of observations and model time steps, on the spatial distance between points of observations and model grid points, on the standard error of observations, etc. The optimality of combining experimental data with numerical models in this case greatly depends on the strength of nudging coefficients. The method of adjoint model variational data assimilation has been used recently in the parameter estimation regime to determine ‘optimal’ forcing (Zou et al., 1992; Stauffer and Bao, 1993). Zou et al. (1992) applied the technique to a three-dimensional spectral adiabatic general circulation model to find constant nudging coefficients appearing in three model equations. Stauffer and Bao (1993) used a onedimensional shallow-water model in their search for three nudging coefficients depending on the horizontal coordinate. It is clear that the stronger the relaxation terms in prognostic equations are, the closer the model will reproduce observations. Therefore, a cost function whose minimum is sought in variational data assimilation, and that usually represents a measure of the distance between simulations and observations, had to be augmented in these studies by a penalty term depending on the deviation of the nudging coefficients from some a priori estimates. This essentially transferred the uncertainty from the nudging coefficients onto these preliminary es-
diagnostic
985
scheme
timates and penalty coefficients. Stauffer and Bao (1993) also studied the sensitivity of the results to the penalty term, and found that, as one would expect, the ‘optimal’ nudging coefficients increase as the weight of the penalty term decreases. Unlike the above mentioned studies, we apply nudging only at the diagnostic step in this study. One can speculate that there must be some finite optimal forcing which will produce best estimates of a,. Indeed, if o( = 0 in equation (3) the model may eventually deviate from the initial empirical climatology in the process of iterative integration towards a steady state. If o( is too large, however, the model climatology will almost entirely be determined by the fictitious forcing in equation (3) and the corresponding estimates of a, may not be sufficient to adequately force the model towards observations at the simulation stage. The purpose of this work is to find the optimal finite nudging and obtain the best estimates of us. Another feature of the present scheme that makes it different from those employed by Zou et al. (1992) and Stauffer and Bao (1993) is that a direct optimization procedure (see Section 3) is used for our comparatively low-dimension problem as opposed to the costly adjoint model technique. Let us now summarize the whole procedure used in this study which differs slightly from that employed by Akmaev (1994). For a given o( = a(z) the model with the modified zonal momentum and thermodynamic equations [equations (1) and (3) respectively] is initialized by the zonal mean temperature and zonal wind distributions from CIRA-86 for each month, and z = 0. The model is then integrated for n = 300 steps, or 12.5 model days, with the zonal forcing being calculated according to equation (2). This integration time interval has been found sufficient to reach a quasi-equilibrium state. The retrieved monthly mean zonal forcing and the vertical diffusion coefficient determined as in Akmaev (1994) are then used in an annual run of the unmodified model, and the following objective or cost function is computed:
F=
2rlog:p,
,p*) ;,
rrn jP,d(logp) P2
n/2 x
( T,,, -
T,”)‘cos 4p d+.
(4)
986
R. A. Akmaev
Here F = u: is the annual altitude global mean temperature deviation squared, pi is set at the model lower boundary to 100 mb, and p2 = 1.8 x low4 mb is the pressure level two model layers below the upper boundary at an altitude of about 110 km, set to filter out the effect of the upper boundary condition (Akmaev, 1994). The length of month m is TV, the length of year is represented by T = C,T,, #J is latitude, T,,,(c$,p) is the monthly averaged simulated temperature, and T,O(@,p) is the ‘observed’ temperature. In equation (4) we only compare temperatures, as the empirical tables contain geostrophic zonal winds simply calculated from temperature distributions. or values obtained with different nudging are then compared to find, as described in the following sections, the optimal nudging that delivers the minimal rms temperature deviation. The ‘observed’ temperatures T,” are calculated from the CIRA-86 tables by interpolation and extrapolation onto the model horizontal Gaussian and vertical grids. In data assimilation schemes, it is usually model results that are projected onto the points of observations in space and/or time. However, in this scheme we use the inverse transformation for the following reasons. As mentioned above, T,” and the corresponding empirical zonal winds are used to initialize the model at the diagnostic stage, and thus need to be transformed onto the model spatial grid anyway. To compute the integrals in equation (4) some kind of interpolation has to be used too. The empirical data are presented on a regular space mesh and have, incidentally, been obtained by interpolation from original experimental results. Essentially, the nudging terms represent a mathematical device, a ‘trick,’ the sole purpose of which is to force a numerical model to produce a desirable result, and it is perhaps pointless to seek any physical meaning in the strength of the nudging coefficients or in their spatial or temporal distribution (see also Section 4). Again, we stress a very significant difference between the traditional applications of the nudging technique referred to earlier in this section and that employed in this work. We apply nudging only at the diagnostic stage. The forcing obtained from the diagnostics is then used in direct simulations performed with the unmodified model. This is done without any fictitious nudging, unlike simulations done in the
traditional approach. As is shown in Section 4, the retrieved zonal drag makes it possible for the model to very closely reproduce the empirical climatology. One can therefore argue that the forcing provides a representative quantitative estimate of the strength of real physical processes taking place in the atmospheric layer under consideration, to the extent, of course, that the empirical chmatology itself can be considered representative.
3. OPTIMIZATION
PROCEDURE
The purpose of optimization in this study is to minimize the objective function F(N) from equation (4) with respect to the nudging coefficients, 01 = a(r), which are assumed to depend on the vertical coordinate only. A series of numerical experiments have been performed with different characters of this height dependence (see Section 4). In a discrete form o( = ((~1,. . , a~)‘, where the prime denotes a transpose, may be considered a vector of a variable dimension L ranging from 0 (no nudging) through to 20, the total number of model layers. Newton-type methods are perhaps most effective for smooth unconstrained optimization problems of moderate dimension (Fletcher, 1980; Gill et al., 1981; Dennis and Schnabel, 1983). QuasiNewton low-memory procedures can even be successfully used in large-scale problems (Navon and Legler, 1987). Basically, a Newton-type method uses the available information on the objective function F(x), its first derivatives, or gradient g(x) =
, (g,...,$- >’ L
1
and its second derivatives, or the Hessian matrix, H(x) = {H;j}, where Hij
= -.
a2F
axi axj
’
i,j=l
, . , L.
(6)
A quadratic model, Q(x), of the objective function is constructed at every kth consecutive search point xtk): Qtk’(x) = Fck’ + (Ax)‘~‘~’ + 0.5(A~)‘H(~)Ax,(7) where Fck) = F(xck)) Ax = X-X(~) , g’k’ = g(x’k’), Hck) = H(x”)), and ;he right-hand side of equation (7) is simply a truncated Taylor series for F.
Optimization of a middle atmosphere diagnostic scheme
A search vector, $4 = s(X’k’~ = _ (H(k))-’ g’k’
(8)
calculated as a product of the inverse Hessian and negative gradient will then exactly point to the minimum ,?ck) of the quadratic approximation if the minimum exists: z.(k)
=
$k)
+ S(k).
(9)
Apart from some fine details, for a general nonlinear function, a linear search is conducted along the direction of search vector ~(~1.The whole process is repeated from the next search point .x(~+” which is determined as a minimum along the search direction:
F(X’k’“)
= t$ F(X’k’ + BP). >
(IO)
If the objective function is available analytically so that its gradient and Hessian matrix can be calculated exactly, perhaps the only computational problem is to invert the Hessian in equation (8). The exact Newton methods are characterized by a high, quadratic, convergence rate which essentially reflects the fact that the closer a particular search point x(k) is to the minimum sought x*, the better the quadratic model Qck)(x) approximates F(x), and the more precisely its minimum ~2~) represents x*. If a search point is far from the minimum, the quadratic model in equation (7) may not be adequate, and a Newton-type procedure can diverge, or the search sck), may not point in the direction of substantial function decrease. Therefore Newton steps have to be combined with some global-search moves to ensure that the search does not terminate, for example, in a local minimum instead of the global one within attainable prespecified limits. Also, some restart procedures can be used to establish a new search direction, if the ‘quadratic’ direction does not provide a sufficient function decrease. Optimization methods utilizing information only on the objective function itself (direct search methods) and/or its first derivatives (steepest descent methods) have been demonstrated to be inferior in their convergence rate. They may often fail entirely for objective functions of a comparatively complex structure.
987
Unfortunately, in many practical applications, as is the case in this study, every objective function value is a result of a complex numerical procedure and its second, and often even first, derivatives are not readily available. A group of so-called quasi-Newton optimization algorithms, closely related to the conjugate direction methods (Hestenes, 1980; Navon and Legler, 1987) exploits the idea of building up, at the expense of a slower convergence rate, the information on the curvature of objective function based on its gradients at consecutive search points. The Hessian matrix is successively updated in the search process starting from some initial approximation. It is also possible, as in some more advanced algorithms, to directly update the Cholesky decomposition factors of the Hessian, and to save a large amount of computation on its inversion. The adjoint model technique essentially provides an effective means to calculate first derivatives of various objective functions which result from some ‘direct’ modeling, particularly for large-scale problems. For problems of moderate dimension, the gradient can be obtained by finite differences. Therein, each gradient component requires at least one additional function call, if forward differences are used, and approximately the same order of additional calls, if central differences have to be used to improve precision. Some important issues have to be addressed in practical applications of optimization algorithms (Fletcher, 1980; Gill et al., 1981; Dennis and Schnabel, 1983). Among them are the proper scaling of the objective function and variables, or variable transformation if necessary. It is highly desirable to prescribe a feasible search region, maximal and minimal search steps, stopping criteria, and some other quantities adequately reflecting the nature of problem at hand. Ill-posed or inappropriately scaled problems may cause substantial loss of accuracy and prevent the algorithm from locating the global minimum, or even cause its total failure. It has been found after some limited experimentation with ‘easy-to-use’ and ‘comprehensive’ optimization routines available in the IMSL and NAG scientific computing libraries, that the latter provides routines which are slightly more efficient for our particular problem in that under equal conditions they generally require a smaller number of function calls to locate a minimum. Of
988
R. A. Akmaev
course, this could be just a matter of more adequate tuning of a particular routine. The results presented in Section 4 have been obtained using the NAG FORTRAN library (Gill et al., 1981) ‘comprehensive’ quasi-Newton routine E04UCE.
4. RESULTS AND DISCUSSION
A series of optimization experiments of different dimensionality L has been conducted. Akmaev (1994) presented the results of diagnostics and simulation, and calculated ar for a second model year. Every calculation of function F from equation (4), reported herein, involves a run of the diagnostic scheme for all twelve months and a run of the model for one year instead of two to save computer time. The differences in or for the first and second model years were found negligible, of the order of 0.003 K or about 0.1%. First we assess the effect of replacement of the parameterization of ozone ultraviolet heating by Strobe1 (1978) and Apruzese ef al. (1982) by that of Shine and Rickaby (1989) without nudging (L = 0) in the thermodynamic equation (see Table 1). This ‘optimization’ of the radiation modTable I, Root-mean-square
L
CT 6)
0 0
3.22 3.06 2.88 2.76 2.71
1 2 20
temperature deviation
Ozone heating Strobe1 (1978); Apruzese et al. (1982) Shine and Rickaby (1989) Shine and Rickaby (1989) Shine and Rickaby (1989) Shine and Rickaby (1989)
ule leads to improved simulations, primarily in the stratosphere, due to slightly enhanced heating throughout the region. The corresponding rms temperature deviation fir decreased from 3.22 K to 3.06 K. Of course, this lower ur is produced by the zonal forcing that slightly differs from that obtained in the previous study (Akmaev, 1994) and corresponds to a slightly different temperature distribution. These minute differences can be seen by comparison of Figs 5 and 7 in this paper with Figs 11 and 15 of Akmaev (1994). It has been found that the optimization routine described in the previous section is stable and efficiently locates the global minimum of the objective function, if F = F(x) is optimized with respect to x instead of a, where
XI = log,, a/,
I = l....,L.
(11)
This transformation also eliminates the only constraint on the variables o(/ 2 0, and the optimization problem becomes truly unconstrained. Optimization has been performed with the optimahty tolerance parameter (Gill et af., 1981) specified at c = 10m4. This approximately corresponds to five correct digits in the final value of F or or. A one-dimensional search (L = I), where o((z) = 01~ is assumed to be constant with altitude, typically converges in 2-3 steps from arbitrary initial conditions to at = 9.93 x 10m4s-i. The corresponding or = 2.88 K. For L = 2, we assume that LXexponentially depends on logp between the bottom and the top model layers, where 01 = (~1 and LX= o(2, respectively. The equivalent assumption is that x in equation (11) changes linearly with logp from xi to x2. First, the result of one-dimensional optimization was used as an initial condition, and the search converged in 3 steps to or = 2.76 K. Each major step in this case required 6-10 evaluations of the objective function to calculate its gradient, to perform a linear search along the search direction as in equation (lo), and global search moves if necessary. To ensure that the minimum found is global, and to explore the objective function behavior in a relatively broad area around the minimum, eight additional experiments have been performed with initial conditions separated from the solution by two orders of magnitude in (Y~and/or 012.The results of these experiments are presented in Fig. 1 where every dot represents a major search step. To get a feel for the accuracy of the search, we notice that, in addition to the above mentioned optimality tolerance, E, another stopping criterion, the so-called line search tolerance (Gill et al., 1981), has been specified at 6 = 0.01. This approximately corresponds to one tenth of the dot diameter in Fig. 1. Roughly speaking, a search is completed if two consecutive major steps fall within the line search tolerance, 6, if the relative difference between corresponding function values does not exceed the optimality tolerance, E, and if some other necessary and sufficient conditions for a minimum are satisfied. The eight initial points are located on the perimeter of the plot. In most cases the search required 4 steps, 3 steps were enough in one case, and 6 steps were required in the two worst cases. The one-dimensional case corresponds to a constrained ((xi = 012) mini-
Optimization
of a middle
atmosphere
6 (K*)
1d-s
,;-4
1d-s a1
10-j
(s-7
Fig. 1. Two-dimensional optimization (contour interval I K2). Dots represent major search steps. The eight starting points separated from the global minimum (big cross in the center) by 2 orders of magnitude in cxt and/or cx2 are located on the perimeter. The onedimensional optimal solution corresponding to a constrained minimum (the constraint at = CQ is represented by a dashed line) is presented in the lower right corner of the plot.
mum. The constraint is represented by a straight dashed line in the lower right comer of the plot, and the corresponding one-dimensional minimum and the two-dimensional search initiated from it are also presented in Fig. 1. The objective function F(a) is clearly convex in the displayed region, with a trough expanding from the upper left to the lower right comer of the plot and reaching a well defined global minimum in the center of the plot. This supports our anticipation of the existence of a finite optimal forcing in the thermodynamic equation, at least in one- and twodimensional cases. Figure 1 also demonstrates how a quasiNewton optimization routine locates the minimum. A unit diagonal matrix is a standard first approximation of the Hessian if no additional information is initially available on the curvature of the objective function. As follows from equation (8), if H(“’ is a unit matrix, the corresponding first step is performed in the direction of negative gradient, the steepest descent direction. In each case presented in Fig. 1, the first step approximately locates the bottom of the valley that is then followed toward the minimum.
diagnostic
scheme
989
It is interesting to notice here that direct search methods, based on function evaluations only, are known to often fail in following slanted valleys similar to that presented in Fig. 1 (e.g. Gottfried and Weisman, 1973). It is also worth noticing that in some cases presented in Fig. 1 comparatively long steps toward the minimum have been made along the valley, almost in parallel to isolines of the objective function. This would not be possible for steepest descent methods that only use local information on the objective function and its gradient. In general, such a method would expend more steps to approach an extremum in the direction of steepest descent. In particular, more steps would be necessary in close proximity to the minimum where the function gradient gets smaller and its exact direction, subject to computational errors and proper scaling of variables, becomes less definite, and if the valley is curved. Finally, a few full-scale (L = 20) experiments have been performed to optimize the nudging coefficients in each model layer independently. In particular, the one- and two-dimensional minima have been used as initial conditions. On average, a full-scale search requires about IO-15 major iterations with approximately 30-35 function evaluations each, so that these experiments are extremely computationally demanding. Thus L = 20 appears to be close to the upper limit of practical applicability of direct optimization algorithms in atmospheric model parameter evaluation for comparatively sophisticated models. In. all cases, the search consistently converges to a minimal or = 2.71 K. Table 1 gives an example of the principle of diminishing returns: a ten-fold increase in the problem dimensionality and in the amount of computations provides only a marginal improvement of the simulation results. Stated differently, these experiments demonstrate that a proper adjustment of a limited number of atmospheric model parameters can be performed using direct optimization techniques. These adjustments may yield results very close to the overall optimum for a particular model. Figure 2 presents optimal nudging coefficients for L=l, 2, and 20. One can only speculate that the increase of LXwith altitude by more than three orders of magnitude for L = 2 is somehow related to the decrease of pressure and density by a factor of about 106, and to the corresponding
R. A. Akmaev
990 OPTIMAL NUDGING COEFFICIENT (s-l) 1 I I I I I ,
Latitude
(deg)
Latitude
(deg)
Fig. 3. Zonal mean temperature for March (contour interval 10 K): (a) simulations for L = 20; (b) CIRA-86.
Fig. 2. Optimal nudging coefficient: L = 1(vertical line); L = 2 (slanted straight line); L = 20 (curved line).
decrease in characteristic times of physical processes. It is interesting that the irregular structure of o((z) for L = 20 appears to be typical for largescale nudging optimization. In their experiments with a low penalty for deviations of the nudging coefficients from a prespecified average, Stauffer and Bao (1993) also obtained highly dispersive nudging coefficients even reaching negative values in some points. This erratic behavior appears to be a manifestation of the fact that in 20 dimensions the objective function F(N) is almost ‘flat’ in a comparatively wide region around the twodimensional optimum, and does not have a well defined global minimum, as opposed to the case presented in Fig. 1. This is also confirmed by the comparatively small gain attained by switching from L = 2 to L = 20 (Table l), and leads again to the conclusion that a solution close to optimal can be reached in atmospheric model parameter evaluation problems by way of low-dimensional optimization at a relatively low computational cost. Figures 3 and 4 compare, for L = 20, the simulated monthly mean temperatures with the empirical model for March and December, respectively. Figures 5 and 6 present the total zonal forcing for the same months obtained with no nudging (L = 0) and with optimal nudging (L = 20), respectively. One can see some noticeable changes in the shape and strength of the forcing, for example, in the upper mesosphere where the zonal drag is comparatively strong. In particular, for the solstice conditions the maximum of forcing located near the mesopause appears to split into two in
Latitude
(cleg)
Latitude
(deg)
Fig. 4. As in Fig. 3 but for December.
Latitude (dag)
Latitude
(deg)
Fig. 5. Total diagnosed zonal forcing without nudging (L = 0): (a) March (contour interval IO ms-’ day-‘); (b) December (contour interval 20 m s-l day-‘).
Latitude
(deg)
Latitude
(deg)
Fig. 6. As in Fig. 5 but for optimal (L = 20) forcing.
Optimization of a middle atmosphere diagnostic scheme
Lditvde
(kg)
Latitude
(kg)
Fig. 7. Temperature difference (contour interval 2 K) between simulations and CIRA-86 without nudging (L = 0): (a) March; (b) December. DEC (b)
Latitude
(deg)
Latitude
(deg)
Fig. 8. As in Fig. 7 but for optimal (L = 20) nudging.
991
scheme provides adequate estimates of the waveinduced zonal drag in the middle atmosphere. It is clear that this work can be extended in the direction of optimization of the simulation model itself. Upper and middle atmospheric models usually contain some poorly known parameters such as, for example, those that describe momentum and heat deposition by dissipating and breaking waves, intensity of turbulence produced by these waves etc. These parameters appear to be the ideal candidates for application of formal optimization procedures. It appears that the results of this work and the methods explored herein can, for example, be applied to optimization of currently available gravity wave parameterizations in two ways. First, the parameterization schemes can be tested and thoroughly adjusted on the diagnosed monthly wave drag distributions. This seems to be a much more efficient way compared to the traditional ‘manual tuning.’ Second, parameterization schemes with a comparatively low number of ‘free’ adjustable parameters can be directly optimized within a simulation model to produce the best agreement with empirical climatologies.
both hemispheres. This feature was first reported by Akmaev et al. (1992). These changes in the
5. CONCLUSION
forcing produce an overall improvement of the simulated temperatures (compare Figs 7 and 8). The intermediate cases (L=l and 2) are not presented here for brevity, and we just mention that they produce the zonal forcing and simulated temperature distributions generally fitting between the two extremes. It is important to notice that the temperatures corresponding to the L = 2 case are very close to those obtained for L = 20, as can also be seen from Table 1. Medvedev and Fomichev (1994) obtained zonal forcing for the same atmospheric region and empirical climatology using a similar radiative code within the traditional framework of computing the so-called diabatic circulation. In a previous paper (Akmaev, 1994) qualitative comparisons of their results with ours were presented and possible sources of apparent differences were discussed. It should be stressed here again that only the present diagnostic-simulation scheme uses objective quantitative criteria based upon which it is now possible to establish the representativeness of diagnostics. The results presented herein and by Akmaev (1994) confirm that the new diagnostic
Results of a series of optimization experiments with an assimilative middle atmospheric diagnostic-simulation scheme based on the use of a spectral model are presented. The scheme uses the annual global mean temperature deviation as an objective criterion of the ability of the simulations to reproduce empirical climatology. This quantitative criterion makes it possible to perform parameter optimization of the diagnostic scheme and/or the simulation model. It is demonstrated, in particular, that a proper adjustment of the strength of nudging forcing in the thermodynamic equation at the diagnostic stage can reduce the prognostic rms temperature deviation from the CIRA-86 empirical model to approximately 2.7 K in the 15-I 10 km layer. It is demonstrated that direct optimization procedures can effectively be used in atmospheric model parameter estimation problems of moderate dimensionality and may be recommended for parameter optimization of simulation models. The ‘optimal nudging’ data assimilation technique used in this study at the diagnostic stage is
R. A. Akmaev
992
apparently subject to fast saturation with increasing dimensionality of the problem, suggesting that a low-cost, low-dimension optimal nudging may yield results close to those obtained by a large-scale optimization. Acknowledgements-Maura Hagan and two anonymous reviewers offered many useful comments on early versions of the manuscript. The National Center for Atmospheric Research ence Foundation.
is sponsored
by the National
Sci-
REFERENCES Akmaev, R. A. (I 994) Diagnostics and simulation of an annual cycle in the middle atmosphere. 1 Geophys. Res. 99, 18,933-18,950. Akmaev, R. A., Fomichev, V. I., Gavrilov, N. M. and Shved, G. M. (1992) Simulation of the zonal mean climatology of the middle atmosphere with a threedimensional spectral model for solstice and equinox conditions, 1 Atmos. Terr Phys. 54, 119-128. Anthes, R. A. (1974) Data assimilation and initialization of hurricane prediction models. 1 Atmos. Sci. 31, 702-719. Apruzese, J. I?, Schoeberl, M. R. and Strobel, D. F. (1982) Parameterization of IR cooling in a middle atmosphere dynamics model. I. Effects on the zonally averaged circulation. 1 Geophys. Res. 87, 8951-8966. Courtier, I?, Derber, J., Errico, R., Louis, J.-F. and VukiCeviC, T. (1993) Important literature on the use of adjoint, variational methods and the Kalman filter in meteorology. Tellus &A, 342-357. Daley, R. (I 99 I) Atmospheric Data Analysis. Cambridge University Press, 457 pp. Dennis, J. E. and Schnabel, R. B. (1983) Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, Englewood Cliffs, NJ, 378 pp. Fleming, E. L., Chandra, S., Barnett, J. J. and Corney, M. (I 990) Zonal mean temperature, pressure, zonal wind
and geopotential height as functions of latitude. Adv. Space Res. 10, (12)l l-(12)59. Fletcher, R. (1980) Practical Methods of Optimization. John Wiley & Sons, Chichester, 120 pp. Ghil, M. and Malanotte-Rizzoli, P (1991) Data assimilation in meteorology and oceanography. Adv. Geophys. 33, 141-266. Gill, P E., Murray, W. and Wright, M. H. (1981) Pructical Optimization, Academic Press, London, 401 pp. Gottfried, B. S. and Weisman, J. (1973) Introduction to Optimization Theory Prentice-Hall, Englewood Cliffs, NJ, 571 pp. Hestenes, M. R. (1980) Conjugate Direction Methods in New York, 325 pp. Optimization. Springer-Verlag, Medvedev, A. S. and Fomichev, V. I. (I 994) Net radiative heating and diagnostics of the diabatic circulation in the 15-l IO km height layer. J Atmos. Terr Phys. 56, 1571-1584. Murgatroyd, R. J. and Singleton, F. (1961) Possible meridional circulations in the stratosphere and mesosphere. Quart. J Roy Met. Sot. 87, 125-135. Navon, I. M. and Legler, D. M. (1987) Conjugategradient methods for large-scale minimization in meteorology. Mon. tiea. Rev. 115, 1479-1502. Shine, K. P. (1989) Sources and sinks of zonal momentum in the middle atmosphere diagnosed using the diabatic circulation. Quart. J. Roy Met. Sot. 115,
265-292. Shine, K. P and Rickaby, J. A. (1989) Solar radiative heating due to absorption by ozone. In Ozone in the Atmosphere, eds R. D. Bojkov and P. Fabian, pp. 597-600. A. Deepak Publ., Hampton, VA. Stauffer, D. R. and Bao, J.-W. (1993) Optimal determination of nudging coefficients using the adjoint Tellus 45A, 358-369. equations. Strobel, D. F. (1978) Parameterization of the atmospheric heating rate from I5 to 120 km due to 02 and 03 absorption of solar radiation. J Geophys. Res. 83, 6225-6230. Zou, X., Navon, I. M. and Le Dimet, E X. (1992) An optimal nudging data assimilation scheme using parameter estimation. Quart. J Roy Met. Sot. 118, 1163-l 186.
Journal
Pergamon PII:
of Amosphrric
and
Optimization of a middle atmosphere diagnostic scheme Rashid A. Akmaev * High Altitude Observatory, National Center for Atmospheric Research, Boulder, Colorado, U.S.A. (Received in final form 2 January 1996; accepted 20 January 1996)
Abstract-A new assimilative diagnostic scheme based on the use of a spectral model was recently tested on the CIRA-86 empirical model. It reproduced the observed climatology with an annual global rms temperature deviation of 3.2 K in the 15-l IO km layer. The most important new component of the scheme is that the zonal forcing necessary to maintain the observed climatology is diagnosed from empirical data and subsequently substituted into the simulation model at the prognostic stage of the calculation in an annual cycle mode. The simulation results are then quantitatively compared with the empirical model, and the above mentioned rms temperature deviation provides an objective measure of the ‘distance’ between the two climatologies. This quantitative criterion makes it possible to apply standard optimization procedures to the whole diagnostic scheme and/or the model itself. The estimates of the zonal drag have been improved in this study by introducing a nudging (Newtonian-cooling) term into the thermodynamic equation at the diagnostic stage. A proper optimal adjustment of the strength of this term makes it possible to further reduce the rms temperature deviation of simulations down to approximately 2.7 K. These results suggest that direct optimization can successfully be applied to atmospheric model parameter identification problems of moderate dimensionality. @ 1997 Elsevier Science Ltd
1.
ing data assimilation technique (e.g. Anthes, 1974; Ghil and Malanotte-Rizzoli, 1991) is iteratively integrated towards a steady state. The final result of these iterations is the total zonal mean zonal drag obtained in such a way as to keep the model climatology in the vicinity of empirical data. The close connection of the diagnostic scheme to the simulation model facilitates an objective assessment of the retrieved zonal drag. At the second stage, the zonal forcing is substituted back into the unmodified model, and the model climatology is then quantitatively compared with the initial empirical climatology. Akmaev (1994) tested this scheme on the CIRA-86 empirical model (Fleming et al., 1990). The simulations using the retrieved zonal forcing were shown to reproduce the observed climatology with an annual global rms temperature deviation of 3.2 K or about 2% in the 15-I IO km layer. The closeness of the simulations to the empirical model suggested for the first time that the diagnosed forcing yielded reliable quantitative estimates of the intensity of real physical processes governing the global-scale circulation in the middle atmosphere.
INTRODUCTION
Since the work of Murgatroyd and Singleton (1961), diagnostics of observed climatologies have been used to study the mechanisms maintaining the gross features of the global-scale circulation and tracer transport in the middle atmosphere. The zonal drag or forcing imposed by dissipating and breaking waves of all scales from global to mesoscale, is of particular importance and interest (e.g. Shine, 1989; Medvedev and Fomichev, 1994). Akmaev (1994) recently proposed a new assimilative middle atmospheric diagnostic scheme and compared it with more traditional diagnostic procedures (see also, Shine, 1989; Medvedev and Fomichev, 1994, and references therein). The scheme employed in this study and detailed by Akmaev (1994) is based on the use of a spectral model developed by Akmaev et al. (1992). At the diagnostic stage the model with the zonal momentum equation modified in the spirit of the nudg-
* Present address: Department of Aerospace Engineering Sciences, Campus Box 429, University of Colorado, Boulder, CO 80309-0429, U.S.A.
983
R. A. Akmaev
984
Both the spectral model and the diagnostic scheme as a whole have some adjustable parameters. The diagnostic scheme is also quite flexible and it may be subjected to further modification via the introduction of additional parameters. This flexibility and the existence of the above mentioned quantitative criterion make it possible to apply thorough optimization procedures to the simulation model, the diagnostic module, or to the whole diagnostic-simulation complex. The first results of such optimization experiments are reported in this paper. It has been proposed, for example, that estimates of the zonal drag could be improved by introducing a nudging, or Newtonian-cooling term into the thermodynamic equation at the diagnostic stage (Akmaev, 1994). As is demonstrated below, a proper adjustment of the strength of this Newtonian-cooling term makes it possible to further reduce the rms temperature deviation down to approximately 2.7 K. These results strongly support the objective quantitative approach to middle atmospheric modeling and diagnostics. The results also clearly demonstrate that direct optimization can successfully be applied to atmospheric model parameter identification problems of moderate dimensionality as opposed, for example, to the more demanding adjoint model technique (e.g. Daley, 1991; Ghil and Malanotte-Rizzoli, 1991) (see also the bibliography collected by Courtier et al., 1993), the method of choice for large-dimension optimization problems. Section 2 briefly describes the diagnosticsimulation scheme and its modifications introduced in this study. Section 3 presents direct optimization routines suitable for the problem. In Section 4 the results of a series of optimization experiments are presented and discussed.
2. DIAGNOSTICSIMUL.ATlON
SCHEME
The diagnostic-simulation scheme employed in this study is based on the middle atmospheric spectral model developed and described in more detail by Akmaev et a/. (1992) and Akmaev (1994). A global three-dimensional spectral model with a triangular truncation T13 corresponding to a moderate latitude-longitude resolution of about 8.6” x 9” is used. The model covers the altitude range between 100 mb (-15 km) and
17 pressure scale heights (- 120 km) with 20 layers of equal thickness in log-pressure coordinates with a vertical resolution of approximately 5 km. Zonally-symmetric diagnostics are performed as predetermined by the zonal symmetry of empirical data, so the model has been run in a two-dimensional mode with the longitudinal dependency of model variables ‘turned off.’ As was suggested by Akmaev (1994) the radiation transfer module has been updated by replacing the parameterization of the ozone ultraviolet absorption of Strobe1 (1978) and Apruzese et al. (1982) with that of Shine and Rickaby (1989), providing better results in the stratosphere due to enhanced heating primarily in the Huggins bands. The diagnostic scheme and the approximations involved are presented and extensively discussed by Akmaev (1994). At the diagnostic stage the model is initialized with the CIRA-86 empirical temperatures and zonal winds (Fleming et af., 1990) for each month. It is then integrated towards a quasi-equilibrium state and the total zonal drag, aZ, necessary to maintain the empirical climatology is determined from the following iterative procedure. The zonal momentum equation can be written in the form:
au = R, -t-a=, at
-
where u is the zonal mean zonal wind and R, represents all other terms in the equation that can be explicitly computed in a zonally symmetric model. At every nth time step aZ is estimated as
z = 0.5 (ZT-” -R:). The averaging of two consecutive iterations in equation (2) is used to prevent the model from spurious oscillations. The initial estimate 2 may be set to zero, or taken from previous diagnostic results. This procedure generally yields reasonable estimates of aZ if the integration continues for IO-15 model days from the initial zero conditions and if the model remains in close vicinity to the empirical climatology during this integration. Akmaev (1994) suggested that these estimates could potentially be improved if an additional nudging term is inserted into the thermodynamic equation to more strongly constrain the model climatology towards observations at the diagnostic stage. Specifically,
Optimization
aT
,t=RT-a(T-T“),
of a middle
atmosphere
(3)
where RT is the right-hand side of the standard thermodynamic equation, 01 is the nudging (Newton-cooling) coefficient, and To is the ‘observed’ temperature. In this study, o( = o((z) is assumed to depend on the vertical coordinate only. A conceptually straightforward and computationally efficient nudging or relaxation technique that may be considered a degraded form of the Kalman filter (Daley, 1991; Ghil and MalanotteRizzoli, 1991) was first introduced for atmospheric data assimilation by Anthes (1974) and continues to be extensively investigated for use in meteorology and oceanography. It is usually applied to prognostic model equations, in a fashion similar to equation (3), to force the model to ‘account’ for available observations, or, in other words, to blend observations and simulations in an ‘optimal’ manner. Generally, o( is a matrix which may depend on the time differential between times of observations and model time steps, on the spatial distance between points of observations and model grid points, on the standard error of observations, etc. The optimality of combining experimental data with numerical models in this case greatly depends on the strength of nudging coefficients. The method of adjoint model variational data assimilation has been used recently in the parameter estimation regime to determine ‘optimal’ forcing (Zou et al., 1992; Stauffer and Bao, 1993). Zou et al. (1992) applied the technique to a three-dimensional spectral adiabatic general circulation model to find constant nudging coefficients appearing in three model equations. Stauffer and Bao (1993) used a onedimensional shallow-water model in their search for three nudging coefficients depending on the horizontal coordinate. It is clear that the stronger the relaxation terms in prognostic equations are, the closer the model will reproduce observations. Therefore, a cost function whose minimum is sought in variational data assimilation, and that usually represents a measure of the distance between simulations and observations, had to be augmented in these studies by a penalty term depending on the deviation of the nudging coefficients from some a priori estimates. This essentially transferred the uncertainty from the nudging coefficients onto these preliminary es-
diagnostic
985
scheme
timates and penalty coefficients. Stauffer and Bao (1993) also studied the sensitivity of the results to the penalty term, and found that, as one would expect, the ‘optimal’ nudging coefficients increase as the weight of the penalty term decreases. Unlike the above mentioned studies, we apply nudging only at the diagnostic step in this study. One can speculate that there must be some finite optimal forcing which will produce best estimates of a,. Indeed, if o( = 0 in equation (3) the model may eventually deviate from the initial empirical climatology in the process of iterative integration towards a steady state. If o( is too large, however, the model climatology will almost entirely be determined by the fictitious forcing in equation (3) and the corresponding estimates of a, may not be sufficient to adequately force the model towards observations at the simulation stage. The purpose of this work is to find the optimal finite nudging and obtain the best estimates of us. Another feature of the present scheme that makes it different from those employed by Zou et al. (1992) and Stauffer and Bao (1993) is that a direct optimization procedure (see Section 3) is used for our comparatively low-dimension problem as opposed to the costly adjoint model technique. Let us now summarize the whole procedure used in this study which differs slightly from that employed by Akmaev (1994). For a given o( = a(z) the model with the modified zonal momentum and thermodynamic equations [equations (1) and (3) respectively] is initialized by the zonal mean temperature and zonal wind distributions from CIRA-86 for each month, and z = 0. The model is then integrated for n = 300 steps, or 12.5 model days, with the zonal forcing being calculated according to equation (2). This integration time interval has been found sufficient to reach a quasi-equilibrium state. The retrieved monthly mean zonal forcing and the vertical diffusion coefficient determined as in Akmaev (1994) are then used in an annual run of the unmodified model, and the following objective or cost function is computed:
F=
2rlog:p,
,p*) ;,
rrn jP,d(logp) P2
n/2 x
( T,,, -
T,”)‘cos 4p d+.
(4)
986
R. A. Akmaev
Here F = u: is the annual altitude global mean temperature deviation squared, pi is set at the model lower boundary to 100 mb, and p2 = 1.8 x low4 mb is the pressure level two model layers below the upper boundary at an altitude of about 110 km, set to filter out the effect of the upper boundary condition (Akmaev, 1994). The length of month m is TV, the length of year is represented by T = C,T,, #J is latitude, T,,,(c$,p) is the monthly averaged simulated temperature, and T,O(@,p) is the ‘observed’ temperature. In equation (4) we only compare temperatures, as the empirical tables contain geostrophic zonal winds simply calculated from temperature distributions. or values obtained with different nudging are then compared to find, as described in the following sections, the optimal nudging that delivers the minimal rms temperature deviation. The ‘observed’ temperatures T,” are calculated from the CIRA-86 tables by interpolation and extrapolation onto the model horizontal Gaussian and vertical grids. In data assimilation schemes, it is usually model results that are projected onto the points of observations in space and/or time. However, in this scheme we use the inverse transformation for the following reasons. As mentioned above, T,” and the corresponding empirical zonal winds are used to initialize the model at the diagnostic stage, and thus need to be transformed onto the model spatial grid anyway. To compute the integrals in equation (4) some kind of interpolation has to be used too. The empirical data are presented on a regular space mesh and have, incidentally, been obtained by interpolation from original experimental results. Essentially, the nudging terms represent a mathematical device, a ‘trick,’ the sole purpose of which is to force a numerical model to produce a desirable result, and it is perhaps pointless to seek any physical meaning in the strength of the nudging coefficients or in their spatial or temporal distribution (see also Section 4). Again, we stress a very significant difference between the traditional applications of the nudging technique referred to earlier in this section and that employed in this work. We apply nudging only at the diagnostic stage. The forcing obtained from the diagnostics is then used in direct simulations performed with the unmodified model. This is done without any fictitious nudging, unlike simulations done in the
traditional approach. As is shown in Section 4, the retrieved zonal drag makes it possible for the model to very closely reproduce the empirical climatology. One can therefore argue that the forcing provides a representative quantitative estimate of the strength of real physical processes taking place in the atmospheric layer under consideration, to the extent, of course, that the empirical chmatology itself can be considered representative.
3. OPTIMIZATION
PROCEDURE
The purpose of optimization in this study is to minimize the objective function F(N) from equation (4) with respect to the nudging coefficients, 01 = a(r), which are assumed to depend on the vertical coordinate only. A series of numerical experiments have been performed with different characters of this height dependence (see Section 4). In a discrete form o( = ((~1,. . , a~)‘, where the prime denotes a transpose, may be considered a vector of a variable dimension L ranging from 0 (no nudging) through to 20, the total number of model layers. Newton-type methods are perhaps most effective for smooth unconstrained optimization problems of moderate dimension (Fletcher, 1980; Gill et al., 1981; Dennis and Schnabel, 1983). QuasiNewton low-memory procedures can even be successfully used in large-scale problems (Navon and Legler, 1987). Basically, a Newton-type method uses the available information on the objective function F(x), its first derivatives, or gradient g(x) =
, (g,...,$- >’ L
1
and its second derivatives, or the Hessian matrix, H(x) = {H;j}, where Hij
= -.
a2F
axi axj
’
i,j=l
, . , L.
(6)
A quadratic model, Q(x), of the objective function is constructed at every kth consecutive search point xtk): Qtk’(x) = Fck’ + (Ax)‘~‘~’ + 0.5(A~)‘H(~)Ax,(7) where Fck) = F(xck)) Ax = X-X(~) , g’k’ = g(x’k’), Hck) = H(x”)), and ;he right-hand side of equation (7) is simply a truncated Taylor series for F.
Optimization of a middle atmosphere diagnostic scheme
A search vector, $4 = s(X’k’~ = _ (H(k))-’ g’k’
(8)
calculated as a product of the inverse Hessian and negative gradient will then exactly point to the minimum ,?ck) of the quadratic approximation if the minimum exists: z.(k)
=
$k)
+ S(k).
(9)
Apart from some fine details, for a general nonlinear function, a linear search is conducted along the direction of search vector ~(~1.The whole process is repeated from the next search point .x(~+” which is determined as a minimum along the search direction:
F(X’k’“)
= t$ F(X’k’ + BP). >
(IO)
If the objective function is available analytically so that its gradient and Hessian matrix can be calculated exactly, perhaps the only computational problem is to invert the Hessian in equation (8). The exact Newton methods are characterized by a high, quadratic, convergence rate which essentially reflects the fact that the closer a particular search point x(k) is to the minimum sought x*, the better the quadratic model Qck)(x) approximates F(x), and the more precisely its minimum ~2~) represents x*. If a search point is far from the minimum, the quadratic model in equation (7) may not be adequate, and a Newton-type procedure can diverge, or the search sck), may not point in the direction of substantial function decrease. Therefore Newton steps have to be combined with some global-search moves to ensure that the search does not terminate, for example, in a local minimum instead of the global one within attainable prespecified limits. Also, some restart procedures can be used to establish a new search direction, if the ‘quadratic’ direction does not provide a sufficient function decrease. Optimization methods utilizing information only on the objective function itself (direct search methods) and/or its first derivatives (steepest descent methods) have been demonstrated to be inferior in their convergence rate. They may often fail entirely for objective functions of a comparatively complex structure.
987
Unfortunately, in many practical applications, as is the case in this study, every objective function value is a result of a complex numerical procedure and its second, and often even first, derivatives are not readily available. A group of so-called quasi-Newton optimization algorithms, closely related to the conjugate direction methods (Hestenes, 1980; Navon and Legler, 1987) exploits the idea of building up, at the expense of a slower convergence rate, the information on the curvature of objective function based on its gradients at consecutive search points. The Hessian matrix is successively updated in the search process starting from some initial approximation. It is also possible, as in some more advanced algorithms, to directly update the Cholesky decomposition factors of the Hessian, and to save a large amount of computation on its inversion. The adjoint model technique essentially provides an effective means to calculate first derivatives of various objective functions which result from some ‘direct’ modeling, particularly for large-scale problems. For problems of moderate dimension, the gradient can be obtained by finite differences. Therein, each gradient component requires at least one additional function call, if forward differences are used, and approximately the same order of additional calls, if central differences have to be used to improve precision. Some important issues have to be addressed in practical applications of optimization algorithms (Fletcher, 1980; Gill et al., 1981; Dennis and Schnabel, 1983). Among them are the proper scaling of the objective function and variables, or variable transformation if necessary. It is highly desirable to prescribe a feasible search region, maximal and minimal search steps, stopping criteria, and some other quantities adequately reflecting the nature of problem at hand. Ill-posed or inappropriately scaled problems may cause substantial loss of accuracy and prevent the algorithm from locating the global minimum, or even cause its total failure. It has been found after some limited experimentation with ‘easy-to-use’ and ‘comprehensive’ optimization routines available in the IMSL and NAG scientific computing libraries, that the latter provides routines which are slightly more efficient for our particular problem in that under equal conditions they generally require a smaller number of function calls to locate a minimum. Of
988
R. A. Akmaev
course, this could be just a matter of more adequate tuning of a particular routine. The results presented in Section 4 have been obtained using the NAG FORTRAN library (Gill et al., 1981) ‘comprehensive’ quasi-Newton routine E04UCE.
4. RESULTS AND DISCUSSION
A series of optimization experiments of different dimensionality L has been conducted. Akmaev (1994) presented the results of diagnostics and simulation, and calculated ar for a second model year. Every calculation of function F from equation (4), reported herein, involves a run of the diagnostic scheme for all twelve months and a run of the model for one year instead of two to save computer time. The differences in or for the first and second model years were found negligible, of the order of 0.003 K or about 0.1%. First we assess the effect of replacement of the parameterization of ozone ultraviolet heating by Strobe1 (1978) and Apruzese ef al. (1982) by that of Shine and Rickaby (1989) without nudging (L = 0) in the thermodynamic equation (see Table 1). This ‘optimization’ of the radiation modTable I, Root-mean-square
L
CT 6)
0 0
3.22 3.06 2.88 2.76 2.71
1 2 20
temperature deviation
Ozone heating Strobe1 (1978); Apruzese et al. (1982) Shine and Rickaby (1989) Shine and Rickaby (1989) Shine and Rickaby (1989) Shine and Rickaby (1989)
ule leads to improved simulations, primarily in the stratosphere, due to slightly enhanced heating throughout the region. The corresponding rms temperature deviation fir decreased from 3.22 K to 3.06 K. Of course, this lower ur is produced by the zonal forcing that slightly differs from that obtained in the previous study (Akmaev, 1994) and corresponds to a slightly different temperature distribution. These minute differences can be seen by comparison of Figs 5 and 7 in this paper with Figs 11 and 15 of Akmaev (1994). It has been found that the optimization routine described in the previous section is stable and efficiently locates the global minimum of the objective function, if F = F(x) is optimized with respect to x instead of a, where
XI = log,, a/,
I = l....,L.
(11)
This transformation also eliminates the only constraint on the variables o(/ 2 0, and the optimization problem becomes truly unconstrained. Optimization has been performed with the optimahty tolerance parameter (Gill et af., 1981) specified at c = 10m4. This approximately corresponds to five correct digits in the final value of F or or. A one-dimensional search (L = I), where o((z) = 01~ is assumed to be constant with altitude, typically converges in 2-3 steps from arbitrary initial conditions to at = 9.93 x 10m4s-i. The corresponding or = 2.88 K. For L = 2, we assume that LXexponentially depends on logp between the bottom and the top model layers, where 01 = (~1 and LX= o(2, respectively. The equivalent assumption is that x in equation (11) changes linearly with logp from xi to x2. First, the result of one-dimensional optimization was used as an initial condition, and the search converged in 3 steps to or = 2.76 K. Each major step in this case required 6-10 evaluations of the objective function to calculate its gradient, to perform a linear search along the search direction as in equation (lo), and global search moves if necessary. To ensure that the minimum found is global, and to explore the objective function behavior in a relatively broad area around the minimum, eight additional experiments have been performed with initial conditions separated from the solution by two orders of magnitude in (Y~and/or 012.The results of these experiments are presented in Fig. 1 where every dot represents a major search step. To get a feel for the accuracy of the search, we notice that, in addition to the above mentioned optimality tolerance, E, another stopping criterion, the so-called line search tolerance (Gill et al., 1981), has been specified at 6 = 0.01. This approximately corresponds to one tenth of the dot diameter in Fig. 1. Roughly speaking, a search is completed if two consecutive major steps fall within the line search tolerance, 6, if the relative difference between corresponding function values does not exceed the optimality tolerance, E, and if some other necessary and sufficient conditions for a minimum are satisfied. The eight initial points are located on the perimeter of the plot. In most cases the search required 4 steps, 3 steps were enough in one case, and 6 steps were required in the two worst cases. The one-dimensional case corresponds to a constrained ((xi = 012) mini-
Optimization
of a middle
atmosphere
6 (K*)
1d-s
,;-4
1d-s a1
10-j
(s-7
Fig. 1. Two-dimensional optimization (contour interval I K2). Dots represent major search steps. The eight starting points separated from the global minimum (big cross in the center) by 2 orders of magnitude in cxt and/or cx2 are located on the perimeter. The onedimensional optimal solution corresponding to a constrained minimum (the constraint at = CQ is represented by a dashed line) is presented in the lower right corner of the plot.
mum. The constraint is represented by a straight dashed line in the lower right comer of the plot, and the corresponding one-dimensional minimum and the two-dimensional search initiated from it are also presented in Fig. 1. The objective function F(a) is clearly convex in the displayed region, with a trough expanding from the upper left to the lower right comer of the plot and reaching a well defined global minimum in the center of the plot. This supports our anticipation of the existence of a finite optimal forcing in the thermodynamic equation, at least in one- and twodimensional cases. Figure 1 also demonstrates how a quasiNewton optimization routine locates the minimum. A unit diagonal matrix is a standard first approximation of the Hessian if no additional information is initially available on the curvature of the objective function. As follows from equation (8), if H(“’ is a unit matrix, the corresponding first step is performed in the direction of negative gradient, the steepest descent direction. In each case presented in Fig. 1, the first step approximately locates the bottom of the valley that is then followed toward the minimum.
diagnostic
scheme
989
It is interesting to notice here that direct search methods, based on function evaluations only, are known to often fail in following slanted valleys similar to that presented in Fig. 1 (e.g. Gottfried and Weisman, 1973). It is also worth noticing that in some cases presented in Fig. 1 comparatively long steps toward the minimum have been made along the valley, almost in parallel to isolines of the objective function. This would not be possible for steepest descent methods that only use local information on the objective function and its gradient. In general, such a method would expend more steps to approach an extremum in the direction of steepest descent. In particular, more steps would be necessary in close proximity to the minimum where the function gradient gets smaller and its exact direction, subject to computational errors and proper scaling of variables, becomes less definite, and if the valley is curved. Finally, a few full-scale (L = 20) experiments have been performed to optimize the nudging coefficients in each model layer independently. In particular, the one- and two-dimensional minima have been used as initial conditions. On average, a full-scale search requires about IO-15 major iterations with approximately 30-35 function evaluations each, so that these experiments are extremely computationally demanding. Thus L = 20 appears to be close to the upper limit of practical applicability of direct optimization algorithms in atmospheric model parameter evaluation for comparatively sophisticated models. In. all cases, the search consistently converges to a minimal or = 2.71 K. Table 1 gives an example of the principle of diminishing returns: a ten-fold increase in the problem dimensionality and in the amount of computations provides only a marginal improvement of the simulation results. Stated differently, these experiments demonstrate that a proper adjustment of a limited number of atmospheric model parameters can be performed using direct optimization techniques. These adjustments may yield results very close to the overall optimum for a particular model. Figure 2 presents optimal nudging coefficients for L=l, 2, and 20. One can only speculate that the increase of LXwith altitude by more than three orders of magnitude for L = 2 is somehow related to the decrease of pressure and density by a factor of about 106, and to the corresponding
R. A. Akmaev
990 OPTIMAL NUDGING COEFFICIENT (s-l) 1 I I I I I ,
Latitude
(deg)
Latitude
(deg)
Fig. 3. Zonal mean temperature for March (contour interval 10 K): (a) simulations for L = 20; (b) CIRA-86.
Fig. 2. Optimal nudging coefficient: L = 1(vertical line); L = 2 (slanted straight line); L = 20 (curved line).
decrease in characteristic times of physical processes. It is interesting that the irregular structure of o((z) for L = 20 appears to be typical for largescale nudging optimization. In their experiments with a low penalty for deviations of the nudging coefficients from a prespecified average, Stauffer and Bao (1993) also obtained highly dispersive nudging coefficients even reaching negative values in some points. This erratic behavior appears to be a manifestation of the fact that in 20 dimensions the objective function F(N) is almost ‘flat’ in a comparatively wide region around the twodimensional optimum, and does not have a well defined global minimum, as opposed to the case presented in Fig. 1. This is also confirmed by the comparatively small gain attained by switching from L = 2 to L = 20 (Table l), and leads again to the conclusion that a solution close to optimal can be reached in atmospheric model parameter evaluation problems by way of low-dimensional optimization at a relatively low computational cost. Figures 3 and 4 compare, for L = 20, the simulated monthly mean temperatures with the empirical model for March and December, respectively. Figures 5 and 6 present the total zonal forcing for the same months obtained with no nudging (L = 0) and with optimal nudging (L = 20), respectively. One can see some noticeable changes in the shape and strength of the forcing, for example, in the upper mesosphere where the zonal drag is comparatively strong. In particular, for the solstice conditions the maximum of forcing located near the mesopause appears to split into two in
Latitude
(cleg)
Latitude
(deg)
Fig. 4. As in Fig. 3 but for December.
Latitude (dag)
Latitude
(deg)
Fig. 5. Total diagnosed zonal forcing without nudging (L = 0): (a) March (contour interval IO ms-’ day-‘); (b) December (contour interval 20 m s-l day-‘).
Latitude
(deg)
Latitude
(deg)
Fig. 6. As in Fig. 5 but for optimal (L = 20) forcing.
Optimization of a middle atmosphere diagnostic scheme
Lditvde
(kg)
Latitude
(kg)
Fig. 7. Temperature difference (contour interval 2 K) between simulations and CIRA-86 without nudging (L = 0): (a) March; (b) December. DEC (b)
Latitude
(deg)
Latitude
(deg)
Fig. 8. As in Fig. 7 but for optimal (L = 20) nudging.
991
scheme provides adequate estimates of the waveinduced zonal drag in the middle atmosphere. It is clear that this work can be extended in the direction of optimization of the simulation model itself. Upper and middle atmospheric models usually contain some poorly known parameters such as, for example, those that describe momentum and heat deposition by dissipating and breaking waves, intensity of turbulence produced by these waves etc. These parameters appear to be the ideal candidates for application of formal optimization procedures. It appears that the results of this work and the methods explored herein can, for example, be applied to optimization of currently available gravity wave parameterizations in two ways. First, the parameterization schemes can be tested and thoroughly adjusted on the diagnosed monthly wave drag distributions. This seems to be a much more efficient way compared to the traditional ‘manual tuning.’ Second, parameterization schemes with a comparatively low number of ‘free’ adjustable parameters can be directly optimized within a simulation model to produce the best agreement with empirical climatologies.
both hemispheres. This feature was first reported by Akmaev et al. (1992). These changes in the
5. CONCLUSION
forcing produce an overall improvement of the simulated temperatures (compare Figs 7 and 8). The intermediate cases (L=l and 2) are not presented here for brevity, and we just mention that they produce the zonal forcing and simulated temperature distributions generally fitting between the two extremes. It is important to notice that the temperatures corresponding to the L = 2 case are very close to those obtained for L = 20, as can also be seen from Table 1. Medvedev and Fomichev (1994) obtained zonal forcing for the same atmospheric region and empirical climatology using a similar radiative code within the traditional framework of computing the so-called diabatic circulation. In a previous paper (Akmaev, 1994) qualitative comparisons of their results with ours were presented and possible sources of apparent differences were discussed. It should be stressed here again that only the present diagnostic-simulation scheme uses objective quantitative criteria based upon which it is now possible to establish the representativeness of diagnostics. The results presented herein and by Akmaev (1994) confirm that the new diagnostic
Results of a series of optimization experiments with an assimilative middle atmospheric diagnostic-simulation scheme based on the use of a spectral model are presented. The scheme uses the annual global mean temperature deviation as an objective criterion of the ability of the simulations to reproduce empirical climatology. This quantitative criterion makes it possible to perform parameter optimization of the diagnostic scheme and/or the simulation model. It is demonstrated, in particular, that a proper adjustment of the strength of nudging forcing in the thermodynamic equation at the diagnostic stage can reduce the prognostic rms temperature deviation from the CIRA-86 empirical model to approximately 2.7 K in the 15-I 10 km layer. It is demonstrated that direct optimization procedures can effectively be used in atmospheric model parameter estimation problems of moderate dimensionality and may be recommended for parameter optimization of simulation models. The ‘optimal nudging’ data assimilation technique used in this study at the diagnostic stage is
R. A. Akmaev
992
apparently subject to fast saturation with increasing dimensionality of the problem, suggesting that a low-cost, low-dimension optimal nudging may yield results close to those obtained by a large-scale optimization. Acknowledgements-Maura Hagan and two anonymous reviewers offered many useful comments on early versions of the manuscript. The National Center for Atmospheric Research ence Foundation.
is sponsored
by the National
Sci-
REFERENCES Akmaev, R. A. (I 994) Diagnostics and simulation of an annual cycle in the middle atmosphere. 1 Geophys. Res. 99, 18,933-18,950. Akmaev, R. A., Fomichev, V. I., Gavrilov, N. M. and Shved, G. M. (1992) Simulation of the zonal mean climatology of the middle atmosphere with a threedimensional spectral model for solstice and equinox conditions, 1 Atmos. Terr Phys. 54, 119-128. Anthes, R. A. (1974) Data assimilation and initialization of hurricane prediction models. 1 Atmos. Sci. 31, 702-719. Apruzese, J. I?, Schoeberl, M. R. and Strobel, D. F. (1982) Parameterization of IR cooling in a middle atmosphere dynamics model. I. Effects on the zonally averaged circulation. 1 Geophys. Res. 87, 8951-8966. Courtier, I?, Derber, J., Errico, R., Louis, J.-F. and VukiCeviC, T. (1993) Important literature on the use of adjoint, variational methods and the Kalman filter in meteorology. Tellus &A, 342-357. Daley, R. (I 99 I) Atmospheric Data Analysis. Cambridge University Press, 457 pp. Dennis, J. E. and Schnabel, R. B. (1983) Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, Englewood Cliffs, NJ, 378 pp. Fleming, E. L., Chandra, S., Barnett, J. J. and Corney, M. (I 990) Zonal mean temperature, pressure, zonal wind
and geopotential height as functions of latitude. Adv. Space Res. 10, (12)l l-(12)59. Fletcher, R. (1980) Practical Methods of Optimization. John Wiley & Sons, Chichester, 120 pp. Ghil, M. and Malanotte-Rizzoli, P (1991) Data assimilation in meteorology and oceanography. Adv. Geophys. 33, 141-266. Gill, P E., Murray, W. and Wright, M. H. (1981) Pructical Optimization, Academic Press, London, 401 pp. Gottfried, B. S. and Weisman, J. (1973) Introduction to Optimization Theory Prentice-Hall, Englewood Cliffs, NJ, 571 pp. Hestenes, M. R. (1980) Conjugate Direction Methods in New York, 325 pp. Optimization. Springer-Verlag, Medvedev, A. S. and Fomichev, V. I. (I 994) Net radiative heating and diagnostics of the diabatic circulation in the 15-l IO km height layer. J Atmos. Terr Phys. 56, 1571-1584. Murgatroyd, R. J. and Singleton, F. (1961) Possible meridional circulations in the stratosphere and mesosphere. Quart. J Roy Met. Sot. 87, 125-135. Navon, I. M. and Legler, D. M. (1987) Conjugategradient methods for large-scale minimization in meteorology. Mon. tiea. Rev. 115, 1479-1502. Shine, K. P. (1989) Sources and sinks of zonal momentum in the middle atmosphere diagnosed using the diabatic circulation. Quart. J. Roy Met. Sot. 115,
265-292. Shine, K. P and Rickaby, J. A. (1989) Solar radiative heating due to absorption by ozone. In Ozone in the Atmosphere, eds R. D. Bojkov and P. Fabian, pp. 597-600. A. Deepak Publ., Hampton, VA. Stauffer, D. R. and Bao, J.-W. (1993) Optimal determination of nudging coefficients using the adjoint Tellus 45A, 358-369. equations. Strobel, D. F. (1978) Parameterization of the atmospheric heating rate from I5 to 120 km due to 02 and 03 absorption of solar radiation. J Geophys. Res. 83, 6225-6230. Zou, X., Navon, I. M. and Le Dimet, E X. (1992) An optimal nudging data assimilation scheme using parameter estimation. Quart. J Roy Met. Sot. 118, 1163-l 186.