# A wavelet representation of synoptic-scale coherent structures

Matthieu Plu, Philippe Arbogast, Alain Joly

► To cite this version:

Matthieu Plu, Philippe Arbogast, Alain Joly. A wavelet representation of synoptic-scale coherent structures. Journal of the Atmospheric Sciences, 2008, 65 (10), pp.3116-3138. 10.1175/2008JAS2618.1. meteo-00330026

HAL Id: meteo-00330026

https://meteofrance.hal.science/meteo-00330026

Submitted on 9 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Wavelet Representation of Synoptic-Scale Coherent Structures

MATTHIEU PLU, PHILIPPE ARBOGAST, AND ALAIN JOLY

GAME/CNRM, Météo-France, CNRS, Toulouse, France

(Manuscript received 7 September 2007, in final form 31 January 2008)

## ABSTRACT

Midlatitude cyclogenesis as interpreted in the framework of either baroclinic development or potential vorticity thinking heavily relies on the concept of synoptic-scale anomaly. Given the existence of potential vorticity inversion and attribution, what is at stake to provide a mathematical definition for this concept is a complete finite-amplitude alternative to the linear-based theory of cyclogenesis. The existence of a reasonably objective way to represent anomalies in both real and idealized flows would not only help understanding cyclogenesis, it would also have many other applications for both theory and in practical forecasts. Inspired by the recent theory of wavelet representation of coherent structures in two-dimensional fluid mechanics, a wavelet representation of three-dimensional potential vorticity anomalies is built. This algorithm relies on the selection of the appropriate two-dimensional wavelet coefficients from the stationary wavelet transform in order to guarantee the critical translation-invariance property. The sensitivity of the algorithm to the position, size, and shape of the structures is assessed.

The wavelet extraction is then applied to the upper-level precursor of a real-case storm of December 1999 and is compared to a basic monopolar extraction. Using potential vorticity inversion and forecasts with a primitive-equation model, it is found that both anomalies have similar implications on the development of the surface cyclone. However, the coherence in time of the extracted wavelet structure in the forecast and analysis sequence is more satisfactory than the extracted monopole: this suggests that the underlying mathematical description of an anomaly proposed here does, indeed, point toward the direction of an actual physical reality of the concept.

## 1. Introduction

The history of modern meteorology can be characterized by the theoretical framework with which it addressed the midlatitude cyclogenesis problem. This depiction appears to have oscillated between “two contrasting viewpoints” (Davies 1997): One is that cyclogenesis evolves from “an intrinsic instability of the flow,” and the other “links cyclogenesis to the influence of an upper-level trough-like feature.” During the twentieth century, each viewpoint has been prominent, at one time influencing the way other phenomena were described, and then less considered at other times. One strength of the linear instability approach in the published literature, including books, reviewed in depth by Pierrehumbert and Swanson (1995), is the elegance of its underlying mathematical model of cyclones. This

plane-wave model is directly borrowed from decades of experience in mathematical physics and its many fields of application. Eady (1949) and Charney (1947) unfold this approach to derive a reasonable prediction of the spatial scales of midlatitude cyclones growing against a steady large-scale jet flow. Thanks to the familiar underlying mathematical representation, something like the wavelength, one of the predictions of the theory, appears like a well-posed output fit for verification. However, it is not so simple to see growing plane waves in the real atmosphere. The alternative approach, on the other hand, has been induced from daily experience with weather maps. Its main points were refined and strengthened by their use in the making of the demanding forecast routine required by World War II operations (Petterssen 2001). It has even been formalized mathematically up to a point, as for example by Sutcliffe (1939, 1947) and Kleinschmidt (1951). One of its key tenets is that the new surface cyclone to be explained, or to be forecast, is preceded by some “anomaly” of vorticity or potential vorticity. While these anomalies are both critical and ubiquitous in that

<footnote>*Corresponding author address:* Matthieu Plu or Phillippe Arbogast, Météo-France, CNRM/GMAP/RECYF, 42 avenue G. Coriolis, 31057 Toulouse CEDEX, France.
E-mail: matthieu.plu@meteo.fr</footnote>
framework, their mathematical representation and the
ability to identify and track them unambiguously in any
given situation does not match the powerful plane-wave
model when it comes to formalizing field analysis.

The century-long contest between these views has
been fruitful in that it has left a common heritage. Part
of it is the quasigeostrophic theory and framework, as
well as the related invertibility of potential vorticity.
The current position is that the comparison between
the predictions of the linear theory and real situations is
confusing (Pierrehumbert and Swanson 1995). The the-
oretically elegant wave model is so unfit for comparison
to real data that alternatives, such as wave packets, are
considered for replacing it in the classical semilinear
framework of instability (Lee and Held 1993). Even
when the generalized baroclinic instability theory of
Farrell (1999) is considered, and although it is able to
predict highly localized structures, namely, implicit
wave packets, a quantitative assessment of its predic-
tion against one real case yields very poor agreement,
including in its description of the cyclone's overall
shape (Descamps et al. 2007). Conceptually, then, the
"interacting anomalies" approach seems closer to ob-
served facts, although it has been realized in the mean-
time that there are much more diverse situations con-
ducive to cyclogenesis than the upper-level precursor
paradigm, such as effective low-level precursors inter-
acting with a frontlike confluent zone (Arbogast 2004)
or jet-flow entrances, maxima, and exits (Uccellini et al.
1984; Shapiro and Keyser 1990; Rivière and Joly 2006).
Yet, it has the potential flexibility to account for this
growing complexity, especially since it can be cast into
the widely accepted framework provided by potential
vorticity and its inversion beyond quasigeostrophy
(Hoskins et al. 1985; Davis 1992; Arbogast et al. 2008).
These existing bases enable one to work both on ide-
alized, theoretical problems, semianalytically or not,
and on real cases. Indeed, it is essential for theory vali-
dation that the same framework and mathematical
tools should be applied to idealized as well as actual
cases. Potential vorticity attribution provides a power-
ful way to link wind and temperature anomalies to a
source or a cause in the potential vorticity distribution.
Wind and temperature fields attributed to distinct po-
tential vorticity features of the flow enable explicit
computation of their interactions, including deforma-
tion effects or energy conversions, everything that is
needed to set up a quantitative understanding of the
various processes operating in a given flow configura-
tion—everything except a somewhat universal way to
identify the anomalies in any given potential vorticity
distribution (Takayabu 1991).

A complete alternative approach to the linear insta-

bility view of cyclogenesis is nearly available, but it
stumbles on the lack of a proper mathematical repre-
sentation of the various structures present in any given
distribution of potential vorticity, if not analytically
then at least numerically. The primary purpose of the
present work is to trigger a search for such a represen-
tation of the ubiquitous but hazy "potential vorticity
anomaly" given that waves do not seem to provide the
proper concept. This representation should enable
identification, extraction, inversion, and tracking in
both real and idealized cases. It must also achieve a
high level of "objectivity"; that is, it must have very few,
possibly none, arbitrary parameters.

Looking around at the evolution of the close fields of
physics, it is tempting to call on the concept of *coherent
structure* and to associate many of the anomalies of
interest to meteorology to such structures: they are lo-
calized features that keep their coherence in time. The
time scales of air particles traveling inside a coherent
structure are shorter than the typical time scales of the
evolution of such a structure. Hakim (2000) computed
a climatology of upper-level vorticity anomalies and
then classified them according to their nonlinearities. A
main result of this study is that high-amplitude distur-
bances are common and that they cannot be interpreted
as linear or even weakly nonlinear waves. Hakim con-
firms that linear theories cannot represent the coherent
structures.

Some methods of extraction of the coherent struc-
tures associated with extratropical cyclones exist al-
ready. Such techniques usually rely on image segmen-
tation (Hodges 1994), temporal tracking of vortices
(Ayrault and Joly 2000), or contour detection of
anomalies of mean sea level pressure (Wernli and
Schwierz 2006). Their purpose is often to perform cli-
matologies of cyclones, in relative vorticity or pressure.
The extraction of upper-level potential vorticity coherent
structures is a difficult task. For instance, Wernli
and Sprengler (2007) identify potential vorticity
streamers and cutoffs using a form criterion, but this
technique still relies on subjective parameters. These
methods give some elements of the coherent structures,
such as their center, size, and intensity, but they do not
provide a complete description of them, particularly the
possible secondary vortices around them. For these rea-
sons, they are not sufficient to describe all features of
synoptic-scale coherent structures.

