Abstract

In Part I of this series on ensemble-based exigent analysis, a Lagrange multiplier minimization technique is used to estimate the exigent damage state (ExDS), the “worst case” with respect to a user-specified damage function and confidence level. Part II estimates the conditions antecedent to the ExDS using ensemble regression (ER), a linear inverse technique that employs an ensemble-estimated mapping matrix to propagate a predictor perturbation state into a predictand perturbation state. By propagating the exigent damage perturbations (ExDPs) from the heating degree days (HDD) and citrus tree case studies of Part I into their respective antecedent forecast state vectors, ER estimates the most probable antecedent perturbations expected to evolve into these ExDPs. Consistent with the physical expectations of a trough that precedes and coincides with the anomalously cold temperatures during the HDD case study, the ER-estimated antecedent 300-hPa geopotential height trough is approximately 59 and 17 m deeper than the ensemble mean at around the time of the ExDP as well as 24 h earlier, respectively. Statistics of the explained variance and from leave-one-out cross-validation runs indicate that the expected errors of these ER-estimated perturbations are smaller for the HDD case study than for the citrus tree case study.

1. Introduction

Gombos and Hoffman (2013, hereafter Part I) described ensemble-based exigent analysis, a technique to estimate the “exigent” or worst-case scenario (WCS), the forecast that maximizes the damage for a particular weather event and specified risk or confidence level. For each ensemble member, a damage function estimates the potential damage or cost (e.g., the heating demand) from the weather parameters (e.g., the near-surface air temperature). The potential damage is then weighted by what is at risk (e.g., the gridpoint number of inhabitants) to evaluate the actual damage. The exigent damage state (ExDS) is then calculated from the damage ensemble using a Lagrange multiplier optimization. The exigent scenario may be useful to emergency planners, insurers, and the general public because it is the unique [for multivariate Gaussian (mG) ensembles] forecast that maximizes the event-wide damage for a specified risk or confidence level.

This article combines exigent analysis with another ensemble-based technique called ensemble regression (ER; Gombos and Hansen 2008; Gombos 2009; Gombos et al. 2012) to estimate the most probable atmospheric model states expected to precede and coincide with the ExDSs, using the two case studies in Part I as examples. These pre-exigent conditions determined by ER are potentially useful to 1) understand the atmospheric dynamics that may lead to the ExDS (Gombos et al. 2012), 2) preemptively update the ExDS probability with incoming observations in advance of relatively slow data assimilations (Gombos 2009), and 3) identify, for purposes including supplementary forecast guidance and adaptive observing, the antecedent atmospheric features to which the ExDS is most sensitive (Gombos et al. 2012).

Ensemble regression is a multivariate linear inverse technique that uses ensemble model output to make inferences about the linear relationships between vector-valued forecast and/or analysis fields, often gridded “maps” of meteorological variables. ER uses the ensemble members of these fields as training samples to estimate a covariance-based matrix that maps an ensemble perturbation in the predictor field(s) to the most probable ensemble perturbation in the predictand field(s).

The multivariate nature of ER distinguishes it from a number of other univariate approaches that use ensemble-based statistics to infer relationships between scalar or vector predictors and scalar predictands, including the ensemble synoptic analysis of Hakim and Torn (2008), the ensemble transform Kalman filter targeting approach of Bishop et al. (2001), and the analysis sensitivity diagnostic of Liu et al. (2009). Whereas univariate techniques use only correlations, ER makes use of the joint distribution between the predictors and predictands to assess sensitivities, thereby ensuring that the predictand perturbation is a statistically feasible member of the predictand distribution. For example, ensemble synoptic analysis might find the statistical relationship between each point in the antecedent geopotential height field and the subsequent surface pressure at the center of a storm, while ER might find the statistical relationship between the entire antecedent geopotential height field and the entire subsequent surface pressure field.

Previously, using the ensemble covariances of potential vorticity, potential temperature, and geopotential height (Z) to approximate an ER operator, Gombos and Hansen (2008) showed that potential vorticity ER and the dynamical piecewise potential vorticity inversion of Davis and Emanuel (1991) yield nearly identical Z perturbations. Gombos et al. (2012) used ER to study the sensitivity of Supertyphoon Sepat's (2007) 1000-hPa potential vorticity track to the position and strength of the antecedent 500-hPa Z field.

Here, ensemble regression is applied to calculate the pre-exigent conditions for the two case studies of Part I—one on heating demand on 8–9 January 2010 and another on freeze-damage to Florida citrus trees on 11 January 2010. Part I estimated the heating degree days (HDD) 90%-WCS (i.e., the WCS at the 90% confidence level; Fig. 8 in Part I) to yield about 1.26% more heating demand than the ensemble average forecast and the citrus tree 90%-WCS (Fig. 15 in Part I) to damage approximately 14.2 million trees, about 4.3 times more than the ensemble average. For each case study, Part II now uses ER to map the exigent perturbation from Part I backward in time to estimate the corresponding pre-exigent conditions.

This article is organized as follows. Section 2a introduces some notation and sections 2b and 2c briefly review the concepts and mathematics of exigent analysis and of ER, respectively. Section 2d introduces the ER predictand Mahalanobis distance probability qy. Section 3 describes two correlation error metrics used to quantify the goodness of the ER. The heating demand and citrus tree case studies are presented in sections 4 and 5, respectively. Section 6 provides a summary and discussion.

2. Methods

a. Notation

Let x be a column vector representing an arbitrary model state or series of such states. Let p and y define transformations of x, which may make a subset of x and/or calculate diagnostic quantities from x at the analysis time and/or at one or more forecast times. In one example that follows, p, the ER predictor, will be a vector of forecast gridpoint HDD, and y, the ER predictand, will be a vector of forecast gridpoint 300-hPa geopotential heights. In this case, the time of p is after the time of y, and the ExDS will be used to predict, in a statistical sense, the antecedent geopotential heights associated with the 90%-WCS heating demand.

