Abstract

Four-dimensional variational data assimilation (4DVAR) is an established data assimilation method that finds the finite-amplitude perturbation that best fits observations consistent with a priori information and model dynamics. The response of a simulated tropical cyclone to specially designed finite perturbations of selected model variables was studied with a modified version of 4DVAR. The usual goal of minimizing data misfits was replaced with a goal of reducing damaging surface winds at the end of six hours of forecast time. For this purpose a property value cost function based on topography was defined. The case studied was a 20-km simulation of a hurricane approaching the Hawaiian Islands. Each prognostic variable in turn—temperature, winds, humidity, vertical velocity, and perturbation pressure—and all prognostic variables at once were used as the control vector for the optimization problem. Of all prognostic variables examined, temperature and the horizontal wind were the most effective at reducing damaging surface winds. The wind-only perturbation was very similar to the wind component of the perturbation calculated when all prognostic variables were used at once. Calculated perturbations had scales of 0.25°C or 1 m s−1, but changes at a few grid points near the center of the storm were an order of magnitude greater. Vertical velocity and humidity perturbations alone were ineffective at reducing damaging winds. The perturbation pressure experiment failed to converge but did substantially reduce the damaging winds.

1. Introduction

Theoretical and model studies have established that the dynamics governing the atmosphere can be extremely sensitive to small changes in initial conditions (e.g., Rabier et al. 1996). This suggests that the earth's atmosphere is chaotic. Chaos implies sensitivity to small perturbations. In a realistic numerical weather prediction (NWP) model, since small differences in initial conditions can grow exponentially, small but correctly chosen perturbations induce large changes in the evolution of the simulated weather. Current operational NWP practices—including data assimilation, generation of ensembles, and targeted observations—illustrate this daily (Hoffman 2002).

The Four-dimensional variational data assimiliation (4DVAR) is a data assimilation method that determines the optimal initial state for a weather forecast. It solves the Bayesian estimation problem of finding the most likely state given the background and observations, which is also consistent with our knowledge of the dynamics of the atmosphere. This method is now widely used in both operational and research applications for numerical weather prediction. For example, the operational European Centre for Medium-Range Weather Forecasts (ECMWF) global forecasts are based on 4DVAR data assimilation. In practice, after suitable approximations, 4DVAR becomes a nonlinear least square problem. In this problem the control vector describes the perturbation or analysis increment at the start of the 4DVAR interval and a nonlinear model evolution provides the analysis at other times that is compared to observations. Essentially the background error covariances are used in the background cost function to measure the size of the analysis increment, while the sum of the error covariances due to model error, representativeness error, and observing system error are used in an observation cost function to measure the data mismatch during the time interval considered. (When a covariance is used to measure a vector, the process can be explained in this way: The vector is rotated to the coordinate system defined by the eigenvectors of the covariance matrix, and then scaled component-wise by the square root of the associated eigenvalues. The squares of the scaled components, which have expected value one and are uncorrelated, are then summed.) In solving the nonlinear optimization problem, 4DVAR iteratively uses an adjoint model linearized about the current nonlinear model simulation. When the process has converged sufficiently, the solution model state is used to initialize a nonlinear model simulation that extends past the end of the 4DVAR time interval.

The 4DVAR problem can be extended in a variety of interesting ways by adding a cost function that measures an alternative aspect of the forecast weather different from the mismatch between simulated and measured observations. One use is for numerical weather control experiments. A second use is to make forecasts of a specified extreme event consistent with the initial state uncertainty. We have added two different types of cost functions to the 4DVAR motivated by our interests in weather control. The first type of cost function, used by Hoffman et al. (2002), measures the difference between the simulated atmospheric state at the end of the 4DVAR interval and a target or desired atmospheric state. The second type of cost function, used here, measures the property damage due to damaging winds. This cost function is defined below in section 2b. In the same way that normal modes or singular vectors have precise structures, a precise structure is also found when we use 4DVAR with one of these cost functions. In fact the calculations presented here determine the leading nonlinear singular vector for an unusual form of the objective function (i.e., for an unusual form of the quantity to be optimized). Here we minimize wind damage, while in “ordinary” linear singular vector calculations the objective might be to maximize the perturbation energy at 24 h.

From the numerical weather control point of view, perturbations determined by 4DVAR with these new cost functions effectively control the evolution of the simulated atmosphere. In this case the size of a perturbation should be measured by its implementation cost. This cost could include the cost of fuel or other expendable, as well as the difference between cash payments made to those harmed and received from those benefited. Note that it may not be cost-effective or even feasible to directly alter some atmospheric quantities. For example, it might be possible to change the temperature of the atmosphere but not the humidity. In the real world the cost function must include many consequences. For a hurricane the costs of high wind and heavy precipitation and flooding on life and property should all contribute. Minimizing just one of these costs might greatly increase one of the others. Furthermore, uncertainty in all aspects of the weather control system—initial state estimate, forecast model, implementation of perturbations—must be quantified to reduce the possibility of obtaining a negative result. From the extreme event point of view the perturbations calculated must be plausible. In this case the size of the perturbation should be measured by the usual background error covariance term. This term is a measure of the likelihood of the perturbation given our estimate of the uncertainty of the current state of the atmosphere. In the experiments reported here we use an intermediate method of measuring the size of the perturbations that is based on the initial tendencies following Zou et al. (1997).

The present set of experiments attempts to determine which variables are most important for the future evolution of a hurricane's damaging winds as it approaches landfall. Similar experiments could be conducted for other meteorological phenomena, but hurricanes are of particular interest for several reasons. Hurricanes are very difficult to forecast, in part because hurricane forecasts are highly sensitive to specified initial and boundary conditions. Therefore relatively small finite-amplitude perturbations can have large effects. (The other difficulty of course is lack of good observations, both for initial conditions and to develop parameterizations of physical processes appropriate for hurricane conditions.) Hurricanes are also very costly in terms of life and infrastructure, both at sea and when making landfall. In addition, the cost of evacuation in affected communities is also large.