Wavelet theory seems to provide a promising frame-
work for the extraction of coherent structures. The
wavelet transform has been used for modeling coherent
structures in turbulent flows for several years now
(Farge 1992). The underlying assumption is that turbu-
lent flows may be separated into two components: a
coherent part, which carries the major part of momentum and energy of the total flow, and an incoherent part that may be represented statistically (Farge et al. 1999). Nonlinear filtering of the two-dimensional orthogonal wavelet representation is an efficient and accurate way to separate these components.

Following these results, it is tempting to investigate how to depict isolated coherent structures from meteorological fields using wavelets. Indeed, a basic assumption of this paper is that the wavelet representation is suitable to define and extract these localized high-amplitude features of the flow because of the local representation of scales that it provides. The purpose of the present work is to go further than Farge et al. (1999) and to determine locally individual coherent structures out of the coherent flow using the wavelet representation.

The first part of the present article opens with a brief presentation and a justification of the wavelet tools that are used: then it presents the extraction of the coherent field from a meteorological case study and the proposed algorithm to extract an individual coherent structure. The second part shows the application of the extraction on the upper-level precursor of an intense mid-latitude storm and, using potential vorticity inversion, depicts the dynamics of the coherent structure that is extracted. A comparison with a structure extracted with another method is performed. The paper finishes with a discussion of the importance and possible applications of this algorithm for a better understanding of synoptic-scale storms and for numerical weather prediction.

## 2. An algorithm to extract isolated coherent structures

### a. Overview of the algorithm

Before exploring the algorithm in detail, it is important to bring up the assumptions and the framework on which it relies. The coherent structures that are sought are three-dimensional synoptic features of potential vorticity. A separation between the horizontal and vertical directions is assumed, consistent with the hydrostatic balance of the synoptic atmosphere. Thus, the vertical dimension of a synoptic-scale coherent structure is obtained by performing a sequence of two-dimensional horizontal extractions and by looking for colocalized features at several levels.

The main hypothesis that underlies this work is that a two-dimensional coherent structure may be represented by a collection of neighboring wavelets [as Yano et al. (2004b), for one-dimensional signals], where “neighbor” refers to spatial and spectral proximity. To meet the requirements of nonredundancy and possible

reconstruction, the wavelets that compose an individual structure must belong to the same orthogonal basis. If there is a choice between several bases, an optimal wavelet basis should be determined. The result of the method must be weakly sensitive to the position, shape, and size of the structure, which will be assessed a posteriori.

The algorithm applied to a complete set of wavelet coefficients of a two-dimensional field follows the scheme:

* thresholding the wavelet coefficients (also called nonlinear filtering) in order to remove the incoherent part of the field,

* simple detection of the possible coherent structures,

* choice of an optimal orthogonal basis for each structure, and

* local selection on this basis of the wavelets that constitute each structure.

After a brief discussion of the wavelet tools that are available, the successive steps of the method are presented in the continuation of this section.

### b. Wavelet discrete representation

A wavelet transform [see Mallat (1998) for theoretical developments] is defined by its mother wavelet $\psi$, which is a function with zero average. The dilatations and translations of $\psi$ in space generate a set of wavelets:

$$ \psi_{u,s}(x) = \frac{1}{\sqrt{s}} \left( \frac{x-u}{s} \right), u \in \mathbb{R}, s \in \mathbb{R}, \quad (1) $$

that defines a one-dimensional continuous wavelet transform. The wavelet representation of a function $f$ is given by its wavelet coefficients:

$$ Wf(u, s) = \langle f, \psi_{u,s} \rangle = \int_{-\infty}^{+\infty} f(x) \psi_{u,s}(x) dx. \quad (2) $$

If the dilatations are dyadic and if the set of translations is discrete {$s = 2^j, u = 2^j p, j \in \mathbb{Z}, p \in \mathbb{Z}$}, then Eq. (2) generates a discrete wavelet transform (DWT). Moreover, if some (restrictive) conditions for the mother wavelet are met, the corresponding set of discrete wavelets may be an orthogonal wavelet basis.

An orthogonal DWT basis on the sphere would be an ideal technique to represent synoptic-scale structures on the earth, but such a basis is not yet known—if it ever will be. The most interesting, but still unsatisfactory, theories in this domain are the spherical biorthogonal wavelet transform (Schröder and Sweldens 1995) and the semiorthogonal multiresolution technique (Rosca 2005), which are not suitable because of
<figure><img src="image_2.png" /><figcaption>FIG. 1. Visualization of the Coiflet with four vanishing moments for the (left) wavelet function $\psi$ and (right) scaling function $\phi$. Each function is represented in (top) physical space and (bottom) Fourier space.</figcaption></figure>

nonorthogonality. As a consequence, the orthogonal wavelet transform of the present study is the two-dimensional DWT that is applied to the fields after their projection on a plane grid. Such a transform relies on successive one-dimensional DWTs; it is also used for the separation of coherent fields in fluid mechanics (Farge et al. 1999).

Some of the most famous orthogonal wavelet transforms are Haar, Daubechies, Meyer, Symlet, and Coiflet. A main utility of the wavelet transform is to define local features in physical and spectral space. The spatial localization depends on the size of the support of the mother wavelet. The spectral localization of the wavelet transform is defined by its number of vanishing moments. A notable theorem states that the size of the support of the mother wavelet is proportional to its number of vanishing moments, so both properties cannot be optimal and their choice must result from a compromise.

The choice of a suitable wavelet for meteorological fields has been examined by some authors (see, e.g., Yano et al. 2004a). For the purpose of extracting localized physical structures from the field, there should be a compromise between a small support, a smooth and

symmetrical physical shape, and a sufficient number of vanishing moments. The need for a well-localized representation favors a compactly supported wavelet. The shape of the wavelet in physical space is important for the regularity of the structures (Yano et al. 2004a); Haar and Daubechies wavelets are thus ill adapted. The number of vanishing moments should be at least around four (K. Schneider 2005, personal communication). The Coiflet with four vanishing moments seems to be optimal regarding these constraints. The mother wavelet $\psi$ (similar to a high-pass filter) and the scaling function $\phi$ (low-pass filter) of this Coiflet are shown in Fig. 1. The corresponding two-dimensional wavelets in physical space are presented in Fig. 2. A two-dimensional wavelet is defined by

* its wavelet scale $j$, which means that its physical scale covers $\alpha 2^j$ grid points, where the coefficient $\alpha$ depends only on the wavelet and the scaling functions (for the Coiflet with four vanishing moments, $\alpha$ is quite small owing to its compact support);

* its direction, which depends on whether the function computed on the $x$ and $y$ axes is the wavelet or the scaling function [following Mallat (1998), they are
<figure><img src="image_3.png" /><figcaption>FIG. 2. Physical representation of some two-dimensional wavelets (Coiflet with four vanishing moments). The (left) vertical (V), (middle) diagonal (D), and (right) horizontal (H) wavelets are (top to bottom) from scale 1 (the smallest) to scale 3. The wavelet coefficients that generate these functions have the same value, so all wavelets here have the same norm.</figcaption></figure>

called, respectively, horizontal ($\psi_x\phi_y$), diagonal ($\psi_x\psi_y$), or vertical ($\phi_x\psi_y$)]—these wavelet directions refer to the two-dimensional plane and they should not be confused with the physical three-dimensional directions (i.e., the vertical wavelet direction does not correspond to the vertical atmospheric levels).

The usual space for discrete wavelet representation is an array in which the wavelet coefficients are ordered according to their scale and direction (Fig. 3). The two-dimensional domain must be obtained after a conformal projection for which the map factor does not vary too much. For all computations herein, a conformal Lambert projection (Table 1 and Fig. 4) will be used. The average resolution is 120 km and the map factor

does not vary much throughout the areas of interest; over the midlatitude Atlantic storm track, its deviation from 1 does not exceed 5%. This value reaches its maximum (30%) in the polar latitudes, which is not a zone of interest for the present study.

The wavelet scales are not rigorously equivalent to a two-dimensional elliptic truncation because the spectra of the wavelet and the scaling functions do not have an exact cutoff frequency (Fig. 1) and the x and y directions of the domain are supposed to be separable. However, since the number of vanishing moments of the wavelet has been controlled, the wavelet scale, which depends on the wavelet direction (Fig. 2), gives a correct indication of its two-dimensional scale. Where the map factor of the projection does not vary too
<figure><img src="image_4.png" /><figcaption>FIG. 3. Usual representation of the wavelet coefficients of a field defined on the 128 × 64 domain (Fig. 4) through the discrete wavelet transform (DWT) up to four scales. The remaining scale (5) is called the approximation. Every subrectangle contains the wavelet coefficients of the same scale and direction. The different points inside a subrectangle are the wavelet coefficients for different translations of the wavelet.</figcaption></figure>