Given the ensemble members , let define an ensemble-averaged quantity and let and p′ define perturbations with respect to for an ensemble member and for an arbitrary model state, respectively. Let be the matrix formed by concatenating the column vectors , let refer to the exigent damage perturbation (ExDP), and let be the covariance of . Let the y variables be defined in the same way as the p variables. Let refer to the ER estimate. Note that the hat accent refers to an ER-estimated quantity, not to principal components (PCs) as in Part I. Let tp and ty denote the forecast lead times of the p and y variables, respectively.

Where subscripts are necessary, the ensemble members are indexed by n = 1, … , Nens; the elements of p are indexed by i = 1, … , I; and the elements of y are indexed by j = 1, … , J—making an I × Nens matrix and a J × Nens matrix. For example, in section 4, the damage ensemble is 2482 × 50. See Table 1 and Table 1 of Part I for a legend of the symbols and abbreviations used in this paper.

Table 1.

Meanings of symbols and abbreviations introduced in this article. See Table 1 in Part I for symbols used in this article that are introduced in Part I.

Meanings of symbols and abbreviations introduced in this article. See Table 1 in Part I for symbols used in this article that are introduced in Part I.
Meanings of symbols and abbreviations introduced in this article. See Table 1 in Part I for symbols used in this article that are introduced in Part I.

b. Exigent analysis

To define the WCS, Part I introduced the Mahalanobis distance quantile (MDQ) and its associated cumulative density function (CDF) value, the Mahalanobis distance probability (MDP), denoted by q, to describe the unusualness of a perturbation p′ with respect to the probability distribution function (PDF) defined by the ensemble mean and covariance . It is assumed here that this PDF is mG and is a reasonable approximation of the true state-dependent PDF. In ensemble phase space, any p′ whose squared statistical (Mahalanobis) distance (Mahalanobis 1936) to the origin (i.e., the mean of the ensemble perturbations) equals is located on the surface of an ellipsoid defined by . In exigent analysis, this ellipsoid is considered a confidence region of the ensemble PDF, an mG generalization of a confidence interval. The size of the ellipsoid is determined by the specified value (here ) of the MDP. The exigent MDQ, referred to as the -MDQ or , is given by the inverse chi-squared cumulative distribution function with ν degrees of freedom, such that . A randomly chosen perturbation from the ensemble PDF, at this Mahalanobis distance, is on the -MDQ ellipsoid. Equivalently, the exigent MDP (i.e., ) equals the integral of the PDF over the volume of the -MDQ ellipsoid.

Part I defined the damage functional, Jd = wTp, the weighted sum of damage over the I grid points. Here, w is a vector of weights of length I that ascribe a user-determined absolute or relative importance to each grid point. The cumulative probability of the damage functional is referred to as the damage functional probability (DFP) or by the symbol β. The DFP is the probability that Jd is less than some constant C and so equals the integral of the damage functional PDF in the half-space bounded by the hyperplane Jd = C in the direction −w.

In Part I, the specified MDP, , and the damage ensemble covariance, , together define the exigent damage perturbation as the perturbation at the -MDQ that maximizes . Equivalently, the ExDP is the perturbation on the specified confidence ellipsoid whose area-integrated weighted damage is maximized. Using the method of Lagrange multipliers, Part I showed that the ExDP at the -MDQ is

 
formula

where and .

Figures 1 and 2, respectively, show the ExDPs for for the HDD and citrus tree case studies. Note that Figs. 8 and 15 in Part I are the w-weighted versions of the ExDPs in Figs. 1 and 2, respectively. In this article, ER is used to estimate contemporaneous and antecedent conditions associated with these ExDPs.

Fig. 1.

The HDD ExDP, for . Figure 8 in Part I is the 90%-WCS anomaly damage map corresponding to this ExDP. Areas with low absolute values are white. Plus signs are city locations given in Fig. 5 in Part I. Units are HDD.

Fig. 1.

The HDD ExDP, for . Figure 8 in Part I is the 90%-WCS anomaly damage map corresponding to this ExDP. Areas with low absolute values are white. Plus signs are city locations given in Fig. 5 in Part I. Units are HDD.

Fig. 2.

The citrus tree ExDP, for . Figure 15 in Part I is the 90%-WCS anomaly damage map corresponding to this ExDP. Areas with low absolute values are green. Units are fraction of citrus trees damaged.

Fig. 2.

The citrus tree ExDP, for . Figure 15 in Part I is the 90%-WCS anomaly damage map corresponding to this ExDP. Areas with low absolute values are green. Units are fraction of citrus trees damaged.

c. ER

Ensemble regression (Gombos and Hansen 2008; Gombos 2009; Gombos et al. 2012) uses the ensemble members and as training samples to compute a multivariate linear regression operator given by

 
formula

Right multiplying Eq. (2) by and then rearranging yields

 
formula

where pp = and yp denotes the cross covariance of and . Equation (3) shows that is a function only of the covariance of the combined predictor and predictand ensemble (Hakim and Torn 2008). The  appendix describes extended exigent analysis, an equivalent alternative to using ER to estimate the preexigent perturbation that embeds the ER into the exigent analysis.

Assuming a sufficiently large and calibrated ensemble and a linear relationship between and , given any prescribed predictor perturbation p′, is used to approximate

 
formula

the most probable perturbation (with respect to the ensemble mean at ty) in the span of when the perturbation (with respect to the ensemble mean at tp) in the span of equals . The quantity is called the effective predictor perturbation and is equal to the projection of p′ onto the subspace spanned by (see Gombos and Hansen 2008 for details). In summary, ER uses the predictor and predictand ensemble covariances to approximate a regression operator , which is used to estimate the most probable perturbation when is the portion of the prescribed perturbation p′ resolved by (i.e., not in the null space of) .

In realistic applications, I and J are typically much greater than Nens, making rank deficient and and not strictly invertible. Moreover, strong spatial dependences due to geographic proximity are likely to render ill-conditioned. To alleviate problems related to multi-collinearity, inflated regression parameter variances, and computational tractability, the ER machinery is regularized by projecting and onto their leading np and ny PCs, respectively. See Gombos et al. (2012) and Part I for details. Also see the conclusion of Part I for a discussion of how the characteristics of the ensemble may affect the exigent analysis. There it is mentioned that neither the PC filtering approach used here, nor an alternative approach to reduce the impact of small sample size by applying the Gaspari and Cohn (1999) covariance localization, resulted in major changes to the ExDPs.