In experiments reported here we calculate perturbations to minimize a cost function that estimates property damage due to surface wind in a simulation of central Pacific Hurricane Iniki of 1992. All calculations make use of the complete nonlinear MM5 model including moist physics, except that the adjoint MM5 model that includes linearized moist physics is used to calculate the gradient of the objective function for use by the minimization algorithm. Hurricanes are very nonlinear phenomena and, while 4DVAR solves the fully nonlinear problem, the use of linearization within the iterative minimization process can result in slow convergence or even nonconvergence. Hurricane Iniki caused extensive damage to property and vegetation on parts of the Hawaiian Islands and killed six people (Chun et al. 1993; Lawrence and Rappaport 1994). The storm made landfall on Kauai at 0130 UTC 12 September 1992, with a central pressure of 945 hPa. Maximum sustained winds over land were estimated at 60 m s−1 with gusts as high as 80 m s−1. Our unperturbed simulation is much weaker due to the limitations of the fifth generation version of the Pennsylvania State University–National Center for Atmospheric Research (NCAR) Mesoscale Model (MM5), which we used, and the initial conditions. Nevertheless, damaging winds are significantly reduced in controlled simulations using relatively small perturbations to all variables together or to just the horizontal wind or temperature alone.

2. The modified 4DVAR method

The implementation of 4DVAR used in this study is described by Zou et al. (1997). It has been applied to assimilate such nonconventional meteorological variables as zenith delay observations from global positioning system (GPS) satellites (De Pondeca and Zou 2001) and cloud-cleared brightness temperatures from geostationary operational environmental satellite (GOES) sounders (Zou et al. 2001).

a. Initial perturbation cost function

To mathematically define the objective function that will be minimized by 4DVAR, we first define the unperturbed simulation U, from time 0 to T, with corresponding states U(0) and U(T). We then use 4DVAR to find an optimal controlled simulation C that simultaneously minimizes the difference from the initial state [i.e., C(0) − U(0)] and the estimated wind damage later in the forecast period. In the experiments reported the model is assumed perfect in the sense that the unperturbed simulation is taken as reality. These experiments by themselves leave open the possibility that response of the atmosphere to the calculated perturbations might be very different than the response of the model. As a partial answer, in similar experiments for Hurricane Andrew (Henderson et al. 2005) perturbations calculated at 20-km resolution with simple physics were successful at limiting wind damage when applied in a 6.67-km simulation that used more complex physics.

In these preliminary experiments, the size of the initial perturbation is measured by a simple quadratic norm:

 
formula

Here Jb is the background cost functional, x defines the control vector variables (e.g., the temperature, or the horizontal wind components, or all variables), and i, j, and k index the grid points in the three spatial dimensions. In Eq. (1), the flux or coupled form of the variables is used since this is the form of the primitive equations in the MM5. For example, p*u is the coupled eastward wind component, where p* is the reference pressure difference between the bottom and top model boundaries. The reference state varies in the vertical only, therefore p* depends only on the model surface topography. Below, in Fig. 6, we present the components of Jb for temperature at different model layers as the square root of the terms in square brackets in Eq. (1), normalized by the number of grid points, and dimensionalized assuming p* = 950 hPa, the value over the ocean.

Fig. 6.

Profiles of rms difference of temperature at 0 h with respect to the unperturbed simulation 𝒰 for experiments 𝒞[T] (dotted line) and 𝒞[X] (solid line). The components of Jb are converted into rms differences as described in the text. For visualization the layer rms differences are connected by lines.

Fig. 6.

Profiles of rms difference of temperature at 0 h with respect to the unperturbed simulation 𝒰 for experiments 𝒞[T] (dotted line) and 𝒞[X] (solid line). The components of Jb are converted into rms differences as described in the text. For visualization the layer rms differences are connected by lines.

The scales Sxk depend only on variable and layer. In the present experiments, following Zou et al. (2001), Sxk is calculated as the maximum absolute difference between two model states for each variable at each layer, in this case between U(0) and U(δt) with δt = 40 × 60 s. Figure 1 shows the vertical profiles of the Sxk used in our experiments, again dimensionalized assuming p* = 950 hPa. In general, these scales vary smoothly in the vertical. Because the storm is moving northward the Sxk for u is greater than for υ. As mentioned in the introduction, to determine a feasible extreme model forecast, an estimate of the initial error covariance should be used to define Jb as is done in operational 4DVAR applications. Note that if flow-dependent background error statistics were used, this would allow relatively large analysis increments near the center of a hurricane because of the larger uncertainty there. The Sxk correspond to a diagonal background error covariance matrix and thus provide only a rough normalization to equalize the contributions of variables of different quantities (e.g., temperature versus humidity). Also, for weather control, the control variables and cost function defining the size of the perturbation do not depend on the atmospheric state, but on the parameters describing the perturbations to the system. For example, if orbiting solar reflectors were the control mechanism, then the control vector would describe the orientation of the reflectors as a function of time and the costs would measure the use of expendables. At best the Sxk may be considered a very crude estimate of the relative costs of introducing perturbations at different levels or in different variables.

Fig. 1.

Profiles of scaling factors for temperature (T), eastward wind component (u), and northward wind component (υ) plotted using solid, dotted, and dashed lines, respectively. Values are calculated from the coupled variables and then dimensionalized for this plot assuming p* = 950 hPa. The horizontal scale is °C for temperature and m s−1 for the wind components.

Fig. 1.

