- Email: [email protected]
Improving the calculation of particle trajectories in the extra-tropical troposphere using standard NCEP fields
Atmospheric Environment 36 (2002) 483–490
Improving the calculation of particle trajectories in the extra-tropical troposphere using standard NCEP fields Nathan Paldora,*, Yona Dvorkina, C. Basdevantb a
Ring Department of Atmospheric Sciences, The Hebrew University of Jerusalem, Jerusalem 91904, Israel b Laboratoire de Me!te!orologie Dynamique, CNRS Ecole Normale Supe!rieure, F-75231 Paris, France Received 9 January 2001; received in revised form 16 May 2001; accepted 20 May 2001
Abstract The calculation of particle trajectories in the extra-tropical troposphere is improved by a hybrid model that employs the temperature and geopotential fields to supplement the velocity field. The hybrid model uses the temperature and geopotential fields to construct the Montgomery Stream function, which, together with a Rayleigh friction force and the Coriolis force determine the evolution of a ‘‘correctional velocity’’ based on Newton’s 2nd law of motion. This velocity, however, is decoupled from the continuity equation so its horizontal divergence does not affect the pressure. The improvement of the trajectory calculation is obtained by integrating a linear combination of National Centers for Environmental Predictions’ (NCEP) velocity field and the ‘‘correctional velocity’’ computed from NCEP’s temperature and geopotential fields. The improvement of the model-generated trajectories over those obtained from a straightforward advection by the velocity field is verified by comparing the calculated trajectories to the observed trajectories of 379 constant-level balloons launched in 1971 as part of the EOLE experiment. For flight times between 2 and 10 weeks the new algorithm generates trajectories that are statistically closer to the observed EOLE trajectories than those obtained from advection by the velocity field only. There are, however, several balloon flights where for certain, isolated, values of its parameters the hybrid model generates trajectories that are actually less accurate than those of straightforward advection. For balloon flights between 2 and 9 weeks the worst hybrid model trajectories are only about 4% less accurate than those of straightforward advection. The model’s improvement over advection by the velocity field reaches a maximum value of over 15% for 4–7 week long trajectories while for 1 week long, or longer than 10 week long trajectories the model offers no improvement over straightforward advection by the velocity field. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Lagrangian dispersion; Constant Level Balloons
1. Introduction A hybrid model was recently developed (Dvorkin et al., 2001) to reconstruct four, about 20-day-long, trajectories of three balloons launched in the tropical stratosphere as part of the EQUATEUR experiment *Corresponding author. Fax: +972-2-566-2581. E-mail addresses: [email protected] (N. Paldor), [email protected] (C. Basdevant).
conducted by the Laboratoire de M!et!eorologie Dynamique for the French Space agencyFCNES. This field experiment is a precursor of the STRATEOLE experiment planned for 2006 in the Antarctic stratosphere to look at the breakup of the Polar Vortex, and the subsequent penetration of CFC molecules into it, during spring time. The 4 trajectories calculated for EQUATEUR from the hybrid model were shown to be up to 5 times more accurate than those obtained by advecting the balloons with the ECMWF horizontal
1352-2310/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 2 - 2 3 1 0 ( 0 1 ) 0 0 3 0 5 - 3
484
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
velocity fields. An average improvement of all four trajectories of about 100% was achieved by fitting the values of the only two model parameters. The intent of the present study is to report on an application of the same hybrid model to a much larger set of historic balloons, launched in the extra-tropical troposphere during the EOLE experiment in 1971. The larger data set used here in the confirmation of the model results, makes the improvement of the hybrid model statistically significant and enables a breakdown of the model improvement into weekly bins. Its verification with tropospheric balloons makes the model relevant to large-scale dispersal of particles spewed into the atmosphere (e.g. pollutants). The 2–3 month long time scale relevant to the upcoming STRATEOLE experiment is dictated by both technological constraints of the balloons’ durability in the harsh environment at the polar stratosphere and the time it takes for the polar vortex to collapse. This ambitious length of balloon trajectory differs markedly from other works that have employed a similar approach but for much shorter time scalesFon the order of a few days only. Stohl (1998) reviews these works and shows that the errors in their O (5-day) long trajectories exceed 60% despite their short duration. Obviously, this error level is too large for 2–3 month long flights. The main differences between the present paper and Dvorkin et al. (2001) are: (1) the use of NCEP, instead of ECMWF, fields; (2) the application of the hybrid model to extra-tropical tropospheric trajectories instead of the tropical stratospheric ones and (3) the use of a much larger set of balloons in the validation of the model’s improvement. The paper is organized as follows: The hybrid model is described in Section 2 and the data used for simulating the trajectories and for verifying the improvement of the model generated trajectories are described in Section 3. The improvement of the hybrid model’s trajectories compared with those of advection by the velocity field is presented and verified in Section 4. The paper ends in Section 5 with concluding remarks.
2. The hybrid model The novel element in the formulation of the hybrid model is the use of more than just the velocity field in the advection of a particle. By definition, the trajectory of a particle is a solution of the advection equation: dXp ðtÞ=dt ¼ Vp ðXp ðtÞ; tÞ;
ð1Þ
where Vp ðXp ðtÞ; tÞ and Xp ðtÞ are the particle velocity and location (coordinates), respectively, and t is time. An accurate reconstruction of the trajectory, Xp ðtÞ; requires that the velocity at the particle location be accurately
known at every moment. This is, of course, unrealistic as no velocity is known with 100% accuracy on large spatial scales and for long time on a continuous basis. A crude estimate of the particle velocity is obtained from standard fields produced by large organizations such as ECMWF (European Centre for Medium Range Weather Forecasting) and the US based, National Center for Environmental Prediction (NCEP). The fields calculated (and calibrated with available observations) by these agencies are available to the general public on a standard spatial and temporal grid. One of these fields is Va Fthe air velocity field. Assuming that a particle is advected by the surrounding air a straightforward scheme for forecasting the particle trajectory (the ‘‘pure }}pureadvection’’schemeÞistosubstitutetheknownVaforVppinEq:ð1Þ; i:e:tointegratetheequation\fleqalignno{dXp(t)/ dt=Va(Xp(t), t). The trajectories, Xp ðtÞ; generated from Eq. (2) are usually highly inaccurate compared with observed trajectories. One reason for this inaccuracy is the large uncertainties in the velocity field and the massive interpolation required when evaluating the velocity Va at the particle position, Xp ðtÞ; from its given values on a standard grid. In this regard the modification to Eq. (2) suggested here can be taken as a form parameterization of ‘‘sub-grid’’ scale processes that are not taken into account in the calculation of Va by NCEP. A more fundamental, rather than technical, source of error in the application of Eq. (2) has to do with the fact that the horizontal air velocity, Va ; is calculated in conjunction with the continuity equation. This equation, which relates the horizontal divergence of Va to changes in the pressure, constitutes the main driving force for the changes in Va itself. Obviously, a particle, or a cluster of particles, moving in the atmosphere does not affect the pressure in the surrounding air, so substituting Va for Vp is synonymous to erroneously imposing the continuity equation on particles. Although, strictly speaking, Va is only calculated at the grid points, so the continuity equation is wrongly imposed only there; this is a technical deficiency of the available computing resources only. The desire (and future intent) of any numerical modeler is to calculate the fields on as many points (in space and time) as possible. One can, therefore, expect the mesh of the fields to become finer as our computing capability increases, so the erroneous application of the continuity equation on particle trajectories will occur on more points in space and time. In addition, since the velocity field is interpolated from the grid points to the particle locationFthe error in the imposition of the continuity equation on the grid is interpolated to the interior of the cells. In the hybrid model promulgated here the particle velocity is a linear combination of the air velocity field, Va ; and a ‘‘correctional velocity’’, Vc ; i.e.
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
dXp ðtÞ=dt Vp ¼ ð1 aÞVa þ aVc ;
ð3Þ
where 0pap1 is a measure of the deviation of the particle velocity from the air velocity i.e. a ¼ 0 implies Vp ¼ Va : J.ager (1992) introduced the concept of weighing the velocity field with a ‘‘dynamical’’ velocity. The idea was subsequently applied in Stohl and Seibert (1998) to 6day long trajectories but the calculated trajectories were never compared with observations. Since the particle and air velocities, Vp and Va ; both satisfy Newton’s 2nd law of motion it is clear that Vc ; too, has to satisfy it. Accordingly, the calculation of the evolution of Vc is carried out in full compliance with Newton’s law of motion where the acceleration results from the (vectorial) sum of the pressure gradient force, the Coriolis force and drag i.e. dVc =dt ¼ rC f k#Vc gVc ;
485
instead of the Cartesian form of advection given by Eqs. (2–3). More details on the hybrid model are given in Dvorkin et al. (2001). For certain values of model parameters the trajectory could not be completed because once the calculated balloon arrives at the pole, the calculation cannot be continued due to the singularity of Eqs. (4–5) at f ¼ 7p=2:
3. Data Two types of data used in the present study are briefly described in this section. The first is NCEP fields, used both for running the hybrid model and for advecting the balloons by the airflow, and the second is the EOLE balloon trajectories used for comparing the two.
ð4Þ 3.1. NCEP fields
where C is the Montgomery stream function (calculated on the balloon’s density level), f is the Coriolis parameter (including the geometric terms due to Earth’s spherical surface), k is a unit vector in the vertical direction, g is the Rayleigh friction coefficient and # designates the cross (vector) product. The Montgomery stream function, C, on the right hand side of Eq. (4) is based on the temperature (T) and geopotential (G) fields that are both provided along with Va : Since the balloons fly on isopycnal (not isentropic) surfaces the stream function is given by C ¼ G þ RT where R is the gas constant for air (287 J/(K kg)). The second term on the RHS of Eq. (4) also contains the geometric correction due to Earth’s spherical surface. Eqs. (3) and (4) contain only two parametersFa and gFthat can be varied so as to yield the best fitting trajectories. Compared with the pure advection model, Eq. (2), the hybrid model makes use of the geopotential and temperature fields in addition to the velocity field, Va : In the practical application, the nondimensional form of systems (3 and 4) is used where time is scaled on 1 day/(4p) and length is scaled on the radius of the Earth R ¼ 6400 km. These scales imply that 1=g; which measures the velocity relaxation time (i.e. the rate of exponential decay once the pressure gradient vanishes) is measured on 24 h/ð4pÞE1:9 h. The commensurate velocity scale is (4pR km)/(24 h)E930 m/s, so a realistic nondimensional velocity is on the order of 0.01. The coordinates used in the global scale problem considered here are the particle’s latitude, f; and longitude, l; that evolve in (nondimensional) time due to the eastward (up ) and northward (vp ) components of Vp according to df=dt ¼ vp ¼ ð1 aÞva þ avc ; dl=dt ¼ up =cosðfÞ ¼ ½ð1 aÞua þ auc =cosðfÞ;
ð5Þ
The United States NCEP issues velocity, temperature and geopotential fields on a global scale (with a resolution of 2.51 longitude by 2.51 latitude) on 17 pressure levels every 6 h. The NCEP operational model assimilates meteorological observations from a network of global stations and the model output is then calibrated so as to minimize the model’s errors relative to the available data. The fields from different times are reanalyzed to ensure homogeneity of accuracy and precision. The first step in using NCEP fields in our hybrid model is to interpolate the fields to the particle’s isopycnal surface using cubic spline interpolation in log(pressure) of their values on isobaric surfaces. The gridded data on the relevant isopycnic surface (0.33 Kg/ m3 for the EOLE balloons, see below) are interpolated horizontally to the particle’s exact location at the prescribed time using 12-point bi-cubic spline interpolation (3 grid points to the North, East, South and West of the particle’s location). Interpolation in time between the embedding 6-h intervals is linear. Experiments done with other, more sophisticated, interpolation scheme (e.g. 24-point in space, cubic in time) showed no significant difference in the resulting trajectory. More details on the reanalyzed NCEP fields are given in Kalnay et al. (1996). 3.2. EOLE balloon trajectories The 480 EOLE balloons were released from three launching sites in the Southern Hemisphere during July– November 1971. The balloons were ballasted to fly along the 0.33 kg/m3 density surface (i.e. nearly the 11,700 m or 205 hPa levels) and had an average lifetime of 103 days. The daily averaged height deviation (i.e. the
486
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
Fig. 1. A histogram of the number of EOLE balloons participating in each weekly bin. The bins attempted here span a trajectory length of order 1–15 weeks.
vertical motion not related to the daily heating) was on the order of 2 m only while 20 m daily vertical excursions were recorded due to the daily heating cycle. Only two balloons left the Southern Hemisphere altogether at some point in the course of their flight and crossed the equator. Several other balloons have visited the latitude band north of 101S (but south of 5 degrees North) while the vast majority of the balloons have remained in the extra-tropical Southern Hemisphere for the entire duration of their flight. More details on the EOLE experiments are given in Morel and Bandeen, (1973) and Morel and Larcheveque, (1974). Of all the EOLE balloons 379 balloons have recorded data for longer than 1 week and all these trajectories are included in the present study. The trajectories considered for comparison between the hybrid model and the advection by the surrounding air were arranged into weekly bins ranging in length from 1 to 15 weeks. As expected, the number of trajectories used for verifying the hybrid model results decreases monotonically with the increase in the trajectory length (i.e. bin number). While 379 balloons were used in the comparison of 1 weeklong trajectories only 164 were available in the 15 weeklong ones. The histogram of the number of balloons used in each bin of weekly spaced trajectory length is shown in Fig. 1.
4. Results For each of the EOLE balloon trajectories a pure advection trajectory was calculated based on the NCEP’s velocity field (i.e. an integration of Eq. (5) with a ¼ 0) starting from the balloon’s initial location and
velocity (estimated from the distance between the 2 first points divided by the time-difference). A set of 567 additional trajectories was calculated from the hybrid model for the 21 a-values (0.01 and 0.05–1.0 every 0.05) by 27 g-values (0.01 and 0.05–1.3 every 0.05). The minimal value of g ¼ 0:01 corresponds to a relaxation time of 8 days, and the maximal value of 1.3 to a relaxation time of less than 1.5 h so values outside this range need not be considered. The value of a ¼ 0:01 was added to the 0.05 grid in order to guarantee that values smaller than 0.05 are not overlooked. The hybrid model equations Eqs. (4) and (5) and the advection (i.e. a ¼ 0) model were both integrated in time using a Runge-Kutta-Verner 5th order scheme with variable time step that adjusts to yield an accuracy of 104. Time series of position (i.e. coordinates) and velocity components were issued every 0.1 nondimensional time units, i.e. 0:6 h=pE11 min. The observed EOLE trajectories were also interpolated from the original time series, that had a variable temporal spacing, to the same, 0.1 nondimensional time units, equidistant time series. The fit between each of these trajectories (be it one of the hybrid model trajectories or the pure advection one) and the observed, EOLE, trajectory was quantified by computing, at each weekly bin (e.g. week J) the mean separation between the observed and calculated trajectories of balloon #i at week J M i ðJ; a; gÞ ¼ ð1=KÞ
K X
d i ðJ; a; gÞ:
ð6Þ
j¼1
In this expression d i ðj; a; gÞ is the length of the arc (on a unit sphere) between the observed and calculated balloons in flight #i at time jdt and K ¼ 280pJE880 J is the number of time steps (i.e. 10 times the number of nondimensional time units) in J weeks. For week J the average of M i ðJ; a; gÞ over all participating balloons (i.e. average over i), MðJ; a; gÞ; was calculated for comparison with MðJ; a ¼ 0Þ; which is the average of the pure advection’s mean separation at week J: For any particular pair of ða; gÞ values, and for each of the 15 weekly bins (i.e. J in Eq. (6) between 1 and 15) the Improvement Factor (IF) (J; a; g), of the hybrid model was calculated as follows: The sum of mean separation after J weeks, over all relevant balloons, of the hybrid model’s trajectories was calculated from Eq. (6) using the particular (a; g) values. This sum was divided into the sum, over all balloons, of the mean separation at the same bin of the pure advection model. The explicit expression of the improvement factor, IF(J; a; g), is IFðJ; a; gÞ ¼ MðJ; a ¼ 0Þ=MðJ; a; gÞ:
ð7Þ
A contour map of the IF at 4 weeks is shown in Fig. 2. The maximal IF value of 1.15 occurs at a ¼ 0:85; g ¼
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
487
Fig. 2. A contour map of the improvement factor (IF) of the hybrid model in the (a, g) plane for J ¼ 4 i.e. 4 weeklong trajectories. Closed contours of IF=1.14 occur in several ‘‘islands’’ in this plane. The maximal IF value of 1.151 occurs at ða; gÞ ¼ ð0:85; 0:85Þ but close IF values occur at other (a; g) pairs (e.g. IF=1.149 near ða; gÞ ¼ ð0:50; 0:25Þ). At other weekly bins the IF contours have the same hyperbola-like general structure but the values of the contours are different. Notice that there are no contours with IFo1.0.
0:85 and very similar IF values also occur at the same (a; g) pairs in weeks 4–7. There are no IF contours smaller than 1.0 in Fig. 2 but there are 5 out of the 567! (a; g) pairs where M(4; a; g) is slightly larger than M(4; a ¼ 0) and min{IF}=0.99. These ‘‘bad’’ mean separations generated by the hybrid model occur at very small values of g (0.01Fa relaxation time of 8 days!Fand 0.05) and a (o0.15). The hyperbola-like structure of the IF contours in Fig. 2 is typical of all weekly bins but the IF values vary with J: The main results regarding the statistics of best (a; g) pairs and (best/worst) IF values at all weekly bins are summarized in Table 1. For 2pJp4 there are no IF values smaller than 1.0 but at J ¼ 1 and 10 the minimal values of IF deviate from 1.0 more than its maximal value. In all bins between 2 and 10 weeks the number of (a; g) pairs with IF>1.0 is greater than the number of pairs with IF o1.0 while for J ¼ 1 only 175 pairs (out of 567) had IF>1.0. For J > 10 no improvement obtains from the application of the hybrid model as is evident from the results shown in Table 1. The dependence of the best IF value on J is shown in Fig. 3. The best IF exceeds 1.0 at all 10 weekly bins but
IF>1.1 only at 2–9 weeks. It reaches a maximal value of 1.15 for 4–7 weeklong trajectories (1.14 for J ¼ 3 weeks) and it decreases to values of about 1.05 for J ¼ 1 and 1.08 for J ¼ 10: These IF values were obtained by searching for the maximal IF values over all 567 (a; g) pairs at weekly bin J; so the (a; g) pairs that yield the IF values shown in Fig. 3 for different J are not necessarily the same. Many pairs of (a; g) yield an improvement over pure advection, i.e. IF(J; a; g)>1.0, for all 1pJp10: One such pair is ða; gÞ ¼ ð0:85; 0:85Þ and the IF (J; 0.85, 0.85) function is shown in Fig. 4. The IF values (though slightly smaller than those in Fig. 3) significantly exceed 1.0 for all weeks considered. The pair ða; gÞ ¼ ð0:85; 0:85Þ also yields the PLlargest sum of weekly improvement factors (i.e. J=1 IF (J; a, g) is maximal for ða; gÞ ¼ ð0:85; 0:85Þ) for any Lp12:
5. Concluding remarks It should be emphasized that despite the ubiquitous improvement that the ‘‘correctional velocity’’ offers
488
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
Table 1 Summary of results of best fitting a (2nd column) and g (3rd column) pairs and best/worst IF factors (4th column) as a function of J (1st column). The 5th column shows the number of (a; g) pairs (out of 567) that yield IF>1.0. The improvement for 11pJp14 is insignificant as (1.0Fworst{IF})>(best{IF}F1.0) and the number of (a; g) pairs with IFo1.0 exceeds that with IF>1.0. No data is available for J ¼ 15 since too many (a; g) pairs did not have even one balloon lasting the entire 15 weeks flight before terminating at the ‘‘sink’’ at the pole J (weekly bin)
Best a
Best g
Best/worst IF
# of (a; g) pairs with IF>1.0 (out of 567)
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.70 0.80 0.50 0.85 0.85 0.85 0.85 1.00 0.85 0.85 0.20 0.20 0.20 0.20
1.20 0.95 0.70 0.85 0.85 0.85 0.85 0.25 0.30 0.30 1.30 1.30 1.30 1.30
1.05/0.95 1.11/1.00 1.14/1.00 1.15/1.00 1.15/0.99 1.15/0.99 1.15/0.99 1.12/0.98 1.10/0.96 1.08/0.91 1.07/0.89 1.06/0.86 1.06/0.84 1.05/0.81
175 561 562 562 549 548 549 538 490 388 255 120 78 47
Fig. 3. The maximal IF values for the different weekly bins. The values of a and g that yield the particular maximal IF(J) vary with J (the weekly bin). No data is presented for J > 10 as the best/worst IF outweigh one another and the improvement is insignificant.
Fig. 4. The improvement factor of the hybrid model as a function of weekly bin for ða; gÞ ¼ ð0:85; 0:85Þ: IF(J; 0.85, 0.85). Other (a; g) pairs exist that yield similar IF(J; a; g). For J > 10 the improvement is insignificant.
when combined with the velocity field, it has no physical meaning in itself and should only be used to properly incorporate the geopotential and temperature fields into the velocity field. To demonstrate this point we need to look at the calculated and observed trajectories of a particular balloon. In Fig. 5 we show the various trajectories relevant to the first 2 weeks of EOLE balloon #183. The IF of the best fitting hybrid model’s trajectory is 9.30 and is obtained for a ¼ 0:35 and g ¼ 1:3: Both the observed trajectory (solid curve, upper panel) and the best fitting one generated by the hybrid model (dotted curve, upper panel) span the observed
B1001 of longitude and B151 of latitude. By comparison, the trajectory that results from pure advection (solid curve, lower panel) spans about 501 of latitude and more than 3001 of longitude! It is clear, therefore, that pure advection by the NCEP velocity field provides a very poor descriptor of the observed trajectory. The trajectory due to the ‘‘correctional velocity’’ alone, too, bears no resemblance, whatsoever, to the observed trajectory (dotted curve, lower panel). However, a linear combination of the two corresponding (notice that the velocities shown are not the ones actually used in Eq. (5) since the balloon’s coordinates determine the actual
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
489
Fig. 5. Top panelFthe observed trajectory of EOLE balloon #183 (solid curve) and its best fitting hybrid model calculated trajectory (dotted curve) after 2 weeks. The improvement factor of the best fitting hybrid model trajectory is 930%! Lower panelFtrajectories due to pure advection by NCEP velocity field (solid curve) and to ‘‘correctional velocity’’ alone (dotted curve). Each of the two trajectories in the lower panel, alone provides a very poor estimate of the real trajectory but a proper combination of the two corresponding velocities yields an accurate one. Note the ‘‘inertial’’ loops in the dynamical model (a ¼ 1:0) trajectory.
velocities used there) velocities with a ¼ 0:35 offers a much better reproduction of the observed trajectory (i.e. IF=9.30). The better fit with the observed trajectory is achieved by properly employing NCEP’s velocity, geopotential and temperature fields instead of the velocity field alone. From the outset, the undertaking of this work of forecasting balloon trajectories for periods of over 2 months is an ambitious one. To explain why the hybrid model performs differently on different time scales (i.e. IF varies with J) we note that both Morel and Bandeen (1973) and in Morel and Larcheveque (1974) report on a dramatic change with time of the dispersal rates of EOLE balloons. During the first week the balloons disperse very fastFthe e-folding dispersal time is about 1.3 dayFafter which time the balloons reach a stage of nearly uniform distribution. The saturated stage occurs when the mean pair separation between balloons is about 2000 km and the balloons are distributed all over the extra-tropical latitudes of the Southern Hemisphere. Given the fast initial dispersal, it can be expected that accurate trajectory reconstruction will be difficult to achieve in the first week. By comparison, in the saturated stage the statistical properties of the individual balloon trajectories
remain constant with time so that a satisfactory reconstruction of these trajectories is possible. The dependence of our results, as measured by the improvement factor, on the length of the trajectories is in accordance with the observations of the EOLE balloons described above. At all weekly bins between 2 and 10 the hybrid model provides improved (statistically) trajectories compared with the pure advection ones. This is manifested both in the number of (a; g) pairs where IF>1.0 and in its maximal value being larger than 1.0 by up to 15% while its minimal value being only a few percent below 1.0. In contrast, for 1 weeklong trajectories (as well as for trajectories shorter than 1 week or longer than 10 weeksFresults not shown here) the hybrid model does not provide any improvement (in a statistically significant way) over the pure advection by the velocity field. The minimal value of IF deviates from 1.0 more than its maximal value and the number of (a; g) pairs where IF>1.0 is smaller than that with IFo1.0. The improvement of the present hybrid model compared with the pure advection model is pronounced at a values that are sufficiently large. When a is too small the hybrid model does not provide any significant improvementFsince, by definition, IF(a ¼ 0)=1.0. This
490
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
expected result is very conspicuous in the numerical computations shown in Fig. 2. In comparison, in the midtropospheric, ‘‘blended’’ model of Stohl and Seibert (1998) an equation, similar to Eq. (3) (suggested by J.ager; 1992), was employed but with time-dependent weightFa: The value of the weight (i.e. a) in their model was determined at each time step by the ratio between the integration time step, which is determined dynamically by the CFL condition, and the time interval of the velocity fields. Since accurate integration requires that the former time steps be small, their results should be compared to the small a limit of the present study. As it turns out, the two models agree with each otherFStohl and Seibert (1998) found that their blended model performs worse than any other 2D kinematic (e.g. pure advection) model they tried while Fig. 2 indicates that no improvement occurs in the hybrid model when a is too small. An additional point of comparison between the present model and the blended model of Stohl and Seibert (1998) is the small J limit of Fig. 3 (see also Table 1). Here it is clear that for J ¼ 1 the present hybrid model offers no improvement over the pure advection model. The Stohl and Seibert (1998) model was applied to 6-day long trajectories and, as in our case, the resulting trajectories were not better than any 2D kinematic model trajectories. In the review of Stohl (1998) an inter-model comparison and error analysis was done but a study of the model performance against observations was not attempted. By comparison, in the present study the hybrid model’s improvement is verified by comparing the calculated trajectories with the observations. This puts the conclusions reached in this study on much firmer grounds. On time scale between 2 and 9 weeks our results indicate that the improvement factor of the hybrid model exceeds 1.1. The reason for the drop in the IF values at longer times is not clear. However, the maximum value of the IF periods around 4–7 weeks makes the hybrid model particularly efficient for the prediction of trajectories on the time scales of 1–2 months, which is relevant to the STRATEOLE experiment.
Acknowledgements Financial support for this study was provided by the US-Israel Binational Science via a research grant to the Hebrew University of Jerusalem and by the European Sciences Foundation, which provided YD with travel funds allowing him to spend a summer at LMD where the EOLE trajectories and NCEP fields were processed. Discussions held with F. Vial were very helpful in clarifying some of the subtle issues involved with EOLE trajectories and M. Pinsky, of the Hebrew University, provided helpful programming advice.
References Dvorkin, Y., Paldor, N., Basdevant, C., 2001. Modeling balloon trajectories in the tropical stratosphere with a hybrid model using analyzed fields. Quarterly Journal of Royal Meteorological Society 127 (573A), 975–988. J.ager, A., 1992. Isentrope trajektorien und ihre Anwendung auf das Konzept der potentiellen Vorticity bei orograpisch induzierten Lee-Zyklogenesen. Diploma Thesis, University of Innsbruck, Austria. Kalnay, E.M., Kanamitsu, R., Kistler, W., Collins, D., Deaven, L., Gandin, M., Iredell, S., Saha, G., White, J., Woollen, Y., Zhu, M., Chelliah, W., Ebisuzaki, W., Higgins, J., Janowiak, K.C., Mo, C., Ropelewski, J., Wang, A., Leetmaa, R., Reynolds, R., Jenne and Joseph, D., 1996. The NCEP/NCAR 40-year reanalysis project. Bulletin of American Meteorological Society 77, 437–471. Morel, P., Bandeen, W., 1973. The EOLE experiment: early results and current objectives. Bulletin of American Meteorological Society 54, 298–306. Morel, P., Larcheveque, M., 1974. Relative dispersion of constant-level balloons in the 200 mb general circulation. Journal of Atmospheric Sciences 31, 2189–2196. Stohl, A., 1998. Computation, accuracy and applications of trajectoriesFa review and bibliography. Atmospheric Environment 32, 947–966. Stohl, A., Seibert, P., 1998. Accuracy of trajectories as determined from the conservation of meteorological tracers. Quarterly Journal of Royal Meteorological Society 124, 1465–1484.
Improving the calculation of particle trajectories in the extra-tropical troposphere using standard NCEP fields Nathan Paldora,*, Yona Dvorkina, C. Basdevantb a
Ring Department of Atmospheric Sciences, The Hebrew University of Jerusalem, Jerusalem 91904, Israel b Laboratoire de Me!te!orologie Dynamique, CNRS Ecole Normale Supe!rieure, F-75231 Paris, France Received 9 January 2001; received in revised form 16 May 2001; accepted 20 May 2001
Abstract The calculation of particle trajectories in the extra-tropical troposphere is improved by a hybrid model that employs the temperature and geopotential fields to supplement the velocity field. The hybrid model uses the temperature and geopotential fields to construct the Montgomery Stream function, which, together with a Rayleigh friction force and the Coriolis force determine the evolution of a ‘‘correctional velocity’’ based on Newton’s 2nd law of motion. This velocity, however, is decoupled from the continuity equation so its horizontal divergence does not affect the pressure. The improvement of the trajectory calculation is obtained by integrating a linear combination of National Centers for Environmental Predictions’ (NCEP) velocity field and the ‘‘correctional velocity’’ computed from NCEP’s temperature and geopotential fields. The improvement of the model-generated trajectories over those obtained from a straightforward advection by the velocity field is verified by comparing the calculated trajectories to the observed trajectories of 379 constant-level balloons launched in 1971 as part of the EOLE experiment. For flight times between 2 and 10 weeks the new algorithm generates trajectories that are statistically closer to the observed EOLE trajectories than those obtained from advection by the velocity field only. There are, however, several balloon flights where for certain, isolated, values of its parameters the hybrid model generates trajectories that are actually less accurate than those of straightforward advection. For balloon flights between 2 and 9 weeks the worst hybrid model trajectories are only about 4% less accurate than those of straightforward advection. The model’s improvement over advection by the velocity field reaches a maximum value of over 15% for 4–7 week long trajectories while for 1 week long, or longer than 10 week long trajectories the model offers no improvement over straightforward advection by the velocity field. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Lagrangian dispersion; Constant Level Balloons
1. Introduction A hybrid model was recently developed (Dvorkin et al., 2001) to reconstruct four, about 20-day-long, trajectories of three balloons launched in the tropical stratosphere as part of the EQUATEUR experiment *Corresponding author. Fax: +972-2-566-2581. E-mail addresses: [email protected] (N. Paldor), [email protected] (C. Basdevant).
conducted by the Laboratoire de M!et!eorologie Dynamique for the French Space agencyFCNES. This field experiment is a precursor of the STRATEOLE experiment planned for 2006 in the Antarctic stratosphere to look at the breakup of the Polar Vortex, and the subsequent penetration of CFC molecules into it, during spring time. The 4 trajectories calculated for EQUATEUR from the hybrid model were shown to be up to 5 times more accurate than those obtained by advecting the balloons with the ECMWF horizontal
1352-2310/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 2 - 2 3 1 0 ( 0 1 ) 0 0 3 0 5 - 3
484
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
velocity fields. An average improvement of all four trajectories of about 100% was achieved by fitting the values of the only two model parameters. The intent of the present study is to report on an application of the same hybrid model to a much larger set of historic balloons, launched in the extra-tropical troposphere during the EOLE experiment in 1971. The larger data set used here in the confirmation of the model results, makes the improvement of the hybrid model statistically significant and enables a breakdown of the model improvement into weekly bins. Its verification with tropospheric balloons makes the model relevant to large-scale dispersal of particles spewed into the atmosphere (e.g. pollutants). The 2–3 month long time scale relevant to the upcoming STRATEOLE experiment is dictated by both technological constraints of the balloons’ durability in the harsh environment at the polar stratosphere and the time it takes for the polar vortex to collapse. This ambitious length of balloon trajectory differs markedly from other works that have employed a similar approach but for much shorter time scalesFon the order of a few days only. Stohl (1998) reviews these works and shows that the errors in their O (5-day) long trajectories exceed 60% despite their short duration. Obviously, this error level is too large for 2–3 month long flights. The main differences between the present paper and Dvorkin et al. (2001) are: (1) the use of NCEP, instead of ECMWF, fields; (2) the application of the hybrid model to extra-tropical tropospheric trajectories instead of the tropical stratospheric ones and (3) the use of a much larger set of balloons in the validation of the model’s improvement. The paper is organized as follows: The hybrid model is described in Section 2 and the data used for simulating the trajectories and for verifying the improvement of the model generated trajectories are described in Section 3. The improvement of the hybrid model’s trajectories compared with those of advection by the velocity field is presented and verified in Section 4. The paper ends in Section 5 with concluding remarks.
2. The hybrid model The novel element in the formulation of the hybrid model is the use of more than just the velocity field in the advection of a particle. By definition, the trajectory of a particle is a solution of the advection equation: dXp ðtÞ=dt ¼ Vp ðXp ðtÞ; tÞ;
ð1Þ
where Vp ðXp ðtÞ; tÞ and Xp ðtÞ are the particle velocity and location (coordinates), respectively, and t is time. An accurate reconstruction of the trajectory, Xp ðtÞ; requires that the velocity at the particle location be accurately
known at every moment. This is, of course, unrealistic as no velocity is known with 100% accuracy on large spatial scales and for long time on a continuous basis. A crude estimate of the particle velocity is obtained from standard fields produced by large organizations such as ECMWF (European Centre for Medium Range Weather Forecasting) and the US based, National Center for Environmental Prediction (NCEP). The fields calculated (and calibrated with available observations) by these agencies are available to the general public on a standard spatial and temporal grid. One of these fields is Va Fthe air velocity field. Assuming that a particle is advected by the surrounding air a straightforward scheme for forecasting the particle trajectory (the ‘‘pure }}pureadvection’’schemeÞistosubstitutetheknownVaforVppinEq:ð1Þ; i:e:tointegratetheequation\fleqalignno{dXp(t)/ dt=Va(Xp(t), t). The trajectories, Xp ðtÞ; generated from Eq. (2) are usually highly inaccurate compared with observed trajectories. One reason for this inaccuracy is the large uncertainties in the velocity field and the massive interpolation required when evaluating the velocity Va at the particle position, Xp ðtÞ; from its given values on a standard grid. In this regard the modification to Eq. (2) suggested here can be taken as a form parameterization of ‘‘sub-grid’’ scale processes that are not taken into account in the calculation of Va by NCEP. A more fundamental, rather than technical, source of error in the application of Eq. (2) has to do with the fact that the horizontal air velocity, Va ; is calculated in conjunction with the continuity equation. This equation, which relates the horizontal divergence of Va to changes in the pressure, constitutes the main driving force for the changes in Va itself. Obviously, a particle, or a cluster of particles, moving in the atmosphere does not affect the pressure in the surrounding air, so substituting Va for Vp is synonymous to erroneously imposing the continuity equation on particles. Although, strictly speaking, Va is only calculated at the grid points, so the continuity equation is wrongly imposed only there; this is a technical deficiency of the available computing resources only. The desire (and future intent) of any numerical modeler is to calculate the fields on as many points (in space and time) as possible. One can, therefore, expect the mesh of the fields to become finer as our computing capability increases, so the erroneous application of the continuity equation on particle trajectories will occur on more points in space and time. In addition, since the velocity field is interpolated from the grid points to the particle locationFthe error in the imposition of the continuity equation on the grid is interpolated to the interior of the cells. In the hybrid model promulgated here the particle velocity is a linear combination of the air velocity field, Va ; and a ‘‘correctional velocity’’, Vc ; i.e.
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
dXp ðtÞ=dt Vp ¼ ð1 aÞVa þ aVc ;
ð3Þ
where 0pap1 is a measure of the deviation of the particle velocity from the air velocity i.e. a ¼ 0 implies Vp ¼ Va : J.ager (1992) introduced the concept of weighing the velocity field with a ‘‘dynamical’’ velocity. The idea was subsequently applied in Stohl and Seibert (1998) to 6day long trajectories but the calculated trajectories were never compared with observations. Since the particle and air velocities, Vp and Va ; both satisfy Newton’s 2nd law of motion it is clear that Vc ; too, has to satisfy it. Accordingly, the calculation of the evolution of Vc is carried out in full compliance with Newton’s law of motion where the acceleration results from the (vectorial) sum of the pressure gradient force, the Coriolis force and drag i.e. dVc =dt ¼ rC f k#Vc gVc ;
485
instead of the Cartesian form of advection given by Eqs. (2–3). More details on the hybrid model are given in Dvorkin et al. (2001). For certain values of model parameters the trajectory could not be completed because once the calculated balloon arrives at the pole, the calculation cannot be continued due to the singularity of Eqs. (4–5) at f ¼ 7p=2:
3. Data Two types of data used in the present study are briefly described in this section. The first is NCEP fields, used both for running the hybrid model and for advecting the balloons by the airflow, and the second is the EOLE balloon trajectories used for comparing the two.
ð4Þ 3.1. NCEP fields
where C is the Montgomery stream function (calculated on the balloon’s density level), f is the Coriolis parameter (including the geometric terms due to Earth’s spherical surface), k is a unit vector in the vertical direction, g is the Rayleigh friction coefficient and # designates the cross (vector) product. The Montgomery stream function, C, on the right hand side of Eq. (4) is based on the temperature (T) and geopotential (G) fields that are both provided along with Va : Since the balloons fly on isopycnal (not isentropic) surfaces the stream function is given by C ¼ G þ RT where R is the gas constant for air (287 J/(K kg)). The second term on the RHS of Eq. (4) also contains the geometric correction due to Earth’s spherical surface. Eqs. (3) and (4) contain only two parametersFa and gFthat can be varied so as to yield the best fitting trajectories. Compared with the pure advection model, Eq. (2), the hybrid model makes use of the geopotential and temperature fields in addition to the velocity field, Va : In the practical application, the nondimensional form of systems (3 and 4) is used where time is scaled on 1 day/(4p) and length is scaled on the radius of the Earth R ¼ 6400 km. These scales imply that 1=g; which measures the velocity relaxation time (i.e. the rate of exponential decay once the pressure gradient vanishes) is measured on 24 h/ð4pÞE1:9 h. The commensurate velocity scale is (4pR km)/(24 h)E930 m/s, so a realistic nondimensional velocity is on the order of 0.01. The coordinates used in the global scale problem considered here are the particle’s latitude, f; and longitude, l; that evolve in (nondimensional) time due to the eastward (up ) and northward (vp ) components of Vp according to df=dt ¼ vp ¼ ð1 aÞva þ avc ; dl=dt ¼ up =cosðfÞ ¼ ½ð1 aÞua þ auc =cosðfÞ;
ð5Þ
The United States NCEP issues velocity, temperature and geopotential fields on a global scale (with a resolution of 2.51 longitude by 2.51 latitude) on 17 pressure levels every 6 h. The NCEP operational model assimilates meteorological observations from a network of global stations and the model output is then calibrated so as to minimize the model’s errors relative to the available data. The fields from different times are reanalyzed to ensure homogeneity of accuracy and precision. The first step in using NCEP fields in our hybrid model is to interpolate the fields to the particle’s isopycnal surface using cubic spline interpolation in log(pressure) of their values on isobaric surfaces. The gridded data on the relevant isopycnic surface (0.33 Kg/ m3 for the EOLE balloons, see below) are interpolated horizontally to the particle’s exact location at the prescribed time using 12-point bi-cubic spline interpolation (3 grid points to the North, East, South and West of the particle’s location). Interpolation in time between the embedding 6-h intervals is linear. Experiments done with other, more sophisticated, interpolation scheme (e.g. 24-point in space, cubic in time) showed no significant difference in the resulting trajectory. More details on the reanalyzed NCEP fields are given in Kalnay et al. (1996). 3.2. EOLE balloon trajectories The 480 EOLE balloons were released from three launching sites in the Southern Hemisphere during July– November 1971. The balloons were ballasted to fly along the 0.33 kg/m3 density surface (i.e. nearly the 11,700 m or 205 hPa levels) and had an average lifetime of 103 days. The daily averaged height deviation (i.e. the
486
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
Fig. 1. A histogram of the number of EOLE balloons participating in each weekly bin. The bins attempted here span a trajectory length of order 1–15 weeks.
vertical motion not related to the daily heating) was on the order of 2 m only while 20 m daily vertical excursions were recorded due to the daily heating cycle. Only two balloons left the Southern Hemisphere altogether at some point in the course of their flight and crossed the equator. Several other balloons have visited the latitude band north of 101S (but south of 5 degrees North) while the vast majority of the balloons have remained in the extra-tropical Southern Hemisphere for the entire duration of their flight. More details on the EOLE experiments are given in Morel and Bandeen, (1973) and Morel and Larcheveque, (1974). Of all the EOLE balloons 379 balloons have recorded data for longer than 1 week and all these trajectories are included in the present study. The trajectories considered for comparison between the hybrid model and the advection by the surrounding air were arranged into weekly bins ranging in length from 1 to 15 weeks. As expected, the number of trajectories used for verifying the hybrid model results decreases monotonically with the increase in the trajectory length (i.e. bin number). While 379 balloons were used in the comparison of 1 weeklong trajectories only 164 were available in the 15 weeklong ones. The histogram of the number of balloons used in each bin of weekly spaced trajectory length is shown in Fig. 1.
4. Results For each of the EOLE balloon trajectories a pure advection trajectory was calculated based on the NCEP’s velocity field (i.e. an integration of Eq. (5) with a ¼ 0) starting from the balloon’s initial location and
velocity (estimated from the distance between the 2 first points divided by the time-difference). A set of 567 additional trajectories was calculated from the hybrid model for the 21 a-values (0.01 and 0.05–1.0 every 0.05) by 27 g-values (0.01 and 0.05–1.3 every 0.05). The minimal value of g ¼ 0:01 corresponds to a relaxation time of 8 days, and the maximal value of 1.3 to a relaxation time of less than 1.5 h so values outside this range need not be considered. The value of a ¼ 0:01 was added to the 0.05 grid in order to guarantee that values smaller than 0.05 are not overlooked. The hybrid model equations Eqs. (4) and (5) and the advection (i.e. a ¼ 0) model were both integrated in time using a Runge-Kutta-Verner 5th order scheme with variable time step that adjusts to yield an accuracy of 104. Time series of position (i.e. coordinates) and velocity components were issued every 0.1 nondimensional time units, i.e. 0:6 h=pE11 min. The observed EOLE trajectories were also interpolated from the original time series, that had a variable temporal spacing, to the same, 0.1 nondimensional time units, equidistant time series. The fit between each of these trajectories (be it one of the hybrid model trajectories or the pure advection one) and the observed, EOLE, trajectory was quantified by computing, at each weekly bin (e.g. week J) the mean separation between the observed and calculated trajectories of balloon #i at week J M i ðJ; a; gÞ ¼ ð1=KÞ
K X
d i ðJ; a; gÞ:
ð6Þ
j¼1
In this expression d i ðj; a; gÞ is the length of the arc (on a unit sphere) between the observed and calculated balloons in flight #i at time jdt and K ¼ 280pJE880 J is the number of time steps (i.e. 10 times the number of nondimensional time units) in J weeks. For week J the average of M i ðJ; a; gÞ over all participating balloons (i.e. average over i), MðJ; a; gÞ; was calculated for comparison with MðJ; a ¼ 0Þ; which is the average of the pure advection’s mean separation at week J: For any particular pair of ða; gÞ values, and for each of the 15 weekly bins (i.e. J in Eq. (6) between 1 and 15) the Improvement Factor (IF) (J; a; g), of the hybrid model was calculated as follows: The sum of mean separation after J weeks, over all relevant balloons, of the hybrid model’s trajectories was calculated from Eq. (6) using the particular (a; g) values. This sum was divided into the sum, over all balloons, of the mean separation at the same bin of the pure advection model. The explicit expression of the improvement factor, IF(J; a; g), is IFðJ; a; gÞ ¼ MðJ; a ¼ 0Þ=MðJ; a; gÞ:
ð7Þ
A contour map of the IF at 4 weeks is shown in Fig. 2. The maximal IF value of 1.15 occurs at a ¼ 0:85; g ¼
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
487
Fig. 2. A contour map of the improvement factor (IF) of the hybrid model in the (a, g) plane for J ¼ 4 i.e. 4 weeklong trajectories. Closed contours of IF=1.14 occur in several ‘‘islands’’ in this plane. The maximal IF value of 1.151 occurs at ða; gÞ ¼ ð0:85; 0:85Þ but close IF values occur at other (a; g) pairs (e.g. IF=1.149 near ða; gÞ ¼ ð0:50; 0:25Þ). At other weekly bins the IF contours have the same hyperbola-like general structure but the values of the contours are different. Notice that there are no contours with IFo1.0.
0:85 and very similar IF values also occur at the same (a; g) pairs in weeks 4–7. There are no IF contours smaller than 1.0 in Fig. 2 but there are 5 out of the 567! (a; g) pairs where M(4; a; g) is slightly larger than M(4; a ¼ 0) and min{IF}=0.99. These ‘‘bad’’ mean separations generated by the hybrid model occur at very small values of g (0.01Fa relaxation time of 8 days!Fand 0.05) and a (o0.15). The hyperbola-like structure of the IF contours in Fig. 2 is typical of all weekly bins but the IF values vary with J: The main results regarding the statistics of best (a; g) pairs and (best/worst) IF values at all weekly bins are summarized in Table 1. For 2pJp4 there are no IF values smaller than 1.0 but at J ¼ 1 and 10 the minimal values of IF deviate from 1.0 more than its maximal value. In all bins between 2 and 10 weeks the number of (a; g) pairs with IF>1.0 is greater than the number of pairs with IF o1.0 while for J ¼ 1 only 175 pairs (out of 567) had IF>1.0. For J > 10 no improvement obtains from the application of the hybrid model as is evident from the results shown in Table 1. The dependence of the best IF value on J is shown in Fig. 3. The best IF exceeds 1.0 at all 10 weekly bins but
IF>1.1 only at 2–9 weeks. It reaches a maximal value of 1.15 for 4–7 weeklong trajectories (1.14 for J ¼ 3 weeks) and it decreases to values of about 1.05 for J ¼ 1 and 1.08 for J ¼ 10: These IF values were obtained by searching for the maximal IF values over all 567 (a; g) pairs at weekly bin J; so the (a; g) pairs that yield the IF values shown in Fig. 3 for different J are not necessarily the same. Many pairs of (a; g) yield an improvement over pure advection, i.e. IF(J; a; g)>1.0, for all 1pJp10: One such pair is ða; gÞ ¼ ð0:85; 0:85Þ and the IF (J; 0.85, 0.85) function is shown in Fig. 4. The IF values (though slightly smaller than those in Fig. 3) significantly exceed 1.0 for all weeks considered. The pair ða; gÞ ¼ ð0:85; 0:85Þ also yields the PLlargest sum of weekly improvement factors (i.e. J=1 IF (J; a, g) is maximal for ða; gÞ ¼ ð0:85; 0:85Þ) for any Lp12:
5. Concluding remarks It should be emphasized that despite the ubiquitous improvement that the ‘‘correctional velocity’’ offers
488
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
Table 1 Summary of results of best fitting a (2nd column) and g (3rd column) pairs and best/worst IF factors (4th column) as a function of J (1st column). The 5th column shows the number of (a; g) pairs (out of 567) that yield IF>1.0. The improvement for 11pJp14 is insignificant as (1.0Fworst{IF})>(best{IF}F1.0) and the number of (a; g) pairs with IFo1.0 exceeds that with IF>1.0. No data is available for J ¼ 15 since too many (a; g) pairs did not have even one balloon lasting the entire 15 weeks flight before terminating at the ‘‘sink’’ at the pole J (weekly bin)
Best a
Best g
Best/worst IF
# of (a; g) pairs with IF>1.0 (out of 567)
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.70 0.80 0.50 0.85 0.85 0.85 0.85 1.00 0.85 0.85 0.20 0.20 0.20 0.20
1.20 0.95 0.70 0.85 0.85 0.85 0.85 0.25 0.30 0.30 1.30 1.30 1.30 1.30
1.05/0.95 1.11/1.00 1.14/1.00 1.15/1.00 1.15/0.99 1.15/0.99 1.15/0.99 1.12/0.98 1.10/0.96 1.08/0.91 1.07/0.89 1.06/0.86 1.06/0.84 1.05/0.81
175 561 562 562 549 548 549 538 490 388 255 120 78 47
Fig. 3. The maximal IF values for the different weekly bins. The values of a and g that yield the particular maximal IF(J) vary with J (the weekly bin). No data is presented for J > 10 as the best/worst IF outweigh one another and the improvement is insignificant.
Fig. 4. The improvement factor of the hybrid model as a function of weekly bin for ða; gÞ ¼ ð0:85; 0:85Þ: IF(J; 0.85, 0.85). Other (a; g) pairs exist that yield similar IF(J; a; g). For J > 10 the improvement is insignificant.
when combined with the velocity field, it has no physical meaning in itself and should only be used to properly incorporate the geopotential and temperature fields into the velocity field. To demonstrate this point we need to look at the calculated and observed trajectories of a particular balloon. In Fig. 5 we show the various trajectories relevant to the first 2 weeks of EOLE balloon #183. The IF of the best fitting hybrid model’s trajectory is 9.30 and is obtained for a ¼ 0:35 and g ¼ 1:3: Both the observed trajectory (solid curve, upper panel) and the best fitting one generated by the hybrid model (dotted curve, upper panel) span the observed
B1001 of longitude and B151 of latitude. By comparison, the trajectory that results from pure advection (solid curve, lower panel) spans about 501 of latitude and more than 3001 of longitude! It is clear, therefore, that pure advection by the NCEP velocity field provides a very poor descriptor of the observed trajectory. The trajectory due to the ‘‘correctional velocity’’ alone, too, bears no resemblance, whatsoever, to the observed trajectory (dotted curve, lower panel). However, a linear combination of the two corresponding (notice that the velocities shown are not the ones actually used in Eq. (5) since the balloon’s coordinates determine the actual
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
489
Fig. 5. Top panelFthe observed trajectory of EOLE balloon #183 (solid curve) and its best fitting hybrid model calculated trajectory (dotted curve) after 2 weeks. The improvement factor of the best fitting hybrid model trajectory is 930%! Lower panelFtrajectories due to pure advection by NCEP velocity field (solid curve) and to ‘‘correctional velocity’’ alone (dotted curve). Each of the two trajectories in the lower panel, alone provides a very poor estimate of the real trajectory but a proper combination of the two corresponding velocities yields an accurate one. Note the ‘‘inertial’’ loops in the dynamical model (a ¼ 1:0) trajectory.
velocities used there) velocities with a ¼ 0:35 offers a much better reproduction of the observed trajectory (i.e. IF=9.30). The better fit with the observed trajectory is achieved by properly employing NCEP’s velocity, geopotential and temperature fields instead of the velocity field alone. From the outset, the undertaking of this work of forecasting balloon trajectories for periods of over 2 months is an ambitious one. To explain why the hybrid model performs differently on different time scales (i.e. IF varies with J) we note that both Morel and Bandeen (1973) and in Morel and Larcheveque (1974) report on a dramatic change with time of the dispersal rates of EOLE balloons. During the first week the balloons disperse very fastFthe e-folding dispersal time is about 1.3 dayFafter which time the balloons reach a stage of nearly uniform distribution. The saturated stage occurs when the mean pair separation between balloons is about 2000 km and the balloons are distributed all over the extra-tropical latitudes of the Southern Hemisphere. Given the fast initial dispersal, it can be expected that accurate trajectory reconstruction will be difficult to achieve in the first week. By comparison, in the saturated stage the statistical properties of the individual balloon trajectories
remain constant with time so that a satisfactory reconstruction of these trajectories is possible. The dependence of our results, as measured by the improvement factor, on the length of the trajectories is in accordance with the observations of the EOLE balloons described above. At all weekly bins between 2 and 10 the hybrid model provides improved (statistically) trajectories compared with the pure advection ones. This is manifested both in the number of (a; g) pairs where IF>1.0 and in its maximal value being larger than 1.0 by up to 15% while its minimal value being only a few percent below 1.0. In contrast, for 1 weeklong trajectories (as well as for trajectories shorter than 1 week or longer than 10 weeksFresults not shown here) the hybrid model does not provide any improvement (in a statistically significant way) over the pure advection by the velocity field. The minimal value of IF deviates from 1.0 more than its maximal value and the number of (a; g) pairs where IF>1.0 is smaller than that with IFo1.0. The improvement of the present hybrid model compared with the pure advection model is pronounced at a values that are sufficiently large. When a is too small the hybrid model does not provide any significant improvementFsince, by definition, IF(a ¼ 0)=1.0. This
490
N. Paldor et al. / Atmospheric Environment 36 (2002) 483–490
expected result is very conspicuous in the numerical computations shown in Fig. 2. In comparison, in the midtropospheric, ‘‘blended’’ model of Stohl and Seibert (1998) an equation, similar to Eq. (3) (suggested by J.ager; 1992), was employed but with time-dependent weightFa: The value of the weight (i.e. a) in their model was determined at each time step by the ratio between the integration time step, which is determined dynamically by the CFL condition, and the time interval of the velocity fields. Since accurate integration requires that the former time steps be small, their results should be compared to the small a limit of the present study. As it turns out, the two models agree with each otherFStohl and Seibert (1998) found that their blended model performs worse than any other 2D kinematic (e.g. pure advection) model they tried while Fig. 2 indicates that no improvement occurs in the hybrid model when a is too small. An additional point of comparison between the present model and the blended model of Stohl and Seibert (1998) is the small J limit of Fig. 3 (see also Table 1). Here it is clear that for J ¼ 1 the present hybrid model offers no improvement over the pure advection model. The Stohl and Seibert (1998) model was applied to 6-day long trajectories and, as in our case, the resulting trajectories were not better than any 2D kinematic model trajectories. In the review of Stohl (1998) an inter-model comparison and error analysis was done but a study of the model performance against observations was not attempted. By comparison, in the present study the hybrid model’s improvement is verified by comparing the calculated trajectories with the observations. This puts the conclusions reached in this study on much firmer grounds. On time scale between 2 and 9 weeks our results indicate that the improvement factor of the hybrid model exceeds 1.1. The reason for the drop in the IF values at longer times is not clear. However, the maximum value of the IF periods around 4–7 weeks makes the hybrid model particularly efficient for the prediction of trajectories on the time scales of 1–2 months, which is relevant to the STRATEOLE experiment.
Acknowledgements Financial support for this study was provided by the US-Israel Binational Science via a research grant to the Hebrew University of Jerusalem and by the European Sciences Foundation, which provided YD with travel funds allowing him to spend a summer at LMD where the EOLE trajectories and NCEP fields were processed. Discussions held with F. Vial were very helpful in clarifying some of the subtle issues involved with EOLE trajectories and M. Pinsky, of the Hebrew University, provided helpful programming advice.
References Dvorkin, Y., Paldor, N., Basdevant, C., 2001. Modeling balloon trajectories in the tropical stratosphere with a hybrid model using analyzed fields. Quarterly Journal of Royal Meteorological Society 127 (573A), 975–988. J.ager, A., 1992. Isentrope trajektorien und ihre Anwendung auf das Konzept der potentiellen Vorticity bei orograpisch induzierten Lee-Zyklogenesen. Diploma Thesis, University of Innsbruck, Austria. Kalnay, E.M., Kanamitsu, R., Kistler, W., Collins, D., Deaven, L., Gandin, M., Iredell, S., Saha, G., White, J., Woollen, Y., Zhu, M., Chelliah, W., Ebisuzaki, W., Higgins, J., Janowiak, K.C., Mo, C., Ropelewski, J., Wang, A., Leetmaa, R., Reynolds, R., Jenne and Joseph, D., 1996. The NCEP/NCAR 40-year reanalysis project. Bulletin of American Meteorological Society 77, 437–471. Morel, P., Bandeen, W., 1973. The EOLE experiment: early results and current objectives. Bulletin of American Meteorological Society 54, 298–306. Morel, P., Larcheveque, M., 1974. Relative dispersion of constant-level balloons in the 200 mb general circulation. Journal of Atmospheric Sciences 31, 2189–2196. Stohl, A., 1998. Computation, accuracy and applications of trajectoriesFa review and bibliography. Atmospheric Environment 32, 947–966. Stohl, A., Seibert, P., 1998. Accuracy of trajectories as determined from the conservation of meteorological tracers. Quarterly Journal of Royal Meteorological Society 124, 1465–1484.