For the particular ER application where and , respectively, equal full model state vector perturbations and tp precedes tf, ER is effectively a linear least squares approximation [subject to sampling and other errors discussed in section 3 and Gombos and Hansen (2008)] to the full nonlinear model used to integrate the ensemble. For these applications, is analogous to the transition matrices used in linear inverse modeling (LIM; Penland 1989) and climate applications of the fluctuation dissipation theorem (FDT; Leith 1975) that have been used to predict low-frequency Northern Hemisphere 700-hPa geopotential height anomalies (Penland and Ghil 1993), seasonal and interannual tropical Atlantic sea surface temperatures (Penland and Matrosova 1998), and the excitation of Atlantic storm tracks by tropical heating (Gritsun and Branstator 2007). The ER operator, however, is approximated from state-dependent synoptic-scale ensemble statistics, whereas LIM and FDT operators are typically estimated from stationary decadal-scale climate statistics.

For more general applications, such as the ones in this article, for which one (transformed) subset of the state vector at tp (e.g., the HDD forecast ensemble) is used to estimate the most probable perturbation (e.g., the ER-estimated preexigent perturbation) of a different subset of the state vector at ty (e.g., the 300-hPa Z forecast ensemble), ER is simply ensemble-based multivariate regression with vector-valued predictors and predictands. For these general applications, ER uses least squares to relate the fields of interest and is not intended to approximate the atmospheric model dynamics.

For the ER applications in this article, in order to estimate atmospheric conditions antecedent to the ExDP, the prescribed perturbation is , the ExDP estimated from Eq. (1). Note that the ExDP is, by construction, always spanned by , so that in this case . When the predictand ensemble is a subset of or the entire atmospheric state vector antecedent to tp, the ER-estimated predictand perturbation approximates the most probable (perturbation) conditions at ty spanned by . The sum of and approximates the projection of the true preexigent perturbation onto the space spanned by .

d. Ensemble regression Mahalanobis distance probability, dy

As summarized in section 2b, q is the value of the CDF of the Mahalanobis distance of an ensemble perturbation with respect to a particular ensemble distribution. Here, and in Part I, this ensemble distribution is that of the damage state vector p. A similar quantity is now defined for the ER predictand distribution. Let qy denote the value of the CDF of the Mahalonbis distance of an ER predictand perturbation with respect to this distribution. In analogy to Eq. (4) of Part I, qy is defined by

 
formula

where is the Mahalanobis distance of with respect to the predictand ensemble and ν denotes the degrees of freedom of . Equation (5) assumes is mG and perfectly calibrated and that ν and the covariance of are known exactly. See the discussion of Eqs. (4) and (5) in Part I for more details.

3. ER error metrics

Throughout this article it is assumed that the ensembles follow mG distributions and that the predictors explain a significant fraction of the total variance of the respective predictands. The validity of the mG assumption is assessed using the Gaussian QQ plot (e.g., Wilks 2006) in Part I. This section outlines two metrics that are used to assess the validity of the covariability assumptions: the composite correlation coefficient (Glahn 1968) and leave-one-out cross validation (LOOCV; e.g., Wahba and Wendelberger 1980; Michaelson 1987; Wilks 2006; Gombos et al. 2012). The and LOOCV error statistics for the HDD and citrus tree case studies are presented in sections 4b and 5b, respectively.

a. Composite correlation coefficient

The composite correlation coefficient (Glahn 1968), , is used to assess the degree to which the np predictor PCs are linearly related to the ny predictand PCs. Here, can be considered a multivariate extension to the ordinary (Pearson) correlation (e.g., Wilks 2006) for the assessment of the joint statistical linear association between vector-valued predictors and vector-valued predictands. Additionally, gives the fraction of the total variance of the ny predictand PCs accounted for by the np predictor PCs:

 
formula

where tr denotes the trace of the indicated argument, Σ is a diagonal matrix of the variances of the ny predictand PCs, and is the diagonal matrix of multiple correlation coefficients [i.e., the coefficients of determination in multiple regression; e.g., Glahn (1968); DelSole and Tippett (2009)]. Using j = 1, 2, … , ny to index the predictand PCs, the jth entry of the diagonal of is equal to the multiple correlation coefficient between the jth predictand PC and the np predictor PCs. Note that, because PCs are defined to be orthogonal, the square of the jth entry of the diagonal of is nothing more than the sum of the squares of the ordinary correlation coefficients between the jth predictand PC and each of the np predictor PCs (or, equivalently, the reduction of the variance of the jth predictand PC by the np predictor PCs) (e.g., Abdi 2007). Note that quantifies the statistical linear relationship between the predictor PCs and predictand PCs, not between the full predictor and predictand fields; however, large values of do suggest a strong statistical association between the full fields as long as the PCs explain a substantial fraction of the variance of their respective full fields.

b. LOOCV

The LOOCV technique discussed in Gombos (2009) and Gombos et al. (2012) is applied to estimate the expected ER accuracy. The technique estimates the error of ER-estimated perturbations by 1) removing a single predictor ensemble perturbation (and its twin ensemble perturbation, which is equal in magnitude and opposite in direction about the ensemble mean at the analysis time) and the corresponding predictand perturbation (and its twin ensemble perturbation) from and , respectively, when calculating ; 2) applying to the left-out predictor perturbation ; and then 3) comparing the ER-estimated predictand perturbation to the actual left-out predictand perturbation . Note that the LOOCV is performed using the leading np and ny PCs.