Profiles of scaling factors for temperature (T), eastward wind component (u), and northward wind component (υ) plotted using solid, dotted, and dashed lines, respectively. Values are calculated from the coupled variables and then dimensionalized for this plot assuming p* = 950 hPa. The horizontal scale is °C for temperature and m s−1 for the wind components.

b. Damage cost function

In our experiments the total cost function is defined as

 
formula

where the subscript d stands for damage. The damage cost function,

 
formula

is written in terms of physical damage estimates based on a relationship between surface wind speeds and economic damage. Here D denotes damage and λ is a weighting factor that is used to equalize the contribution of the damage cost function to the total cost function. The damage at each grid point is the product of the fractional wind damage and the property value (Pij). The fractional damage (Unanwa et al. 2000) depends upon two threshold wind speeds: the lower threshold (U0) is the wind speed at which damage to property first occurs, while the second (U1) is the wind speed at which complete destruction occurs. Between these two threshold values, we model the increase in damage using a cosine curve. Thus,

 
formula

where U is the simulated horizontal surface wind speed. The property value P is unitless and varies with location. In the experiments described here, Eqs. (3) and (4) are evaluated only at time t = 6 h (i.e., at the end of the 4DVAR interval). In all cases U0 = 25 m s−1 and U1 = 90 m s−1 and the surface wind is taken to be the lowest model layer wind for this calculation. In general, the results are sensitive to the choice of λ, but experimentation showed that large changes in this parameter are needed to obtain a clearly visible effect. Also a large value is needed since only a very small fraction of model grid points experience wind damage, while most grid points contribute to the background term. The particular value used of λ = 400 000 was found to provide solutions with reduced wind damage while retaining all of the basic characteristics of the original simulated storm at the initial time.

3. Experimental procedures

We use the MM5 system to simulate a hurricane based on the meteorological situation observed during the time period leading up to the landfall of Hurricane Iniki on the Hawaiian island Kauai in 1992. In our experiments t = 0, the beginning of the 4DVAR interval, corresponds to 0600 UTC 11 September 1992. The current prototype experiments use a relatively coarse, 20-km resolution. As a result and due to the simple parameterized physics used here, our simulations are only crude representations of observed track and intensity. However we note that MM5 produces very detailed and accurate simulations of tropical cyclones when high-resolution and advanced physical parameterizations are used (e.g., Liu et al. 1999).

a. Mesoscale model

The MM5 used in our experiments is described by Grell et al. (1994) and by Dudhia (1993). Here, the MM5 computational grid covers an approximately 3000 × 4000 km2 horizontal domain with ten sigma layers in the vertical from the surface to 50 hPa. Hurricane Iniki remains far enough from the domain edges that boundary effects are small during the course of the experiments (with the possible exception of the experiment using perturbation pressure as the control variable). The sigma vertical coordinate is a terrain-following normalized pressure coordinate (Holton 1992). All experiments described here use nonhydrostatic dynamics, a 60-s time step and a 20-km Mercator grid. The parameterizations of surface fluxes, radiative transfer, and cumulus convection are currently limited to relatively simple schemes in the MM5 4DVAR system. For our experiments, the physical parameterizations include the Medium-Range Forecast Model (MRF) planetary boundary layer (PBL) and the Anthes Kuo convection scheme. Large-scale stable (i.e., nonconvective) precipitation occurs whenever a layer reaches saturation. Excess moisture rains out immediately with no reevaporation as it falls. The longwave radiation computation of simple cooling, including cloud effects, is applied every 30 time steps.

b. Hurricane initialization in the gridded datasets

Gridded National Centers for Environmental Prediction (NCEP) reanalysis fields (Kalnay et al. 1996) from the NCAR archives, and operational NCEP sea surface temperature (SST) analyses, were used to initialize the MM5 model and provide boundary conditions during the 4DVAR and forecast periods. The boundary conditions define both the model state along the lateral boundaries and surface parameters, such as SST. Other surface parameters were derived from the databases included in the MM5 distribution. Figure 2 shows the SST over the entire computational domain. The very warm surface waters south of Hawaii provide the energy needed for the formation and intensification of tropical cyclones. While the model grid is 158 × 194 for the experiments reported here, the smaller domain of 92 × 95 grid points centered over the storm is used in all subsequent map figures except Fig. 4.

Fig. 2.

Computation domain showing sea surface temperature at the start of the 4DVAR interval, 0600 UTC 11 Sep 1992. The contour interval is 1°C. The inner rectangle indicates the area shown in all other map figures except Fig. 4.

Fig. 2.

Computation domain showing sea surface temperature at the start of the 4DVAR interval, 0600 UTC 11 Sep 1992. The contour interval is 1°C. The inner rectangle indicates the area shown in all other map figures except Fig. 4.

Fig. 4.

Property values based on smoothed topography. Property values are unitless.

Fig. 4.

Property values based on smoothed topography. Property values are unitless.

The NCEP reanalyses have only a hint of the actual strength and structure of Hurricane Iniki. In the future it may be possible to use improved observations and data assimilation to accurately initialize a mesoscale hurricane forecast model. For the present experiments we added an analytic representation of a hurricane vortex using the method of Davis and Low-Nam (2001) 6 h before the start of the 4DVAR minimization interval (i.e., at t = −6). The Davis and Low-Nam method is part of the MM5 system. The bogus storm is axisymmetric and is based on specifying the storm position and the radius and magnitude of the maximum wind in the lowest model layer. Given these parameters an analytical Rankine wind vortex is added to the existing large-scale wind field to provide the bogus wind field. From the bogus wind field a temperature field is calculated to be in nonlinear gradient balance at all levels. After the bogus procedure, we let the model representation of the hurricane equilibrate during the 6 h leading up to the start of the 4DVAR interval at 0600 UTC 11 September 1992 (Fig. 3). Even after the bogus procedure and the 6-h forecast, however, the simulated storm is too weak compared to the best track estimate valid at this time of almost 60 m s−1.