much, the concept of wavelet scale is therefore meaningful as a measure on the sphere of the scale or of the typical size of a structure.

The DWT suffers from a lack of translation invariance, which is definitely a serious drawback (Fournier 2000). Figure 5 shows the main reason why translation invariance is critical for the extraction of coherent structures. By just modifying the translation of the DWT basis, a monopolar structure has a very different wavelet representation. For the first basis 13 wavelet coefficients are sufficient to reproduce the structure, and for the second basis 36 coefficients are needed. If the wavelet representation is not compact enough, then the structures may be incorrectly reconstructed and there may be some interference with any other neighboring structure. An important consequence follows: an arbitrary choice or a global optimal choice (such as presented by Fournier 2000) for the wavelet basis is unsatisfactory; one needs to adapt the translated wavelet basis to each structure in the flow.

<table><caption>TABLE 1. Parameters of the Lambert conformal projected domain.</caption><tr><td>Latitude where the cone is tangent to the sphere</td><td>45°N</td></tr><tr><td>Vertical longitude</td><td>30°W</td></tr><tr><td>Coordinates of the southwest corner of the domain</td><td>5°N, 85°W</td></tr><tr><td>Longitude of the northeast corner of the domain (its latitude is deduced from isotropy of the grid)</td><td>85°E</td></tr><tr><td>Number of grid points of the domain (x × y)</td><td>(128 × 64)</td></tr></table>

The solution that is proposed is to compute all of the DWT obtained by translation of the grid. The resulting transform is called the stationary wavelet transform (SWT) [also called undecimated wavelet transform, introduced by Coifman and Donoho (1995)]. It gives a redundant representation of the field. However, as the SWT is equivalent to a set of translated DWT, some wavelets from the SWT define orthogonal bases. More formally, the SWT consist in applying to Eq. (1) the set of dilatations and translations {$s = 2^j, u, j \in \mathbb{Z}, u \in \mathbb{Z}$}.

Since the DWT generates an orthogonal basis of the two-dimensional ($x \times y$) domain, there are ($x \times y$) DWT coefficients (Fig. 3), whatever the final scale of the wavelet transform. However, for a SWT with maximum scale $j$, there are [$x \times y \times (3j + 1)$] coefficients in all (Fig. 6).

The coherent structures will be sought in the whole set of SWT coefficients. First, a separation between incoherent and coherent parts is performed.

### c. Nonlinear filtering of the field

The aim of this section is to show a method to filter out the weakest wavelet coefficients and to de-noise the original field in order to detect the coherent structures more easily. Indeed, the nonlinear filtering should help to separate the coherent and the incoherent part of a meteorological field, in a similar way as Farge (1992) does in fluid mechanics. Wavelet thresholding techniques are commonly used to compress meteorological fields (Yano et al. 2004a).

The nonlinear filtering algorithm is presented in the appendix. It is expected to be suitable for the purpose of extracting the coherent field, which is composed of the coherent structures (Farge et al. 1999). The assumption of turbulence—that the incoherent field is isotropic Gaussian noise—must be verified a posteriori for the synoptic-scale meteorological fields.

An example of wavelet thresholding on DWT and on SWT coefficients is shown in Fig. 7 for the potential vorticity field on the isentropic surface 315 K. The original field has been extracted from the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) database (Uppala et al. 2005). The accuracy and the efficiency of the nonlinear filtering methods are tested by assessing the coherence of the incoherent part, which must be weak. The coherent part of the DWT (10.5% of the total number of coefficients) represents 99.7% of the square of the norm of the total field; the corresponding reconstructed norm for SWT is 99.5%. The skewness of the incoherent part is 0.18 for SWT and -0.03 for DWT. The kurtosis is 4.3 for SWT and 4.6 for DWT. These scores show that the incoherent part is not far from a
<figure><img src="image_5.png" /><figcaption>FIG. 4. Two-dimensional domain (size: 128 × 64 grid points) on which the wavelet transforms are computed. It results from a Lambert projection with the parameters given in Table 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Wavelet coefficients from the scales 2 to 5 (zoom of Fig. 3 because the scale-1 coefficients are negligible here) of the same monopolar structure for two translated DWT bases. For the (top) first basis, the structure has a more compact representation than that for the (bottom) second basis.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Wavelet coefficients of a field defined on the 128 × 64 domain (Fig. 4) through the stationary wavelet transform (SWT) up to four scales. The redundancy of the transform is illustrated by the number of coefficients, which is several times the size of the domain.</figcaption></figure>

Gaussian distribution (skewness 0 and kurtosis 3). Moreover, a local diagnosis on the structure marked by **T** in Fig. 7 has been performed. This structure is an upper-level precursor that will be described later. The reconstruction of its amplitude in the coherent part is 94.6% of the original amplitude for DWT and 96.2% for SWT. This intense local feature of the flow is well kept by both methods of nonlinear filtering, which confirms their ability to represent the local high-amplitude patterns of the field.

As expected, the coherent field seems to be smoother for the SWT reconstruction. The better performance of SWT at de-noising images is a known property (Coifman and Donoho 1995) and is due to its translation-invariant property. Although the SWT is not useful for compressing the field due to its redundant representation, the nonlinear filtering allows one to reduce the number of useful SWT coefficients. In the given example, around 55% of the SWT coefficients are kept after thresholding.

### d. Searching for the individual coherent structures

The reconstructed field after SWT nonlinear filtering highlights the cores of vorticity (e.g., **T** in Fig. 7) because it does not have the noisy oscillations that may be present in the original field. These cores of positive large amplitude are the coherent structures that are sought by the algorithm. Their centers are simply detected in the SWT-filtered field as the grid points where the field reaches a local maximum.

Once the center of a coherent structure is detected, the next step is to select the local wavelets that build the structure at best.

### e. Principles for optimizing the local wavelet representation

An idealized study for a one-dimensional signal should help to illustrate some general principles. SWT representations in wavelet space (u, j) are computed, which consists of plotting the values of each wavelet coefficient of position u and scale j in a (u, j) plane. A Gaussian monopole in wavelet space (Fig. 8) has a local maximum at every scale that is located at the same point as the peak of the monopole. Moreover, along this line of maxima, an absolute maximum in wavelet space is reached at scale 4. The question now is to choose, inside the SWT representation, the DWT basis that would best represent the monopole with the smallest number of coefficients. In Fig. 8, the position of the wavelets from two different bases are shown. One has a wavelet located at the peak of the structure at scale 4, and the other at scale 5. Keeping the wavelet whose coefficient has the maximum value is consistent with the purpose of optimizing the compactness of the wavelet representation. As a consequence, it is reasonable to choose a DWT basis that owns the maximal wavelet coefficient located at the peak of the structure. This coefficient is called the *main coefficient* of the structure.

In this DWT basis, the structure projects preferably on a limited band of wavelet scales between 1 and 4. Indeed, the coefficients at scale 5 for this basis are quite weak. Therefore, it is acceptable to state that there is a maximal wavelet scale for representing the structure, that is, the scale of the main coefficient. Such a property is essential to create a compact structure since involving too large scales would give birth to a global structure in space, which is one of the main drawbacks of the plane-wave representation.

### f. Determination of DWT basis to represent the structure

In two dimensions, the determination of the DWT basis needed for the structure is inspired by these general principles. The (u, j) plane of Fig. 8 becomes a parallelepiped (x, y, j), similar to the one presented in Fig. 6. Contrary to the one-dimensional case, there are three types of directional wavelets. In addition to the scale of the main coefficient, it is also necessary to establish its direction.

It is notable that the vertical and horizontal wavelets have an anisotropic shape, whereas the diagonal wavelets are nearly isotropic. To detect and extract structures independently from their rotation, it is preferable to assume the main coefficient to be diagonal. This guarantees a rotation-invariant detection and extrac-
<figure><img src="image_8.png" /><figcaption>FIG. 7. Nonlinear wavelet filtering of the potential vorticity field on the isentropic level 315 K at 0600 UTC 27 Dec 1999. The (top) original field is decomposed into (left) a coherent and (right) an incoherent part (negative contours are dashed) through the filtering of (middle) DWT coefficients or (bottom) the SWT coefficients. Contour intervals are 1 pvu (potential vorticity unit equivalent to $10^{-6} \text{ m}^2 \text{ s}^{-1} \text{ K kg}^{-1}$) for the original and coherent fields, 0.1 pvu for the incoherent fields.</figcaption></figure>