After repeating the LOOCV for each ensemble member, the median of the nens anomaly correlation coefficients (ACCs; e.g., Wilks 2006) between each and pair is a measure of the ER error. The median is used rather than the mean to reduce the influence of outliers on the expected value. For the particular ER application where and equal full state vector perturbations and tp precedes ty, the median ACC is a measure of how closely approximates the full nonlinear model used to integrate the ensemble. For more general applications, LOOCV only approximates the expected ACCs for the specific ER application of interest. For the backcasting LOOCV applications in this article, the median ACC approximates the expected value of the ACC between the ER-estimated preexigent perturbation and the true perturbation antecedent to the ExDP. The variance of the ACCs is a measure of the spread of the expected ACC value.

4. HDD case study

This section uses ER to estimate the atmospheric conditions that precede and coincide with the HDD ExDP (Fig. 1). Section 4a describes the HDD data and defines the predictor ensemble, predictor perturbation, and predictand ensembles used for the ER. Section 4b presents and LOOCV errors statistics. Section 4c presents and analyzes the ER-estimated antecedent perturbations, depicted in Figs. 3 and 4, associated with the HDD ExDP.

a. HDD data and synoptic overview

The overall synoptic scenario during the time period of this HDD case study was marked by a strong upper-level trough that swept eastward through the central and eastern United States. The magenta lines in Fig. 3 contour the ensemble mean forecast 300-hPa Z at forecast lead times of tf = 0, 12, 24, and 36 h. Northerly cold-air advection associated with the strong trough brought anomalously cold temperatures toward the domain of interest. This set the stage for a potentially extreme bout of cold temperatures, for which the HDD ExDP (measured in HDD units) is depicted in Fig. 1. The remainder of this section uses ER to estimate the evolution of the position and strength of this trough throughout the exigent scenario.

Fig. 3.

The 300-hPa Z-only portion of the ER predictand perturbation at lead times ty = (a) 0, (b) 12, (c) 24, and (d) 36 h, associated with the HDD ExDP (Fig. 1). Filled contours (m) depict the perturbation with respect to the ensemble mean 300-hPa Z at the respective lead times. Black lines contour the 300-hPa Z of the full ER predictand state (i.e., the ensemble perturbation plus the ensemble mean) and magenta lines (labeled in green, m) contour the 300-hPa Z ensemble mean. The title of each panel states the respective value of qy, the MDP of the ER-predicted state. Areas with low absolute values are white. Plus signs are city locations given in Fig. 5 in Part I.

Fig. 3.

The 300-hPa Z-only portion of the ER predictand perturbation at lead times ty = (a) 0, (b) 12, (c) 24, and (d) 36 h, associated with the HDD ExDP (Fig. 1). Filled contours (m) depict the perturbation with respect to the ensemble mean 300-hPa Z at the respective lead times. Black lines contour the 300-hPa Z of the full ER predictand state (i.e., the ensemble perturbation plus the ensemble mean) and magenta lines (labeled in green, m) contour the 300-hPa Z ensemble mean. The title of each panel states the respective value of qy, the MDP of the ER-predicted state. Areas with low absolute values are white. Plus signs are city locations given in Fig. 5 in Part I.

The predictor ensemble used here is the HDD damage ensemble from Part I. That is, the predictor ensemble is the forecast daily average HDD for 8–9 January 2010 computed using the Nens = 50 ECMWF ensemble forecast T data initialized at 0000 UTC 8 January 2010, retrieved from the THORPEX Interactive Grand Global Ensemble (TIGGE) dataset (Bougeault et al. 2010), and linearly interpolated to 0.25° resolution at points between 31.75° and 40°N and 98° and 80°W. The forecast HDD is computed by applying the HDD damage function [Eq. (9) in Part I] to the average forecast T from the four forecast times tp centered on 2100 local (central standard) time on 8 January 2010 (i.e., the , , , and forecasts). The notation in black on the timeline in Fig. 5 depicts the sequence of events for the HDD case study. The HDD damage ensemble is projected onto its leading np = 7 PCs (see appendix C and Fig. 4 in Part I).

Fig. 4.

As in Fig. 3, but for the 2-m-temperature-only portion. Filled contours (contour interval 0.5°C) depict the perturbation with respect to the ensemble mean 2-m temperature at the respective lead times, but note that the perturbation for t = 0 h has been multiplied by 5 and that for t = 12 h by 2. For reference, magenta lines repeat the contours of the ExDP from Fig. 1.

Fig. 4.

As in Fig. 3, but for the 2-m-temperature-only portion. Filled contours (contour interval 0.5°C) depict the perturbation with respect to the ensemble mean 2-m temperature at the respective lead times, but note that the perturbation for t = 0 h has been multiplied by 5 and that for t = 12 h by 2. For reference, magenta lines repeat the contours of the ExDP from Fig. 1.

Fig. 5.

Timeline depicting the chronology of events for the ERs described in the text. As in Fig. 3 in Part I, but with ty added. Times in black correspond to the HDD case study and times in gray correspond to the citrus tree case study. Dots denote the initialization time t0, circles the forecast times tf, asterisks the times that define the damage ensemble tp, and squares the ER predictand lead times ty. Each time in the first line along the x axis represents HHDD, where HH is the UTC hour and DD is the day in January 2010. Each time in the second line represents LLDD, where LL refers to the local CST hour.

Fig. 5.

Timeline depicting the chronology of events for the ERs described in the text. As in Fig. 3 in Part I, but with ty added. Times in black correspond to the HDD case study and times in gray correspond to the citrus tree case study. Dots denote the initialization time t0, circles the forecast times tf, asterisks the times that define the damage ensemble tp, and squares the ER predictand lead times ty. Each time in the first line along the x axis represents HHDD, where HH is the UTC hour and DD is the day in January 2010. Each time in the second line represents LLDD, where LL refers to the local CST hour.

Because the goal here is to estimate the antecedent perturbations associated with the HDD ExDP, the predictor perturbation in this case is the HDD ExDP depicted in Fig. 1. Note that Fig. 8 in Part I is the corresponding 90%-WCS damage map, the population-density-weighted version of the ExDP in Fig. 1.