Fig. 3.

Lowest model layer wind field (hereafter the surface wind field) at the start of the 4DVAR interval (i.e., t = 0 at 0600 UTC 11 Sep 1992) for Hurricane Iniki 6 h after application of the bogus procedure. Wind speed thresholds that correspond to increasing Saffir–Simpson category are contoured with increasingly thick line widths. The plotted contours correspond to the lowest wind speeds associated with tropical depressions (12 m s−1), tropical storms (17 m s−1), category 1 hurricanes (33 m s−1), and category 2 hurricanes (43 m s−1). Shading denotes the extent of wind speeds above 25 m s−1—the lower threshold of damaging winds. Wind symbols are in m s−1 with a full barb indicating 5 m s−1.

Fig. 3.

Lowest model layer wind field (hereafter the surface wind field) at the start of the 4DVAR interval (i.e., t = 0 at 0600 UTC 11 Sep 1992) for Hurricane Iniki 6 h after application of the bogus procedure. Wind speed thresholds that correspond to increasing Saffir–Simpson category are contoured with increasingly thick line widths. The plotted contours correspond to the lowest wind speeds associated with tropical depressions (12 m s−1), tropical storms (17 m s−1), category 1 hurricanes (33 m s−1), and category 2 hurricanes (43 m s−1). Shading denotes the extent of wind speeds above 25 m s−1—the lower threshold of damaging winds. Wind symbols are in m s−1 with a full barb indicating 5 m s−1.

c. Generation of property values

The two-dimensional property value field (Fig. 4) was generated by smoothing topography. The smoother averages all points within 400 km resulting in a gradient of property values for nearshore water points. This approach is taken since we found that localizing property values only over the islands permitted damaging winds directly offshore. Also, spreading out property values results in a smooth, well-behaved cost function. Property values are unitless, but might be considered tens of thousands of dollars per square kilometer.

4. Results

a. Summary of experiments

The experiments described here seek a perturbation that minimizes damaging surface winds at the end of the 6-h 4DVAR interval at t = 6 h, or 1200 UTC 11 September 1992. In the first experiment, denoted 𝒞[T], the control vector is the field of coupled temperature p*T only. In general, the control vector is a list of all the quantities that may be changed directly by the minimization. That is, we minimize J with respect to the control vector. The experiment name 𝒞[T] refers to both the 4DVAR analysis and to the subsequent fully nonlinear controlled forecast that evolves from the analyses. The forecast using unperturbed initial conditions is denoted 𝒰. Here 𝒞 stands for controlled and 𝒰 for unperturbed. In addition to 𝒞[T], several other experiments with different control vectors are described and key results are reported. In experiment 𝒞[X] the control vector is the entire model state. Analogs to 𝒞[T] with variations on the control vector are experiments 𝒞[V], 𝒞[w], 𝒞[q], and 𝒞[p′]—in which the control vector is composed of the coupled horizontal wind components p*u and p*υ, vertical velocity p*w, specific humidity p*q, and perturbation pressure p*p′, respectively. A maximum of 10 iterations were allowed in the conjugate gradient algorithm used in the 4DVAR minimization. To restrict the control vector to a subset of X, each calculated gradient is projected onto that subset. For example, in experiment 𝒞[T] the complete gradient is calculated at each iteration, but then all entries not associated with the temperature field are set to zero before the gradient is used in the minimization procedure. The minimization ended early for the 𝒞[w] and 𝒞[q] cases because the minimizer could not determine a clearly defined direction in which to proceed after five or six iterations.

b. Baseline experiment: 𝒞 [T]

In Table 1 we report the final values of Jd, the components of Jb, and J for all experiments. As the minimum is approached, J should decrease and the gradient should approach zero. Note that the initial values of J and Jd are 80 727 in all cases and the initial values for Jb are all zero. In the table we see that for 𝒞[T] the value of J is reduced to 5% of its initial value after 10 iterations of the minimization algorithm. Figure 5 shows the cost function versus iteration for experiment 𝒞[T]. The asymptotic leveling off of the total and component cost functions indicates convergence of the minimization, and the much lower final value of J indicates success at finding a small perturbation that reduces damaging surface winds.

Table 1.

Hurricane Iniki experiments and cost functions. The experiment names in the first column indicate the control vector in use in the topography-based wind damage cost function Iniki experiments. Columns 2 and 3 give the number of iterations and function evaluations required by the optimization. Columns 4 and 5 list the final total cost function J and the ratio of the final to initial value expressed as a percent. The initial total cost function is 80 727 in all experiments. Next, the final values of the components of the cost function due to the perturbation at the initial time are listed for each model variable. The final value of Jd follows. The final three columns list the magnitude of the gradient of the cost function at the start and end of the minimization along with the ratio of the final to initial value expressed as a percent.

Hurricane Iniki experiments and cost functions. The experiment names in the first column indicate the control vector in use in the topography-based wind damage cost function Iniki experiments. Columns 2 and 3 give the number of iterations and function evaluations required by the optimization. Columns 4 and 5 list the final total cost function J and the ratio of the final to initial value expressed as a percent. The initial total cost function is 80 727 in all experiments. Next, the final values of the components of the cost function due to the perturbation at the initial time are listed for each model variable. The final value of Jd follows. The final three columns list the magnitude of the gradient of the cost function at the start and end of the minimization along with the ratio of the final to initial value expressed as a percent.
Hurricane Iniki experiments and cost functions. The experiment names in the first column indicate the control vector in use in the topography-based wind damage cost function Iniki experiments. Columns 2 and 3 give the number of iterations and function evaluations required by the optimization. Columns 4 and 5 list the final total cost function J and the ratio of the final to initial value expressed as a percent. The initial total cost function is 80 727 in all experiments. Next, the final values of the components of the cost function due to the perturbation at the initial time are listed for each model variable. The final value of Jd follows. The final three columns list the magnitude of the gradient of the cost function at the start and end of the minimization along with the ratio of the final to initial value expressed as a percent.
Fig. 5.