tion. Similarly to the (u, j) plane, maps of diagonal coefficients at every scale j (Fig. 9) are obtained by plotting the diagonal wavelet coefficients of scale j and position (x, y) in a plane. These maps are used to determine the scale of the main coefficient of every coherent structure. Directly inspired by the one-dimensional principles, two rules are followed:

* if a scale is involved in a coherent structure, then there is a local maximum in the diagonal coefficients at this scale at the point (± a little tolerance) of the peak of the coherent structure;

* the scale of the main coefficient is the highest one for which the value of the diagonal coefficient at the point deduced from the preceding rule is stronger than for the scales below.

Figure 9 illustrates these two constraints on the maps of diagonal coefficients at scales 3, 4, and 5 for a particular

coherent structure T whose peak is indicated by a cross, from the field shown in Fig. 7. For scales 3 and 4, this maximum is collocated with a local maximum in the SWT coefficients, which is no longer the case for scale 5: the first rule is only fulfilled at scales 3 and 4. Moreover, the value of the diagonal coefficient is higher for scale 4 than for scale 3. Therefore, the diagonal wavelet at scale 4 is the main wavelet for building the structure T. The reconstruction process will validate this choice afterward. Reconstruction at scale 3 (not shown) creates a too small structure, while at scale 5 (not shown either) it does not keep its compactness and connectedness.

g. *Building the coherent structure locally*

From the above section, a DWT basis and an upper-bound scale are deduced; what remains is to select the local wavelet coefficients that are in this basis and be-
<figure><img src="image_9.png" /><figcaption>FIG. 8. Projection of a one-dimensional Gaussian monopole (black bold function) in the SWT wavelet space ($u, j$), where $u$ is the position and $j$ the scale of each wavelet. The values of the wavelet coefficients are indicated in shadings (positive for black contours, negative for white contours). The positions of the wavelets for two DWT bases are shown: one has a wavelet at the peak of the structure at scale 4 (crosses) and the other at scale 5 (dots). Below scale 3, these two bases are the same. At scale 5 for the cross-marked basis, there are two exclusive symmetrical possibilities: one is black, the other is gray.</figcaption></figure>

low this scale in order to describe the coherent structure. At the same time this problem is solved, some choices that have been made will receive a supplementary justification. No assumption must be made about the shape of the structure and the possible secondary poles around the central monopole.

The example of a Gaussian monopole confirms that the choice of an optimal DWT basis is critical: in Fig. 5 the detected wavelet scale of the monopole is 4, and the optimal representation is obtained by the DWT basis that has a diagonal wavelet of scale 4 at the same grid point as the monopole. This illustrates again the good choice of a diagonal coefficient as the main coefficient for representing the structure. Besides, there are only four other nonnegligible coefficients at scale 4: two vertical coefficients and two horizontal coefficients that surround the central diagonal coefficient.

The generalization of this scheme (one diagonal + two horizontal + two vertical) for the description of any coherent structure is now justified. The position of the diagonal wavelets is such that they are at the center of the patterns of symmetry of the DWT basis (Fig. 10); therefore, taking a diagonal wavelet as the main coefficient leads to the most compact representation. It is closely flanked by two horizontal wavelets above and below and by two vertical wavelets at the left and at the right of it at the same scale. A horizontal or vertical main wavelet would not offer such simplicity and symmetry.

At lower scales, the reconstruction must involve all of the wavelets that are colocalized with these five large-scale wavelets. By analogy with what is known in

the literature as the cone of influence of a physical point (Mallat 1998), the wavelets that are colocalized with a large-scale wavelet may be called its cylinder of influence (Fig. 11). To build the coherent structure, it is sufficient to define the cylinder of influence from the main coefficient, between scale 1 and the scale of the main coefficient. For the Coiflet with four vanishing moments, the diameter for the cylinder should be $2^j$, where $j$ is the wavelet scale of the structure (Fig. 11).

Figure 12 illustrates the process of construction of structure **T** from Fig. 7 following these principles. The role of small-scale wavelet coefficients is to tighten or relax the local gradients in a more realistic way: the reconstructed local features (Fig. 12) look like the local features of **T** in Fig. 7. At the maximal scale, the four wavelet coefficients that surround the central diagonal coefficient allow some freedom for the existence of secondary poles (in Fig. 12, top left panel, the negative poles appear due to the “horizontal” or “x wavelets”). Then, the small-scale coefficients change the size and the shape of the three poles of the coherent structure.

### h. General properties of the wavelet extraction

An ideal method of extraction should be able to keep the main characteristics of the coherent structures. In other words, the result should be nearly insensitive to the method itself. To conclude with the presentation of the algorithm, a short review of the quality of the reconstruction of the main characteristics of a coherent structure is follows:

* its *position*—the SWT representation is designed to have the optimal translated orthogonal basis, so
<figure><img src="image_10.png" /><figcaption>FIG. 9. Maps of diagonal coefficients in the plane (x, y) at (top) scale *j* = 3, (middle) scale *j* = 4, and (bottom) scale *j* = 5. See text for more details. The solid (dashed) lines are positive (negative) coefficients. The contour interval at scale *j* is twice as large as the one at scale *j* − 1. The black crosses indicate the point where the coherent structure **T** in physical space reaches its maximum.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Position of the wavelets of a DWT basis at a given scale *j* (*j* = 3 here) on the grid. The peaks of the diagonal (d), horizontal (h), and vertical (v) wavelets are regularly spaced following the scheme presented on this figure. The diagonal wavelets play a central role: they are closely flanked by two horizontal and two vertical wavelets. Note that in the SWT representation, there is a wavelet peak for every direction and every scale at each gridpoint.</figcaption></figure>

the method is insensitive to the position of the struc-
ture;

* its *size* (or *scale*)—the wavelet scale is defined as the one that leads to the most compact wavelet representation of the structure (section 2f), and the dyadic scale invariance of the wavelet representation guarantees the low sensitivity of the extraction to the size of the structure;

* its *aspect ratio*—if the structure is not too far from isotropy (i.e., if its aspect ratio is not larger than 2), the algorithm extracts and reproduces the structure well;

* its *rotation*—some tests of sensitivity (not shown) reveal that the algorithm is able to reconstruct an idealized elliptic (with an aspect ratio around 2) rotated structure well;

* its *shape*—the small-scale wavelet coefficients (the cylinder of influence of the main coefficient) modulates the shape of the structure so that it resembles the original field locally;

* its *multipolarity*—the five-coefficients representation at the large wavelet scale enables the existence of secondary poles around the central monopole, but does not impose them.

The extraction can be considered objective since only
a few input parameters have been set. These are the
mother wavelet and the diameter of the cylinder of

<figure><img src="image_12.png" /><figcaption>FIG. 11. (left) Cone of influence (gray) of size *C* associated with the grid point *p* and (right) cylinder of influence (gray) linked to a wavelet of scale *j*; *j* = 3 in this illustration. The DWT coefficients in the cone and in the cylinder are indicated by the black crosses. The diameter of the cylinder is 2^1, which covers the negative secondary poles of the Coiflet at scale *j* (drawn below the cylinder).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12. Progressive reconstruction of the potential vorticity (PV) structure named T in Fig. 7 on the isentropic level 315 K (contour interval 1 pvu; positive PV is solid, negative PV is dashed). The DWT basis is the one for which the diagonal wavelet at scale 4 has its peak at the same point as the structure. The coherent structure is progressively built from the (top left) main diagonal coefficient, on which all the (top right) coefficients of scale 4 have been added, and then the (middle left) coefficients of scale 3 and the (middle right) complete structure. The bottom panel shows the complete structure on the entire domain in order to reveal its compactness.</figcaption></figure>

influence for reconstruction. However, the values of these parameters have been objectively justified and they should not be changed by users. It follows that this algorithm may be used for the extraction of all kinds of two-dimensional coherent structures, from any meteorological or other scientific field.

The only deficiency of the algorithm in its current stage may appear when the structure to be extracted is a filament with a high anisotropy. A different exploitation of the wavelet coefficients will have to be sought to correct this problem in future work. For instance, the adaptation of the spatial extension of the cylinder of influence may be a solution.

From now on, it will be assumed that the coherent structures of the two-dimensional field are the structures resulting from the algorithm. The collocation between isobaric levels allows the computation of three-dimensional structures. An important aspect of this algorithm is that the temporal evolution is not involved: as a result, temporal consistency or coherence offers a way to check the behavior of the algorithm. The objective definition of the coherent structures provides a new basis for interpreting their dynamics in the potential vorticity inversion and the attribution framework; an example is presented hereafter.