The predictand ensembles use ECMWF ensemble forecasts of 17 fields: 2-m temperature (T), and the zonal wind (u), meridional wind (υ), air temperature (Ta), and geopotential height (Z) at the 1000-, 850-, 500-, and 300-hPa pressure levels. The ensemble forecasts were initialized at 0000 UTC 8 January 2010 and also linearly interpolated to the same grid as the HDD. For each forecast lead time ty = 0, 12, 24, and 36 h, a separate ER is performed to estimate the most probable perturbation. For each lead time, the Nens = 50 ensemble members of the 17 fields are concatenated into a single predictand ensemble matrix that is projected onto its leading ny PCs. The values for ny for lead times ty = 0, 12, 24, and 36 h are 7, 9, 11, and 11, respectively, based on the respective integer rounded median of the four metrics described in appendix C of Part I (Maaten 2010). See the scree graphs (e.g., Wilks 2006) depicted in Fig. 6. Also note that temperatures are cool enough throughout the domain for all ensemble members so that T and HDD are linearly related.

Fig. 6.

Scree graph for the four HDD predictand ensembles described in section 4a at the ty = 0, 12, 24, and 36-h lead times. The number of meaningful singular vectors to retain and the number of effective degrees of freedom are approximated as np = ν = 7 for ty = 0 h, as np = ν = 11 for ty = 12 h, and as np = ν = 9 for ty = 24 and ty = 36 h (denoted by circles at the respective lead times). Only the leading 12 singular values are shown.

Fig. 6.

Scree graph for the four HDD predictand ensembles described in section 4a at the ty = 0, 12, 24, and 36-h lead times. The number of meaningful singular vectors to retain and the number of effective degrees of freedom are approximated as np = ν = 7 for ty = 0 h, as np = ν = 11 for ty = 12 h, and as np = ν = 9 for ty = 24 and ty = 36 h (denoted by circles at the respective lead times). Only the leading 12 singular values are shown.

b. Assessing the statistical linear relationship of the HDD predictor and predictand ensembles

The quality of the ER-estimated preexigent perturbations for the HDD case are assessed using the two metrics outlined in section 3. The solid line in Fig. 7 depicts the variation of with lead time. At ty = 0, the lead time most distant from the lead times that define the HDD damage ensemble (tp = 18, 24, 30, and 36 h) = 0.57, increases to = 0.82 at ty = 12 h, and then remains nearly constant through ty = 36 h. Note that the maximum of occurs between ty = 12 and 36 h because these lead times are the most contemporaneous with those that define the HDD damage ensemble. An value of 0.82 at ty = 12 h implies that 67% of the variance of the 300-hPa Z predictand ensemble PCs is explained by the HDD damage ensemble PCs and suggests that a linear model may be appropriate to model the relationship between perturbations of these two fields.

Fig. 7.

Boxplots of the ACC coefficients from the HDD LOOCV calculation described in section 3b. The dashed line connects the median ACCs at different lead times. The solid line depicts the composite correlation coefficient of the leading np = 7 PCs of the HDD damage ensemble and the leading ny = 8, 9, or 11 PCs of the predictand ensemble described in section 4a. Note that correlations and boxplots are shown for predictand lead times ty = 6, 18, and 30 h, in addition to the ER predictand lead times ty = 0, 12, 24, and 36 h labeled in Fig. 5.

Fig. 7.

Boxplots of the ACC coefficients from the HDD LOOCV calculation described in section 3b. The dashed line connects the median ACCs at different lead times. The solid line depicts the composite correlation coefficient of the leading np = 7 PCs of the HDD damage ensemble and the leading ny = 8, 9, or 11 PCs of the predictand ensemble described in section 4a. Note that correlations and boxplots are shown for predictand lead times ty = 6, 18, and 30 h, in addition to the ER predictand lead times ty = 0, 12, 24, and 36 h labeled in Fig. 5.

Figure 7 also displays a boxplot (e.g., Wilks 2006) at each lead time illustrating the variability of the ACC at that lead time. The ensemble median ACC between the ER-predictand perturbation and left-out ensemble member equals 0.40 at ty = 0 h, quickly increases to 0.78 at ty = 18 h, and then remains nearly constant through ty = 36 h. These ensemble median ACC values approximate the expected ACC values between the ER-estimated antecedent state vector perturbation and the actual antecedent perturbation, implying considerable confidence in the pattern of the ER predictions presented in the following. Note that the choices for np and ny are based on the four metrics described in appendix C of Part I and are not optimized to maximize the ACC; other choices for np and ny yield even higher median ACCs (not shown).

c. ER-estimated perturbations associated with the HDD ExDP

Using ExDP as the predictor perturbation, ER is employed to investigate the atmospheric conditions that precede and coincide with the HDD ExDP illustrated in Fig. 1. The filled contours in Fig. 3 depict the 300-hPa Z-only portion of the ER predictand perturbations from the four separate ERs with respective lead times ty = 0, 12, 24, and 36 h. Magenta contours depict the ensemble mean 300-hPa Z and black contours illustrate the preexigent 300-hPa trough state, the sum of the ensemble mean and the predictand perturbation. The 300-hPa-only portion is displayed here because it illustrates upper-tropospheric dynamics relevant to the anomalously cold temperatures at tp.

In Fig. 3, the ensemble mean 300-hPa Z represents a 300-hPa trough that deepens as it progresses eastward between ty = 0 h (Fig. 3a) and ty = 36 h (Fig. 3d). Compared to the mean trough (magenta), the preexigent 300-hPa trough state (black) is significantly stronger, with a −17-m perturbation at ty = 0 h that rapidly deepens to approximately −37 and −59 m at ty = 12 and 24 h, respectively. The maximum of the ty = 36 h Z perturbation is approximately 24 m, making the trough weaker than the ensemble mean at that time.

The preexigent 300-hPa Z perturbation (Fig. 3) is consistent with the physical expectations of a trough that precedes and coincides with anomalously cold temperatures in the domain of interest between lead times ty = 18 and 36 h. A deeper and stronger trough would strengthen northerly cold-air advection on its western side, ushering anomalously cold arctic air southward into the domain throughout this time window. At around ty = 36 h, cold air continues to advect southward on the western side of the exigent trough as an upstream ridge progresses eastward toward the center of the domain. See Gombos et al. (2012) for a more detailed example of ER perturbation patterns and their physical interpretations for the case of a tropical cyclone.