Cost function vs iteration for experiment 𝒞[T]. The total cost function J (in thousands) and the individual parts of the cost function for Jb at 0 h and Jd at 6 h are shown as solid, dotted, and dashed lines, respectively.

Fig. 5.

Cost function vs iteration for experiment 𝒞[T]. The total cost function J (in thousands) and the individual parts of the cost function for Jb at 0 h and Jd at 6 h are shown as solid, dotted, and dashed lines, respectively.

Figure 6 shows the profiles of rms perturbation temperature for experiments 𝒞[T] and 𝒞[X]. The components of J are converted into rms differences as described in the discussion of Eq. (1). The 𝒞[T] temperature perturbations are larger in magnitude in the upper troposphere, where perturbations are typically 0.25°C in magnitude. The 𝒞[X] temperature increments are much smaller since in this case the desired effect is obtained principally by altering the wind field. (See the discussion in the next section.)

Although the perturbation magnitudes are small, the domain is large and the energies involved are huge. Henderson et al. (2005) discuss the petajoule energy requirements for controlling a simulation of Hurricane Andrew. Here we estimate the perturbation energy directly from the rms perturbation profiles, such as those shown in Fig. 6, as the integral over the entire mass of the atmosphere of the squares of the perturbation quantities times appropriate constants (Ehrendorfer and Errico 1995; Thépaut et al. 1993). For example, for wind the perturbation energy is given by one-half the mean square perturbation wind speed times the mass of the atmosphere in the model domain. (For temperature, a similar formula is used in which the mean squared perturbation temperature is multiplied by cp/Tr, the ratio of the specific heat for dry air to the reference temperature.) Thus we find for experiment 𝒞[V] the perturbation energy is 20.9 × 1015 J, while for experiment 𝒞[T] the perturbation energy is 17.5 × 1015 J. For comparison in experiment 𝒞[X] the wind perturbation energy is 15.1 × 1015 J while the temperature perturbation energy is 0.1 × 1015 J. As discussed by Ehrendorfer and Errico (1995) there are many caveats associated with the definition of perturbation energy used here. In addition, our estimates ignore variations in p*, map factor, and deviations of density from the reference density. However the effect of these variations is small compared to the effect of the assumed reference temperature, pressure, and specific humidity used in the definition of the perturbation energy.

Figure 7 shows the structure of the perturbation temperature at σ = 0.350, 0.650, and 0.950 for experiment 𝒞[T]. Over the ocean these layers correspond to 382.5, 667.5, and 952.5 hPa. At σ = 0.950, the lowest model level, there is cooling close to the center of the storm and heating to the west. At mid and upper levels near the center there is a complex pattern of stronger heating and cooling. These patterns are not correlated between σ = 0.650 and σ = 0.350. Patterns in the intervening layers show that these features twist and amplify with altitude. At σ = 0.350, a large cold temperature increment is present directly over the center of Iniki. Cold temperature increments over the eye are in direct opposition to the warm core thermodynamic structure of the hurricane and act to destroy the hurricane in place. Although the largest magnitude changes are near the storm center, there are small changes required throughout the rest of the domain. These slowly decay as distance from the storm center increases. In particular, there is evidence of a concentric banded structure away from the hurricane center, the amplitude of which increases with altitude. A concentric banded structure is found to varying degrees in all of our experiments. In preliminary experiments at 40-km resolution a similar structure was found but with larger horizontal scales. These concentric patterns seem to be fundamental to the solution. In an experiment similar to 𝒞[T], but where perturbations were not allowed near the center of the storm, the solution found had only this banded structure but of much larger magnitude in those regions where perturbations were allowed.

Fig. 7.

Structure of the initial perturbation for experiment 𝒞[T]. Horizontal slices of δT are shown at σ = 0.350, 0.650, and 0.950 (from top to bottom) at 0 h. Negative perturbations are shaded, and contours indicate temperature magnitudes of differences exceeding 1° (thin lines) and 3°C (thick lines). For reference, the crosshairs, here and in succeeding perturbation plots, indicate the positions of the storm in the 𝒰 simulation at the appropriate times.

Fig. 7.

Structure of the initial perturbation for experiment 𝒞[T]. Horizontal slices of δT are shown at σ = 0.350, 0.650, and 0.950 (from top to bottom) at 0 h. Negative perturbations are shaded, and contours indicate temperature magnitudes of differences exceeding 1° (thin lines) and 3°C (thick lines). For reference, the crosshairs, here and in succeeding perturbation plots, indicate the positions of the storm in the 𝒰 simulation at the appropriate times.

Figure 8 shows the time evolution of the perturbation for experiment 𝒞[T] calculated as the difference between the two nonlinear simulations, 𝒞[T] minus 𝒰. Horizontal slices of δT, δw, and δV at selected levels are shown at 2-h intervals. Note the largest changes in all three fields are close to the hurricane center at each time. This finding is exaggerated at 6 h as the required reduction in surface wind speed is realized. Note that the perturbations grow quickly only until the end of the 4DVAR interval. Figure 9 shows the evolution of the lowest model layer wind field (hereafter referred to as the surface wind field) for the unperturbed simulation 𝒰 and for experiment 𝒞[T]. The effect of the perturbations on the full wind fields is to effectively suppress the winds to near or below the critical damaging wind level of 25 m s−1 at 6 h, and at 6 h only. The areal extent of the damaging winds (>25 m s−1, shaded) at 6 h has been reduced. The maximum wind speeds are much lower at 6 h, and to a lesser extent at 8 h in experiment 𝒞[T]. Figure 9 also shows that in the controlled case, the hurricane has slightly accelerated along its track relative to the unperturbed case.