### 3. Application: Extraction of the upper-level structure of the T2 storm

The exceptional storm T2 of December 1999, also called Martin by some mass media, hit the coasts of Europe around 1800 UTC 27 December 1999.

The purpose of this section is to show that T2 involves coherent structures of potential vorticity and that the upper-level structure may be extracted using the wavelet algorithm. Such an extraction may help to show the interactions of this coherent structure with its environment during the development stage of the cyclone.

#### a. The life cycle of the upper-level precursor of T2

The dynamical fields that gave birth to T2 deserve attention before the relevant upper-level precursor is extracted. The study will concentrate on the day before the storm hit the coasts of Europe, which corresponds to the phase of cyclonic development. To keep a Lagrangian conservative view (under the adiabatic assumption) of the upper-level dynamics, the field of potential vorticity from ERA-40 on an isentropic surface (315 K) is shown in Fig. 13. These isentropic potential vorticity maps (Hoskins et al. 1985) allow one to diagnose the coherence of the upper-level structure associated with T2. In the same figure, the surface cyclone is indicated through the field of relative vorticity on the isobaric level 850 hPa. As the vorticity in the surface cyclone V strengthens with time, the upper-level structure T lags behind it in the main flow, which suggests an interaction between both structures. They may both be tracked with time, which is consistent with the definition of a coherent structure. These structures evolve in a strong baroclinic region, indicated by a narrow band of strong gradient of potential vorticity on the isentropic level (Schwierz et al. 2004), correlated with an upper-level wind speed maximum (shown hereafter).

The extractions will be done on the fields from 27 December. They will help to show the role of this upper-level coherent structure for the development of the
<figure><img src="image_14.png" /><figcaption>FIG. 13. Time evolution of the potential vorticity field on the isentropic level 315 K (solid lines, interval 1 pvu) and relative vorticity on the isobaric level 850 hPa (gray shadings, min level 0, interval $5 \times 10^{-5} \text{ s}^{-1}$) during the growing phase of T2 in December 1999. Source: ERA-40. The relative position of the upper-level structure **T** and the low-level vortex **V** suggest a significant interaction.</figcaption></figure>

cyclone and to diagnose the interactions with the surface cyclone in the baroclinic region.

b. *Extraction of the coherent structure*

The data were provided by the Action de Recherche Petite Echelle Grande Echelle (ARPEGE)/Integrated Forecast System (IFS) operational analysis, three-dimensional variational data assimilation (3DVAR) in 1999 (Courtier et al. 1991) with truncation T199 on a stretched spherical grid. The average resolution over the Atlantic Ocean is about 40 km. The fields are projected on the same plane as explained above. The switch of data source compared to the previous part is due to the need for sufficient vertical and horizontal resolutions and because of the need to use the ARPEGE framework, which allows the computation of potential vorticity inversion and numerical forecasts.

The coherent structure is sought in the Ertel potential vorticity field on the isobaric levels in steps of 25 hPa between 100 and 850 hPa. The algorithm detects a signature of the coherent structure between levels 200 and 650 hPa: on the higher vertical levels and on the levels below there is no signature of a local maximum of potential vorticity.

The coherent structure is extracted at each of these levels at 0000, 0600, and 1200 UTC 27 December

<figure><img src="image_15.png" /><figcaption>FIG. 14. Potential vorticity on the level of 350 hPa at (top) 0000, (middle) 0600, and (bottom) 1200 UTC 27 Dec, projected on the Lambert conformal projection used for the wavelet representation. The dashed lines 1 and 2 represent the vertical planes of the cross sections of Fig. 15, directed by the arrowheads. The double-headed arrows indicate the dilatation axes calculated for the environment of the anomaly (see text for more details). The length of the axis is proportional to the norm of the deformation vector. Only the axes near the anomaly are drawn. The field after SWT nonlinear filtering is shaded (intervals 1.5 pvu) and the structure that is extracted is represented by black lines (solid for positive, dashed for negative; intervals of 1.5 pvu).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 15. Potential vorticity cross sections at (top) 0000, (middle) 0600, and (bottom) 1200 UTC 27 Dec. Cross sections (left) 1 and (right) 2 of Fig. 14. Shown are: the original field after SWT nonlinear filtering (shaded; isolines from 1.5 to 3 pvu), the structure extracted (black solid lines; intervals 1.5 pvu) and the total field after removal of the structure (dashed lines; intervals from 1.5 to 3 pvu).</figcaption></figure>

(Figs. 14 and 15). The horizontal fields and vertical cross sections reveal a compact structure that evolves in time. The residual field is modified only locally and no spurious discontinuities are created horizontally. The extracted anomaly undergoes a time evolution that may be interpreted by dynamical diagnoses. It is obvious from Fig. 14 that the shape of the structure changes with time. At first elongated in a west-northwest direction, it evolves into a nearly isotropic shape before becoming elongated in a southwest direction. The dynamical consistency of this behavior may be assessed by considering the dilatation axes associated with the deformation field of the environment. Here, the environment must be understood as the residual field after removal of the extracted structure. Figure 14 reveals that the evolution of the shape of the anomaly is consistent with the dilatation axes. There is a persistent

orientation of the axes in a southwest direction, which is the one toward which the elongation of the structure gradually converges.

The cross sections (Fig. 15) reveal the ability of the algorithm of extraction to represent the tropopause trough associated with the upper-level precursor of T2. It moderately gains some amplitude and the vertical extension of the structure increases. Moreover, the coherent structure extends downward and southward between 0000 and 1200 UTC, typical of tropopause folding (Uccellini et al. 1985).

Another strong point of our approach is that it is sufficient to extract the coherent structure in the potential vorticity field. Its signature in all other dynamical fields is then derived after inversion in a dynamically consistent way—in the ARPEGE primitive-equation framework using an implicit balance condition
(Arbogast et al. 2008)—that prevents spurious gravity wave emission in the very short range forecast. For the inversion of the three-dimensional field of potential vorticity, the low-level boundary condition assumes that the potential temperature at the hybrid model level closest to 850 hPa is unchanged. This level belongs to the free atmosphere, where the balance condition is valid. It has been verified that the thermal structure and the vortex associated with the low-level preexisting cyclone are unchanged after inversion.

The attribution may be ambiguous in the context of a nonlinear balance (Davis 1992). To compute the dynamical fields of the structure, the field of potential vorticity after SWT filtering and the field without the coherent structure are both inverted. The attribution of the anomaly fields is obtained as the difference between the results of these two inversions.

### c. Description of the structure and of the large-scale flow

The upper-level (350 hPa) anomaly in terms of relative vorticity and wind (Fig. 16) indicates a main cyclonic structure that is collocated with the structure of potential vorticity. The wind strengthens with time as the cyclone develops. The wind flow attributed to the anomaly (barbs in Fig. 16) and the geopotential field (not shown) reveal two secondary anticyclonic poles around the central cyclonic monopole, particularly at 0000 and 0600 UTC. A pole is upstream of the main cyclone, and the other pole is downstream with regard to the baroclinic zone. This is consistent with the tripo lar configuration on the isentropic potential vorticity maps (Fig. 12) and on the cross sections of the structure (Fig. 15).

Another remarkable feature of the structure is its local contribution to baroclinicity. The wind attributed to the anomaly is maximum (around $30 \text{ m s}^{-1}$) in its southern part, where it has the same eastward direction as the jet stream. Therefore, it plays the role of a jet streak. The potential vorticity coherent disturbance and the dipolar structure of relative vorticity across the jet are consistent with the literature (Pyle et al. 2004; Cuningham and Keyser 2004). Figure 17 shows the effect of removing the structure in terms of upper-level geopotential and wind velocity. At the three times considered, removing the coherent structure is equivalent to deleting the small local trough. After removing the coherent structure, the flow is locally zonal and the maximum wind changes roughly from 90 to $60 \text{ m s}^{-1}$. Therefore, deleting the coherent structure makes the jet stream straighter and more uniform. It is then reasonable

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 16. The 350-hPa relative vorticity (black lines, solid: positive, dashed: negative; contour intervals $5 \times 10^{-5} \text{ s}^{-1}$), wind barbs (minimum speed $3 \text{ m s}^{-1}$), and wind speed (shading; minimum speed $10 \text{ m s}^{-1}$, intervals $10 \text{ m s}^{-1}$) attributed to the extracted and inverted coherent structure at (top) 0000, (middle) 0600, and (bottom) 1200 UTC 27 Dec.</figcaption></figure>

to consider it as the background flow in which the structure evolves.

### d. A comparison with temporal filtering and gridpoint extraction

