Sensitivity of ensemble Lagrangian reconstructions to assimilated wind time step resolution

Abstract. We study the impact of temporal and spatial resolution and changes in modelled meteorological winds in the context of diffusive ensemble Lagrangian reconstructions. In situ tracer measurements are modelled based on coarse resolution global 3-D tracer distributions from a chemistry-transport model and on different time series of meteorological wind fields including a special set of 1-hourly analysed winds which is compared with 3 and 6-hourly operational analysed winds and with 3-hourly ERA-interim reanalysis. Increasing the time resolution of the advecting winds from three to one hour using the operational winds provides an improvement on diffusive reconstructions in the period studied but smaller than that obtained from six to three hours. The positive impact of using 1-hourly winds is similar to that obtained using ERA-Interim 3-hourly winds instead of the 3-hourly ECMWF operational analysis for the same period. This study sets out a technique to quantify differences in time series of meteorological wind fields here applied to assess the optimal space and time resolutions for ensemble Lagrangian reconstructions in the lower stratosphere.


Introduction
In recent years the description of transport by atmospheric flows has seen significant progress. The calculation of Lagrangian backward trajectories is a widely used tool for the study of transport, both in the troposphere and the stratosphere (Wilson et al., 1994;Stohl et al., 2005). In the latter, Correspondence to: I. Pisso (ignacio@jamstec.go.jp) since the seminal works of Sutton (1994) and Norton (1994) the tools and the advecting fields have consistently improved.
The small-scale structure of chemical tracer distributions can be reconstructed based on coarse 3-D tracer fields and Lagrangian trajectories under the assumption that stirring is performed by the resolved scales of motion. In practice, this approach is performing very well where and when motion is basically layerwise. This is true in the stratosphere but also in a large portion of the troposphere as well (Newell et al., 1999). However, individual parcel trajectories are subject to chaotic sensitivity to the meteorological winds and purely advective reconstructions generate a large number of spurious structures as the reconstruction time increases. This effect is circumvented by adding (vertical) stochastic diffusion to ensembles of trajectories (Legras et al., 2005;Pisso and Legras, 2008) based on Green's function properties which significantly improves the quality of tracer reconstructions and eliminates the dependency on the reconstruction time. Nevertheless, the choice of temporal and spatial resolution of wind fields used to calculate trajectories remains a major problem. It has been shown that the standard 6-hourly sampled winds may produce large overestimates of the tracer concentration fluctuations (Stohl et al., 2004;Legras et al., 2005).
The combined use of analyses and forecasts yields 3hourly sampled analysed wind fields (Stohl et al., 2005;Bregman et al., 2006) which provide more accurate results for trajectory calculations and reconstructions of chemical compounds. The question arises whether further improvement can be obtained at higher temporal resolutions. An additional question is how modifications introduced in the model used for assimilation and weather forecast may modify the reconstructions.
In this work we study the impact of changes in the temporal resolution by using a special data set of analysed wind fields at 1-hourly resolution and we compare the results obtained with two versions of the ECMWF model and assimilation system. We use as a reference for comparisons a O 3 balloon profile already studied by Pisso and Legras (2008) which exhibits an interesting case of extratropical intrusion in the subtropics that can be used as a benchmark for various tests. In order to obtain more robust statistics, we also study the sensitivity of reconstructing a series of ozone profiles from the SHADOZ database (Thompson, 2003) for the period over which the special wind dataset is available. To our knowledge, this work had never been undertaken before since no operational center provides forecast fields with a temporal resolution shorter than 3 h.