Fig. 8.

Evolution of the perturbation for experiment 𝒞[T]. Horizontal slices are shown (from left to right) for δT at σ = 0.350, δw at σ = 0.650, and δV at σ = 0.950 each at (from top to bottom) 2, 4, 6, and 8 h. Negative perturbations are shaded. Thin and thick contours indicate magnitudes of differences greater than 1° and 3°C for temperature, 30 and 60 × 10−2 m s−1 for vertical velocity, and 3 and 10 m s−1 for wind speed. Wind symbols are in m s−1 with a full barb indicating 5 m s−1

Fig. 8.

Evolution of the perturbation for experiment 𝒞[T]. Horizontal slices are shown (from left to right) for δT at σ = 0.350, δw at σ = 0.650, and δV at σ = 0.950 each at (from top to bottom) 2, 4, 6, and 8 h. Negative perturbations are shaded. Thin and thick contours indicate magnitudes of differences greater than 1° and 3°C for temperature, 30 and 60 × 10−2 m s−1 for vertical velocity, and 3 and 10 m s−1 for wind speed. Wind symbols are in m s−1 with a full barb indicating 5 m s−1

Fig. 9.

The surface wind field for the (left) unperturbed simulation 𝒰 and (right) experiment 𝒞[T] at (from top to bottom) 4, 6, and 8 h. Plotted as in Fig. 3. Note that the difference of (b) − (a) is plotted in Figs. 8f; (d) − (c) in Fig. 8i; and (f) − (e) in Fig. 8l.

Fig. 9.

The surface wind field for the (left) unperturbed simulation 𝒰 and (right) experiment 𝒞[T] at (from top to bottom) 4, 6, and 8 h. Plotted as in Fig. 3. Note that the difference of (b) − (a) is plotted in Figs. 8f; (d) − (c) in Fig. 8i; and (f) − (e) in Fig. 8l.

c. Parameter sensitivity experiments

Figure 10 shows the surface wind field at 6 h for the 𝒞[X] experiment. As seen in the small extent of the shaded area in Fig. 10, damaging winds were reduced the most in this experiment. As in all of the successful analyses (i.e., in all cases in which the minimization converged), the extent and intensity of damaging winds (winds >25 m s−1) were drastically reduced in the hours near the evaluation time of the wind damage cost function. It should be noted that the extent and intensity of the winds increased rapidly in the hours following this time in these experiments. Also, in experiments 𝒞[X] and 𝒞[V] the storm remains nearly stationary from 6 to 9 h. This was followed by a northwesterly motion of the simulated vortex that resulted in the storm passing the latitude of Hawaii much later than in the unperturbed case and much farther west. In some experiments, the sea level pressure field (which is related to the overall temperature through the depth of the atmosphere) appeared to be temporarily mispositioned with respect to the wind field (compared to gradient balance) close to the time of the wind damage cost function evaluation.

Fig. 10.

Surface wind field at 6 h for experiment 𝒞[X]. Plotted as in Fig. 3.

Fig. 10.

Surface wind field at 6 h for experiment 𝒞[X]. Plotted as in Fig. 3.

Figure 11 shows the vertical structure of the wind perturbation δV for experiments 𝒞[X]. The largest calculated changes are at σ = 0.650. At all levels the changes are most pronounced within a radius of ∼250 km from the center of the storm. Figure 12 shows the pattern of perturbation temperature δT at σ = 0.350 for experiment 𝒞[X], which should be compared to the top panel of Fig. 7. Note that the 𝒞[T] temperature increments are 2–4 times larger than those of 𝒞[X]. Plots at other levels are qualitatively similar.

Fig. 11.

Structure of the initial wind perturbation δV for experiment 𝒞[X]. Horizontal slices are shown (from top to bottom) at σ = 0.350, 0.650, and 0.950 at 0 h. Plotted as in Fig. 8c.

Fig. 11.

Structure of the initial wind perturbation δV for experiment 𝒞[X]. Horizontal slices are shown (from top to bottom) at σ = 0.350, 0.650, and 0.950 at 0 h. Plotted as in Fig. 8c.

Fig. 12.

Initial perturbation temperature δT for experiment 𝒞[X]. Plotted at σ = 0.350 at 0 h as in Fig. 7.

Fig. 12.

Initial perturbation temperature δT for experiment 𝒞[X]. Plotted at σ = 0.350 at 0 h as in Fig. 7.

Based on the experiment using X as the control vector we expect that wind is the most effective control variable. To explore this further we conducted additional sensitivity experiments in which the control vector is restricted to each of V, w, q, and p′ in turn. The huge reduction in J seen in Table 1 when the control vector is either V or X shows that 4DVAR is more successful in preventing wind damage at 6 h when the horizontal velocity field is allowed to vary. It should be noted that the choice of scales (i.e., Sxk) affect the numeric results in Table 1. However in experiments like 𝒞[T] in which the control vector comprises a single variable, the relative magnitudes of variables with different dimensions is immaterial since only one of these (e.g., temperature) enters the calculation of Jb. Figure 13 shows the perturbations for these parameter sensitivity experiments at key levels at the start of the 4DVAR interval t = 0 h. The largest perturbations in each case are near the center of Iniki. Increments of smaller magnitude with a concentric pattern of alternating signs are found at larger distances from the storm's center. That the larger increments occur at the center of the storm is expected since hurricanes are largely sustained by physical processes in the storm's eyewall. The concentric patterns appear wavelike, propagating both inward and outward with respect to the hurricane center. The presence of strong horizontal wind increments near the center of the storm that are similar in size and form in experiments 𝒞[X] and 𝒞[V] is consistent with the fact that 4DVAR reduces wind damage most efficiently via changes in the initial horizontal wind field. (Compare the panels for V at σ = 0.950 in Figs. 11 and 13.)