A comparison is proposed between the wavelet extraction and a subjective method of extraction that may be found in the literature nowadays. As suggested by Davis and Emanuel (1991), a temporal filtering is applied to the field in order to separate the anomaly from the basic state. A Lanczos time filter is used (Duchon 1979), with an approximate cutoff period of 6 days. Following Chaigne and Arbogast (2000), the anomaly is determined in the high-frequency field as the set of grid points that have positive potential vorticity and are low
<figure><img src="image_20.png" /><figcaption>FIG. 17. Geopotential on the 350-hPa level [contour intervals 10 decameter geopotential (damgp)] and wind velocity at 350 hPa (minimum shading of 40 m s⁻¹, intervals of 20 m s⁻¹) at (top) 0000, (middle) 0600, and (bottom) 1200 UTC 27 Dec. The original fields (before extraction) are plotted on the left and the fields after removal of the coherent structure are on the right.</figcaption></figure>

cated inside a circle typical of the size of the structure. This method has been applied, 0600 UTC 27 December, on the same isobaric levels as for the wavelet extraction (200–650 hPa). The determination of the radius of the monopole is ambiguous since the structure on the high-frequency field is not well separated from the neighboring structures. Therefore, a mean radius, which appears in Fig. 18, has been chosen subjectively within reasonable geographical bounds.

The structure **T** in Fig. 13 is extracted by these two different techniques. The one extracted with the wavelet method will be called **W** (for wavelet), and the one obtained by the extraction of the monopole after temporal filtering will be called **F** (for temporal filter).

The structure **F** (Fig. 18) reproduces correctly the three-dimensional tropopause folding. The main difference with **W** is the absence of secondary poles around the central monopole, which may be observed on the vertical cross sections (Fig. 18).

The inversion and the result of attribution of dynamical fields reveal other differences (Fig. 19). The relative vorticity anomaly is roughly monopolar, but this is a

consequence of a built-in assumption. Though there is a weak anticyclonic vortex south of the main cyclone, it does not have strong wind acceleration in its southern part; it is not a jet streak. Furthermore, the change in the local curvature of the height contours after removing the structure is not as complete as for **W**.

To assess the performance of both extractions, some short-range simulations with the primitive-equation operational model ARPEGE/IFS (truncation T358, stretched grid C2.4, full physics) were run for 0600 UTC 27 December with different initial conditions. The diagnosis is simply the mean sea level pressure field at 1800 UTC 27 December (Fig. 20). Three forecasts are computed, starting from the initial states with the total inverted field, without **W** at 0600 UTC, or without **F**. Although the modifications of potential vorticity have been done at the middle and upper levels of the troposphere (above 650 hPa), the structures have a strong impact on the forecast of the surface cyclone. The removal of any initial structures **W** or **F** completely deletes the surface cyclone at final time.

A further diagnosis helps to better understand the
<figure><img src="image_21.png" /><figcaption>FIG. 18. Alternative extraction of the upper-level anomaly involved in the development of storm T2. This simpler technique uses temporal filtering and imposes a monopolar structure (see text for details); it gives an alternative anomaly called **F**. The fields are shown (top) on the 350-hPa isobaric level (Lambert conformal projection), and on both cross sections (bottom left) 1 and (bottom right) 2 indicated by the arrows. The background field (gray shadings) is the total potential vorticity field, and the solid black line is the extracted structure. Same legend as in Figs. 14 and 15. The circle in black-dashed contour demarcates the region inside which the structure is defined.</figcaption></figure>

dynamical role of both structures on the development of the storm. Some intermediate terms in the previous short-range forecasts are used to diagnose the fields attributed to each anomaly by means of differences between the forecasted fields. Figure 21 shows the vertical velocity attributed to the anomalies in the model evolution after both extractions at 0600 UTC. Vertical velocity is a common signature of the baroclinic interaction between synoptic-scale vortices. For both **W** and **F** extractions, the anomaly fields show the existence of a downward motion upstream and an upward motion downstream of the upper-level anomaly, consistent with vertical motion in baroclinic systems (Holton 1992) and with idealized cyclonic development (Schär and Wernli 1993). Moreover, the vertical motion attributed to the anomaly strengthens with time. The core of the evolved vertical velocity ascending zone related to **W** keeps a slightly simpler shape than the ascent linked to the evolving anomaly **F**.

Hence, both structures have similar impacts on the

dynamics of the storm. These results suggest that the upper-level monopolar cyclonic vortex, common to **W** and **F**, is the most important upper-level feature for the development of the surface cyclone (Fig. 20). However, further dynamical differences in **W** and **F** will be investigated in the following.

### e. Diagnosis of coherence

The purpose of this section is to diagnose the coherence of the extracted structures by looking at their model evolutions between instant 0 and *t*. Let $\mathbf{x}_0$ and $\mathbf{x}_t$ be the state vectors of the atmosphere at these two times. With the perfect model assumption, the relation $\mathbf{x}_t = \mathcal{M}(\mathbf{x}_0)$ is guaranteed, where the operator $\mathcal{M}$ refers to the nonlinear integration of the model. At time 0, an anomaly $\Delta\mathbf{x}_0$ is extracted, as well as $\Delta\mathbf{x}_t$ at time *t*, with a technique such as a wavelet or temporally filtered extraction. The residual from $\mathbf{x}$ after removing the anomaly $\Delta\mathbf{x}$ is called its environment ($\mathbf{x} - \Delta\mathbf{x}$).
<figure><img src="image_22.png" /><figcaption>FIG. 19. Result of the inversion for the extraction of the structure F, 0600 UTC 27 Dec. Legend as in (top) Fig. 16 and (bottom) Fig. 17.</figcaption></figure>

It is hypothesized that an anomaly is perfectly coherent regarding the model if

$$ \Delta \mathbf{x}_t = \mathcal{M}(\mathbf{x}_0) - \mathcal{M}(\mathbf{x}_0 - \Delta \mathbf{x}_0). \quad (3) $$

Equation (3) directly states that the extracted anomaly $\Delta \mathbf{x}_t$ at time $t$ is equivalent to the propagated anomaly from time 0 by the nonlinear model. It may also be seen as a generalization to nonlinear models of the classical linear equation of perturbations used in stability analysis. If the model is linear or linearized, the operator $\mathcal{M}$ becomes $\mathbf{M}$ and Eq. (3)

$$ \Delta \mathbf{x}_t = \mathbf{M} \Delta \mathbf{x}_0, \quad (4) $$

which is the classical resolvent equation for a linear problem. The classical mathematical definition of an anomaly comes from the solutions of Eq. (4), in particular, in the special cases when $\mathbf{M}$ is separable in space.

Another interpretation of Eq. (3) comes after rewriting it:

$$ \mathbf{x}_t - \Delta \mathbf{x}_t = \mathcal{M}(\mathbf{x}_0 - \Delta \mathbf{x}_0), \quad (5) $$

which means that the environment of the extracted anomaly and the environment of the propagated anomaly are the same. In other words, the environment

<figure><img src="image_23.png" /><figcaption>FIG. 20. Mean sea level pressure at 1800 UTC 27 Dec (contour intervals 5 hPa; bold references 1000 hPa) and 850-hPa wind speed (shading from 20 m s⁻¹; intervals 10 m s⁻¹) at the time the cyclone T2 strikes the European coasts. (top left) The ARPEGE 3DVAR analysis is presented. (bottom left) The ARPEGE 12-h simulation starting from the result of the inversion from the total distribution of initial potential vorticity is consistent with the analysis. The 12-h simulation after removal of the initial structure (top right: W, wavelet extracted; bottom right: F, monopolar extraction) both delete the surface cyclone.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIG. 21. Illustration of the vertical motion attributable to the upper-level anomalies using the ARPEGE simulations from 0600 UTC (left: initial state; right: 6-h simulation), for both methods of extraction: (top) **W** and (bottom) **F**. The large-scale flow is taken as the geopotential at 350 hPa (thin solid lines; intervals 20 damgp) after removing the structure. The anomaly fields are the relative vorticity at 350 hPa (bold lines, solid: positive and dashed: negative; intervals $5 \times 10^{-5} \text{ s}^{-1}$) and the vertical velocity at 600 hPa (shaded for upward velocity, dotted for subsidence, from $2 \times 10^{-1} \text{ Pa s}^{-1}$, intervals $4 \times 10^{-1} \text{ Pa s}^{-1}$).</figcaption></figure>

may be propagated by the model: it has its own existence, independently from the anomaly.