Meanwhile, surface temperatures decrease with lead time and pockets of cold-air anomalies evolve in a manner consistent with the HDD ExDP in Fig. 1. Figure 4 depicts the HDD ExDP from Fig. 1 (magenta lines) overlaid on top of the preexigent 2-m-temperature-only portion of the ER predictand perturbations from the four separate ERs with respective lead times ty = 0, 12, 24, and 36 h (filled contours). In the first two panels, the ER predictions are of small amplitude and, for plotting purposes, have been multiplied by 5 and 2 at ty = 0 and 12 h, respectively. At ty = 0 h, the preexigent T perturbation has a cold anomaly in the northern central portion of the domain. By ty = 12 h, this pocket begins to spread zonally, and a separate cold pocket forms to the southeast. At ty = 24 h and particularly at ty = 36 h, these cold pockets become collocated with the HDD ExDP (magenta lines). The ACCs between these ER predictions and the HDD ExDP are −0.58, −0.42, −0.82, and −0.84 for ty = 0, 12, 24, and 36 h, respectively, and the ACC between the sum of the 24- and 36-h ER predictions and the ExDP is −0.94. Considering that the weights, w, are population estimates, Fig. 4 depicts the evolution of a cold outbreak that targets the pattern of the HDD ExDP and is consistent with the ensemble statistics from lead time to lead time.

The title of each panel in Figs. 3 and 4 states the value of qy [Eq. (5)] for the respective predictand perturbations. Note that in each case. Since the Mahalanobis distance is invariant for linear affine transformations (e.g., Delsole and Tippett 2008), is expected to equal qy for “perfect” ER applications for which = 1 and is well fit for out-of-sample perturbations, invertible, and does not contain a null space. However, perfect ER operators are not to be expected because of null spaces associated with rank deficiency and PC truncation and because of the nonlinearity and complexity of atmospheric dynamics that results in < 1. Therefore, , even for well-chosen, quasi-linearly evolving mG ensembles.

5. Citrus tree case study

In contrast to the HDD case study, the statistical relationship between the citrus tree damage ensemble and the relevant upstream antecedent state vector predictand ensemble is weak for reasons explained below in section 5b.

a. Citrus tree data

The predictor ensemble for ER applications in this section is the citrus tree damage ensemble used in Part I. That is, the citrus freeze damage function [Eq. (10) from Part I] is applied to the Nens = 50 ECMWF ensemble forecast T data initialized at 1200 UTC 10 January 2010 and retrieved from the TIGGE dataset (Bougeault et al. 2010) to estimate the citrus freeze damage ensemble. This damage ensemble is used to approximate the forecast covariance of the fraction of trees damaged at tp = 24 h (near dawn local time). Note that a single forecast time tp is used to define the citrus damage ensemble, whereas four forecast times are used to define the HDD damage ensemble. The notation in gray on the timeline in Fig. 5 depicts the sequence of events for the citrus tree case study. These ECMWF T ensemble data are linearly interpolated to 0.25° resolution at points between 25.75° and 29.5°N and 82.75° and 80°W. The citrus tree damage ensemble is projected onto its leading np = ν = 5 PCs (see appendix C and Fig. 4 in Part I).

The predictor perturbation for this section equals the citrus tree ExDP depicted in Fig. 2, the forecast perturbation of the fraction of citrus trees expected to be damaged during the 90%-WCS freeze event at tp. Again, note that Fig. 15 in Part I is the corresponding 90%-WCS anomaly damage map, the citrus-tree-density-weighted version of the ExDP in Fig. 2.

The predictand ensembles used to assess the predictor–predictand relationship in the following subsection are the same as those described in section 4a, except that the ensemble data are initialized at 1200 UTC 10 January 2010 and interpolated to 0.25° resolution between 22° and 40°N and 103° and 80°W, which is a somewhat larger area than the area used for the HDD case. Predictand ensembles are projected onto their leading ny = 7, 10, 11, 10, or 11 PCs for ty = 0, 6, 12, 18, or 24 h, respectively, based on the four metrics (Maaten 2010) described in appendix C of Part I.

b. Citrus predictor–predictand relationship and expected errors

The values for the citrus tree case study (solid line in Fig. 8) are markedly worse than are those for the HDD case study; = 0.38 at ty = 0 h and remains nearly constant through ty = 24 h. Correspondingly, the median cross-validated ACCs (dashed line in Fig. 8) are also poor for the citrus tree case study; the ensemble median ACC between the ER-predictand perturbation and left-out ensemble member equals 0.12 at ty = 0 h and remains nearly constant through ty = 24 h.

Fig. 8.

As in Fig. 7, but for the citrus tree case study described in section 5a. Note that correlations and boxplots are shown for predictand lead times ty = 0, 6, 12, and 18 h, in addition to the ER predictand lead time ty = 24 h labeled in Fig. 5.

Fig. 8.

As in Fig. 7, but for the citrus tree case study described in section 5a. Note that correlations and boxplots are shown for predictand lead times ty = 0, 6, 12, and 18 h, in addition to the ER predictand lead time ty = 24 h labeled in Fig. 5.

The poor expected quality of the estimate of the preexigent conditions implied by the low and LOOCV values is in part attributable to the citrus damage function [Eq. (10) from Part I] being nonlinear. The citrus damage fraction is defined to equal zero over a large portion of the domain for most ensemble members (8969 times out of I × Nens = 195 × 50 = 9600) where the temperature is above −4.2°C, and to equal one over a small part of the domain for a few ensemble members (14 times out of 9600) where the temperature is below −6.7°C. However, as discussed below, the citrus damage ensemble is at least very strongly linearly related to the contemporaneous T ensemble.