Fig. 13.

Initial perturbations for experiments with different control vectors. (a) Wind at σ = 0.950 is plotted for experiment 𝒞[V], (b) vertical velocity at σ = 0.650 for 𝒞[w], (c) specific humidity at σ = 0.650 for 𝒞[q], and (d) perturbation pressure at σ = 0.950 for 𝒞[p′]: all at 0 h. Negative perturbations are shaded. Thin and thick contours indicate magnitudes of differences greater than 3 and 10 m s−1 for wind speed, 30 and 60 × 10−2 m s−1 for vertical velocity, 5 and 10 g kg−1 for specific humidity, and 5 and 10 hPa for perturbation pressure. Wind symbols are in m s−1 with a full barb indicating 5 m s−1.

Fig. 13.

Initial perturbations for experiments with different control vectors. (a) Wind at σ = 0.950 is plotted for experiment 𝒞[V], (b) vertical velocity at σ = 0.650 for 𝒞[w], (c) specific humidity at σ = 0.650 for 𝒞[q], and (d) perturbation pressure at σ = 0.950 for 𝒞[p′]: all at 0 h. Negative perturbations are shaded. Thin and thick contours indicate magnitudes of differences greater than 3 and 10 m s−1 for wind speed, 30 and 60 × 10−2 m s−1 for vertical velocity, 5 and 10 g kg−1 for specific humidity, and 5 and 10 hPa for perturbation pressure. Wind symbols are in m s−1 with a full barb indicating 5 m s−1.

Figure 14 shows the surface wind field at 6 h for the sensitivity experiments. Because 4DVAR minimizes an integrated quantity (Jd), different patterns of wind speed are expected at the end of the 4DVAR interval for different control vectors. A comparison with the surface wind speed field in the unperturbed simulation (Fig. 9) shows that, in general, in the controlled simulations the position of Iniki is farther north and west and that damaging winds are substantially reduced at 6 h, although less so in 𝒞[w] and more so in 𝒞[V]. In experiment 𝒞[p′] it is noteworthy that due to initial imbalances gravity waves advance from the lateral boundaries at 70 m s−1 and interact with the inner environment of the storm at 6 h (not shown). Apparently this is the most efficient way for 4DVAR to weaken the storm under the constraints of experiment 𝒞[p′]. In experiment 𝒞[q] we note that there are isolated regions of supersaturation (relative humidity greater than 100%) early in the simulation that are a consequence of large positive perturbations close to the center of the storm and over the large island of Hawaii (not shown).

Fig. 14.

Surface wind field at 6 h for experiments (a) 𝒞[V], (b) 𝒞[w], (c) 𝒞[q], and (d) 𝒞[p′]. Plotted as in Fig. 3.

Fig. 14.

Surface wind field at 6 h for experiments (a) 𝒞[V], (b) 𝒞[w], (c) 𝒞[q], and (d) 𝒞[p′]. Plotted as in Fig. 3.

5. Summary and concluding remarks

We use the 4DVAR system developed for the fifth- generation Pennsylvania State University–NCAR Mesoscale Model (MM5) to determine which variables are most important for the future evolution of a hurricane's damaging winds as it approaches landfall. We find that relatively large changes to humidity or vertical velocity do not have a large influence on damaging winds at the end of the 4DVAR period. Moderate changes to temperature and relatively small changes to the horizontal winds can have a substantial effect. These results are consistent with the high value placed on wind observations in conventional applications of 4DVAR in the Tropics. These experiments should be considered prototypes: our simulations of hurricanes could be improved with more sophisticated physical parameterizations and higher resolution. Further, in the context of controlling the weather, we introduce perturbations to the model atmosphere as instantaneous changes. In a more realistic simulation the vector of control parameters that is optimized might be a description of the temporal and spatial patterns of feasible forcing. For example, these parameters might describe the pattern and strength of additional heating supplied to the atmosphere by a space solar power downlink in the 183-GHz water vapor spectral region (Hoffman et al. 2005). To do this one would calculate heating rate profiles for the current state of the atmosphere including the effect of cloud and then calculate the optimal timing and amplitudes of the allowable heating. Of course, the optimal heating might not be precisely attainable, and experiments should be conducted to examine the sensitivity of the result to variations between the perturbation calculated and that actually applied. In spite of the simplifications made, our experiments demonstrate that changes to atmospheric temperature or horizontal winds alone are sufficient to control simulated hurricanes. Model resolutions allowing more realistic simulations should be used in future experiments. Preliminary experiments at higher (6.67 km) and lower (40 km) resolutions suggest that our results are robust and also that the magnitudes of the calculated perturbations decrease with increasing resolution. However, it is possible that entirely different behavior will be found at the even higher resolution that is required for more realistic hurricane forecasting and that resolves internal hurricane dynamics.

Clearly it will be a long time before it is possible to control a tropical cyclone in reality, but our approach may be useful in the near term in assessing extreme, though possible, forecast outcomes for various interests. While the present experiment minimizes property damage, we can also estimate the worst case scenario if we simply multiply the property damage cost function by −1, so as to instead maximize property damage. Such a worst case forecast might help planning evacuations. For different user communities the most critical aspect of a weather forecast will generally be different, but in many cases it may be possible to quantify weather effects on profit and loss. The method described here could then be modified to calculate the best and worst case scenarios under different assumptions of how resources are deployed. For example, the cost of buying electric power on the spot market can be very much higher than buying power a few days ahead of time. An electric utility might therefore make good use of the worst case weather forecast (i.e., worst case in terms of maximum cooling degree days on a hot summer day). Experiments along these lines are currently being planned.

