Ensemble model output statistics for wind vectors

1 MB·English

by  Nina Schuhen

Ensemble model output statistics for wind vectors
Nina Schuhen, Thordis L. Thorarinsdottir and Tilmann Gneiting
Institute of Applied Mathematics, University of Heidelberg, Germany
2 January 13, 2012
1
0
2

n Abstract
a
J A bivariate ensemble model output statistics (EMOS) technique for the postprocessing of

2 ensemble forecasts of two-dimensional wind vectors is proposed, where the postprocessed
1
probabilistic forecast takes the form of a bivariate normal probability density function. The

postprocessed means and variances of the wind vector components are linearly bias-corrected
]
P
versions oftheensemble meansandensemblevariances, respectively, andtheconditional cor-
A
relation between the wind components is represented by a trigonometric function of the en-
.
t semble mean wind direction. In a case study on 48-hour forecasts of wind vectors over the
a
t North American Pacific Northwest with the University of Washington Mesoscale Ensemble,
s
[ the bivariate EMOSdensity forecasts werecalibrated and sharp, andshowed considerable im-

provementovertherawensembleandreferenceforecasts,includingensemblecopulacoupling.
1
v
2
1 1 Introduction
6
2
1. The past two decades have seen a change of paradigms in weather forecasting, in that ensemble
0 prediction systems have been developed and implemented operationally (Leutbecherand Palmer,
2
2008). Ensemble systems seek to reflect and quantify sources of uncertainty in numerical weather
1
: forecasts, such as imperfections in initial conditions and incomplete mathematical representations
v
i oftheatmosphere. Despitetheubiquitouspositivespread-skillrelationship(Whitakerand Loughe,
X
1998; Grimitand Mass, 2002), ensemble forecasts tend to be biased, and typically they are un-
r
a derdispersed (HamillandColucci, 1997), in that the ensemble spread is too small to be realistic.
Furthermore, differing spatial resolutions of the forecast grid and the observation network may
need tobereconciled.
Toaddresstheseshortcomings,varioustechniquesforthestatisticalpostprocessingofensemble
modeloutputhavebeendeveloped(WilksandHamill,2007),withensemblemodeloutputstatistics
(EMOS)ornonhomogeneousGaussianregression(Gneitinget al.,2005;ThorarinsdottirandGneiting,
2010) being a state of the art method. The EMOS technique transforms a raw ensemble forecast
intoapredictiveprobabilitydensityfunction,andsimultaneouslycorrectsforbiasesanddispersion
errors. EMOSmethodshave been developed fortemperature and surface pressure(Gneitinget al.,
2005; Hagedorn etal., 2008; Kannet al., 2009), where the predictive density is normal and the
methodis often referred to as nonhomogeneousGaussian regression, for quantitativeprecipitation
(Wilks,2009),andforwindspeed(Thorarinsdottirand Gneiting,2010;Thorarinsdottirand Johnson,
2011). Inalltheseimplementations,thepredictivedensityappliestoaunivariateweatherquantity.
1
Local EMOS
3 m/s
Figure 1: Raw ensemble and EMOS postprocessed forecasts of surface wind vectors at stations in the
Olympic Peninsula and Puget Sound area in the US state of Washington, valid October 20, 2008 at 00
UTC, at a prediction horizon of 48 hours. The eight members of the University of Washington Mesoscale
Ensemble(UWME;EckelandMass2005)areshownasgraydots. The75%predictionellipseandthemean
vector for the EMOS density forecast are shown in dark gray, and the verifying wind vector is represented
byasmallblackcircle.
In this current paper, we propose and develop an EMOS technique for a bivariate weather
quantity, namely surface wind vectors, comprising both zonal and meridional components or, in
an alternativebut mathematicallyequivalent representation, wind speed and wind direction. Prob-
abilistic forecasts of wind conditions are critical in a wide range of applications, including air
traffic control, ship routing, recreational and competitivesailing, and wind energy, where their so-
cietaland monetaryvalueis huge(Marquiset al., 2011). Untilveryrecently, windspeedand wind
direction have been addressed independently in statistical postprocessing, without taking depen-
dencies into account (Bao et al., 2010; Sloughteret al., 2010; ThorarinsdottirandGneiting, 2010;
ThorarinsdottirandJohnson,2011). However,inmanyoftheaforementionedapplicationsitisim-
portant to honor the full information about the bivariate structure of the future wind vector that is
provided by the ensemble. Thus, our EMOS postprocessed forecasts take the form of elliptically
symmetric bivariate normal densities, as illustrated in Figure 1 in an application to the University
ofWashingtonMesoscaleEnsemble(UWME;Eckel and Mass2005).
The remainder of the paper is organized as follows. In Section 2 we provide the details of
2
the bivariate EMOS technique. A case study is presented in Section 3, where we consider 48-
hourahead forecasts ofwindvectorsovertheNorth AmericanPacific Northwestin2008based on
the eight-member UWME. The paper closes in Section 4, where we hint at future developments
and discuss the similarities and differences between our EMOS technique, the BMA approach of
Sloughter (2009); Sloughteret al. (2011), ensemble copula coupling (ECC; Schefzik 2011) and
the postprocessing method proposed by Pinson (2011), all of which are directed at the bivariate
postprocessing of ensemble forecasts of wind vectors. The Appendix describes our verification
methods.
2 Ensemble model output statistics for wind vectors
A wind vector is determined by wind speed and wind direction, or by its zonal (west-east) and
meridional(north-south)components, which we denote by u and v, respectively. We now develop
an ensemble model output statistics (EMOS) method for wind vectors, where the postprocessed
probabilistic forecast takes the form of a bivariate probability density function. The method is
tailored to ensembles with relatively few members, such as the eight-memberUniversityof Wash-
ington Mesoscale Ensemble (UWME; Eckel and Mass 2005), and we illustrate it using forecast
andobservationdatafrom thisensemble.
2.1 Bivariate normal distribution
Our EMOS postprocessed forecast takes the form of an elliptically symmetric, bivariate normal
probabilitydensityfunctionforthewindvector(u,v),withtheparametersofthisdistributionbeing
specified in terms of the ensemble forecast. The analytic form of the bivariate normal probability
densityfunctionis
1
f(u,v) = (1)
2πσ σ 1−ρ2
u v uv
p1 (u−µ )2 (u−µ )(v −µ ) (v−µ )2
u u v v
×exp − −2ρ + .
uv
(cid:18) 2(1−ρ2 ) (cid:18) σ2 σ σ σ2 (cid:19)(cid:19)
uv u u v v
Thetasknowistospecifythefiveparametersin(1),namelythemeanvalues,µ andµ ,ofthewind
u v
vectorcomponentsuand v,thecorrespondingmarginalvariances,σ2 andσ2, respectively,and the
u v
correlation coefficient, ρ , between the wind components, in their dependence on the ensemble
uv
forecast.
Bivariate normal density forecasts for wind vectors have also been proposed by Gneitinget al.
(2008), though in very crude form. The ingenious method of Pinson (2011) estimates a dilation
and translation of an ensemble forecast of wind vectors based on bivariate normal densities, and
ourapproach is verysimilarin itstreatmentofthemean andvarianceparameters. However,major
differences between the approach of Pinson (2011) and our method lie in the form of the post-
processed forecast, which is a probability density function in our case, rather than a dilated and
translated ensemble, and in the explicit modeling of the correlation coefficient ρ in our method.
uv
Foramoredetailed comparisonand recommendationsforajudiciouschoiceofthemostappropri-
atemethod,givenanyparticularensembleandtask at hand,werefer to Section4.
3
In the description that follows, we consider a general ensemble with m members, and denote
theindividualwindvectorforecasts by (u ,v ),...,(u ,v ), respectively.
1 1 m m
2.2 Means
In our standard implementation of the bivariate EMOS technique, the means µ and µ are bias-
u v
corrected versionsoftherespectiveensemblemeans,in that
µ = a +b u¯ and µ = a +b v¯, (2)
u u u v v v
whereu¯ = 1 m u andv¯ = 1 m v . Thebiascorrection parameters a , b ,a andb arees-
m i=1 i m i=1 i u u v v
timatedfrom tPrainingdata,in wayPsdescribed below. In aslightlymoreambitiousimplementation,
themean componentsµ andµ areaffine functionsoftheindividualensemblememberforecasts,
u v
namely
µ = a +b u +···+b u and µ = a +b v +···+b v ,
u u u,1 1 u,m m v v v,1 1 v,m m
wheretheregressionparametersa ,b ,...,b ,a andb ,...,b areestimatedfromtraining
u u,1 u,m v v,1 v,m
data. Thisgeneral versionappliestoensembleswithnon-exchangeablemembersonlyandreduces
tothestandardversionwhenb = ··· = b = 1 andb = ··· = b = 1. Inourexperiences
u,1 u,m m v,1 v,m m
withtheUWME,whichhasnon-exchangeablemembers,thegeneralversiongaveonlyveryslightly
improvedpredictiveperformance,andsowereportresultsforthestandardimplementation(2)only.
2.3 Variances
We specify the marginal variances σ2 and σ2 of the bivariate normal density forecast (1) as affine
u v
functionsoftherespectiveensemblevariances,in that
σ2 = c +d s2 and σ2 = c +d s2, (3)
u u u u v v v v
wheres2 = 1 m (u −u¯)2 ands2 = 1 m (v −v¯)2. Thedispersioncorrectionparametersc ,
u m i=1 i v m i=1 i u
du, cv and dv aPre estimated from training dPata, as described below. To guarantee the nonnegativity
of the variances, we constrain the parameters to be nonnegative, using the technique described by
ThorarinsdottirandGneiting(2010).
2.4 Correlation coefficient
Thekey characteristicand majorinnovationin ourwork is theexplicitmodelingof thecorrelation
coefficientρ ofthepostprocessedbivariatenormaldensityforecast (1).
uv
To motivateour specification of ρ , we consider ensembleforecast and observation data from
uv
the UWME in calendar year 2007, comprising a total of 23,250 forecasts cases at 79 meteorolog-
ical stations in the Pacific Northwest. Figure 2 distinguishes nine sectors for (u,v) wind vector
forecasts and observations. Figure 3 shows wind vectorobservationsin 2007, conditionallyon the
ensemble mean wind vector falling into any given sector. The conditional distributions are partly
elliptically contoured, particularly in the first sector, and partly skewed, with an orientation that
depends strongly on the sector, and they reflect the discretized nature of wind observations, as
discussedin moredetail inSection 3.
4
Figure 2: Sectors for wind vector forecasts and observations in the (u,v) plane. Sector 1 is the circular
region that is centered at the origin and corresponds to wind speeds less than or equal to two meters per
second. Sectors 2–9 are assigned clockwise, with Sectors 2–3 corresponding to a south-westerly, Sectors
4–5anorth-westerly, Sectors6–7anorth-easterly, andSectors8–9asouth-easterly wind.
The left-hand panel in Figure 4 plots the correlation coefficient between the wind components
in the scatter plots in Figure 3 as a function of the wind directions that correspond to the centers
of sectors 2–9. The panel demonstrates that the ensemble mean wind direction ought to have a
profound influence on the correlation coefficient ρ in the postprocessed EMOS density forecast
uv
(1). Thus, we model the correlation coefficient ρ as a trigonometric function of the ensemble
uv
meanwinddirection,θ, measured indegrees, inthat
2π
ρ = rcos (kθ+ϕ) +s, (4)
uv
(cid:18)360 (cid:19)
wheretheparameters r, s, k and ϕ are estimatedfrom trainingdata, in ways described below. The
coefficientsr and s concern theoverallmagnitudeofthecorrelationcoefficient andneed tosatisfy
|r|+|s| ≤ 1. Theparameterk correspondstothenumberofperiodsofthetrigonometricfunction,
whichweconstrainto beeitherk = 1,2 or3,and thedirectionϕ encodesphaseinformation.
We apply the correlation model (4) to forecasts with ensemble mean vectors in all nine wind
sectors, including the first sector. Alternatively, if the ensemble mean wind vector falls into the
central first sector, one can take ρ to be equal to the empirical correlation coefficient in the cor-
uv
responding scatterplot. In our experiences with the UWME, the two approaches resulted in nearly
identicalpredictiveperformance.
2.5 Estimation
OurEMOSpostprocesseddensityforecast forawindvectortakes theformofthebivariatenormal
probabilitydensity (1), where the relationshipsof µ ,µ ,σ2,σ2 and ρ to the raw ensemble fore-
u v u v uv
cast are specified in equations (2), (3) and (4). It remains to estimate the parameters that govern
theseequations,and weaddressthistask inthreephases.
Before describing the phases of our estimation scheme, it is worth noting that we follow
ThorarinsdottirandGneiting(2010)andconsidertwodistinctvariants,whichwecalltheRegional
EMOSandtheLocalEMOSmethod,respectively. In theRegionalEMOSmethod,onlyonesetof
5
Sector 9 Sector 2 Sector 3
5 5 5
1 1 1
0 0 0
1 1 1
5 5 5
V 0 V 0 V 0
5 5 5
− − −
0 0 0
1 1 1
− − −
5 5 5
1 1 1
− − −
−15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15
U U U
Sector 8 Sector 1 Sector 4
5 5 5
1 1 1
0 0 0
1 1 1
5 5 5
V 0 V 0 V 0
5 5 5
− − −
0 0 0
1 1 1
− − −
5 5 5
1 1 1
− − −
−15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15
U U U
Sector 7 Sector 6 Sector 5
5 5 5
1 1 1
0 0 0
1 1 1
5 5 5
V 0 V 0 V 0
5 5 5
− − −
0 0 0
1 1 1
− − −
5 5 5
1 1 1
− − −
−15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15
U U U
Figure 3: Wind vector observations over the Pacific Northwest in 2007, conditional on the ensemble mean
forecast falling into one of the sectors defined in Figure 2. The unit for the wind components is meters per
second.
6
Pacific Northwest Sea−Tac Airport (Station KSEA)
0 0
1. 1.
5 5
0. 0. 39
48
n n 97
o o
orrelati0.0 1075 490 4119 5555 2545 orrelati0.0 24 1 1 1
C C
1794 3512
5 5
0. 829 0. 5
− −
0 0
1. 1.
− −
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350
Wind Direction in Degrees Wind Direction in Degrees
Figure 4: Left: Correlation coefficient between the wind components in the scatter plots in Figure 3 as a
function of the wind direction that corresponds to the center of the sector. The correlation coefficients and
observation counts correspond to sectors 2–9, respectively. The dashed curve shows the correlation model
(4) as fitted by the weighted least squares method. Right: Same as left, but considering data from Sea-Tac
Airportonly.
parametersisestimatedandusedtoproduceforecastsovertheentireensembledomain,suchasthe
Pacific Northwest for the UWME. Training sets thus comprise data from all stations. In contrast,
theLocalEMOSmethodusestrainingdatafromthestationathandonly,andthusobtainsadistinct
setofparameters foreach station.
We now turn to the first phase of our estimation scheme, in which we fit the correlation model
(4). We do this offline, once and for all, based on historic, out-of-sample forecast and observation
data. Specifically, for the UWME, we use data from calendar year 2007 to form the conditional
scatterplotsin Figure 3 and computeand plot thecorresponding empirical correlation coefficients,
asillustratedinFigure4. Wethendecideaboutasuitablevalueforthenumberofcycles,k,andfit
the remaining parameters, r, s and ϕ, of the correlation model (4) by a weighted non-linear least
squares technique, using the R function NLS with the weights being proportionalto the numberof
observationsinthesectors(RDevelopmentCoreTeam,2011). Theuseofaweightedleastsquares
techniqueis critical, particularly for the Local EMOS technique, as local wind patterns may result
in very few observations being available in any given sector. This first phase of the estimation is
done once and for all, using historic data from 2007, and the fitted correlation model is applied
throughoutcalendaryear2008,which wetookas ourtestperiod.
The left-hand panel in Figure 4 shows the fitted correlation model for the Regional EMOS
method, where we chose k = 2 and obtained weighted least squares estimates of r = 0.20,
s = −0.15 and ϕ = −61.9 degrees, respectively. As noted, the Regional EMOS method uses
these parameters throughout the UWME domain and throughout the test period in calendar year
2008. Theright-handpanelshowstheLocalEMOScorrelationmodelatSeattle-Tacoma(Sea-Tac)
Airport,wherewetookk = 1andobtainedweightedleastsquaresestimatesofr = 0.24,s = 0.07
and ϕ = 70.5 degrees, respectively. The fitted Local EMOS correlation models at the remaining
stationsintheUWMEdomainare providedand illustratedin theappendix ofSchuhen (2011).
In contrast to the first phaseof our estimationscheme, the second and third stages proceed on-
7
line, that is, they userolling trainingperiods consistingof data from therecent past. The Regional
EMOS method uses all available data from the Pacific Northwest from the last n days prior to the
forecast being made. Missing data are simply omitted from the training set. For the Local EMOS
method, the training period comprises data from the station at hand from the n most recent days
whereforecasts and observationswere available. In eithercase, we talkof aslidingn-day training
period.
In the second phase of our estimation scheme, the parameters a , b , a and b in the specifi-
u u v v
cation (2) of the mean vector are estimated from the training data by standard linear least squares
regression. In the third phase, the parameters c , d , c and d in the specification (3) of the
u u v v
marginal variances are estimated on the same rolling training period by the maximum likelihood
technique, with all other parameters being held fixed. In other words, we maximize the logarithm
ofthelikelihoodfunction,namely
l(c ,d ,c ,d ) = logf(x,t)(c ,d ,c ,d ), (5)
u u v v (x,t) u u v v
P
asafunctionoftheparametersc ,d ,c andd inthevariancemodel(3),whichareallconstrained
u u v v
to be nonnegative. The sum in the log likelihood function extends over all locations x and times t
forwhichtherearedatainthetrainingset. Anysingletermoftheformf(x,t)(c ,d ,c ,d )refersto
u u v v
thebivariatenormaldensity(1)evaluatedat theverifyingvaluesu = u(x,t) andv = v(x,t),withµ ,
u
µ andρ setatthenumericalvaluesimpliedbythemeanmodel(2)andthecorrelationmodel(4),
v uv
based on the parameter estimates from the first two phases of our estimation scheme, and putting
σ2 = c +d s2|(x,t) andσ2 = c +d s2|(x,t),wheres2|(x,t) ands2|(x,t) denotetheensemblevariances
u u u u v v v v u v
at location x and valid time t, respectively. The optimization is performed numerically with the
OPTIM function in R, using the Broyden-Fletcher-Goldfarb-Shanno algorithm with initial values
providedby thepreviousday’sestimates.
The interpretation of the right-hand side of (5) as a log likelihood function is valid only if the
forecast errors are independent between times and locations. While this is usually not the case, an
alternativeinterpretation as a mean logarithmicscore permits us to view theestimates as optimum
scoreestimates,which aretailoredto theestimationofforecasting models(Gneitingetal., 2005).
2.6 Example
Figure5andTable1illustratethepostprocessedLocalandRegionalEMOSdensityforecastsofthe
surface wind vector at Sea-Tac Airport, valid October 20, 2008 at 00 UTC, at a prediction horizon
of 48 hours. Thus, the forecast concerns the same valid time and the same prediction horizon as
the Local EMOS forecasts illustrated in Figure 1, where the station at Sea-Tac Airport is located
in the south-east Puget Sound area. The postprocessed density forecasts correct for the biases
and underdispersion in the raw UWME. The Local EMOS forecast retains the positivecorrelation
structureintheraw ensembleand issharperthan theRegional EMOSforecast.
3 Case study: Forecasting surface wind vectors over the Pacific
Northwest
We now consider the out-of-sample predictive performance of our two-dimensional EMOS tech-
niqueinacasestudyforwindvectorforecastsovertheNorthAmericanPacificNorthwestin2008,
8
Local EMOS Regional EMOS
4 4
2 90% 2 90%
0 0
V V
25%
−2 −2 25%
50%
50%
75%
4 4 75%
− −
90%
90%
6 6
− −
−4 −2 0 2 4 6 −4 −2 0 2 4 6
U U
Figure5: ContourplotofthepostprocessedLocalEMOS(left)andRegionalEMOS(right)densityforecasts,
alongwiththerawensembleforecast,ofthesurfacewindvectoratSea-TacAirport,validOctober20,2008at
00UTC,atapredictionhorizonof48hours. TheeightmembersoftheUniversityofWashingtonMesoscale
Ensemble (UWME; Eckel and Mass 2005) are shown as gray dots, along with the 90% prediction ellipse,
whichisbased onabivariate Gaussian fittothe ensemble values. The25%, 50%, 75% and90% prediction
ellipses and the mean vector for the postprocessed EMOS density forecast are shown in dark gray. The
verifying wind vector is represented by the small black circle at (u,v) = (1.45,−0.53). The units are in
meterspersecond.
Table1: Predictivemeans,variancesandcorrelationforthepostprocessedLocalandRegionalEMOSdensity
forecasts in Figure 5 of the surface wind vector at Sea-TacAirport, valid October 20, 2008 at00 UTC,at a
prediction horizon of48hours. Theparameters refertothegeneral bivariate normalprobability density (1),
withthewindcomponents represented inmeterspersecond.
Method Local EMOS RegionalEMOS
µ 0.84 0.11
u
µ 0.05 −0.19
v
σ2 1.99 3.87
u
σ2 4.00 4.31
v
ρ 0.33 −0.02
uv
9
based on the UniversityofWashington MesoscaleEnsemble(UWME; Eckel and Mass 2005). We
compare to the raw ensemble forecast and various reference techniques, such as ensemble copula
coupling, and assess the performance of the bivariate EMOS technique when forecasts of wind
speedare desiredonly.
3.1 University of Washington Mesoscale Ensemble
In our case study, the test set consists of forecasts of surface (10 meter) wind vectors based on the
UniversityofWashingtonMesoscaleEnsemble(UWME;Eckel and Mass2005)withvaliddatein
calendar year 2008, at a prediction horizon of 48 hours. The UWME is an eight-member multi-
analysisensemblethenbasedontheFifth-GenerationPennState/NCARMesoscaleModel(MM5)
with initial and lateral boundary conditions obtained from operational centers around the world.
The forecasts are made on a 12 km grid and the region covered is the Pacific Northwest region
of Western North America, including the US states of Washington, Oregon and Idaho, as well as
the southern part of the Canadian province of British Columbia. The forecasts were bilinearly
interpolatedfromthefoursurroundinggridpointstotheobservationlocationsandrotatedtomatch
thetruedirectionat each station.
Surface wind vector observations were provided by the weather observation stations in the
AutomatedSurfaceObservingSystemnetwork(National WeatherService,1998). Thevectorwind
quantity studied here is horizontal instantaneous surface wind, where ‘instantaneous’ means that
the wind was measured and averaged over the last two minutes before the valid time at 00 UTC.
The wind vector observations were recorded as wind speed and wind direction, where wind speed
was rounded to the nearest whole knot, where a knot equals 0.5144 meters per second, while
values below two knots were recorded as zero. The observations are thus discretized, as is easily
recognizable in Figure 3. Quality control procedures as described in Baars (2005) were applied to
theentiredataset,removingdatesand locationswithanymissingforecasts orobservations.
Forcalendar year2008, 19,282pairs ofensembleforecasts and observationswereavailableon
291distinctdaysand at79 distinctobservationlocations. Additionaldatafromtheyears 2006and
2007 were used to provide an appropriate rolling training period for all days in 2008, for the first
phaseestimationofthecorrelationmodel(4)andtoestablishtheoptimallengthoftherollingtrain-
ingperiod. FurtherinformationabouttheUWME,nowusingtheWRFmesoscalemodel,aswellas
realtimeforecastsandobservations,canbefoundonlineathttp://www.atmos.washington
.edu/~ens/uwme.cgi.
3.2 Training periods for Regional and Local EMOS
As noted, we distinguishRegional and Local EMOS forecasts. The Regional EMOS method uses
all available training data to estimate a single set of parameters that is used throughout the Pacific
Northwestdomain,andcanbeuseddirectlyonthemodelgridaswell. TheLocalEMOStechnique
usestrainingdatafromthestationathandonlytoobtainadistinctsetofparametersateachstation.
Thus, the method applies at observation stations only, and can not be used directly on the model
grid.
To determine a suitable value of the length n of the rolling training period for phases two and
three of our estimation scheme, as described in Section 2, we considered experiments with wind
vector forecasts based on the UWME in calendar year 2007. In these experiments, phase one of
10