These interpretations of Eq. (3) justify its definition for the coherence of an anomaly in the context of a nonlinear model. Hence, Eq. (3) may be used to assess the coherence of the extracted structures at time 0 and *t* by comparing its right- and left-hand sides.

Some forecasts from 0600 UTC 27 December with ARPEGE/IFS are used to diagnose the coherence of the extracted upper-level anomalies **W** and **F** after a 6-h evolution. The propagated anomalies after extraction in the initial conditions [rhs of Eq. (3)] are compared with the extracted anomalies from the verifying analysis [lhs of Eq. (3)], at 1200 UTC 27 December. The results are presented in Fig. 22 for potential vorticity on the 350-hPa level. The extracted structure **W** from the verifying analysis at 1200 UTC reveals a similar pattern as its 6-h propagation from 0600 UTC. A similar comparison for **F** at 1200 UTC does not show such good agreement because of its too strong built-in assumptions: the extracted anomaly from the verifying analysis is monopolar and roughly circular, whereas the propagated anomaly has a multipolar complex shape. Therefore, Eq. (3) is nearly valid for the anomaly **W**, contrary to the anomaly **F**. More objective evidence of these state-

ments is provided by the correlation of potential vorticity between the anomaly fields that are extracted from the verifying analysis and the propagated anomalies. Correlations are computed over a limited domain that covers the anomalies ($35^\circ-65^\circ\text{N}, 40^\circ\text{W}-10^\circ\text{E}$), between levels 250 and 600 hPa (Fig. 23), for both techniques of extractions, thus creating four cases. The anomaly **F** generates low correlations at every level. The best scores for the anomaly are obtained when the wavelet technique is used for both extractions in the initial conditions and in the verifying analysis.

As a consequence, the time evolution of **W** and **F** through the model reveals that the wavelet extraction is more coherent with the dynamics than a monopolar extraction, with respect to Eq. (3).

Several reasons may explain the improved coherence of the wavelet extracted structure **W** over **F**. The wavelet extraction is able to produce multipolar configurations, which may be more stable than monopoles. Indeed, a monopole induces secondary circulations in the surroundings that may accelerate dispersion. Moreover, the suppression of a monopole generates a field that is not continuously derivable, contrary to the deletion of the structure extracted with wavelets. The resulting discontinuities in the dynamical fields may also enhance
<figure><img src="image_25.png" /></figure>

FIG. 22. Coherence of the extracted anomalies: (left) **W** and (right) **F**. Fields are the potential vorticity of the anomaly (bold lines, solid: positive and dotted: negative; interval 1 pvu) and the geopotential of the large-scale flow (thin solid lines, intervals 20 damgp). The anomalies extracted from the (top) verifying analysis are the direct results of the extraction from the analysis at 1200 UTC 27 December. The (bottom) propagated anomalies are obtained as the difference between the 6-h simulations (from 0600 UTC) initialized by the unmodified initial analysis and the initial analysis without the initial anomaly. The comparison between the top and bottom panels gives an indication of the coherence of the anomalies (see text).

dispersion in the simulation. For instance, the upper-level relative vorticity of the propagated anomaly (Fig. 21) remains quite compact in space for **W**, whereas for **F** it spreads along the jet. Another difference between the two methods of extraction that may explain their coherence is the method of computation of the large scale (wavelet scale or temporal filtering).

The coherence of an anomaly with the model is a central property if the purpose is to find the relevant degrees of freedom in a field. A structure that does not show such a coherence may stimulate irrelevant modes and modify the overall dynamics of the system.

## 4. Conclusions

The main results of this paper are summarized here:

* a new algorithm based on the stationary wavelet transform that is able to extract the coherent structures of a three-dimensional flow has been presented;

* this method of extraction is more objective than existing methods in that it has less built-in features and thresholds/parameters and it is weakly sensitive to the properties of a coherent structure, that is, position, size, and shape; moreover, the method allows one to extract multipolar vortices;

* the algorithm applied to the upper-level potential vorticity precursor of a midlatitude storm shows its ability to extract relevant structures of the dynamics; the evolution of the wavelet derived anomaly can be understood from the action from its environment, in particular from deformation;

* the comparison of this extraction with a method based on a monopolar extraction after a temporal filtering reveals the same impact on the forecast of the surface cyclone;

* however, the assessment of the coherence of each structure through the equations of the model reveals a far better coherence of the wavelet-extracted structure than the monopolar structure with the primitive-equations model; the wavelet-derived anomaly does seem to capture something akin to an actual coherent structure, a real, meaningful organization of the flow.

The capability to extract and properly handle synoptic-scale coherent structures of potential vorticity opens some new tracks to study synoptic-scale dynamics and its predictability. The present approach makes some
<figure><img src="image_26.png" /><figcaption>FIG. 23. Correlation of the potential vorticity between the extracted anomaly at 1200 UTC and the propagated anomaly at 1200 UTC (simulations from 0600 UTC). Each curve is named after two letters, the first being the technique used for the extraction in the verifying analysis (at 1200 UTC), and the second for the one used in the initial conditions (at 0600 UTC) to compute the propagated anomalies (at 1200 UTC). The correlation is the best when the wavelet method is used for both extractions.</figcaption></figure>

progress on the way toward a description of cyclogenesis that does not invoke linear stability arguments and on defining “anomalies.”

The dynamics and the interactions of coherent structures could be studied using such an objective method of extraction. All kinds of coherent structures are concerned, from large-scale jet streaks to mesoscale structures, but different meteorological fields and dynamical processes are also considered. In the context of baroclinic development, the role of vertical velocity has been emphasized, but the energetics of extratropical cyclones may also be described. Additionally, as the method reproduces correctly the shape of the coherent structures, the interactions determined by the deformation (Rivière and Joly 2006) may be diagnosed. Energetics and deformation should receive more attention in forthcoming work.

At least two important applications are foreseeable for predictability issues. First, the automatic extraction of coherent structures gives a new method to initialize an ensemble prediction system. Plu and Arbogast (2005) showed that, in an idealized framework, perturbing the initial coherent structures is a reliable method for obtaining the probability density function for the short-range forecast of cyclogenesis. That article concluded with the statement that there is no algorithm to extract the relevant structures that give birth to a storm; the present paper provides one. The present observation that the wavelet extraction provides anomalies that are coherent with the dynamics is an additional argument for using this algorithm.

Second, the assimilation of coherent structures from images could be achieved using this extraction. For instance, it might be expected that the upper-level precursors that can be detected in the water vapor channel from satellites would be assimilated with more accuracy using the wavelet representation rather than a gridpoint radiance representation.

*Acknowledgments.* The authors thank Pr. K. Schneider from the University of Aix-Marseilles for very stimulating discussions about the use of wavelets for extracting coherent structures in fluid mechanics. The useful comments from G. Rivière are also sincerely acknowledged. This paper has gained much in clarity after the suggestions of two anonymous reviewers and after a careful proofreading by J. Maziejewski. This work was supported by the French Ministry of Research through the Cyprim project (ACI Aléas et changements globaux).

# APPENDIX

## Nonlinear Wavelet Filtering

The nonlinear wavelet filtering is an optimal method to estimate a signal that is corrupted by some stochastic noise (Donoho 1993). It consists of separating the DWT coefficients above and below an optimal threshold $s$ that may be obtained through an iterative method (Azzalini et al. 2005):

$\tilde{f}_>$: wavelet signal to be split into $\tilde{f}_>$, the coherent part and $\tilde{f}_<$, the incoherent part.

$N = x \times y$: size of the domain (and also number of DWT coefficients); $\sigma[\cdot]$: standard deviation of the two-dimensional function in the argument over the computational domain.

Initialization:

$\tilde{f}_>= 0$, the coherent part (above threshold)

$\tilde{f}_<= \tilde{f}$, the incoherent part (below threshold)

$s = (2\sigma^2[f_<]\ln N)^{1/2}$ and $s_0 = 0$: initial thresholds.

While $s > s_0$, $s_0$ takes the value of $s$. Apply threshold $s$ on $|\tilde{f}_<|$:

* the wavelet coefficients that are above $s$ are added to $\tilde{f}_>$

* the wavelet coefficients that are below $s$ remain in $\tilde{f}_<$
  $s = (2\sigma^2[f_<]\ln N)^{1/2}$

**End**

As the threshold decreases at each step and remains positive, this algorithm always stops at a value for $s$ that is either positive or zero.
Another filtering method consists in thresholding the SWT coefficients with the threshold given by the DWT iterative algorithm. The thresholding of SWT coefficients has been proposed by Coifman and Donoho (1995) to avoid the spurious effects owing to the lack of translation invariance. After filtering all DWT fields, the SWT-reconstructed field is obtained by averaging the reconstructed DWT fields.