Tracer data
The ozone profile ( Fig. 1) used for the case study in the first part of this work was collected by the Solid State Ozone Sensor instrument (Hansford et al., 2005) during the HIBIS-CUS campaign (Pommereau et al., 2007) over Bauru,Brazil (49.0 W,22.4 S) on 13 February 2004. The sampling frequency is about 1 Hz, yielding more than 10 000 measurement points. The ozone peak observed around 17 km altitude in this profile has been associated with a synoptic extratropical intrusion and is accurately reproduced by ensembles of diffusive trajectories (Pisso and Legras, 2008).
In order to provide a statistical analysis for different locations and times during the analysed period (1 January 2000 to 1 March 2004), available profiles from the SHADOZ database (Thompson (2003), http://croc.gsfc.nasa. gov/shadoz/) have been used to create a set of reconstructions. The data are originating from five stations (Ascension, Fiji, Irene, Reunion and Samoa) in the tropics and subtropics. The dataset is made of 23 profiles containing a total of 17 743 measurement points.

Back trajectories
Diffusive back-trajectories are calculated with TRACZILLA (Legras et al., 2003(Legras et al., , 2005Pisso and Legras, 2008;Pisso et al., 2009), a modified version of FLEXPART (Stohl et al., 2005) which uses a direct vertical interpolation of ECMWF winds, from hybrid levels, with a time step of 15 min. In this work, different diffusive reconstructions of the HIBISCUS O 3 profile are performed using different time/spatial resolutions.
A special wind dataset has been produced for this study, using the version Cy28r1 of ECMWF model operational from 1 January 2004 00:00 UT to 1 March 2004 00:00 UT (corresponding to the time of the HIBISCUS campaign), and has been archived on the ECMWF MARS database. The originality of this dataset is that the first 12-h forecast fields are archived every hour instead of the standard 3-hourly time interval. The 12-h short term forecasts start from the 4d-Var analysis at 00:00 UTC and 12:00 UTC and the version of the ECMWF model used has a T511 spectral truncation corresponding to a horizontal resolution of the order of 40 km with 60 vertical levels.
The reconstructions based on 6-hourly winds use the ECMWF standard operational analyses at 00:00 UT, 06:00 UT, 12:00 UT and 18:00 UT interpolated on a 1 • ×1 • grid. Those based on 3-hourly winds use the same analyses and at the intermediate times the 3-h and 9-h forecasts started from the 00:00 UT and 12:00 UT analyses interpolated on a 1 • ×1 • grid. Unlike previous studies, we have performed a set of reconstructions based on 1-hourly winds which extend the 3-hourly experiment by using the specially archived 1hourly forecasts described above. In order to test the consistency with the finer time resolution, we not only use analyses interpolated on a 1 • ×1 • grid, but also refine the spatial resolution interpolating the wind fields on a horizontal 0.5 • ×0.5 • grid. In all the cases the vertical resolution is the original one (60 levels). In addition, ERA-Interim fields have been used in order to estimate the effect of changes in the model on the diffusive reconstructions. In this case the time resolution is set to 3 h using at the intermediate times the 3-h and 9-h forecasts started from the 00:00 UT and 12:00 UT reanalyses interpolated on a 1 • ×1 • grid. ERA-Interim is the latest ECMWF reanalysis (Simmons et al., 2007) which is based on the version Cy31r2 of the model with a spectral truncation of T255L.
Three values of diffusivity D (0.1 m 2 s −1 , 0.5 m 2 s −1 , 1 m 2 s −1 ) have been used to reconstruct the measured profiles for the temporal and spatial resolution of each wind field. The experimental design for the high resolution profile is based on the study by Pisso and Legras (2008) who investigated the effect of changes in diffusivity using the standard time resolution of 3 hours. In the reconstructions presented here, unlike in previous studies, the seed of the pseudo random perturbations are designed to be the same for every timestep and wind configuration to ensure the changes are due to differences in the advective field only. The same analysis has been reproduced also for the selected ozone profiles in the SHADOZ database in order to provide statistical information. In the latter case all the fields are interpolated on a 1 • ×1 • grid. The back-trajectories were calculated for 300 h for the HIBISCUS profile (reconstructions shown at 192 h in Fig. 1) and for 504 h for the profiles in the SHADOZ database (reconstructions ranging from 48 h to 504 h for the statistical analysis). No sources or sinks of ozone are considered along the trajectories in this study.

3-D chemical fields
Ozone mixing ratios are interpolated from REPROBUS three-dimensional is a 3-D chemical-transport model with a comprehensive treatment of gas phase and heterogeneous chemical processes in the stratosphere (Lefèvre et al., 1998). Longlived species, including ozone, undergo advection by a semi-Lagrangian scheme forced by 3-hourly ECMWF meteorological data constructed from analyses and forecasts.
The model is integrated on 42 hybrid pressure levels that extend from the ground up to 0.1 hPa, with a horizontal resolution of 2 • latitude by 2 • longitude. The simulation was initiated on 1 April 2002 using the 3-D ECMWF analysis for ozone and an April zonal mean computed from a previous long-term simulation of REPROBUS for other species.

HIBISCUS ozone profile
As expected, the diffusive reconstructions displayed in Fig. 1 present fine scale structures which are absent from the REPROBUS dataset interpolated at the time and location of the measurements. This illustrates how spatial information about the fine scale structure can be extracted by the recon-struction method from the initial chemical field and a time series of wind fields, both at coarser resolution. This is seen more clearly with small values of diffusivity (D=0.1 m 2 s −1 ) because high values tend to smooth the small-scale variability. The comparison between first and second rows in Fig. 1 shows that inserting 3-hourly forecasts between the 6-hourly analyses improves substantially the reconstructions. This is consistent with the results obtained by Legras et al. (2005) for the polar stratosphere and with the findings of Bregman et al. (2006) and Berthet et al. (2006) who demonstrated better quality chemical fields in Eulerian and semi-Lagrangian CTMs.
For D=0.1 m 2 s −1 (panel a) the 6-hourly reconstruction exhibits spurious small fluctuations and the ozone peak is shifted downwards by almost 1000 m.
For D=0.5 m 2 s −1 (panel b) and D=1 m 2 s −1 (panel c) the fluctuations disappear because of the increase of diffusivity but the peak remains shifted with respect to the measurements.
As shown in Pisso and Legras (2008), for the 3hourly reconstructions (second row in Fig. 1  The comparison between the second and the third rows in Fig. 1 shows small differences between reconstructions carried out with a 3-hourly temporal resolution at 1 • ×1 • spatial resolution compared to those with a 1-hourly temporal resolution at 0.5 • ×0.5 • spatial resolution. Finer spatial resolution has been performed aiming to balance the higher temporal resolution. The same reconstructions performed with 1-hourly winds interpolated on a 1 • ×1 • grid yield similar results (not shown).
To explain the small differences obtained between the 1hourly and the 3-hourly reconstructions, we consider the locations of the trajectory points at t=−192 h (the backward time used for the reconstructions) before the measurements to analyse the influence of the resolution on the advective transport. Figure 2 displays, at t=−192 h, the location of a cloud of particles produced from a single point on the profile and advected backwards using different values of D and of the sampling time of the advecting winds. This ensemble of particles is chosen because it corresponds to the most interesting feature of the measured profile, that is the ozone peak at 17 km altitude. Figure 2, first column (D=0.1 m 2 s −1 ), clearly shows that the clouds of particles for the 1-hourly and 3-hourly reconstructions are close and largely overlap while they are disjoint from the cloud of particles provided by the 6-hourly reconstruction. For D=0.5 m 2 s −1 , the cloud of particles provided by the 6hourly reconstructions is close to the eastern part of the cloud of points provided by the 1-hourly and 3-hourly reconstructions. For D=1 m 2 s −1 the differences tend to be smaller, because the competing diffusion hinders the effect of the global perturbation in the wind field. See also additional material for an animated representation of the effect of backward advection, (http://www.atmos-chem-phys.net/10/3155/ 2010/acp-10-3155-2010-supplement.zip).
Refinement of the time resolution can be interpreted as a global perturbation on the flow associated with the time series of wind fields. The schematic in Fig. 3 illustrates the relationship between sensitivity to initial conditions and turbulent diffusion. Consider a point lying on either side of a separatrix associated with the stable manifold of an hyperbolic point which is displaced by the global perturbation of the advecting field. A neighborhood of this point (associated Atmos. Chem. Phys., 10, 3155-3162, 2010 www.atmos-chem-phys.net/10/3155/2010/ Fig. 3. Schematic on structural stability: qualitative different behaviors occur when the ratio between dissipative (initial) blob radius (associated with D) and the scale of wind field perturbation is smaller (blue, lower case 'c') or larger (red, capital 'C') than 1. The perturbation is here the displacement of the separatrix associated with the stable manifold of an hyperbolic point. The solid lines are the unperturbed case while the dashed lines correspond to the perturbed field, that we assume resulting from the increase in time resolution. If the relative size of the blob with respect to the perturbation is small, the separatrix passes across the blob and the evolution changes significantly, whereas if the relative size is large, the change is only moderate.
19 Fig. 3. Schematic on structural stability: qualitative different behaviors occur when the ratio between dissipative (initial) blob radius (associated with D) and the scale of wind field perturbation is smaller (blue, lower case "c") or larger (red, capital "C") than 1. The perturbation is here the displacement of the separatrix associated with the stable manifold of an hyperbolic point. The solid lines are the unperturbed case while the dashed lines correspond to the perturbed field, that we assume resulting from the increase in time resolution. If the relative size of the blob with respect to the perturbation is small, the separatrix passes across the blob and the evolution changes significantly, whereas if the relative size is large, the change is only moderate.
with an in situ measurement) corresponds to a cloud of particles, whose dissipative radius (the radius of the cloud after the initial stage dominated by diffusion, Legras et al., 2005) depends on the stochastic noise related to the diffusivity coefficient D. The distribution resulting from backward integration of an advected cloud (c 0 ) may change in the unperturbed (c 1 ) and perturbed case (c 2 ) (on either side of the separatrix) if the initial radius is smaller than the displacement of the separatrix (lower case indicates lower diffusion). In that case, the cloud of particles (c 0 ) visits different locations and samples different air masses. On the other hand, if the initial radius of the cloud (C 0 , upper case indicates higher diffusion) is large with respect to the displacement of the separatrix, the domains occupied by the clouds advected by the original and the perturbed flows overlap and, as backward integration evolves, the difference tends to be small. The schematic in Fig. 3 points out why, for large enough values of D, the change in wind resolution has a reduced effect, whereas for small diffusivity the change in resolution may be more significant. The effect on the reconstructions is large if the deviation of the flow due to increased resolution exceeds the diffusive radius or the resolution of the chemical field used for interpolation at the trajectory location at the end of reconstruction time.
To make an objective comparison between the backward positions associated with different reconstructions, we have calculated the distance between clouds in the sense of point sets weighted by the particle density. The distance between two clouds (here C indicates an arbitrary cloud) is then defined as dist (C 1 ,C 2 ;t)= |f (C 1 ;x,t)−f (C 2 ;x,t)|dx where f (C;x,t) is the probability distribution in space of the cloud C at time t. A histogram of the trajectory positions consistent with the global model grid was used to estimate the probability density f (.;.). The function dist (.,.;.) is a distance over the classes of clouds defined by their histogram and has a value between 0 (the clouds occupy exactly the same grid boxes) and 2 (the clouds are totally disjoint). Figure 4 shows the evolution in time of these distances (C 6h − C 3h , C 6h − C 1h , C 3h − C 1h ) for different values of D. For all D values, the distance between the 1-hourly and 3-hourly reconstructions is smaller than between 1-hourly and 6-hourly and between 3-hourly and 6hourly reconstructions. This indicates that the 1-hourly and 3-hourly wind fields provide similar trajectories and therefore similar dynamical fields which are both different from the 6-hourly wind fields. Both Figs. 2 and 4 show large differences between the location of the cloud of particles for the same time/spatial resolution but different diffusivity values. Using 3-hourly and 1-hourly wind fields, the trajectories for D=0.1 m 2 s −1 mainly come from the subtropics 192 h before the measurements while they originate from a large range of latitudes from the subtropics to the polar regions for D=0.5 m 2 s −1 . This means that the diffusivity is a key parameter in the reconstruction method significantly affecting the location of the trajectories. This is consistent with the results in Fig. 1 showing larger differences in the reconstructed profiles when varying D values than when using 1-hourly instead of 3-hourly time resolution.
Consequently, our results suggest that when D is 0.5 m 2 s −1 (i.e., the optimal diffusivity parameter for this case as shown in Pisso and Legras, 2008), increasing the time resolution from 3 to 1 h does not significantly modify the large scale (synoptic) structures in tracer reconstructions. Such behaviour is qualitatively different from the single trajectory perspective, in which sensitivity to initial conditions has to be taken always into account.

SHADOZ profiles
The study of SHADOZ ozone profiles in this section provides statistical information which is complementary to the detailed study of the HIBISCUS ozone profile. Contrarily to the analysis of the HIBISCUS profile, no particular features are analysed on SHADOZ data but statistical diagnostics are calculated.
An objective metric of the quality of a reconstruction can be calculated as the average integrated difference (Pisso et al., 2009) from a set of reconstructions performed using operational analysis winds at 1, 3 and 6 h resolution and from ERA interim analysis winds with 3 hour resolution. A 1 • ×1 • spatial resolution was used for all these analysed winds.   Fig. 4. Distance between trajectory clouds (average for ensembles initialised at 16, 17 and 18 km along the vertical HIBISCUS profile) as a function of backwards advection time for different advecting winds and turbulent diffusivities. Blue: distance between 6 hourly and 3 hourly winds (C 6h − C 3h ), Green: distance between 6 hourly and 1 hourly winds (C 6h − C 1h ), Red: distance between 3 hourly and 1 hourly winds (C 3h − C 1h ). Dotted lines: D=0.1 m 2 s −1 . Dashed lines: D=0.5 m 2 s −1 . Continuous lines: D=1 m 2 s −1 . The abscissa represents time of backward advection in hours. The ordinate correspond to the function distance (dist (C 1 ,C 2 ;t)= |f (C 1 ;x,t)−f (C 2 ;x,t)|dx, where f (C;x,t) is the probability distribution in space of the cloud C at time t) between pairs of particle clouds, dimensionless. Table 1. Mean values and standard deviation along the studied period of the integrated difference between reconstructions and observations for each D and wind set: EI, 1 h, 3 h, and 6 h correspond to ERA Interim, one hourly, three hourly and six hourly winds respectively. Units in ppmv of Ozone.
where N d is the number of measurement points in the SHADOZ database for the day d, T d is the number of matches with REPROBUS output fields for back trajectories initialised at day d, C meas i is the ozone mixing ratio at the ith measurement point for day d and C rec it are the corresponding mixing ratios of the diffusive reconstructions, depending on the diffusivity D and the advecting windfield W . is a function of d, the date, D, the diffusivity, and W , the wind fields used for advection. Figure 5 shows the values of this metric over the period for which 1-hourly winds are available. Table 1 summarises the mean values of the scoring metric along the analysed period for different wind fields and diffusivity values. Standard deviations are included as well but it has to be taken into account that this estimator may be biased as the distribution of is not expected to be Gaussian. The profiles reach a height of 30 km and the maximum values of ozone are of the order of 10 ppmv. N d , the number of measurement points per day, reaches the order of 10 3 and T d , the time of back integration, reaches 504 h with up to 4 REPROBUS matches.
In most of the cases, the 1-hourly winds perform better than the 3-hourly winds but the improvement is lower than that of 3-hourly winds with respect to the 6-hourly winds. There are some points toward the end of the analysed period in which the 3-hourly winds perform better than the 1-hourly winds. It can be noticed that the performance of the ERA-Interim 3-hourly winds is comparable to that of the 1-hourly winds from the ECMWF operational model in 2004 and that the standard deviation obtained with ERA-Interim 3-hourly winds is the smallest of the set (Table 1), illustrating how changes in the model yield significant improvements.
Although this is not in the main focus of this work, the analysis of the Hibiscus profile is consistent with the results in Pisso and Legras (2008) that a diffusivity of 0.1 m 2 s −1 is too low and 1 m 2 s −1 is too high to provide accurate reconstructions of ozone profiles. In the case of D lower than 0.1 m 2 s −1 , the signal is hindered by spurious fluctuations caused by the undersampling of gravity waves. On the other hand, too large D flattens the laminated structures identified in the profile and the integrated difference in mixing ratio along the profile becomes larger. In the case of the statistical analysis along whole SHADOZ profiles, the distance between reconstructions with 1 m 2 s −1 and the dataset is smaller, but is not excluded that too much artificial diffusion is needed to overcome spurious fluctuations due possibly to undersampling of gravity waves (cf. Legras et al., 2005).
It is interesting to remark that for the purposes of this study, the addition of diffusion removes the sensitivity to initial perturbations smaller than the dissipative radius, suggesting a way to quantify the information gained during the process of assimilation.

Conclusions
We have studied the influence of temporal and spatial resolution of advective transport on Lagrangian tracer profile reconstructions. We have focused on the subtropics in January and February 2004 (period for which a special dataset of 1hourly global wind fields has been prepared), comparing high resolution O 3 in-situ measurements with diffusive ensemble reconstructions calculated with 6, 3 and 1-hourly wind fields provided by the ECMWF meteorological analysis and ERA Interim reanalysis.
Increased time resolution improves the quality of the reconstructions in a statistical sense, likely to be related to the reduction of small and mesoscale inaccuracies , for example related to the undersampling of gravity waves. The improvement in using 1-hourly instead of 3-hourly is statistically significant but it is smaller with respect to using 3-hourly instead of 6-hourly winds.
In some cases of presence of large-scale structures governed by synoptic scale motions such as in the HIBISCUS profile the use of 1-hourly winds instead of 3-hourly winds does not significantly change the quality of the reconstructions as illustrated by the geometry of the advected cloud of particles used to approximate the Green's functions needed for the reconstruction.
Although the results obtained in this study are limited to the times and locations of the measurements used to calculate the reconstructions, they support the use of a 3-hourly time resolution and a 1 • ×1 • spatial resolution. One could argue that the reconstructions done in this study are limited by the REPROBUS output fields having a coarse 2 • ×2 • resolution. However, previous studies (Legras et al., 2005) suggest that the coarse resolution of initial tracer fields does not play a major role provided reconstructions are integrated over a long enough time interval. The results obtained from the Hibiscus profile show that using 1-hourly wind fields in the reconstructions does not add significant information on the large scale atmospheric dynamics. This suggests that in the ECMWF analyses and 12 h forecasts for 2004, the information contained at the scale of 1 h is not significantly different from the information contained at 3 h. This does not necessarily apply to more recent versions of the analysis performed at higher resolution, but can be taken as valid also for analyses or reanalyses performed at a lower resolution.
It is also suggested that the model improvement represented by the ERA-Interim reanalysis, in spite of being performed at a horizontal resolution lower than the 2004 operational analysis, makes the quality of the corresponding reconstructions with 3-hourly winds comparable to the one obtained with the 1-hourly winds of this operational analysis. This is a hint that model improvement can be as rewarding as resolution improvement.
This study sets out a technique to assess the optimal space and time resolutions for ensemble Lagrangian reconstructions in the lower stratosphere. The robustness of the reconstructions is associated with the structural stability of atmospheric flows under perturbations (Demazure, 2000) like those introduced by changes in resolution and parametrization tunings. Including diffusion in the Lagrangian calculations, a time resolution can be reached such that further refinement does not add significant information, in contrast to a single trajectory approach. It is interesting to note that the results regarding the optimal diffusion that represents unresolved turbulent scales are robust to the change of resolution and temporal sampling. As Lagrangian calculations can achieve arbitrary resolutions, the described scheme is suited to study such changes regardless of the resolution of the fields involved. The optimal temporal/spatial resolution for Lagrangian off-line calculations is in practice modeldependent. The kind of experiments described here are likely to be useful in the context of the assessment of assimilation schemes and the effect of changes in the assimilation parameters on the final analysed field. It would be interesting to apply the tools described here in case of ensemble assimilation in order to assess the behaviour of different members of the ensemble.