The usefulness of any simulation of atmospheric tracers using low-resolution winds relies on both the dominance of large spatial scales in the strain and time dependence that results in a cascade in tracer scales. Here, a quantitative study on the accuracy of such tracer studies is made using the contour advection technique. It is shown that, although contour stretching rates are very insensitive to the spatial truncation of the wind field, the displacement errors in filament position are sensitive. A knowledge of displacement characteristics is essential if Lagrangian simulations are to be used for the inference of airmass origin. A quantitative lower estimate is obtained for the tracer scale factor (TSF): the ratio of the smallest resolved scale in the advecting wind field to the smallest “trustworthy” scale in the tracer field. For a baroclinic wave life cycle the TSF = 6.1 ± 0.3 while for the Northern Hemisphere wintertime lower stratosphere the TSF = 5.5 ± 0.5, when using the most stringent definition of the trustworthy scale. The similarity in the TSF for the two flows is striking and an explanation is discussed in terms of the activity of potential vorticity (PV) filaments.
Uncertainty in contour initialization is investigated for the stratospheric case. The effect of smoothing initial contours is to introduce a spinup time, after which wind field truncation errors take over from initialization errors (2–3 days). It is also shown that false detail from the proliferation of finescale filaments limits the useful lifetime of such contour advection simulations to ∼3σ−1 days, where σ is the filament thinning rate, unless filaments narrower than the trustworthy scale are removed by contour surgery. In addition, PV analysis error and diabatic effects are so strong that only PV filaments wider than 50 km are at all believable, even for very high-resolution winds. The minimum wind field resolution required to accurately simulate filaments down to the erosion scale in the stratosphere (given an initial contour) is estimated and the implications for the modeling of atmospheric chemistry are briefly discussed.
Advection by atmospheric flows results in the stretching and folding of material contours and thence to an exponential rate of decrease in the scale of passive tracer features and a corresponding increase in tracer gradients that is frequently referred to as a scale cascade (Cocke 1969; Juckes and McIntyre 1987; Haynes and Anglade 1997). This occurs even though the wind field is typically dominated by large-scale features such as synoptic-scale weather systems. On this basis, analogies have been drawn with dynamical systems theories for structurally simple steady flows with time-periodic perturbations that result in complex trajectory behavior known as chaotic advection (e.g., Aref 1984; Ottino 1989). Exponential stretching of line elements occurs along the axis of dilatation in a steady deformation flow. In chaotic advection studies, stretching and folding of contours occurs in the vicinity of hyperbolic fixed points [e.g., the stagnation points in the comoving frame of a perturbed traveling wave, as studied by Pierrehumbert (1991)]. In flows with more general time dependence the hyperbolic points move and filamentation events can occur when a hyperbolic point crosses a tracer contour. For example, Polvani et al. (1989) estimated the onset of potential vorticity (PV) filamentation at the edge of linearly unstable vortices using the time when the hyperbolic points in corotating streamfunction and the vortex edge cross. Through the chaotic advection process, coarse-grained representations of atmospheric winds can generate finescale tracer filaments that may qualitatively resemble those simulated by high-resolution winds. However, displacement errors in the filament positioning are inevitable. The magnitude of such displacement errors should be known if a Lagrangian simulation is to be used to infer, with any certainty, the origins and history of the air within a filament whose edges have been detected by high gradients in measurements of chemical tracers.
Waugh and Plumb (1994, hereafter WP94) made a preliminary investigation of this problem using the contour advection (CA) technique. In CA calculations particles are positioned along chosen contours and then advected by 2D velocity fields by numerically integrating dx/dt = u(x, t). In order to maintain accuracy, particles are redistributed along the contour (defined by a local cubic spline), at regular time intervals, adding more particles in sections of high curvature. Usually, isentropic winds are used for advection since isentropic surfaces are material surfaces in the absence of diabatic heating or friction. Contours of constant tracer mass mixing ratio, as determined by the interpolation of a gridded tracer field onto the isentropic surface, are used to initialize the particle positions. The contours at a later time give a representation of the tracer isopleths, in the absence of nonconservative processes, which has extremely high spatial resolution, although the concentration is discretized. The technique was adapted from the contour surgery algorithm of Dritschel (1989a) and used by Norton (1994) and WP94 to study the polar stratospheric vortex. WP94 compared PV contours that were passively advected using wind fields that had been truncated spatially and temporally to the same contours from an accurate contour dynamics calculation in which the PV filaments were dynamically active. They found that coarse wind fields (on a 5° grid with records every 12 h) produced qualitatively realistic results.
Here a quantitative study is made of the relationship between the advecting wind field resolution and the smallest “trustworthy” scale in a tracer simulated using CA. It is hypothesized that smaller filaments are more likely to be displaced by a perpendicular distance that is some factor (εt) greater than their width; such filaments are defined here to be “bad.” Some time into a CA calculation the average displacement of the bad filaments is expected to reach a steady value that depends only upon the advecting wind field truncation and the nature of the flow. Thus the width of the widest bad filaments is expected to be roughly constant, although all filaments’ widths are decreasing with time; the proportion of bad filaments therefore increases with time. This information can subsequently be used to decide the minimum wind field truncation that should be used to accurately simulate tracer filaments of a chosen scale or larger.
The investigation is conducted for flow in two distinct atmospheric regimes: midlatitude weather systems (section 2) and the polar stratospheric vortex (section 3). The first case is simplified by using a baroclinic wave life cycle as a proxy for generic weather system behavior. In particular, the flow is adiabatic and frictionless and therefore layer-wise two-dimensional along isentropic surfaces. The Ertel PV field is used to initialize the contours. In this experiment the PV field is initially smooth and is therefore identical in the coarse-grained representations. In contrast, PV in the stratospheric case is strongly dependent upon truncation of the atmospheric analyses, introducing the added complexity of contour initialization.
Section 4 discusses the role of the dynamics in ensuring that the larger scales in the wind field have most influence on the positioning and formation of tracer filaments. In particular, the steep slope of the kinetic energy spectrum remains approximately constant even though the widths of PV filaments continually decrease with time. However, in the vicinity of unstable PV filaments that roll up nonlinearly, the strain is dominated by local rather than far-field PV anomalies. In this situation it is shown that the amplitude of undulations of unstable filaments is severely underestimated by CA with coarse-grained winds. There is a brief discussion on the halt of the tracer scale diminution by nonconservative processes in section 5. Starting from an estimate of the smallest scales existing in a stratospheric tracer field, the minimum-resolution wind field required for an accurate simulation of the tracer advection problem is estimated.
2. Advection by midlatitude weather systems
a. Contour advection simulations
The simulation of a tracer field evolving through advection by midlatitude weather systems is performed using the baroclinic wave life cycle experiments of Thorncroft et al. (1993, hereafter THM). These experiments isolate the essential dry (primitive equation) dynamics of weather system development on a baroclinically unstable jet. Here the LC1 experiment of THM has been repeated at much higher resolution (spectral truncation T341 with 30 levels in the vertical, concentrated near the tropopause). The evolution of surface temperature and interior potential vorticity closely resembles that observed in real cases of cyclogenesis. In particular, the wind fields, when wave development becomes nonlinear, exhibit most features characteristic of maturing weather systems (e.g., fronts). These experiments have the added advantage that isentropic surfaces are material surfaces and potential vorticity acts as a conserved tracer (aside from the effects of spatial truncation and numerical hyperdiffusion) so that PV contours on isentropic surfaces are material contours.
The greatest contour lengthening is expected to occur on isentropic surfaces that pass through the steering level of the baroclinic wave because the system relative meridional displacements will be greatest there. The 300-K potential temperature surface was chosen because it is the lowest isentropic surface that does not intersect the ground, thus avoiding any problems with the interpolation of velocity close to the ground. This surface passes through the steering level of the wave to the poleward side of the jet. Figure 1 shows the evolution of PV on the 300-K surface. By day 6 the baroclinic wave has begun to wrap up nonlinearly. There is cyclonic rollup to the poleward side of the jet and signs of anticyclonic turning in the subtropics. By day 8 a very long PV trough has been produced by the anticyclonic wave-breaking event. This trough is unstable and rolls up forming a cutoff cyclonic vortex over the next day. The rapid development here involves a strong interaction with the surface potential temperature gradients and is described as a case of secondary cyclone development by Thorncroft and Hoskins (1990). Unsteady straining also breaks up the cyclonic spiral and the PV field mixes through a combination of stirring and numerical dissipation. The strongest feature to remain is the subsynoptic-scale vortex that formed during the secondary cyclogenesis. This persists through to the end of the experiment at day 12. The different aspects of the PV evolution in this life cycle are characteristic of the extratropical storm-track regions, and the life cycle is therefore a useful tool for examining the tracer advection problem.
The effects of wind field truncation on filament displacement were then investigated in a set of CA simulations. Particles were initialized along the 2.0-PVU1 PV contour on the 300-K surface at day 4 of the life cycle. A control simulation of the PV-like contour was then obtained by advecting the contour with the T341 winds (interpolated onto the 300-K surface), linearly interpolating between records at 6-hourly intervals. Hereafter the temporal resolution of the advecting wind field will be denoted by Δt. The resulting contour is shown in Figs. 2a–d. The simulated contour is virtually identical to the 2.0-PVU contour of the PV field until day 8 (cf. Fig. 1). After this time the dissipation of the PV field structure results in substantial differences in detail although many features may be clearly identified (e.g., the vortex associated with secondary cyclogenesis). The complexity of the simulated contour increases rapidly. By day 12 the contour has become so convoluted that a single line can intersect the contour at over 200 points!
Accuracy experiments for 3D trajectory calculations indicated that the most severe trajectory errors occur through wind field truncation, especially in the vertical (Methven 1997). In order to keep the rms deviation of trajectories from their initial isentropic surface to less than 1 K in the LC1 experiment (at T85, L30 resolution) it was necessary to use cubic vertical interpolation of the winds. This cubic method was also used here to interpolate the 3D model winds onto the 300-K surface for the CA simulations. Note that, in the 3D trajectory calculations, temporal truncation errors incurred by increasing Δt from 3 to 6 h were found to be far smaller than typical spatial truncation errors. The errors introduced by using a fourth-order Runge–Kutta scheme, with a time step of 0.5 h, to integrate dx/dt = u(x, t) were found to be two orders of magnitude lower than from reducing horizontal truncation from T85 to T42 or the number of model levels from 60 to 30 (Methven 1997).
The accuracy parameter for the CA calculations, μ, was set to 0.02. This means that the initial spacing of particles around a contour lying on a great circle is μa (≈127 km). The cubic spline interpolation errors are then minimized by ensuring that particles are redistributed with a density e = μ2a/8 + μa(R/a)2/3, where a is the earth’s radius and R is a nonlocal curvature measure weighted toward closer sections of contour and places where velocity is changing most quickly (see Dritschel 1989a). Note that particle redistribution occurred every 6 h and contour surgery was not used so there was no limitation to the filament widths and the contour was not split.
Comparison PV contour simulations were obtained by temporally and spatially truncating the high-resolution life cycle winds to a greater degree. At day 4, when the CA simulations were initiated, the PV field was still wavelike and smooth so that it was virtually unchanged by spatial truncation. Nevertheless, the same initial contour (from T341 PV) was used for all simulations. Temporal truncation was achieved by taking wind records less frequently. Spatial truncation was achieved by chopping off the high wavenumbers in the control record to lower truncation limits and then transforming them onto a Gaussian grid at sufficient resolution to prevent the aliasing of quadratic terms, before interpolation from model levels onto the 300-K surface. Triangular truncation of the spherical harmonics was used with highest total wavenumbers N equal to 213, 170, 127, 85, 73, 63, 57, 52, 47, 42, 37, 31, 26, and 21. Quantitative estimates of the accuracy of the contour positioning by course-grained winds can be drawn by taking intersections through the contour and using these to define filament positions and widths. This procedure will be described in section 2b and the results presented in sections 2c and 2d.
Note that the smallest resolved length scale Lr for structures in a spectrally represented field cannot be defined uniquely because it is a function of the shape of the structure as well as the truncation limit. Lander and Hoskins (1996) studied the appearance of small-scale features after transformation to truncated spectral space and back to gridpoint space, in the manner of a spectral transform model using a “quadratic grid.” For instance, a δ function becomes a point spread function upon “truncation.” This function is “Bessel-like” with a central peak of half-width 1.6πa/N, where a is the earth’s radius. They also found that an isolated, narrow filament with a δ-function cross section smears to a width of 2.8πa/N upon truncation. In contrast, a filament that can be described as one half-wavelength of an encompassing wavelike structure with wavenumber N has a width of πa/N. Thus we can say that the smallest resolved feature in the wind field has a width sπa/N where the shape factor s lies between 1 and 3.2. For convenience in the following sections we will take Lr = πa/N although it is important to remember that a larger value may be appropriate.
b. Definition of filament widths and displacements
Advection by the baroclinic wave results in filamentary structures in PV contours on isentropic surfaces. The aim here is to find the widths and midpoints of filaments from the control and comparison advection calculations in order to examine the accuracy of filament simulation. The coordinates of the particles describing a PV contour were transformed so that the center of the region of interest (e.g., a spiral) lay on the pole in the new coordinate system. These coordinates were then projected onto a stereographic projection centered on the new pole (as in Fig. 2b). Lines were chosen to cut the majority of the filaments almost at right angles so that the great circle distances between consecutive points of intersection were approximately the filament widths, and the interval midpoints were defined as the filament positions. The points of intersection were determined by piecewise linear interpolation between the contour nodes. Since there were typically many filaments intersecting each cut, care was taken in developing an algorithm that successfully matched all filaments, in the control and comparison, that shared the same origins.
The high-PV intervals from comparison runs were matched to the nearest high-PV intervals in the control run.
This process was then repeated for the low-PV intervals.
Checks were made that the matched filaments alternated between high and low PV. Where this was not the case an attempt was made to restore the ordering using the second closest matches. If this failed, then ordering was restored by canceling the violating matches.
Any intervals in the comparison runs that were not matched in the control were described as “spurious” filaments.
Any intervals in the control run that had no match in the comparison runs were described as “lost” filaments.
Figure 3 shows the results of matching filaments intersecting the third line from the left in Fig. 2c. This example illustrates the remarkable correspondence of filaments achieved with very coarse-grained winds and the occurrence of lost and spurious filaments. The displacement of a matched filament was defined as the difference in its midpoint position between the comparison and control simulations. An error estimate, ε, for a matched filament’s position was then defined as the filament displacement divided by its width. A variable threshold error, εt, was then defined and filaments for which ε > εt were classified as bad. All matched filaments were associated with their width in the comparison simulation because it was the error associated with filaments of particular scales in the comparison simulation that was under investigation.
For each comparison run a list of spurious filaments was made. Typically spurious filaments occur in pairs where a very thin erroneous filament divides a real filament into two. Obviously it is the thin filament that is in error, not the space neighboring it. These spurious filaments typically arise when coarse graining the winds decreases the separation between hyperbolic points and the comparison contour, enabling an extra filamentation event to occur. The bad filaments and the smallest filament from each spurious pair were described collectively as “unbelievable” filaments.
c. Tracer scale cascade and displacement errors
During the baroclinic wave development there are time-dependent straining regions that are expected to result in the stretching and folding of material contours. Although segments of the contour may wrap up in eddies where contours only lengthen at an algebraic rate, the lengthening rate of the whole contour will be exponential. Since the isentropic flow is approximately nondivergent, exponential filament stretching suggests filament thinning at the same rate. However, it is found that the rate of decrease of the median of the filament width distribution [σ = −d(lnδmed)/dt] is faster by a factor that depends upon the details of filamentation events (discussed below). Figures 4a and 4b show the exponential lengthening and thinning rates as a function of wind field truncation. The rates are obtained by a least squares fit over days 8.5–12.0 and the standard deviations are used for the error bars. Here the smallest resolved scale in the wind field is taken to be πa/N as discussed at the end of section 2a.
It is clear that contour lengthening is indeed exponential to a high degree of accuracy. Interestingly, the rate λ hardly varies with the spatial truncation of the wind field over a huge range (T341 down to T26). However, the rate is decreased substantially by sampling the analyses at a coarser temporal resolution (interpolating between records at 12-h intervals rather than 6-h intervals). The same behavior is seen with σ, although it is about 1.5 times higher than λ. It is also expected that the average2 number of filaments intersected by a cut through a contour, n, will increase as the median filament width δmed decreases. Figure 4c shows that n is related to δmed by an approximate power law:
where n0 and δ0 represent n and δmed at t = 0. Both the control and comparison simulations form a compact line of points on a log–log plot, with slope m = 0.59 ± 0.01. The compactness of the line arises partly because the contours in all experiments start from the same initial condition. Since m < 1, n increases at a slower rate than δmed decreases. This again is the result of filamentation events.
If only one filament formed and experienced a similar stretching history along its length, then σ ≈ λ and m ≈ 1. However, each new filamentation event introduces two extra intervals along a line of intersection. If both of these intervals are much narrower than the parent filament (this is the distinction between folding an old filament and a filamentation event), then the median filament width will be decreased. However, the length of the new filament increases at the same rate as the old; therefore, σ/λ is increased by these filamentation events. Similarly, the increase in n is not affected by the distinction between folding and filamentation and we expect m < 1. The ratio σ/λ and m depend on the frequency of filamentation events and their details of early formation. Polvani and Plumb (1992) showed that the width of an ejected filament increases with the distance that a hyperbolic point crosses over into the parent filament. Note that the value of m was obtained using all CA simulations and all time frames and thus filamentation events typically result in a similar daughter-to-parent width ratio, introducing a scale invariant structure to the filament distribution. Indeed, the value of Kolmogorov capacity, a scale invariant measure like fractal dimension, of the set of points of intersection along each cut through the contour is roughly constant in time. The Kolmogorov capacity along a line is defined within the range of filament widths and varies from 0 to 1 as the distribution of intersection points becomes more space filling. When calculated, as in Methven and Hoskins (1998), for the control contour its average and standard deviation were 0.22 ± 0.03. This implies that although the contour becomes increasingly convoluted it does not become space filling.
The histogram of filament widths compiled from all lines of intersection during days 10–12 is strongly peaked at small filaments and has a roughly power-law form (Fig. 4d) with a slope of −1.6 ± 0.1. Without new filamentation events (m = 1) all widths would decrease at the same rate, keeping the range of widths constant, and the histogram would keep roughly the same slope, −α. When m < 1 and the histogram slope is assumed3 to be approximately constant in time, it is possible to show that σn/σ ≈ (m + 1 − α)/(2 − α), where σn is the rate of decrease of the widest filament’s width. In the LC1 case, m ≈ 0.6 and α ≈ 1.6, giving σn/σ ≪ 1. In other words, the median filament width decreases much more rapidly than the maximum, which is often reset by new folding events near the surf zone edges.
While δmed decreases at an exponential rate (see Fig. 5a), the median of displacements initially increases rapidly and then becomes approximately steady before decreasing with time (Fig. 5b). The decrease occurs when the average filament width decreases to a comparable magnitude and more and more filaments are packed into a limited area. At this stage the comparison filaments may be erroneously matched with control filaments that do not originate from the same section of the initial contour; there is so much false detail that the CA simulation is rendered useless. The half-life of a useful CA simulation can be defined as the time taken for half of the filaments to become narrower than the average displacement errors incurred by wind field truncation:
where σ is the rate of decrease of median width, δ0 is the initial median width, and D is the time-averaged median displacement. Its standard deviation is then obtained by combining the standard deviations of σ, δ0, and D in quadrature. Figure 5c shows that D increases as both temporal and spatial resolution decrease. Since the coarse-grained simulations start with the same contour, thus sharing the same δ0, and σ is quite insensitive to spatial truncation (see Fig. 4b), τ generally decreases as wind field resolution decreases (Fig. 5d). Note that the values of τ correspond to the time when the median displacement starts to decrease. At this time it was also found that the fraction of filaments that are unbelievable stops increasing. For all simulations this fraction never exceeds 0.25.
d. The tracer scale factor
The characteristics of displacements and the tracer scale cascade investigated above are now used to define a scale in the tracer distribution above which almost all features are trustworthy. The aim is to relate the trustworthy tracer scale to the wind field resolution. Since displacement errors are found to be approximately constant (for t < τ) we also expect the widths of the widest unbelievable filaments to be roughly time independent. However, this widest width will fluctuate greatly, so the width corresponding to the pth percentile on the histogram of unbelievable filament widths is used to provide a more stable statistic for the smallest trustworthy scale. We denote this scale by δup where u indicates that all unbelievable filaments are considered, but b and s will be used to indicate the sets of bad or spurious filaments separately (see section 2b for these classifications). Here δu100 and δu50 denote the widest and median of the unbelievable filament widths, respectively.
Figure 6a shows that δb100 is roughly constant although clearly quite variable. For smaller values of p the trustworthy scale decreases with time due to the influence of the tracer scale cascade, and by p = 50% the decrease is almost exponential, although at a much slower rate than σ. Figure 6b shows the time average of δb100 versus the smallest resolved scale in the wind field, Lr. The solid line is for a temporal resolution of Δt = 6 h and it is clear that wider filaments are incorrectly positioned when using coarser spatial resolution for the wind field.
At this stage we introduce the notion of the tracer scale factor (hereafter TSF): the ratio of the smallest resolved scale in the wind field to the smallest trustworthy tracer scale. The TSF is estimated from the reciprocal of the slope of the solid line in Fig. 6b by means of a weighted least squares fit. Note that reducing the temporal resolution from 6 to 12 h does not affect the slope of the plot until small Lr when the effects of temporal truncation take over from spatial truncation of the wind field (the dotted line in Fig. 6b). The TSF is clearly a function of the error threshold εt and the parameter p. As εt is increased the criterion for believable filaments becomes less selective, δbp decreases, and the TSF increases. Figure 7a shows the variation of the TSF with εt and p when using only matched filaments that have failed the error criterion ε > εt (i.e., the trustworthy scale is δbp). It is clear that the εt dependence does not vary greatly with p, indicating that the fluctuations inherent when p approaches 100% do not affect conclusions about the TSF. Fitting the dependence
it was found that c = 1.0 ± 0.1 for all values of p. Thus the TSF is proportional to εt, meaning that T1(p) is approximately the smallest resolved scale in the wind field divided by the typical displacement error for bad filaments. The error in the weighted fit used to obtain the TSF is shown in Fig. 7b. As εt increases the number of bad filaments decreases and the estimate becomes less reliable. However, the fractional error is always less than 20% and T1(100) = 11 ± 1.
It is desirable to estimate the percentage q of filaments with widths greater than the trustworthy scale (δ > δbp) that are bad. Since δbp decreases very slowly with time, q will increase with time as the total fraction of bad filaments rises. Figure 6c shows the variation of q with p on day 11.5 for all simulations. Here q increases as p and therefore δbp are decreased, indicating that the fraction of bad filaments, fbi, in a width range δi ± Δ increases as δi decreases. Figure 6d shows that this is indeed the case for the fraction of filaments in each width bin falling into the spurious and bad categories (fsi and fbi, respectively). Interestingly, although δbp increases with Lr, for fixed p, the q–p relationship does not depend strongly upon Lr. When using δb90 to define the smallest trustworthy scale, the percentage of bad filaments in the trustworthy range is only 5%, even at late times in the CA simulations (from Fig. 6c). Note that the trustworthy scale was defined using δbp, rather than a width δ(q) corresponding to a specified percentage of unbelievable filaments in the range δ < δ(q), because δ(q) varies more rapidly with time than δbp. Also, δ(q) drops to zero as q exceeds its maximum value, which is less than 100% and not known a priori, while δbp only decreases to zero as p tends to zero.
The fraction of spurious filaments, fsi, is almost constant below 100 km, although 25% of filaments with δ < 20 km are spurious (see Fig. 6d). However, the fraction of bad filaments, fbi, varies with εt; for εt < 1 the bad filaments largely determine δup but, for εt ≥ 1 the spurious filaments play an important role. Indeed, Fig. 7c clearly shows that as εt increases spurious filaments dominate δup and the TSF value obtained using bad and spurious filaments levels off. The fractional errors in the TSF are always less than 20% as before (Fig. 7d). The net impact of the spurious filaments is to reduce the TSF for large εt.
In summary, almost all filaments seen in a CA simulation of tracers in midlatitude weather systems are believable when wider than the “smallest trustworthy tracer scale,” given by Lr/TSF. When the set of unbelievable filaments includes all those that are displaced by a distance greater than their width or spurious, then the TSF = 6.1 ± 0.3 (εt = 1, p = 100%). This is a conservative estimate: the TSF should be multiplied by the shape factor arising from the definition of the smallest resolved scale in the wind field (see section 2a). It should also be noted that many filaments narrower than suggested by the TSF and the wind field truncation will be accurately simulated.
3. Advection by the stratospheric polar vortex
a. Contour advection simulations
In the previous section we presented results concerning the accuracy of tracer simulations using coarse-grained representations of the wind field. These results will obviously depend upon the characteristics of the flow, such as the dominant timescales and power spectrum of the wind field. Therefore, to broaden the scope of the results a second set of CA simulations was performed for the wintertime Northern Hemisphere (NH) stratosphere. This flow is very different from the baroclinic wave life cycle; the evolution of the wind field is much slower and the flow is dominated by a single vortex. The winds on the 450-K potential temperature surface were obtained from European Centre for Medium-Range Weather Forecasts (ECMWF) analyses for 16–28 January 1992. Advection during this period has been extensively examined in the literature; Plumb et al. (1994, hereafter P94) have described in detail the PV field simulated by CA. They showed that the vortex is greatly disturbed and is almost split into two by an intrusion of midlatitude air. Stirring during this period is particularly active. Schoeberl and Newman (1995) found that contour lengthening rates were as much as three times higher than when the vortex was relatively undisturbed. Also, the transport of air across the vortex edge, as defined by coarse graining the contour from a CA simulation, was over twice as great as any other fortnight in the winter of 1991/92 (Waugh et al. 1994). As a consequence of the enhanced activity, many finescale filaments are formed providing suitable contours for a detailed investigation of the tracer advection problem. Moreover the correspondence between the wider filaments in PV simulated by CA, and filaments in chemical tracers (N2O and ClO) observed by the National Aeronautics and Space Administration ER-2 aircraft, was verified during this period by Waugh et al. (1994). P94 also found that wide filaments and the vortex edge detected by lidar were collocated with high-PV gradients from CA simulations. The broader features (Δx > 150 km) were indeed simulated accurately by CA using coarse-grained winds to the extent that their displacement was less than their width. This was true even though the effects of filament erosion and diabatic motion are ignored in CA simulations.
The use of atmospheric analyses for CA simulations of PV introduces the added complication of contour initialization. The analyses in the stratosphere are constrained by very few observations with the result that smaller scales in winds and temperature can be rather unreliable. The PV field derived from winds and temperature is particularly noisy and does not capture the majority of finescale filaments observed by aircraft or found in CA simulations. However, during a CA simulation parts of the contour that lie at a large angle to the wind vectors are rapidly sheared over to form narrow filaments. In this way, initialization errors are projected onto ever smaller scales while wider filaments formed during the CA simulations tend to be accurately positioned. The issue of the spinup of CA simulations from uncertain initial contours will be discussed in section 3c.
The methodology for the stratospheric investigation was very similar to that for the life cycle. The control CA simulation was obtained by advecting with winds derived from ECMWF uninitialized spectral analyses at a resolution of T213, L31, and Δt = 6 h. As for the life cycles, winds, potential temperature, and Ertel PV were obtained on a quadratic Gaussian grid by transforming from the basic prognostic variables of the ECMWF model: vorticity, divergence, temperature, and surface pressure. Finally, winds and PV were interpolated onto the 450-K surface using the cubic interpolation method described in section 2a and Methven (1997). As noted above, the PV field in the analyses is rather noisy and the shape of contours is therefore sensitive to spatial truncation. Unless otherwise stated, the initial PV contour was obtained from analyses truncated to T42 that were quite smooth. Higher-resolution PV contours were not used because the added complexity of the contours very rapidly resulted in very finescale filaments that made the CA calculations prohibitively expensive before time-dependent straining could generate many new filaments. Attention was focused on the 28-PVU contour, both because it was topologically simple on 16 January 1992 and because it was the contour that formed the greatest variety of filamentary structures over the next 12 days. Figure 8 shows the evolution of the control contour over the entire period. The evolution of other PV contours over this period is shown in Fig. 4 of P94. Over the first three days filaments largely dependent upon initial contour positioning were formed while the whole vortex folded in the middle. As the folding progressed new filaments were formed that thinned as they were stretched and folded again. After day 8 the contour was sufficiently convoluted to contain a wide enough range of filament widths for the advection accuracy investigation.
Comparison contours were then obtained by advecting the same 28-PVU contour using coarse-grained wind fields. Temporal resolutions of 6, 12, and 24 h were used and spectral truncation was employed to smooth the wind fields spatially (highest retained total wavenumber N = 159, 106, 85, 73, 63, 57, 52, 47, 42, 37, 31, 26, and 21). Figure 9a shows that the contour simulated using T42 winds bears a striking similarity to the control contour (Fig. 8d). However, the coarser temporal resolution of 12 h results in an artificially smooth contour (see Fig. 9b). In particular the wavy edges of the wide filament in the 45°–90°E sector are straightened out. These features may be propagating as Rossby waves on the high-PV gradients or may arise through “bulls-eye” noise in the analyzed wind field. Their characteristic timescale must be less than 12 h; the difficulty of simulating fast growing and propagating waves on PV filaments is discussed in section 4b.
b. Estimation of the tracer scale factor
The positions and widths of filaments were measured using cuts across the contours as explained in section 2b. This information was then used to assess the accuracy of advection by coarse-grained stratospheric winds by the same method as for the midlatitude weather systems.
The exponential contour lengthening rate is depicted as a function of wind field truncation in Fig. 10a. As in the weather system experiment, the lengthening rate varies very little with spatial truncation, although there is slightly more sensitivity for truncations coarser than T42. As demonstrated before the sensitivity to temporal resolution is greater than to spatial smoothing. Also, the rate of decrease in the median filament width (Fig. 10b) is greater than the contour lengthening rate, again reflecting the appearance of new, thin filaments as well as stretching. Note that the stretching rates are in accord with the estimates by Schoeberl and Newman (1995) for the same period and approximately half those in the baroclinic wave due to the slower wind field evolution.
Figure 10c demonstrates that the relationship between the number of filaments intersecting a line and median filament width given by (1) also holds for the stratospheric experiment. The best-fit slope (m) using the control and all comparison experiments (with Δt = 6 h) is 0.51 ± 0.01, illustrating that the rate of increase in the number of filaments is about half of the rate of decrease of median filament width. Again, m < 1 highlights the importance of new filamentation events in stratospheric flows and the likelihood of scale invariance in tracer distributions. The Kolmogorov capacity has no time dependence, as for the LC1 experiment, but a greater average value of 0.26 ± 0.05. This hints at a slightly more space-filling arrangement of filaments in stratospheric tracers. The histogram of filament widths (not shown here) has an approximate power-law shape but with a smaller exponent than the LC1 experiment (α = 1.1 ± 0.1). The steady cascade model (section 2c) predicts that the maximum filament width decreases at approximately half the rate of the median.
As for the baroclinic wave, the displacement errors are approximately constant in time. The time-averaged median displacement is shown as a function of wind field truncation in Fig. 11a. These displacements in conjunction with the filament thinning rate define the useful half-life τ of the stratospheric CA simulations. As the wind field is smoothed or temporal resolution is degraded, τ decreases (Fig. 11b). The value obtained here for τ is shorter than for the baroclinic wave experiment (only 8 days) and much shorter than the length of stratospheric CA simulations in other studies. Of course this period is a particularly active one; normal stretching rates maybe three times lower giving a useful half-life of roughly 24 days [see Eq. (2)]. The value is also small because the median of filament widths is biased to very low values by the many spurious filaments in the stratospheric simulations that have no match with the control and therefore no contribution to the displacement estimate. Indeed, for δ < 20 km only 40% of filaments are believable (ε < 1) while 40% are spurious and 20% are matched but bad (cf. Fig. 6d). Since most spurious filaments are much narrower than their parent filaments, erroneous matching rarely occurs before day 12 and the median of displacements does not decrease for t > τ, unlike in the baroclinic wave experiment. The CA simulations could readily be extended beyond τ by removing filaments smaller than the smallest trustworthy scale using contour surgery (Dritschel 1989a); the vast majority of the filaments removed would be those arising from spurious filamentation events.
The width of the widest unbelievable filament reaches an approximately steady value, δu100, which is a function of the spatial and temporal truncation of the stratospheric analyses. The TSF is obtained by the same method as in section 2d from the reciprocal of the slope of the solid line in Fig. 10d. As before, the TSF is a function of the displacement error threshold εt and the position on the bad filament width histogram, p, chosen to represent the trustworthy tracer scale δbp. Figure 11c shows that the TSF roughly fits dependence [(3)] and as before the TSF is approximately proportional to εt (c = 0.9 ± 0.1), but now T1(100) = 8 ± 1. In other words, the typical displacement of bad filaments is approximately eight times smaller than the smallest resolved scale in the wind field.
When filaments arising through spurious filamentation are included, the TSF decreases for p = 100% and εt > 1, as seen in the baroclinic wave experiment. Note, however, that the TSF increases for small εt and p because there are very many, narrow spurious filaments that bias the unbelievable filament width histogram to small scales so that δup < δbp. The stratospheric TSF increases more rapidly as p decreases because the histogram of filament widths is shallower than for the LC1 experiment. Reasons for the far greater influence of spurious filamentation in the stratospheric case are discussed in section 5.
When insisting that the smallest trustworthy scale is the width of the widest filament that has been displaced by a distance greater than its width or that has no match in the control, then the TSF = 5.5 ± 0.5 (εt = 1, p = 100%). Surprisingly, the values of TSF (for p ≥ 90%) for the stratospheric flow and the baroclinic wave are almost equal even though the stretching rates and contour topology are very different. Therefore, it is not reasonable to say that one can achieve better accuracy in a CA simulation of a tracer, using winds of given spatial resolution, in the stratosphere than for flow in the midlatitude storm-track regions. Moreover, temporal truncation introduces significant displacement errors for the stratospheric flow even though the stretching rates are only half those in the weather systems. However, it is likely to be the case that nonconservative processes resulting in the erosion of tracer filaments act more vigorously in the troposphere (see section 5).
c. Uncertainties in contour initialization
The results of a CA simulation of a tracer isopleth will be dependent upon the initial shape of the contour. When the initial position is obtained by the interpolation of a 3D tracer field onto the chosen isentropic surface, there will be errors arising from interpolation as well as uncertainties in the analyzed tracer field. For example, in the previous section the contour was initialized following a contour of PV as determined from ECMWF analyses. The PV field itself was noisy, so a smoothed (T42) representation of PV was taken for the contour definition. If the same value PV contour had been taken from a less-smoothed PV representation, or from analyses by another meteorological center, then the contour would inevitably retain differences in detail throughout the CA simulation. For example, P94 show the dramatic differences between the ECMWF and National Meteorological Center (NMC, now the National Centers for Environmental Prediction) analyses of PV on the 450-K surface. However, if the same winds were used for advection, then these initialization errors would tend to propagate to ever decreasing length scales, and larger filaments formed by advection could be accurately simulated. The goal here is to estimate the spinup time taken for initialization errors to become less significant than the effects of wind field truncation.
An extra CA simulation, C106, was performed to investigate the effects of contour smoothing; the 28-PVU contour was obtained from ECMWF analyses at T106 resolution, rather than T42 resolution as used for the earlier simulations, and then advected by winds also at T106 resolution. The T106 contour shows much more structure (cf. Figs. 8a and 9e). Here it is not possible to determine whether these features are realistic or the result of a truncated spectral representation of a field with finescale filaments, although P94 suggest that the bulls-eye features may reflect the impact of assimilating isolated radiosonde observations. However, we will treat the “control” run, which starts with the T42 contour and advects with T213 winds, as ground truth and treat the finescale filaments, which appear rapidly as the bumps in the T106 contour are sheared out, as spurious. The filament matching procedure of section 2b is then used to associate filaments from the two simulations and find their widths and displacements.
Figure 12 shows the same measures of filament width and displacement as discussed in section 2c. The CA simulations obtained by advecting the T42 PV contour with T106 and T42 winds are also shown for comparison (labeled as W106 and W42). The median of filament widths (Fig. 12a) clearly illustrates that the many small filaments arise from advecting the T106 contour that are not present with the T42 contour. The added complexity is retained since all filaments are stretched and thinned at the same exponential rate. The median of unbelievable filament widths (Fig. 12b) shows a very rapid decrease in the C106 simulation until day 2 or 3. This is also the time at which unbelievable filaments begin to appear in the simulations using coarse-grained winds (W106 and W42). After this time the median unbelievable width decreases at a similar exponential rate for all simulations. However, the maximum of unbelievable filament widths (Fig. 12c) is approximately steady, as discussed previously, with the W106 simulation being significantly more accurate than W42 and C106.
The most interesting measurement is the median of displacements, D (Fig. 12d). Until day 3, D is much greater for the C106 simulation than for W106; this can only arise through differences in the initial contour. After this period D corresponds very closely for these two simulations, indicating that the contour initialization errors have been swamped by the effects of wind field truncation. Note that D is only defined for filaments that are matched in the control. As much of the C106 contour aligns itself with the dilatation axis of local strain, the spurious filaments thin rapidly, and the edges of the matched filaments approach the positions that would have been attained by smoothing the contour before the CA simulation (W106). To some extent the importance of the initial conditions diminishes and the positions of these matched filaments are then subject to the displacement errors arising from wind field truncation.
The evidence presented above suggests that it is possible to detect a spinup time for CA simulations; for initialization uncertainties arising from smoothing stratospheric contours, its length is 2–3 days. This timescale is related more to the strength of the shear at the polar vortex edge than to the characteristic filament thinning timescale (σ−1) resulting from stretching and folding in the surf zone. It is hypothesized that the spinup time is less sensitive to the exceptional flow activity in January 1992 than the useful half-life of the CA calculation [see Eq. (2)].
Since the displacement errors decrease rapidly during the spinup period, it would seem sensible to start the CA simulation as early as possible. Sutton et al. (1994) and Schoeberl and Newman (1995) have advocated the use of 5-day trajectories in order to have time to develop a desirable degree of contour complexity without introducing too much false detail. However, a contour obtained from analyses on a particular day is not simply a coarse-grained representation of a contour from a mature CA simulation due to diabatic effects and large uncertainties in the analyses. For example, Figs. 5b and 11 of Waugh et al. (1994) illustrate substantial differences between the NMC PV analysis and PV obtained by CA simulation and then coarse graining. Therefore, it is anticipated that, at time t3, the similarity between a contour from a CA simulation initialized at time t2 from analyzed PV and an older contour (initialized at time t1 < t2 < t3) will be less than between the old contour and a contour that is obtained by smoothing the old contour at time t2 and then restarting the CA simulation. This additional problem was studied with another set of simulations conducted by initializing contours at dates staggered over 5 days, from 16 January 1992 to 21 January 1992, and comparing their positions from day 8 to 12 (24–28 January). Comparison of Figs. 9c and 9d illustrates the lower degree of complexity in younger contours. They have fewer finescale filaments, although the “envelope” contour shape is very similar in all cases. Figure 13a shows the time-averaged median of filament displacements relative to the oldest contour. The next oldest contour, started on 17 January 1992, is obviously closest to the control contour but the four younger contours have similar displacements that are much larger than those incurred by wind field truncation alone (cf. Fig. 11a).
The triangles in Fig. 13b mark the smallest trustworthy scale δu90 averaged from day 8 onward by which time even the youngest contour has spun up. When εt = 1, δu90 for the youngest contour is much greater than for the older contours, but when εt is doubled the values are quite similar for all contour ages (see the diamonds in Fig. 13b). The widest unbelievable filaments arise through spurious filamentation because the relative positioning of the contour and hyperbolic points changes with each analysis. While smoothing an old contour simply aligns the contour with streamlines, the PV contours from the latest analysis do not share this alignment even when smoothed. This indicates that the analyses do not evolve as close to a balanced flow as the ECMWF model would do without the assimilation of sparse stratospheric observations. Diabatic effects would also result in the nonconservation of PV between analyses.
The insensitivity of the smallest trustworthy scale to contour age is supported by calculating displacements when using the youngest contour as the control. In this instance, many of the extra filaments present in the oldest contour that were lost when the oldest contour was defined as the control will be classed as spurious and fall into the unbelievable category. Therefore, only the bad filaments will be considered to contribute to the smallest trustworthy scale (δb90) for this comparison. This scale is again around 50 km and varies weakly with contour age (see the squares in Fig. 13b). The evidence suggests that the use of a contour that is much older than the spinup time does not result in a significant improvement in the simulation of tracer isopleths. This rather negative result is supported by Fairlie et al. (1997) who compared the prediction of along-flight variations in N2O using PV analyses and reverse-domain-filled (RDF) PV fields.4 They found that although the use of PV information gave skill in predicting along-flight N2O concentrations, no significant improvement was found with the Lagrangian (RDF PV) forecast. They also pointed out that sideways filament displacements would result in large rms differences between observed and predicted N2O concentrations and significant reductions in correlation coefficients. However, the use of RDF PV, rather than analyzed PV, did not even result in a significant improvement in the cumulative frequency distribution of predicted N2O values, which is insensitive to filament displacements. Although the lack of skill improvement may arise partly through the approximate nature of the linear relationship between PV and N2O assumed in their study, it appears that current PV analyses contain too much error to simulate accurately most PV filaments (δ < 100 km) using Lagrangian techniques such as CA. Indeed, only the polar vortex edge and filaments with widths of 150 km (P94) and 400 km (Waugh et al. 1994) were clearly collocated in observations by the ER-2 aircraft and CA simulations of PV (advecting with winds on a 5° grid with Δt = 24 h). More accurate simulation would be possible if the initialization uncertainties amounted simply to smoothing the true PV field and if diabatic effects were weak.
In summary, when uncertainties in the initial tracer contours are sufficiently small, CA simulations yield useful results in the interval bounded by the spinup time and the useful half-life τ. For stratospheric simulations the spinup time should be at least 3 days while τ depends on the smoothness of the initial contour and filament thinning rates [see Eq. (2)]. The latter limit can be circumvented by the use of contour surgery. Unfortunately, when the tracer analyses are more uncertain, the accuracy of CA simulations can be severely impeded. When using the ECMWF analyses it was found that PV filaments smaller than about 50 km could not be trusted for a wide range of contour ages. In this situation the errors from initialization outweigh those from wind field truncation in the period before most filaments have thinned to unbelievable scales.
4. Dynamics of PV filaments
a. Arguments based on wavenumber spectra
So far we have examined the advection of tracers by prescribed wind fields and found that tracer filaments can be simulated accurately with widths that are smaller than the smallest resolved scales in the wind field by a factor (TSF) of about 10. Such behavior is only expected to occur when the larger scales in the wind field have a dominant influence on the positioning of filaments so that the neglect of the smaller scales results in small displacement errors.
On a sphere, a suitable measure of the power of a field as an isotropic function of scale is given by the n spectrum (Boer 1983; Lander and Hoskins 1996). Figure 14a shows the kinetic energy (KE) n spectrum (here defined as the velocity covariance spectrum) for the LC1 experiment at day 10. There is no obvious cutoff that would suggest that large scales dominate the wind field but there is an obvious change in spectral slope at N ≈ 100. The range 16 < N < 106 has a spectral slope5 of 3.0 ± 0.2 while the higher wavenumbers fit a steeper slope of 5.6 ± 0.1. After day 10 these slopes remain approximately steady. The same features are also seen in the KE spectrum of the stratospheric analyses, although the most obvious change in slope occurs for N ≈ 50 (see Fig. 14b). Over the 12-day experiment the slope in the range 10 < N < 50 lies between 2.5 and 3.1, and in the range 50 < N < 110 lies between 6.5 and 7.1. There is also a shallower tail end to the spectrum for N > 130 that is related to the noise in the uninitialized ECMWF analyses. Again, it is clear that the power in the wind field remains in the large scales throughout the flow evolution.
Kinetic energy spectra are expected to relate to the spectra of PV covariance (or potential enstrophy). Indeed, the ranges of power-law behavior seen in the KE spectra are also seen in the PV covariance spectra (not shown here). The power-law ranges may suggest that there are a number of flow features of different scales and that a smooth change in coarse graining of the wind field might result in a smooth change in the accuracy of tracer simulation. This would support the notion of a global value for the TSF for a given flow. However, a power-law power spectrum can be obtained for a single step or high gradient interface in a field. The power-law ranges in the PV covariance spectrum might be thought to suggest the existence of an ensemble of eddies when it is more likely to be indicative of PV filaments with sharp edges. The KE spectrum would also have a corresponding power-law behavior with a steeper slope due to the smoothing effect of PV inversion. Even when there are many PV filaments with a range of scales and one has prior knowledge of the topology of filaments, it is extremely difficult to infer any geometrical information from spectra. For instance, Methven and Hoskins (1998) showed that spectra were too crude of a measure to distinguish between the spirals formed in tracers by weather systems and regularly spaced sets of parallel filaments. The structure of the wind field, which requires detailed knowledge of the phases between wave components, is crucial to the appearance of hyperbolic points and therefore to the tracer advection problem. For instance, the temporal resolution of the wind field was found to have a strong influence on contour stretching rates (see Figs. 4a and 10a), presumably through the time dependence of hyperbolic points in the flow. Thus it is not possible to estimate the TSF from KE spectra alone.
Nevertheless, it is clear that large-scale strain continues to dominate the advection even though the widths of PV filaments decrease rapidly. This is somewhat surprising since in regions where the growth rate of shear instability is large, the small-scale structures in PV and velocity may play the dominant role in tracer advection. In isolation a PV filament is unstable to lateral perturbations; the simplest example of this mechanism is the Rayleigh instability of a strip of vorticity in 2D incompressible flow. The fastest growing wavelength scales with the width of the filament, and as wave slope reaches order one, the nonlinear dynamics become important and the strip rolls up into a chain of spirals (e.g., Dritschel 1989b). These spiral-shaped vortices have substantial circulations on scales similar to the wavelength. If all PV filaments were active the strain scales would be expected to decrease with the scales in PV. This does not happen due to three features of PV inversion.
First, a large-scale PV anomaly induces a greater circulation than a smaller PV anomaly of the same magnitude; this is termed the scale effect of inversion by Hoskins et al. (1985). Thus narrower filaments in isolation have a lesser influence on the flow. Second, PV inversion has a smoothing effect so that the streamfunction resembles a coarse-grained representation of arrangements of PV filaments. For instance, when a filament has wrapped up into a spiral it will possess a large-scale circulation even though the filaments themselves are thinning at an algebraic rate. The extent of the smoothing and scale effect depends upon the nature of the PV inversion operator for a balanced flow that closely approximates the true flow. Pierrehumbert et al. (1994) have analyzed α-turbulence models where streamfunction ψ is related to PV, q, in Fourier space by ψ̂(k) = −|k|−αq̂(k), where k is the 2D wavenumber vector. As the parameter α increases, the streamfunction becomes smoother and vortex interactions become more nonlocal. In these models the slope of the velocity covariance spectrum, p′, is related to the slope of the PV covariance spectrum, p, by p′ = p + 2(α − 1). For 2D nondivergent barotropic flow α = 2, while α = 1 applies to surface quasigeostrophic dynamics (Held et al. 1995).
Finally, large-scale strain due to far-field PV anomalies can suppress the instability of PV filaments. For the case where a straight filament is placed in a uniform strain, Dritschel et al. (1991) have shown that no amplification of perturbations is possible when the ratio of strain rate to strip vorticity exceeds 0.25 and that the barotropic growth rate is significantly suppressed for much weaker strain. Pierrehumbert et al. (1994) have demonstrated the change in filament “activity” with the interaction range of PV anomalies by using their α-turbulence model. For α < 2 filaments succumb to secondary instability, while for α > 2 the strain due to distant PV anomalies always suppresses instability. Here α = 2 represents the marginal case that is sensitive to the initial PV distribution. Waugh and Dritschel (1991) have also demonstrated that roll-up is less inhibited when the interaction range of PV anomalies is shorter.
In the lower wavenumber range of the stratospheric experiment p′ − p = 2.0 ± 0.2, while in the lower wavenumber range of the LC1 experiment p′ − p decreases from 1.6 to 1.2 between days 8 and 12. However, in the higher wavenumber ranges of both experiments p′ ≈ p. These observations are indicative of a shorter interaction range (α → 1) for PV anomalies at the smaller scales. We tentatively suggest that the close correspondence of TSF values for the weather system and lower stratosphere arises through a similarity in the spatial properties of PV inversion for the two flows even though the timescale for flow evolution is twice as long in the stratospheric case.
b. Simulation of unstable PV filaments
The accuracy of CA simulations of unstable PV filaments is now investigated and the influence on the global value of TSF is discussed. Spatial and temporal truncation of the winds affect the simulated amplitude of contour undulations in the vicinity of a propagating and growing wave in two ways. Its amplitude is underestimated by sampling winds infrequently in space and time when the wave propagates rapidly in a frame moving with the flow [i.e., k(U − c)Δt ∼ 1] because the peaks in the velocity are missed. However, if the disturbance growth rate is sufficiently high, then the amplitude may actually be overestimated by linearly interpolating between wind records. Spatial smoothing also has the effect of reducing maxima in velocity and therefore the simulated amplitude.
A striking example of instability occurs on a spiral-shaped PV filament in the LC2 baroclinic wave life cycle (Thorncroft et al. 1993) when performed at very high resolution (Methven and Hoskins 1998). Its growth rate was shown to be in close accord with barotropic instability theory when accounting for the radial shear in angular velocity in the vortex (Methven 1998). Here, the same experiment is used to investigate the accuracy of simulation of PV contours using CA when instability has reached large amplitude. Figure 15a shows Ertel PV on the 300-K isentropic surface at day 9 when the spiral has wrapped up four turns and the wave slope of shear instability has reached order one. Figures 15b,c show the 2.2-PVU contour simulated by CA using winds at T341 and T57 spatial resolutions, respectively, at a temporal resolution of Δt = 6 h. Although many features of the spiral are very similar, two features are poorly represented: the spiral center and the amplitude of the instability.
The center of the spiral is sensitive to spatial truncation. For coarse winds (N < 60) the number of turns simulated by CA decreases as N decreases. This behavior reflects the sharpness of the maximum in angular velocity around the spiral center, ω(r), which decreases monotonically with radius. Away from the center, ω has an r−q radial dependence where q < 2. The value of q depends upon spiral shape (Methven 1998). However, near the center ω(r) flattens off. The radius of this central patch must be close to the smallest resolved scale at T60 (≈330s km).
Spatial truncation also reduces the amplitude of the simulated disturbances so that they are barely perceptible using winds at T57 resolution (Fig. 15c). Since the azimuthal velocity around the LC2 vortex at the radius of maximum disturbance is about 100 km h−1, a temporal truncation of 6 h has comparable effects to a spatial truncation of 600 km (or T30–T40) through the sampling effect. Indeed, the disturbance amplitude has been severely reduced by using Δt = 6 h and spatial truncation is seen to enhance the effect.
When repeating the LC2 experiment at much lower resolution (T85) the PV spiral is not resolved at all and therefore neither is the shear instability. A CA simulation using these winds produces a very smooth spiral-shaped filament (not shown here). However, the CA simulation of PV using the high-resolution winds truncated to T85 does capture the signature of the instability at a late stage. This is because the narrow filament rolls up into mesoscale vortices that are marginally resolved when the high-resolution winds are truncated. When using high-resolution atmospheric analyses there will be instances when a CA simulation will detect the effects of nonlinear roll-up even though the PV filament itself may have been unresolved in the analyses at an earlier stage [e.g., Appenzeller et al. (1996), Figs. 12 and 15 of Waugh et al. (1994)]. However, in order to accurately position the filament when the wave slope reaches order one it is necessary to use a wind field with a spatial and temporal resolution that is sufficiently high to resolve the filament dynamics; the effect of instability is to reduce the TSF toward unity.
The difficulty of simulating unstable filaments means that many CA simulations may produce artificially smooth PV contours. This problem is particularly acute when using coarse temporal resolution. For instance, the undulations in the wide tongue of PV in Fig. 8d are completely smoothed out by halving the temporal resolution of the wind field (see Fig. 9b). Filaments in the stratosphere are unlikely to behave as passively as a CA simulation with coarse-grained winds would suggest. For example, the high-resolution barotropic model of the stratosphere by Juckes and McIntyre (1987), where filament dynamics was resolved, showed many examples of vortex roll-up.
Rapid changes in chemical tracers, observed at a point or along a flight track, often correspond spatially to high gradients at the edge of tracer filaments. If tracer simulations are to be used to infer the origin and history of air contained within such filaments, it is essential to characterize their accuracy. Waugh and Plumb (1994) have demonstrated that qualitatively realistic tracer filaments can be simulated that are several times narrower than the smallest scales retained in the coarse-grained wind field used for advection. In this paper a quantitative measure of the ratio between the wind field resolution and the smallest trustworthy tracer scale, termed the tracer scale factor (TSF), was investigated. The TSF also gives interesting information on the properties of the flow itself that have not been obtained by any other means. The investigation was performed for two very different atmospheric flows: an extratropical weather system and the lower-stratospheric polar vortex. Flow on an isentropic surface that passed through the steering level of a very high-resolution baroclinic wave life cycle was used as a proxy for flow in a typical weather system. ECMWF analyses of wind on the 450-K surface during January 1992 were used for the stratospheric investigation. The highest possible resolution for each flow was used to advect a tracer contour using the contour advection (CA) technique. The evolution of this contour defined the control simulation. Contours simulated using coarse-grained representations of the control wind field were then compared with the control contour.
The errors in the comparison tracer simulations were assessed by matching filaments in common with the control and then using the ratio ε of each filament’s displacement to its width. A variable threshold error, εt, was used to classify filaments for which ε < εt as “believable” and the remaining matched filaments as bad. As εt increases, more filaments fall into the believable category. Although all filament widths decrease rapidly, the typical displacements and the width of the widest bad filaments are tied to the wind field truncation and are approximately constant. The smallest trustworthy scale was defined to be the scale corresponding to the pth percentile on the histogram of bad filament widths (δbp). The aim was to provide a measure that was less susceptible to fluctuations than the width of the widest bad filament (p = 100%). The TSF was then defined through the slope of the graph of the smallest resolved scale in the wind field versus δbp. It was found that the TSF is approximately proportional to εt, meaning that the TSF value for εt = 1 is equal to the smallest resolved scale in the wind field divided by the typical displacement of bad filaments. The TSF increases as p decreases, but for each value of p the dependence on εt is the same. Therefore, the error characteristics are not strongly affected by fluctuations in the tail of the histogram of bad filaments. For the baroclinic wave the TSF was found to be 11s ± 1s (εt = 1, p = 100%), where the shape factor takes on a value 1 < s < 3.2 depending upon the smallest features in the wind field. In the stratospheric experiment a very similar value of TSF was obtained: 8s ± 1s.
In a snapshot of a comparison contour some filaments may have no match with the control and these are classified as spurious. It is necessary to include the effects of spurious filaments on the TSF. In particular, in the stratospheric simulations there is a much higher fraction of spurious filaments than for the baroclinic wave, although they are typically extremely narrow. This is attributed to the sensitivity of weak hyperbolic points near the vortex edge to wind field smoothing. Some of these hyperbolic points may exist because of the noise incurred by the assimilation of sparse observations. When using the histogram of bad and spurious filaments to redefine the smallest trustworthy scale, the TSF for the baroclinic wave was found to be 6.1s ± 0.3s (εt = 1, p = 100%). In the stratospheric experiment a very similar value of TSF was obtained: 5.5s ± 0.5s. These are the most pessimistic estimates for the TSF.
It was shown that the position of an unstable PV filament can only be simulated accurately, when its wave slope exceeds 1, by winds that have sufficient resolution to resolve the dynamics of the PV filament itself. Locally, the TSF tends toward unity, reducing the usefulness of Lagrangian simulation of PV or any other tracer that crosses the unstable region. In contrast, chaotic advection studies represent the opposite extreme where a large-scale wind field is specified and passive tracers cascade to arbitrarily small scales (TSF ≫ 1). The atmosphere lies at an intermediate point (TSF ≈ 10) due to the scale effect of PV inversion and because far-field PV anomalies suppress the instability of most PV filaments. It was suggested that the similarity in the spatial properties of PV inversion, crudely characterized by the relationship between the power spectra of velocity and PV, partly explains the correspondence between TSF values for the two flows studied. Nevertheless, this result is rather surprising since contour stretching rates are twice as large in the baroclinic wave and the contour topology is very different.
Contour stretching rates have been found to be very insensitive to the spatial truncation of the wind but rather sensitive to temporal truncation. Ngan and Shepherd (1998) have pointed out that this result is consistent with the chaotic advection picture where the autocorrelation timescale for velocity is much shorter following trajectories (τL) than in the Eulerian framework (τE). They found τE ≈ 2.5 days for flow in the lower stratosphere while τL was 4 ± 2 times shorter. Haynes and Anglade (1997) investigated the scale cascade in passive tracers and found that a random straining model of the winds, where the deformation tensor is parameterized using a decorrelation timescale for strain and an rms strain rate, gave quantitatively similar stretching statistics to those obtained using stratospheric winds even though details of the spatial structure of the strain field were not specified. Here, contour stretching rates are found to be underestimated unless Δt ≪ τL. The insensitivity to spatial truncation indicates that smoothing does not radically alter the basic structure of the strain field.
If one is not interested in the determination of the origin of air within individual filaments, then the requirements for the wind field are not so strict. The insensitivity of stretching and filament thinning rates to the spatial truncation of the wind field mean that the typical scales and gradients in simulated tracers will be well represented even if the details of filament positions are inaccurate. Moreover, it was found that the parent-to-daughter width ratio in filamentation events, and thence scale invariance in filament distributions, also appears to be insensitive to wind field truncation. In this case offline chemical transport models that use coarse-grained winds to advect higher-resolution tracers should be able to reproduce the gross statistics of the chemical distributions. The importance of resolving microstructure in chemicals that are stirred together when the volume-averaged reaction rate is limited by the mixing efficiency of the flow has been highlighted by a number of recent studies (e.g., Edouard et al. 1996; Thuburn and Tan 1997).
The filamentary structure in the stratospheric experiment was found to be sensitive to the smoothing of PV before contour initialization. Contours that were obtained from unsmoothed analyses retained greater complexity for all time because the median of filament widths decreased at a rate, σ, that was independent of initial contour shape. However, a spinup period of 2–3 days was found, after which the initial details in the unsmoothed contour had been sheared to such small scales that the displacement errors accumulating through wind field truncation had become more important. In addition, a useful half-life, τ, for CA simulations was defined to be the time when half of the filaments had become thinner than the average displacement. Its value depends upon the smoothness of the initial contour and the filament thinning rate (τ ≈ 3σ−1 for the stratospheric case). It was suggested that CA simulations should only be extended past this time if contour surgery (Dritschel 1989a) is used to remove the filaments that are narrower than the smallest trustworthy scale.
Although stirring results in an exponential rate of thinning of tracer filaments on isentropic surfaces, persistent erosion by nonconservative processes will eventually result in the mixing away of the oldest, thinnest filaments. Every filament on an isentropic surface represents a slice through a long, thin tracer anomaly. Typically these anomalies are much shallower than they are wide. However, the erosion processes often act in a more isotropic manner with the result that mixing in the vertical limits the lifetime of tracer anomalies rather than layer-wise 2D mixing. For the lower stratosphere, erosion of tracers can occur through sporadic episodes of clear air turbulence caused by breaking gravity waves and also through scale selective radiative damping (Haynes and Ward 1993). Balluch and Haynes (1997) obtained an upper limit to the effective vertical diffusion coefficient of 0.01 m2 s−1 by comparing the profiles of tracer filaments from airborne observational campaigns with profiles obtained by modeling sheets of tracer eroded by vertical diffusion while thinning due to horizontal strain and vertical shear. Near the tropopause the inhomogeneity in tracer concentrations and turbulent activity make quantitative estimates of mixing very difficult. Shapiro (1980) has deduced a value as large as 20 m2 s−1 in the upper troposphere just below tropopause folds; this is a region of high Richardson number where Kelvin–Helmholtz instabilities are expected to occur. Tracer filaments on isentropic surfaces that pass through the upper troposphere and lower stratosphere may therefore experience an effective eddy diffusivity that varies by three orders of magnitude!
Assuming that 3D turbulence dominates erosion in the lower stratosphere, the depth of the smallest-scale filaments surviving erosion, H, is given by the balance scale between vertical diffusion (κυ) and the typical horizontal stretching rate (σ), which is (κυ/σ)1/2 (Haynes and Anglade 1997). Taking κυ = 0.01 m2 s−1 and σ = 0.4 days−1 gives H = 50 m. Assuming6 L = NH/f, the smallest width is found to be about 10 km. Contour advection simulations of tracers could be prolonged using contour surgery to eliminate filaments that fall below this lower scale. Assuming that an appropriate value of TSF is ten, the coarsest wind field that could accurately simulate all the retained filaments would be T200. This corresponds to a smallest resolved scale in the wind field of about 100s km.
Although ECMWF analyses, for example, have a nominal resolution of T213, the observations in the stratosphere are extremely sparse introducing significant large-scale uncertainties into the analyses. Since CA simulations of PV contours look qualitatively more realistic than the PV derived from the stratospheric ECMWF analyses it might be supposed that the accuracy of PV filament positioning improves as contours become older (i.e., the contours have been advected over longer intervals). However, it was found that, once older than the spinup period, there was almost no gain in accuracy by increasing contour age. The smallest trustworthy scale was found to be about 50 km. This scale is several times larger than arising from wind field truncation at the resolution used (T106) and differences must be due either to nonconservation of PV in the stratosphere or errors in the analyses. Since the departure of the analyses from the picture obtained by coarse-graining CA-simulated contours occurs over 1 day, it is not likely to be due to diabatic processes alone and therefore it appears that analysis error is significant. The many spurious filamentation events in CA simulations arise when the movement of contours relative to hyperbolic points changes, as described by Polvani and Plumb (1992). The magnitude of deformation at many of these hyperbolic points is small so that wind field smoothing can readily move or eradicate them. It seems likely that some of these points are introduced by the assimilation of sparse observations rather than through the evolution of a balanced flow. For flight planning purposes in stratospheric chemistry campaigns it would appear to be important to locate regions containing finescale filaments using CA simulations generated by wind analyses from a number of different meteorological centers (e.g., Waugh et al. 1994).
Recent work by N. Nakamura (1998, personal communication) demonstrates that even contours for tracer that are initially uncorrelated gradually align with each other in steady flows at a rate that is dependent on transport and diffusion timescales. In unsteady flows like the stratospheric example it appears that small-scale differences between tracer contours are rapidly sheared to scales below the smallest trustworthy scale (related to the tracer scale factor) within a spinup period. However, if the initial differences are large scale and possess sufficient amplitude for the two contours to be affected by slightly different sets of hyperbolic points in the flow, then alignment does not occur rapidly, to the extent that the resolution of the wind field can become secondary consideration in the accuracy of tracer simulations.
Latent heat release impacts on the dynamics of weather systems, tending to enhance their growth rates and to enhance and localize vertical motions, for example, Emanuel et al. (1987). The gross features of cyclogenesis are similar on the inclusion of latent heat release with the main impact arising through a reduction in effective static stability. However, latent heat release also directly rearranges the PV distribution. For instance, positive PV anomalies near the level of diabatic heating produced by “steady condensation” in moist ascending conveyor belts can have an influence on cyclone development as great as upper-level troughs and surface temperature gradients (Wernli and Davies 1997). Low-level strips of PV along the trailing fronts of large-scale cyclones can roll up through shear instability forming frontal wave secondary cyclones (Bishop and Thorpe 1994). In some cases, latent heat release and surface fluxes result in mesoscale vortices that otherwise may not have occurred. For example, Craig (1993) has examined an outbreak of cold, polar air over the warm surface of the North Atlantic that encouraged the development of an amazing variety of polar lows. These examples suggest that small-scale PV anomalies play a more active role than they would in a dry atmosphere and therefore that the TSF would be lower.
Nevertheless, the validity of isentropic tracer simulations is supported by the relationship between water vapor imagery and PV. For example, Appenzeller et al. (1996) observed that features in PV, simulated by CA on the 320-K surface, matched those in contemporaneous water vapor images to such an extent that small-scale vortices arising from the roll-up of PV streamers were collocated and spiral structures were seen within these vortices on both images. The success of their simulations indicates that sequences of PV analyses on 320 K are more continuous than on 450 K, in the sense that they bear a closer coarse-grained resemblance to PV contours from mature CA simulations, and that motion in the dry intrusions must be nearly isentropic. The initial contour obtained from the smoothed PV field must also be collocated with the high gradients in water vapor across the tropopause to such a degree that the usefulness of the tracer simulation is not limited entirely by the initial conditions.
The contour advection code was adapted from David Dritschel’s contour surgery algorithm and made available for the use of the U.K. Universities’ Global Atmospheric Modelling Programme (UGAMP) by Warwick Norton. Thanks to Jeff Cole, Lois Steenman-Clark, and Mike Blackburn for parallelizing the Reading spectral model for use on the CRAY T3D at the Edinburgh Parallel Computing Centre. This enabled integration of the very high-resolution baroclinic wave life cycles. The ECMWF analyses used in the stratospheric experiment were obtained by Paul Berrisford and Kate Searle (for UGAMP). The reviewers’ comments provided useful motivation for exploring the sensitivity to the definition of accuracy. John Methven was funded by a Natural Environment Research Council studentship, GT4/92/284/P, and more recently through a grant from the U.K. Meteorological Office.
Corresponding author address: John Methven, Massachusetts Institute of Technology, Room 54-1721, 77 Massachusetts Ave., Cambridge, MA 02139.
1 PVU = 10−6 m2 s−1 K kg−1.
For each snapshot of a contour, δmed is calculated from the filaments on all intersecting cuts and n is the total number of filaments divided by the number of cuts.
This steady cascade model assumes a histogram of the form ni = (δi/δn)−α with a wide range of filament widths (δ1/δn ≪ 1) in a surf zone whose size is constant (L = Σni = 1 niδi). It is only valid for 1 < α < 2, 0 ⩽ σn/σ ⩽ 1, and therefore α − 1 ⩽ m ⩽ 1. If σn/σ = 1, one finds m = 1 but α is not determined.
RDF PV fields are constructed by assigning the analyzed PV at the start of 5-day back trajectories to each arrival point on a domain-filling grid. This technique was first used by Sutton et al. (1994).
Power-law spectra have the form k−p where p is the spectral slope.
L = NH/f is appropriate to passive tracers when the decorrelation timescales of strain and vertical shear are longer than the inverse rms strain rate and quasigeostrophic scaling applies to the ratio of shear to strain (N/f) (Haynes and Anglade 1997).