REFERENCES

Arbogast, P., 2004: Frontal wave development by interaction between a front and a cyclone: Application to the FASTEX IOP 17. Quart. J. Roy. Meteor. Soc., 130, 1675-1696.

—, K. Maynard, and F. Crépin, 2008: Ertel potential vorticity inversion using a digital-filter-initialization method. Quart. J. Roy. Meteor. Soc., 134, 1287-1296.

Ayrault, F., and A. Joly, 2000: Une nouvelle typologie des dépressions météorologiques: Classification des phases de maturation. Compt. Rend. Acad. Sci. Ser. II: Sci. Terre Planetes, 330, 167-172.

Azzalini, A., M. Farge, and K. Schneider, 2005: Nonlinear wavelet thresholding: A recursive method to determine the optimal denoising threshold. Appl. Comput. Harmonic Anal., 18, 177-185.

Chaigne, E., and P. Arbogast, 2000: Multiple potential vorticity inversions in two FASTEX cyclones. Quart. J. Roy. Meteor. Soc., 126, 1711-1734.

Charney, J. G., 1947: The dynamics of long waves in a baroclinic westerly current. J. Meteor., 4, 135-162.

Coifman, R. R., and D. Donoho, 1995: Translation-invariant denoising. Wavelets and Statistics, A. Antoniadis and G. Oppenheim, Eds., Springer-Verlag, 125-150.

Courtier, P., C. Freydier, J. Geleyn, F. Rabier, and M. Rochas, 1991: The ARPEGE project at Météo-France. ECMWF Seminar Proc., Vol. II, Reading, United Kingdom, 193-231.

Cunningham, P., and D. Keyser, 2004: Dynamics of jet-streaks in a stratified quasi-geostrophic atmosphere: Steady-state representations. Quart. J. Roy. Meteor. Soc., 130, 1579-1609.

Davies, H. C., 1997: Emergence of the mainstream cyclogenesis theory. Meteor. Z., 6, 261-274.

Davis, C. A., 1992: Piecewise potential vorticity inversion. J. Atmos. Sci., 49, 1397-1411.

—, and K. Emanuel, 1991: Potential vorticity diagnostics of cyclogenesis. Mon. Wea. Rev., 119, 1929-1953.

Descamps, L., D. Ricard, A. Joly, and P. Arbogast, 2007: Is a real cyclogenesis case explained by generalized linear baroclinic instability? J. Atmos. Sci., 64, 4287-4308.

Donoho, D. L., 1993: Unconditional bases are optimal bases for data compression and for statistical estimation. Appl. Comput. Harmonic Anal., 1, 100-115.

Duchon, C., 1979: Lanczos filtering in one and two dimensions. J. Appl. Meteor., 18, 1016-1022.

Eady, E., 1949: Long-waves and cyclone waves. Tellus, 1, 33-52.

Farge, M., 1992: Wavelet transforms and their applications to turbulence. Annu. Rev. Fluid Mech., 24, 395-457.

—, K. Schneider, and N. Kevlahan, 1999: Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthogonal wavelet basis. Phys. Fluids, 11, 2187-2201.

Farrell, B. F., 1999: Advances in cyclogenesis theory: Toward a generalized theory of baroclinic development. The Life Cycles of Extratropical Cyclones, M. A. Shapiro and S. Grønås, Eds., Amer. Meteor. Soc., 111-122.

Fournier, A., 2000: Introduction to orthonormal wavelet analysis with shift invariance: Application to observed atmospheric blocking spatial structures. J. Atmos. Sci., 57, 3856-3880.

Hakim, G. J., 2000: Climatology of coherent structures on the extratropical tropopause. Mon. Wea. Rev., 128, 385-406.

Hodges, K., 1994: A general method for tracking analysis and its application to meteorological data. Mon. Wea. Rev., 122, 2573-2586.

Holton, J., 1992: An Introduction to Dynamic Meteorology. 3rd ed. International Geophysics Series, Vol. 48, Academic Press, 511 pp.

Hoskins, B. J., M. E. M. Intyre, and R. W. Robertson, 1985: On the use and significance of isentropic potential vorticity maps. Quart. J. Roy. Meteor. Soc., 111, 877-946.

Kleinschmidt, E., 1951: Über Aufbau und Entstehung von Zyklonen, III Teil. Meteor. Rundsch., 4, 89-96.

Lee, S., and I. M. Held, 1993: Baroclinic wave packets in models and observations. J. Atmos. Sci., 50, 1413-1428.

Mallat, S., 1998: A Wavelet Tour of Signal Processing. Academic Press, 637 pp.

Petterssen, S., 2001: Weathering the Storm: Sverre Petterssen, the D-Day Forecast, and the Rise of Modern Meteorology. Amer. Meteor. Soc., 329 pp.

Pierrehumbert, R., and K. Swanson, 1995: Baroclinic instability. Annu. Rev. Fluid Mech., 27, 419-467.

Plu, M., and P. Arbogast, 2005: A cyclogenesis evolving into two distinct scenarios and its implications for short-term ensemble forecasting. Mon. Wea. Rev., 133, 2016-2029.

Pyle, M. E., D. Keyser, and L. F. Bosart, 2004: A diagnostic study of jet streaks: Kinematic signatures and relationship to coherent tropopause disturbances. Mon. Wea. Rev., 132, 297-319.

Rivière, G., and A. Joly, 2006: Role of the low-frequency deformation field on the explosive growth of extratropical cyclones at the jet exit. Part I: Barotropic critical region. J. Atmos. Sci., 63, 1965-1981.

Rosca, D., 2005: Locally supported rational spline wavelets on a sphere. Math. Comput., 74, 1803-1829.

Schär, C., and H. Wernli, 1993: Structure and evolution of an isolated semi-geostrophic cyclone. Quart. J. Roy. Meteor. Soc., 119, 57-90.

Schröder, P., and W. Sweldens, 1995: Spherical wavelets: Efficiently representing functions on the sphere. Computer Graphics Proceedings (SIGGRAPH 95), 161-172.

Schwierz, C., S. Dirren, and H. C. Davies, 2004: Forced waves on a zonally aligned jet stream. J. Atmos. Sci., 61, 73-87.

Shapiro, M. A., and D. Keyser, 1990: Fronts, jet streams and the tropopause. Extratropical Cyclones, The Erik Palmén Memorial Volume, C. W. Newton and E. O. Holopainen, Eds., Amer. Meteor. Soc., 167-191.

Sutcliffe, R. C., 1939: Cyclonic and anticyclonic development. Quart. J. Roy. Meteor. Soc., 65, 518-524.

—, 1947: A contribution to the problem of development. Quart. J. Roy. Meteor. Soc., 73, 370-383.

Takayabu, I., 1991: Coupling development: An efficient mechanism for the development of extratropical cyclones. J. Meteor. Soc. Japan, 69, 609-628.

Uccellini, L., P. Kocin, R. Petersen, C. Wash, and K. Brill, 1984: The President's Day cyclone of 18-19 February 1979: Synop-
tic overview and analysis of the subtropical jet streak influ-
encing the precyclogenetic period. *Mon. Wea. Rev.*, **112**, 31–
55.

—, D. Keyser, K. Brill, and C. Wash, 1985: The President's Day cyclone of 18–19 February 1979: Influence of upstream trough amplification and associated tropopause folding on rapid cyclogenesis. *Mon. Wea. Rev.*, **113**, 962–988.

Uppala, S. M., and Coauthors, 2005: The ERA-40 re-analysis.
*Quart. J. Roy. Meteor. Soc.*, **131**, 2961–3012.

Wernli, H., and C. Schwierz, 2006: Surface cyclones in the ERA-

40 dataset (1958–2001). Part I: Novel identification method and global climatology. *J. Atmos. Sci.*, **63**, 2486–2507.

—, and T. Sprengler, 2007: Identification and ERA-15 clima-
tology of potential vorticity streamers and cutoffs near the
extratropical tropopause. *J. Atmos. Sci.*, **64**, 1569–1586.

Yano, J. I., P. Bechtold, J.-L. Redelsperger, and F. Guichard, 2004a: Wavelet-compressed representation of deep moist convection. *Mon. Wea. Rev.*, **132**, 1472–1486.

—, R. Blender, C. Zhang, and K. Fraedrich, 2004b: 1/f noise and pulse-like events in the tropical atmospheric surface variabilities. *Quart. J. Roy. Meteor. Soc.*, **130**, 1697–1722.