Another factor contributing to the low and ACC values is the large size of the predictand domain and the small size of the predictor domain; the area covered by the citrus tree damage ensemble is only 2.6% of the area covered by the citrus tree predictand ensemble. Even for cases when the nonlinearity is negligible and the predictors and predictands are indeed physically related, the covariability of the atmospheric state over such a large area is unlikely to be strongly constrained by that of a relatively small region, especially one that is distant in both space and time.

c. ER-estimated perturbations associated with the citrus tree ExDP

Although the entire upstream state vector for this citrus tree case study cannot be accurately estimated, ER nevertheless can accurately estimate certain portions of this state vector. This section displays the ER-estimated contemporaneous T perturbation associated with the citrus tree ExDP and, given the low and ACC values, omits the ER estimate of the antecedent state vector over the much larger upstream domain.

The leading np = 5 PCs of the citrus damage ensemble explains a high fraction (90%) of the variance of the ty = 24 h forecast T over Florida (i.e., for the same domain and time as the damage ensemble) with = 0.95. This is unsurprising given that the damage ensemble is derived from this T ensemble and that the T values over a portion of the domain (approximately 6% of the grid points) are in the linear interval between −6.7° and −4.2°C.

The filled contours of Fig. 9 display the most probable contemporaneous tp = ty = 24 h T perturbation associated with the citrus tree ExDP. As is to be expected, given = 0.95, qy = 0.88 approximately equals and applying this temperature perturbation (Fig. 9) to the citrus damage equation [Eq. (10) in Part I] does in fact very closely approximate the ExDP in Fig. 2 (not shown).

Fig. 9.

The ER-estimated ty = 24 h T perturbation (shading). The MDP of the ER-predicted state is qy = 0.88. Magenta lines contour the variance of T (°C2) for lead time ty = 24 h. Areas with low absolute values are white.

Fig. 9.

The ER-estimated ty = 24 h T perturbation (shading). The MDP of the ER-predicted state is qy = 0.88. Magenta lines contour the variance of T (°C2) for lead time ty = 24 h. Areas with low absolute values are white.

Perhaps counterintuitively, the minimum of this contemporaneous T perturbation (Fig. 9) is located well to the northeast of the maximum citrus tree density (Fig. 12 in Part I). As portrayed by the ensemble variance map of the forecast T ensemble (magenta line contours in Fig. 9), this location mismatch is attributable to the significantly greater variability at the location of the temperature minimum than at the location of high citrus tree density. At the location of the variance maximum, temperatures are free to vary significantly without being extreme relative to the ensemble PDF, allowing for anomalously cold, yet plausible, temperatures at that location. On the other hand, the ensemble temperature variations are significantly smaller at the location of the high citrus tree density (Fig. 12 in Part I), so temperatures any lower at that location would correspond to the MDP exceeding 0.9. The positioning of the temperature perturbation in Fig. 9 is a result of having sufficiently cold temperatures at the locations of the citrus trees to maximize damage and strong positive temperature correlations at neighboring locations with high ensemble variability.

6. Summary and discussion

This is the second part of a two-part series on ensemble-based exigent analysis. In Part I, Gombos and Hoffman (2013) used a Lagrange multiplier technique to derive the equation for the exigent damage perturbation (ExDP) and then estimated the ExDP for two case studies from a cold outbreak in January 2010: one of a HDD 90%-WCS and another of a citrus tree 90%-WCS. This article combines exigent analysis with ensemble regression (ER) to predict the most probable perturbations expected to precede and/or coincide with the ExDPs from Part I. For the HDD case study, the ER results are consistent with physical expectations; the trough that precedes and coincides with the anomalously cold temperatures during the HDD case study associated with the ExDP (i.e., the ER-estimated preexigent 300-hPa geopotential height trough) is approximately 59 and 17 m deeper than the ensemble mean at around the time of the ExDP and 24 h earlier, respectively.

For the HDD case study, leave-one-out cross-validation (LOOCV) statistics suggest that the anomaly correlation coefficient (ACC) between the ER-estimated preexigent perturbation and the true (assuming perfect forecasts) preexigent perturbation varies between 0.40 and 0.78 depending on the lead time. For the citrus tree case study, and LOOCV statistics imply no useful ability to specify the preexigent conditions. This is attributable to the nonlinearity of the damage function and the relatively small size of the predictor ensemble domain causing a weak linear relationship between the damage ensemble and the antecedent state vector. Skillful estimates of the preexigent conditions can generally be expected when using large, mG, calibrated ensembles with strong linear covariance between the damage ensemble and the antecedent predictand ensemble.

Estimating preexigent conditions using ER has many potential applications. For example, the methods described in this article can be used to 1) gain insights into the atmospheric dynamics associated with a particular -WCS using the techniques outlined in Gombos et al. (2012), 2) preemptively forecast (Etherton 2007; Gombos 2009) worst-case scenarios in advance of slow data assimilation systems using supplementary forecast guidance, 3) target locations for adaptive observing in order to reduce the variance of high-impact forecasts (e.g., Ancell and Hakim 2007), and 4) perform damage-minimization weather modification (e.g., Hoffman 2002), where one could potentially use exigent analysis to find the ExDP that minimizes a damage functional, at an MDP level that reduces the damage to a reasonable degree dictated by resources and other requirements, and apply ER to find the associated preexigent perturbation that will minimize the damage. These and other possible ER applications have potential uses for forecasting and planning in future extreme weather situations that pose risks to life and property.

Acknowledgments

The authors gratefully acknowledge funding provided by National Science Foundation Grant 0838196.

APPENDIX

Embedding ER into Exigent Analysis

An equivalent alternative to the two-step algorithm of approximating the preexigent perturbation (by first estimating the ExDP and then using it as the predictor perturbation for ER) is extended exigent analysis, a one-step procedure that embeds ER directly into exigent analysis.

In extended exigent analysis, the predictor and predictand ensembles are concatenated so that is replaced by the (I + J) × nens matrix . Also, w is replaced by the (I + J) × 1 vector , where wp equals the original weighting matrix w and wy is a J × 1 vector of zeros. Then, if q and ν are chosen as in the original exigent analysis, the values of Qp and Qw are unchanged and the solution of the extended problem, obtained from Eq. (1), may be written as

 
formula