Acknowledgments

The NASA Institute for Advanced Concepts supported this research. NCEP analyses were obtained from NCAR. The series of experiments described here are based on a suggestion of M. Halem of Goddard Space Flight Center.

REFERENCES

REFERENCES
Chun
,
A. K. T.
,
R. T.
Martin
,
H. E.
Rosendal
, and
G. H.
Trapp
,
1993
:
1992 tropical cyclones—Central North Pacific. NOAA Tech. Memo. NWSTM PR-38, Central Pacific Hurricane Center, Honolulu, HI, 74 pp. [Available online at http://www.prh.noaa.gov/cphc/summaries/1992.php.]
.
Davis
,
C.
, and
S.
Low-Nam
,
2001
:
The NCAR–AFWA tropical cyclone bogussing scheme. Tech. Memo., Air Force Weather Agency (AFWA), Omaha, NE, 13 pp. [Available online at http://www.mmm.ucar.edu/mm5/mm5v3/tc-report.pdf.]
.
De Pondeca
,
M. F. V.
, and
X.
Zou
,
2001
:
A case study of the variational assimilation of GPS zenith delay observations into a mesoscale model.
J. Appl. Meteor
,
40
,
1559
1576
.
Dudhia
,
J.
,
1993
:
A nonhydrostatic version of the Penn State–NCAR mesoscale model: Validation tests and simulation of an Atlantic cyclone and cold front.
Mon. Wea. Rev
,
121
,
1493
1513
.
Ehrendorfer
,
M.
, and
R. M.
Errico
,
1995
:
Mesoscale predictability and the spectrum of optimal perturbations.
J. Atmos. Sci
,
52
,
3475
3500
.
Grell
,
G. A.
,
J.
Dudhia
, and
D. R.
Stauffer
,
1994
:
A description of the fifth-generation Penn State/NCAR mesoscale model (MM5). Tech. Note 398+1A, NCAR, 122 pp
.
Henderson
,
J. M.
,
R. N.
Hoffman
,
S. M.
Leidner
,
T.
Nehrkorn
, and
C.
Grassotti
,
2005
:
A 4DVAR study on the potential of weather control and exigent weather forecasting.
Quart. J. Roy. Meteor. Soc
,
131
,
3037
3052
.
Hoffman
,
R. N.
,
2002
:
Controlling the global weather.
Bull. Amer. Meteor. Soc
,
83
,
241
248
.
Hoffman
,
R. N.
,
J. M.
Henderson
, and
S. M.
Leidner
,
2002
:
Using 4DVAR to move a simulated hurricane in a mesoscale model. Preprints, 19th Conf. on Weather Analysis and Forecasting/15th Conf. on Numerical Weather Prediction, San Antonio, TX, Amer. Meteor. Soc., J137–J140
.
Hoffman
,
R. N.
,
C.
Grassotti
,
J. M.
Henderson
,
S. M.
Leidner
,
G.
Modica
, and
T.
Nehrkorn
,
2005
:
Controlling the evolution of a simulated hurricane through optimal perturbations: Initial experiments using a 4D variational analysis system. Preprints, 16th Conf. on Planned and Inadvertent Weather Modification, San Diego, CA, Amer. Meteor. Soc., CD-ROM, 3.6
.
Holton
,
J. R.
,
1992
:
An Introduction to Dynamic Meteorology.
3d ed. Academic Press, 511 pp
.
Kalnay
,
E.
, and
Coauthors
,
1996
:
The NCEP/NCAR 40-Year Reanalysis Project.
Bull. Amer. Meteor. Soc
,
77
,
437
471
.
Lawrence
,
M. B.
, and
E. N.
Rappaport
,
1994
:
Eastern North Pacific hurricane season of 1992.
Mon. Wea. Rev
,
122
,
549
558
.
Liu
,
Y.
,
D-L.
Zhang
, and
M. K.
Yau
,
1999
:
A multiscale numerical study of Hurricane Andrew 1992. Part II: Kinematics and inner-core structures.
Mon. Wea. Rev
,
127
,
2597
2616
.
Rabier
,
R.
,
E.
Klinker
,
P.
Courtier
, and
A.
Hollingsworth
,
1996
:
Sensitivity of forecast errors to initial conditions.
Quart. J. Roy. Meteor. Soc
,
122
,
121
150
.
Thépaut
,
J-N.
,
D.
Vasiljevic
,
P.
Courtier
, and
J.
Pailleux
,
1993
:
Variational assimilation of conventional meteorological observations with a multilevel primitive-equation model.
Quart. J. Roy. Meteor. Soc
,
119
,
153
186
.
Unanwa
,
C. O.
,
J. R.
McDonald
,
K. C.
Mehta
, and
D. A.
Smith
,
2000
:
The development of wind damage bands for building.
J. Wind Eng. Ind. Aerodyn
,
84
,
119
149
.
Zou
,
X.
,
F.
Vandenberghe
,
M.
Pondeca
, and
Y-H.
Kuo
,
1997
:
Introduction to adjoint techniques and the MM5 adjoint modeling system. Tech. Note 435-STR, NCAR, Boulder, CO, 110 pp
.
Zou
,
X.
,
Q.
Xiao
,
A. E.
Lipton
, and
G. D.
Modica
,
2001
:
A numerical study of the effect of GOES sounder cloud-cleared brightness temperatures on the prediction of Hurricane Felix.
J. Appl. Meteor
,
40
,
34
55
.

Footnotes

Corresponding author address: Ross N. Hoffman, Atmospheric and Environmental Research, Inc., 131 Hartwell Avenue, Lexington, MA 02421. Email: rhoffman@aer.com