where pp = is the covariance of , yy is the covariance of , and is the covariance of and .

The upper blocks of the right and left forms in Eq. (A1) give the original exigent analysis [Eq. (1)]. The lower blocks give the ER result [Eq. (4)] for the case when the ER predictor is the ExDP , since

 
formula

using Eqs. (1) and (3).

REFERENCES

REFERENCES
Abdi
,
H.
,
2007
: Multiple correlation coefficient. Encyclopedia of Measurement and Statistics, N. Salkind, Ed., Sage Publications, 648–651.
Ancell
,
B.
, and
G.
Hakim
,
2007
:
Comparing adjoint- and ensemble-sensitivity analysis with application to observation targeting
.
Mon. Wea. Rev.
,
135
,
4117
4134
.
Bishop
,
C. H.
,
B. J.
Etherton
, and
S. J.
Majumdar
,
2001
:
Adaptive sampling with the ensemble transform Kalman filter. Part I: Theoretical aspects
.
Mon. Wea. Rev.
,
129
,
420
436
.
Bougeault
,
P.
, and
Coauthors
,
2010
:
The THORPEX Interactive Grand Global Ensemble
.
Bull. Amer. Meteor. Soc.
,
91
,
1059
1072
.
Davis
,
C. A.
, and
K.
Emanuel
,
1991
:
Potential vorticity diagnostics of cyclogenesis
.
Mon. Wea. Rev.
,
119
,
1929
1953
.
Delsole
,
T.
, and
M. K.
Tippett
,
2008
:
Predictable components and singular vectors
.
J. Atmos. Sci.
,
65
,
1666
1678
.
DelSole
,
T.
, and
M. K.
Tippett
,
2009
:
Average predictability time. Part II: Seamless diagnoses of predictability on multiple time scales
.
J. Atmos. Sci.
,
66
,
1188
1204
.
Etherton
,
B. J.
,
2007
:
Preemptive forecasts using an ensemble Kalman filter
.
Mon. Wea. Rev.
,
135
,
3484
3495
.
Gaspari
,
G.
, and
S. E.
Cohn
,
1999
:
Construction of correlation functions in two and three dimensions
.
Quart. J. Roy. Meteor. Soc.
,
125
,
723
757
.
Glahn
,
H. R.
,
1968
:
Canonical correlation and its relationship to discriminant analysis and multiple regression
.
J. Atmos. Sci.
,
25
,
23
31
.
Gombos
,
D.
,
2009
: Ensemble regression: Using ensemble model output for atmospheric dynamics and prediction. Ph.D. thesis, Massachusetts Institute of Technology, 190 pp.
Gombos
,
D.
, and
J.
Hansen
,
2008
:
Potential vorticity regression and its relationship to dynamical piecewise inversion
.
Mon. Wea. Rev.
,
136
,
2668
2682
.
Gombos
,
D.
, and
R. N.
Hoffman
,
2013
:
Ensemble-based exigent analysis. Part I: Estimating worst-case weather-related forecast damage scenarios
.
Wea. Forecasting
,
28
,
537
556
.
Gombos
,
D.
,
R. N.
Hoffman
, and
J. A.
Hansen
,
2012
:
Ensemble statistics for diagnosing dynamics: Tropical cyclone track forecast sensitivities revealed by ensemble regression
.
Mon. Wea. Rev.
,
140
,
2657
2669
.
Gritsun
,
A.
, and
G.
Branstator
,
2007
:
Climate response using a three-dimensional operator based on the fluctuation-dissipation theorem
.
J. Atmos. Sci.
,
64
,
2558
2575
.
Hakim
,
G.
, and
R. D.
Torn
,
2008
: Ensemble synoptic analysis. Synoptic–Dynamic Meteorology and Weather Analysis and Forecasting: A Tribute to Fred Sanders, L. Bosart and H. Bluestein, Eds., Amer. Meteor. Soc., 147–161.
Hoffman
,
R. N.
,
2002
:
Controlling the global weather
.
Bull. Amer. Meteor. Soc.
,
83
,
241
248
.
Leith
,
C. E.
,
1975
:
Climate response and fluctuation dissipation
.
J. Atmos. Sci.
,
32
,
2022
2026
.
Liu
,
J.
,
E.
Kalnay
,
T.
Miyoshi
, and
C.
Cardinali
,
2009
:
Analysis sensitivity calculation in an ensemble Kalman filter
.
Quart. J. Roy. Meteor. Soc.
,
135
,
1842
1851
.
Maaten
,
L.
, cited
2010
: Matlab toolbox for dimensionality reduction. University of California, San Diego. [Available online at http://homepage.tudelft.nl/19j49/Matlab_Toolbox_for_Dimensionality_Reduction.html.]
Mahalanobis
,
P. C.
,
1936
:
On the generalised distance in statistics
.
Proc. Natl. Inst. Sci. India
,
2
,
49
55
.
Michaelson
,
J.
,
1987
:
Cross-validation in statistical climate forecast models
.
J. Climate Appl. Meteor.
,
26
,
1589
1600
.
Penland
,
C.
,
1989
:
Random forcing and forecasting using principal oscillation pattern analysis
.
Mon. Wea. Rev.
,
117
,
2165
2185
.
Penland
,
C.
, and
M.
Ghil
,
1993
:
Forecasting Northern Hemisphere 700-mb geopotential height anomalies using empirical normal modes
.
Mon. Wea. Rev.
,
121
,
2355
2372
.
Penland
,
C.
, and
L.
Matrosova
,
1998
:
Prediction of tropical Atlantic sea surface temperatures using linear inverse modeling
.
J. Climate
,
11
,
483
496
.
Wahba
,
G.
, and
J.
Wendelberger
,
1980
:
Some new mathematical methods for variational objective analysis using splines and cross validation
.
Mon. Wea. Rev.
,
108
,
1122
1143
.
Wilks
,
D. S.
,
2006
: Statistical Methods in the Atmospheric Sciences. 2nd ed. Elsevier, 672 pp.