Ivane Pairaud

Francis Auclair

Combined wavelet and principal component analysis (WEof) of a scale oriented model of coastal ocean gravity waves

Keywords: Wavelet analysis, Empirical orthogonal function, Internal gravity waves, External gravity waves, Coastal oceanography

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.

Coastal ocean numerical modeling is basically the representation of the dynamics of the coastal ocean in a chosen range of length scales and over an associated frequency band, including the modeling of both coherent processes and associated transient processes. The ocean dynamical features can be individually identified by combining wavelet analysis for time and frequency localization and principal component analysis to "decorrelate" physically consistent structures. In the present paper, the so-called WEof analysis is applied for the extraction of external gravity waves and internal gravity wave lower modes in a simple case of a flat bottom, constant Brunt-Väisälä ocean. It is shown that, with some well known restrictive assumptions, WEof analysis is an efficient candidate for the recognition of frequency localized dynamical processes.

Introduction

One factor common to all geophysical sciences is the spread of the physical processes over a large frequency band. We could insist at length upon the fact that any study of the 1 ocean or atmosphere must necessarily begin by a clear specification of the scales of the processes at play. Both spatial and time scales need to be chosen in a way that makes the problem tractable. The physics of capillarity waves propagating on top of the ocean has very little in common with the meandering of the Gulf Stream for instance. This is true for observational, theoretical or numerical studies.

The simplest answer to such difficulties is to focus on a single process at equilibrium. This leads to "process-oriented" studies, which assume that all other processes can be neglected and that the whole study can be based on some strong assumptions of linearity. Most academic studies in geophysical sciences use such an approach, mainly because no alternative is available. It is evident, however, that no physical process is totally independent of processes on the same scale or even of processes appearing on different scales. A linear combination of the various interactions is only valid as long as induced deviations are "small". Therefore, any linearity assumption is to be considered as a serious limitation to such studies, even if this is not always clearly stated. A second serious limitation of process-oriented studies is connected with the representation of the physical phenomenon independently of any transient process, although these transients may be responsible for the largest energy transfers.

As a consequence, such a restrictive focus diverges somehow from the fundamental idea it is based on: extracting and studying the physics at a given scale and at a given time. In coastal ocean studies, as in many geophysical fields, these limitations are clearly unacceptable. Indeed coastal oceanography is devoted to the study of processes spreading over length scales that range from a few meters to a few hundreds of kilometers with associated time scales that cover a few seconds (freak waves), some months (coastal trapped waves) or even longer (in the case of dynamical studies of shoreline and beach topographic morphology). In such a region of the ocean and over such scales, adjustment processes through waves can be as long as the duration of the forcing at the origin of the process. Several forcing sources can also be active over the same scales. This is the case, for instance, when the wind is competing with tidal forcing [START_REF] Van Der Molen | The influence of tides, wind and waves on the net sand transport in the North Sea[END_REF]. Although waves, and transient processes in general, are of fundamental importance in coastal regions, oceanographers cannot study them independently of their causes or their environment and must extract them from a complex non-linear background.

Moving from process-oriented modeling to scale-oriented modeling raises a large number of new problems including the recognition and reconstruction of coherent dynamical structures. Following [START_REF] Wang | The wavelet empirical orthogonal function and its application to the analysis of internal tides[END_REF] and [START_REF] Ouergli | Low frequency modes and their link with summer monsoon activity described by wavelet analysis throughout India[END_REF], we propose to combine Principal Component Analysis (PCA) [START_REF] Preisendorfer | Foundations of Principal Component Selection Rules[END_REF] and time-frequency wavelet analysis [START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF] to study the generation and propagation of fundamental ocean transient processes such as gravity waves. The so-called WEof technique is thus a combination of spatial PCA and time-frequency analysis. Nevertheless, the objective in this paper is not to automatically select the frequency band where energy can mostly be found [START_REF] Bakshi | Multiscale PCA with application to multivariate statistical processes monitoring[END_REF], but rather to show, as a first step, how WEof analysis can be integrated in a study of coastal ocean transient processes. The WEof method presented in this study is an extension of the analysis technique used by [START_REF] Wang | The wavelet empirical orthogonal function and its application to the analysis of internal tides[END_REF] where the computation of dominant modes from the vertical temperature profile gives the internal wave characteristics. [START_REF] Ouergli | Low frequency modes and their link with summer monsoon activity described by wavelet analysis throughout India[END_REF] further analyzes daily sea-level pressure maps to study low-frequency modes of the summer Indian monsoon. Our approach is 3D in space and includes multivariate analysis of both surface elevation and vertical velocities.

PCA, commonly known in geophysical studies as Empirical Orthogonal Functions (Eofs), allows a multivariate, multidimensional field to be decomposed. It provides an orthogonal basis whose vectors successively account for the maximum explained variance of the multivariate physical fields [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF][START_REF] Preisendorfer | Foundations of Principal Component Selection Rules[END_REF]. This tool is of great interest because of the physical meaning of each of the basis vectors [START_REF] Kutzbach | Empirical eigenvectors of sea-level pressure, surface temperature and precipitation complexes over North America[END_REF]. On the other hand, one of the main drawbacks of Eof based methods is directly related to their definition insofar as Eofs are eigenvectors of the covariance matrix of the fields. In consequence, they only provide information about second-order moments. If the statistics of the fields are far from being normal, the meaning of Eofs is not trivial [START_REF] Auclair | Space-time structure and dynamics of the forecast error in a coastal circulation model of the Gulf of Lions[END_REF].

As far as frequency localization is concerned and in the context of a multiscale approach, PCA provides no clear information. [START_REF] Auclair | Space-time structure and dynamics of the forecast error in a coastal circulation model of the Gulf of Lions[END_REF] show, for instance, that a mode corresponding to the free oscillation (a seiche) of a particular gulf can be recovered by PCA, and is thus well associated with the frequency band of such oscillations. However, for the remaining Eofs of the orthogonal basis, no further indication of frequency localization can be obtained.

A fundamental requirement of the necessary time-frequency analysis is the possibility to investigate transient, time localized processes. Neither the classical Fourier transform, which is well suited to stationary and linear processes, nor its windowed transform generalization can consequently satisfy such requirements [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. A wavelet transform, if applied with care, appears as an original, mathematically attractive answer to the frequency-time localization of processes.

WEof is used in the present study to analyze and separate internal and external gravity waves in a flat-bottomed ocean. Such waves are probably the most fundamental and basic transient processes which can be encountered in the ocean adjustment problem. They remain complex processes, however, as they involve both horizontal and vertical propagations, and can be generated at any location of the ocean basin and over a broad frequency band. The interest they have received in large-scale ocean modeling is, to a great extent, associated with the modeling of the energy transfer they are responsible for (Le [START_REF] Provost | Energetics of the M2 barotropic ocean tides: an estimate of bottom friction dissipation from a hydrodynamic model[END_REF][START_REF] Munk | Abyssal recipes II: energetics of tidal and wind mixing[END_REF][START_REF] Gustafsson | Computations of the energy flux to mixing processes via baroclinic wave drag on barotropic tides[END_REF]. Such general and fundamental characteristics motivated this study: any further application of the technique in a geophysical context is necessarily subordinated to a clear understanding of WEof performance in analyzing gravity waves.

The next section briefly introduces wavelet analysis and PCA. Section 3 describes the main properties of internal and external gravity waves together with a presentation of the modeling strategy. In that section, three experiments (Experiments A-C) are carefully analyzed to be used later as Reference simulations. Each of them is more particularly dedicated to the modeling of a few modes of a simple gravity wave propagating over a flatbottomed ocean with constant Brünt-Vaïsälä frequency. In Section 4, the ability of WEof to analyze several configurations of gravity waves, originating from at least two sources, is tested. Three different experiments (Experiments D-F) are designed to study more complicated cases where gravity waves are simultaneously generated over different frequency bands, different time periods or by sources oscillating at different locations. The results are discussed in Section 5.

WEof analysis

The analysis state vector chosen for the study of gravity waves is multivariate and 3D in space. It includes the surface elevation anomaly together with the vertical velocity at each location considered in the analysis. The surface elevation and vertical velocity are not necessarily analyzed at each model grid point but rather over an appropriately chosen "analysis" grid.

Wavelet decomposition

A 1D Continuous Wavelet Transform (CWT) is performed over the time series of each state vector variable [START_REF] Lee | Wave profile measurement by wavelet transform[END_REF]. The mother wavelet is stretched or tightened depending on the scale "dilation" factor (s) and is translated in time by a translation factor (u t ). The continuous wavelet transform is thus a convolution of the wavelet atoms ψ with the raw input signal f:

Wf (u t ,s) = f, Ψ u t ,s = +∞ -∞ f (t) 1 √ s Ψ * t -u t s dt (1)
where the superscript (*) stands for the complex conjugate, this convolution of f with ψ being computed in the Fourier space [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Based on this wavelet transform, the time-frequency energy density or power spectrum can be defined by:

P Wf (u t ,s) =|Wf (u t ,s)| 2 .
(2)

In the present study, the mother wavelet chosen for the wavelet transform is the complex Morlet wavelet [START_REF] Meyers | Spatial and temporal 26-day SST variations in the equatorial Indian Ocean using wavelet analysis[END_REF][START_REF] Torrence | A practical guide to wavelet analysis[END_REF]:

Ψ (t) = 1 π 1/4 e i6t e -t 2 2 (3) with +∞ -∞ |Ψ (t)| 2 dt = 1.
The chosen wavelet centre frequency is F c = 0.94 Hz. The Morlet wavelet is basically a sine wave modulated by a Gaussian function. It has been used quite often in the analysis of geophysical processes such as atmospheric and ocean waves (Foufoula-Georgiou and Kumar, 1994). The "period" (or "pseudo-period") T of the wavelet can be defined based on the scale "s" by:

T = s F c . ( 4
)
The most important limitations of the analysis concern the definition of the Cone of Influence (COI) and the Heisenberg boxes. Wavelet analysis shows large errors at the beginning and at the end of the finite length time series. The COI is associated with edge effects and is defined so that any process located inside the COI is correctly analyzed. Following [START_REF] Torrence | A practical guide to wavelet analysis[END_REF], an e-folding time of √ 2s is chosen. Everything outside the COI is consequently considered as dubious and only the wavelet coefficients inside the COI will be taken into account hereafter. Another limitation comes from the fact that the resolution of the wavelet transform is limited in both space and time. Unlike in Fourier analysis, the frequency resolution depends on the scale "s". The time-frequency window associated with the wavelet is defined by its centre or mean value (t 0 , ω 0 ):

         t 0 = +∞ -∞ t |Ψ (t)| 2 ||Ψ || 2 dt ω 0 = +∞ -∞ ω |Ψ (ω)| 2 ||Ψ || 2 dω (5)
where

||Ψ || 2 = +∞ -∞ |Ψ (t)
| 2 dt, and by its spread over t in time and ω in frequency space:

               t = +∞ -∞ (t -t 0 ) 2 |Ψ (t)| 2 ||Ψ || 2 dt 1 2 ω = +∞ -∞ (ω -ω 0 ) 2 |Ψ (t)| 2 ||Ψ || 2 dω 1 2 (6)
The time-frequency window is eventually given by:

[u t + st 0 -s t, u t + st 0 + s t] × ω 0 s - ω s , ω 0 s + ω s . (7)
Wavelet coefficients are consequently computed over a particular frequency band and not "at a particular frequency". Moreover, two components with close frequencies will be correctly separated only if the two associated boxes do not overlap. Finally, we note that boxes are centered in frequency space but not in the s-scale space. The time-frequency windows are often called Heisenberg boxes. Inside the box, information is spread and the area of the box satisfies:

ω t ≥ 1 2 . ( 8
)
In consequence, any decrease of the frequency resolution results in an increase of the time resolution.

Complex Eof decomposition

The wavelet analysis being performed only over time, it gives no indication of space correlations between several time series or even between the state vector variables. Complex PCA analysis is performed in wavelet space, i.e. over complex wavelet coefficients leading to the computation of the dominant singular vectors of the wavelet coefficient matrix.

PCA is performed through "Incomplete" Singular Value Decomposition (SVD), i.e. by considering only the k dominant singular values λ i (i =1, 2, ..., k). In practice, SVD is performed over the normalized complex wavelet coefficients whose temporal mean has been removed [START_REF] Bjornsson | A Manual for EOF and SVD Analyses of Climate Data[END_REF] and the normalization factor is based on the variance of each variable of the state vector. The variance explained by the i-th mode is defined by λ i / k j=1 λ j . WEofs are finally obtained by transforming back the left-singular vectors into time space.

Reconstruction

An admissibility condition is required for the analyzing function to allow reconstitution of the original signal from its wavelet transform [START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]:

C ψ = +∞ 0 | Ψ (ω)| 2 ω dω<+∞ ( 9
)
where Ψ is the Fourier transform of Ψ which must then have zero mean and be localized in time and frequency space. From the admissibility condition, the wavelet transform is complete and energy is conserved. This admissibility condition is satisfied by the Morlet wavelet leading to:

f (t) = 2 C Ψ real ∞ 0 ∞ -∞ Wf (u t ,s)Ψ (t -u t )du t ds s 2 . ( 10
)
The reconstructed physical field at the "s" scale is thus computed numerically from:

f (t) = 2 C Ψ real s+ s 2 s-s 2 ∞ -∞ Wf (u t ,s) Ψ (t -u t )du t ds s 2 (11)
where s is the size of the Heisenberg box in s-scale space. A similar reconstruction is performed to transform the WEof back into time space. In this case, Wf corresponds to the WEof and not to the complete field.

Gravity waves

Gravity waves are probably the most fundamental transient processes that appear in ocean adjustments. As far as coastal ocean scales are concerned, they are the first to come into play after a perturbation of the potential energy of the ocean. Such waves have been widely observed, and modeled both theoretically and numerically. In the present study, the gravity wave propagation over a linearly stratified, non-rotating, flat-bottomed ocean is investigated both analytically and numerically. In the present section, three Reference experiments (Experiments A-C) are carefully carried out, each of them being dedicated to modeling specific modes of this gravity wave. One objective of the preliminary investigation of these Reference simulations is to compare numerical modeling and analytical solutions. The reason for the choice of three Reference modeling experiments instead of only one single Reference is due to the high cost of modeling processes occurring simultaneously over differing characteristic scales. In relation with the scale-oriented modeling presented in Model configuration is shown together with the analyzed area (darken box) and the location of the sources S 1 and S 2 .

Section 1, Experiments A-C can be seen as three versions of the same Reference simulation, each being associated with a particular space-time modeling window. Table 1 summarizes the main properties of each experiment.

Gravity wave modal decomposition

A local disturbance of the pressure field can be propagated through the ocean without any mass transport by gravity waves (also known as buoyancy waves). Indeed, buoyancy acts as a restoring force for fluid parcels displaced from hydrostatic equilibrium. External gravity waves propagate horizontally at the interface between the ocean and the atmosphere, while internal waves can transport energy obliquely through the ocean density stratification. In the simple case of the propagation of gravity waves over a flat-bottomed, non-rotating, constant Brunt-Väisälä frequency ocean, analytical solutions can easily be found [START_REF] Gill | Atmosphere-Ocean dynamics[END_REF]. The "long" gravity waves studied in the present paper have a rather small vertical acceleration and can consistently be studied under the "quasi-hydrostatic approximation". With such assumptions, their horizontal wave number is small with respect to the vertical one and the wave propagation can further be expressed as a sum of normal modes. It can also be shown that there is only one barotropic mode but there is an infinite discrete set of baroclinic modes. Mode 0 is the external mode with a horizontal phase velocity:

c 0 = c h,0 = gH (12)
where g =9.81ms -2 is the acceleration due to gravity and H is the ocean depth. The remaining internal modes (Mode 1 and higher) propagate both vertically and horizontally and their wave number is thus given by k

2 = k 2 v + k 2 h
where k v and k h are the vertical and horizontal wave numbers, respectively. Each mode satisfies the dispersion relation:

ω 2 n = k 2 h N 2 H 2 n 2 π 2 (13)
where n ∈ N * is the mode number and k v = nπ H . Waves generated at a specific frequency can only propagate at a well-determined angle, and the maximum allowed frequency is the Brunt-Väisälä frequency (N): 0 ≤ ω 2 ≤ N 2 . The wave crests or "lines of constant phase" sweep at right angles to this direction: the phase velocity c n is perpendicular to the group velocity. The speed at which a surface of constant phase moves [START_REF] Bretherton | The general linearised theory of wave propagation[END_REF] horizontally is given by the phase speed. The horizontal phase speed for Mode n is:

c h,n = ω n k h = NH nπ (14)
and the vertical is:

c v,n = ω n k v = NH nπ 1 | tan θ| (15) with k h = k v | tan θ| ,λ h,v = 2π k h,v
the horizontal or vertical wavelength and θ is the angle of energy propagation with respect to the vertical. At a given frequency, the higher the mode, the shorter the wavelength and the smaller the phase speed.

These internal waves have been observed experimentally to lead to the well-known St. Andrew's Cross (Fig. 1) [START_REF] Leblond | Waves in the Ocean[END_REF]. The magnitude of the group velocity is:

c g = dω dk = d dk k h N k = N sinθ k 2 h + k 2 v . (16)
In the case of homogeneous Brunt-Väisälä frequency, energy propagates along rays in the direction of the group velocity c g making an angle θ with the vertical, given by:

cos θ = ω N . ( 17
)
Under the long-wave assumption, θ is close to π/2. Each mode has a particular group velocity so that the wave front of each mode will move at a different speed from the generation area. Within a short distance (so also after a short time), the wave fronts of the different modes are dispersed. In consequence, in this study we will observe different modal wave compositions depending on the model resolution.

Description of the numerical experiments

Before investigating complex configurations of gravity waves, we carefully consider, in the present section, the modeling of a full set of modes associated with a single wave, insisting in particular on the strategy in terms of scale modeling. Several Reference experiments (Experiments A-C) are presented to study the transient response to a given forcing of the ocean in terms of gravity waves.

The numerical studies are carried out using the time-splitting primitive equation coastal ocean model SYMPHONIE described in detail in [START_REF] Estournel | The plume of the Rhone: numerical simulation and remote sensing[END_REF] and [START_REF] Estournel | Observation and modeling of the winter coastal oceanic circulation in the Gulf of Lion under wind conditions influenced by the continental orography (FETCH experiment)[END_REF]. This is a 3D, free surface model based on a generalized sigma coordinate. Primitive equations are integrated over an Arakawa C-grid [START_REF] Arakawa | Vertical differencing of the primitive equations in sigma coordinates[END_REF] under incompressibility and hydrostatic assumptions using a Boussinesq approximation. Following [START_REF] Gaspar | A simple eddy-kinetic-energy model for simulations of the oceanic vertical mixing: tests at station Papa and long-term upper ocean study site[END_REF], processes smaller than the grid scale are modeled using a 1.5 level turbulent closure scheme.

A characteristic radiation condition is applied for the external mode at the open lateral boundaries [START_REF] Oey | A model simulation of circulation in the northeast Atlantic shelves and seas[END_REF]. For the internal mode, a radiative boundary condition for a simple wave equation is used:

du dt =-C d(u -u b ) dx . ( 18
)
According to [START_REF] Palma | On the implementation of passive open boundary conditions for a general circulation model: the barotropic mode[END_REF] and [START_REF] Raymond | A radiation boundary condition for multi-dimensional flows[END_REF], C should depend both on the internal wave phase speed and on its direction. In the present case, C is chosen constant and set to 1 m s -1 for the coastal ocean, where u and x refer to the direction perpendicular to the open boundary and the subscript b refers to the forcing value at the boundary.

The same ocean domain and bathymetry are chosen for all the experiments shown in this paper: in particular, depth is uniform and equal to 250 m and no coast is present. The Brunt-Vaïsälä frequency of the ocean background stratification is also the same and leads to an exponential density profile (z ≥ 0):

ρ(z) = ρ 0 exp N 2 z g
which can be locally assumed to be quasi-linear.

Gravity waves are generated in the Reference experiments by sources oscillating at a period T 1 = 2.65 h corresponding to an angle of propagation θ =87 • , close to the horizontal. This period is chosen so that Mode 1 has a horizontal phase speed of 1 m s -1 . The remaining characteristics of the theoretical modes are given in Table 2. Two kinds of forcing sources are investigated. The first is a perturbation of the surface elevation and the second represents an inner perturbation of the density field through the generation of oscillating vertical velocities. The first one corresponds, for instance, to the generation of gravity waves by Ekman pumping in a stratified ocean [START_REF] Thorpe | The excitation, dissipation, and interaction of internal waves in the deep Ocean[END_REF], while the second one induces an oscillation of the thermocline due to tidal forcing [START_REF] Baines | On internal tide generation models[END_REF][START_REF] New | Internal tidal mixing in the Bay of Biscay[END_REF][START_REF] Jezequel | Interaction of a semidiurnal tide with a continental slope in a continuously stratified ocean[END_REF]. Both forcing sources are described in Appendix A and are located at the center of the domain.

The computation of the wavelength is a rather fastidious matter in a model output field. We recall that wavelet analysis is used in time, i.e. for time series at a given point, the spatial correlations being obtained by PCA over the wavelet coefficients. It is thus interesting to have an idea of the type of gravity wave modes that are propagating at a given time and in a given area of the domain. To reach this objective, two approaches are used in the present study. The wavelength is first objectively measured directly from the model output (as is the case for the Reference experiments) or from the WEof output. To confirm this measured value, the wavelength is computed locally based on the gradients of wave phase and statistics of the computed values are shown for the chosen area. A complete description of the method is given in Appendix B.

Scale-oriented modeling

As external and internal waves exhibit very different time and space scales and can be generated in distinct frequency bands, they can hardly be studied simultaneously on the basis of a single configuration. This would indeed be computationally expensive and methodologically inconsistent. The model grid and domain size are thus carefully adjusted to the scales of the processes of interest. Three configurations have been designed and include 41-161 points analyzed in both horizontal directions, while 71 levels are regularly distributed over the uniform depth so that a given level is at the same depth at all points (Table 3). Waves are forced at a period of T = 2.6 h (and 6.2 h in Experiment D) in a 250-m deep ocean. These waves propagate and make an angle θ of 87 • (and 88.7 • in Experiment D) with the vertical plane. The total number of grid points and the horizontal and vertical grid resolutions are given before and after the analysis, together with the time step and the imposed horizontal viscosity coefficient.

The time needed for the different modes to encompass the entire domain is obviously longer for higher modes and varies linearly with mode number. Moreover, only the lower modes can accurately be studied in the numerical study and the wave energy associated with higher modes is dissipated by horizontal diffusivity and vertical turbulence. Aliasing problems that are associated with frequency folding also occur, resulting in a false increase in lower mode energy. Smaller scale processes associated with higher wave-modes are thus mentioned with care in the following.

In this way, the spatial scales of the propagating processes change with time. In consequence, a given mode is observed by considering the physical field in the right space-time window. Three different grid and parameterization configurations are consequently proposed to focus more particularly on Mode 0 (Experiment A), Modes 1 and 2 (Experiment B) and Mode 5 and higher (Experiment C). Each experiment is dedicated to the study of a particular transient response.

In Experiment A, the horizontal grid resolution is 50 km in both x and y directions over 41 points and the analysis is performed over one fourth of the domain, corresponding to 1485 km (22 points over the diagonal of the hatched area on Table 1, each of the points being separated from the next one by 50 √ 2 = 70.7 km). Experiments B and C are based on the same model configuration but with a horizontal grid resolution of 300 m and a domain diagonal length of 68 km. Only the time scale and the space resolution of the analysis differ: data are analyzed every 1200 m in Experiment B (over a diagonal length of 34 km) and every grid point, i.e. every 300 m in Experiment C (over a diagonal length of 8.5 km). In the first experiment, the wave has to cross a maximum distance of 1485 km to reach the open boundary on the diagonal direction. As theoretical Mode 0 group velocity is 49.5 m s -1 , 8.3 h are needed for this mode to reach the boundary. As a comparison, it would take more than 17 days for Mode 1 and 34 days for Mode 2. However, none of these modes can be modeled as their horizontal wavelength is five times shorter than the horizontal grid scale.

In Experiment B, several modes can be studied: Mode 0 arrives at the boundary after 0.2 h whereas it takes 9.4 h for Mode 1 and 18.9 h for Mode 2. Furthermore, Mode 1 (λ h = 9.55 km) and Mode 2 (λ h = 4.77 km) can both be represented. However, as the length of the simulation is reduced to 13 h, Experiment B is dedicated more specifically to the study of Mode 1.

In the third and last experiment (Experiment C), the resolution of the extracted data points for analysis is smaller (300 m) which allows a smaller domain than in Experiment B. However, because of the sub-grid process problems mentioned before, modes higher than Mode 6 cannot be well separated. The total duration of the experiment is 19 h, less than the time needed by Mode 10 to reach the open boundary (23.6 h) but large enough to allow Mode 5 to reach this boundary (11.8 h). Consequently, after 19 h, the higher modes are only visible near the generation area.

Reference experiments

Fig. 2 shows a diagonal cross section of the vertical velocity (Fig. 2a) and the surface elevation (Fig. 2b) for Reference Experiment A after 23 h. The propagation being symmetrical, only one-quarter of the modeling domain is shown. The vertical velocity cross section clearly displays a wave propagating nearly horizontally. Fig. 2b exhibits an isotropic horizontal propagation. A direct measurement gives λ h = 465.2 ± 35 km for the wavelength and c h = 48.7 ± 3.6ms -1 for its phase speed. The statistical approach confirms this measurement with 80% of the computed length scale in the area in the range 473 ± 35 km. It can further be noted that all the computed wavelengths are between 370 and 550 km (the horizontal model resolution is 70 km over the diagonal). These results are in good agreement with the theoretical results for Mode 0 given in Table 2, The treatment of Experiment B is complicated by the simultaneous presence of several modes with the same order of magnitude (Fig. 3). The distance over which Mode 1 propagates is twice the distance reached by Mode 2 but they are both present in the plot of the vertical velocity. Mode 1 vertical wavelength is larger than the total depth of 250 m and as a consequence only the horizontal wavelength can be measured on Fig. 3: λ h = 8.66 ± 1.2 km with a phase speed of c h = 0.91 ± 0.12 m s -1 . After 8 h, 65% of the points located over the diagonal between 48.8 and 62.4 km from the origin and between 39 and 217 m (i.e. the area where the direct measurement is made) have a wavelength close to the Mode 1 theoretical wavelength (λ h = 9.55 ± 1.2 km). This confirms the presence of Mode 1 wave propagation in this area. The mismatch can be explained by the simultaneous propagation of Mode 0 in this area and the resulting interactions, but also by the presence of a Mode 1 wave front in the same area. Some evidence for the propagation of Mode 2 can also be found in the area at the time.

Fig. 4 shows the same diagonal section of the vertical velocity after 19 h (Experiment C). Unlike in Experiment B, the vertical velocity is extracted at every grid point. Higher modes with non-negligible vertical components can clearly be seen, knowing that at this time, Mode 10 has crossed nearly all the distance to the open boundary. In consequence, even though the energy is especially contained in the central ray, this figure exhibits more than one mode. The main ray in dark grey makes a propagation angle of θ = 86.8 ± 0.6 • , which agrees with the expected 87 • for an internal wave propagating under such assumptions. Estimations of the wavelength are made at points M 1 and M 2 (Fig. 4). A direct measurement from the model field leads to λ h = 2.04 ± 0.21 km, λ v = 100.4 ± 8 m and c h = 0.21 ± 2.2 × 10 -2 ms -1 near the central ray (case M 1 ), which agrees well with the theoretical values for Mode 5 (λ h = 1.9 km, λ v = 100 m and c h,5 =0.2ms -1 ). At the second point (M 2 ), the same approach leads to λ h = 948 ± 210 m, λ v = 50.2 ± 4 m and c h = 0.10 ± 2.2 × 10 -2 ms -1 which agree with the theoretical values for Mode 10 (λ h = 955 m, λ v =50m and c h,10 =0.1ms -1 ). The statistical confirmation is less obvious than for Experiments A and B. On the diagonal section between 34.8 and 38.2 km from the origin and from 39 to 204 m deep (which corresponds to the neighborhood of M 1 ), 26% of the domain contains wavelengths consistent with the theoretical wavelength for Mode 5 (i.e. satisfies λ h = 1.91 ± 0.21 km). In the neighborhood of M 2 (i.e. in the diagonal section between 34.8 and 36 km, from 140 to 240 m deep) only 23% of the points correspond to the propagation of Mode 10. As in Experiment B, this discrepancy is related to the simultaneous presence of the other modes and their interferences. In Experiment C, modes with high wave numbers and small wavelengths propagate inside the domain as time evolves, the downscaling process being limited by the horizontal diffusion and the turbulence scheme. For instance, Mode 4 with λ h = 2.39 ± 0.21 km represents 24% of points in this area and the propagations are spread over the shorter wavelengths. As far as higher modes are concerned, Mode 10 can be clearly observed only near the source (Fig. 4).

Discussion

Numerical transient processes associated with the perturbation of the mass field are in agreement with gravity wave modal decomposition. Both their modeling and their identification are, however, complicated matters, as both time and space scales change quickly in the ocean dynamical field. Classical numerical modeling approaches are not suited to such a representation if neither their time step nor their grid scale can evolve. As a consequence, several models of the same process must be made to study the complete modal decomposition of the gravity waves without drastically increasing computational cost.

The three Reference numerical Experiments A-C validate the ability of the coastal model to excite and propagate external and internal gravity waves and provide us with a numerical 3D modal decomposition of a gravity wave generated by a simple surface oscillating source. This Reference modal decomposition will be used in the next section as a comparison for the WEof decomposition for more complicated configurations.

The analysis of the Reference experiments is in any case restricted to the first 18 h. Indeed, gravity waves radiate energy out of the domain and induced turbulence mixes the area, giving a smoother and smoother vertical density profile as time goes on. In consequence, the characteristics of the gravity modes also change with time. As far as boundary reflections are concerned, sensitivity experiments carried out over twice-as-large domains showed no clear improvement: the sources are far enough from the open boundary to avoid any strong reflection.

Gravity wave modal reconstruction

In Experiments A-C, the complete modal representation of a gravity wave generated by a single source was studied in detail. When several sources generate a few gravity waves propagating in different directions, the identification of the transient processes, and the determination of their origin, phase or even amplitude is much more complicated. In the present section, three experiments are conducted and gravity waves are studied using a multivariate WEof decomposition of model vertical velocity and surface elevation anomaly in a given analyzed region.

Experiment D is more precisely dedicated to frequency localization by wavelet analysis in the case of two gravity waves excited at the same location but at differing frequencies. Experiments E and F are two examples of waves generated at the same frequency (or at least in the same frequency band) but at different times and locations. In the latter experiment, one of the forcing sources is localized in time whereas the other lasts for the whole period. In the former, forcing sources are not imposed at the same location (details can be found in Table 1).

Frequency localization

In Experiment D, two external perturbations are applied for 56 h. The first one (Source 1) generates waves with a period T 1 equal to 2.65 h. It is imposed at the domain center and the surface elevation amplitude is forced to reach 5 cm (see Appendix A). The second forcing source (Source 2) oscillates at the same location but with a period T 2 = 6.12 h and an amplitude of 10 cm. Both external (Mode 0) and internal waves (Mode 1 and higher) are generated with different wavelengths and different phase velocities. The WEof analysis is based on vertical velocity and surface elevation fields over a grid with a horizontal resolution of 300 m over 13.1 km.

Several characteristics of the waves are investigated: frequency localization through a "scalogram" representation, and propagation angles and wavelengths (or phase velocities) associated with each frequency band (Appendix B). PCA analysis was carried out for the first three empirical modes. Propagation angles depend on wave frequency from Eq. ( 17). The theoretical angles for internal waves generated by Source 1 (T 1 = 2.65 h) and Source 2( T 2 = 6.12 h) are 87 • and 88.7 • , respectively. As the data analysis grid used for WEof decomposition is close to the one of Experiment C, we can expect to find similar wavelength values. The scalogram is shown in Fig. 5 for Sources 1 and 2 forcing point (point S D on Fig. 6a). Both periods are present at the center of the domain where surface wave energy is concentrated at approximately 2.6 and 6.1 h, respectively. We can further note the larger amplitudes of the wavelet coefficients in the 6.1 h band. This is due to the doubling of the forcing amplitude. The gap between periods T 1 and T 2 is large enough for the resulting gravity waves to be clearly identified and discriminated by wavelet analysis. Actually a Heisenberg box centered at T 1 is 0.7 h wide, whereas it is 1.8 h wide for T 2 . To be studied independently, the two periods must be separated at least by 3 h in so far as the boxes are not centered in period but in frequency. Two WEof analyses can thus be carried out, one over each frequency band.

Within the context of internal wave propagation, the first WEof explains more than 80% of the field variance in both cases. The results for surface elevation reconstruction are not shown here but they reflect the amplitude values of 5 and 10 cm for waves at approximate periods of T 1 and T 2 , respectively. The vertical velocity reconstructions for WEof 1 are given in Fig. 6a over the Source 1 frequency band and in Fig. 6b over the Source 2 frequency band.

In the case of reconstructed WEof 1 at approximately 2.6 h (Fig. 6a), the angle of propagation is very close to the expected one (87 ± 0.5 • for a theoretical value of 87 • ). As in Experiment C, Modes 5-10 are present. Wavelengths are found to be λ h = 1.9 ± 0.2 km and λ v = 98.4 ± 6 m (Mode 5). This is confirmed by the statistical approach, which shows that 29% of the points over the diagonal between 13.6 and 15.7 km from the origin and between 21.4 and 200 m deep have a horizontal wavelength of 1.91 ± 0.21 km. Moreover, the same statistical approach shows that 31% of the points in the previous region exhibit wavelengths λ h = 2.39 ± 0.21 km that correspond to Mode 4. Modes higher than 6 can eventually be shown to propagate in the region over the diagonal between 13.6 and 14.8 km from the origin and for depths between 154 and 246 m. Mode 10 can be found more particularly in this region with 14% of the points in the range λ h = 0.95 ± 0.21 km. Thus, the wavelengths agree well with theoretical values given in Table 2 except for the estimate of the horizontal wavelength for Mode 10. Fig. 6a suggests λ h = 0.89 ± 0.21 km, which is also the value found for most of the points in the area with statistics, leading to an error of 6.3%. As mentioned before, one reason for such a shift is connected with the consequences of grid scale, mode mixing and interactions. The smoothing of the density profile as time evolves is also a serious candidate to explain such a mismatch. As Mode 11 would be λ h = 0.86 ± 0.21 km, it is clear that the distinction between Modes 10 and 11 cannot be made with a grid spacing of 300 m.

As far as the reconstruction of WEof 1 at approximately 6.1 h (Fig. 6b) is concerned, the wave propagates from the forcing point following a ray which makes an angle of 88.6 ± 0.3 • with the vertical, in good agreement with the expected theoretical value (88.7 • ). The wave mode around the most energetic ray (white space on Fig. 6b) is found to have wavelengths λ h = 2.20 ± 0.21 km and λ v = 49.2 ± 4 m. The statistical approach over a diagonal section located between 13.6 and 16.1 km from the origin and between 32 and 217 m deep associates 21% of the points with these wavelengths. They correspond to Mode 10 for a wave propagating with a period of 6.1 h (theoretical results lead to λ h = 2.2 km and λ v =50m). The most important modes in terms of space coverage are found to be Mode 9 and higher.

The remaining WEofs computed at T 1 or T 2 account for interactions between the modes. WEof 2 computed at approximately 2.6 h gives a reconstructed vertical velocity map (Fig. 6c) that contains waves propagating from the domain center with an angle of 86.8 ± 0.5 • and others with an angle of 88.6 ± 0.4 • , which are close to gravity wave propagation at periods of 2.6 h (87 • for Source 1) and 6.1 h (88.7 • for Source 2).

Thus, we can conclude that the WEof method is able to separate two transient processes at the same location but with differing frequencies. Moreover, the different modes associated with each frequency, in this case of multi frequency fields, are similar to the modes found for single frequency modeling in the Reference experiments. The first WEof for each frequency is similar to the corresponding mode studied in Section 3.4. The higher WEofs clearly show interaction pictures between the two gravity waves propagating at T 1 and T 2 .

However, wavelet analysis is not able to discriminate all processes on the basis of their frequency. The smallest frequency scale that can be identified is given by the width of the corresponding Heisenberg box, i.e. ω = 8.3 × 10 -2 h -1 for Source 1 [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Two processes separated by a frequency difference smaller than ω cannot be efficiently identified in this frequency band (see also [START_REF] Ouergli | Low frequency modes and their link with summer monsoon activity described by wavelet analysis throughout India[END_REF]. For the present experiment, the Heisenberg box criterion associated with the analysis frequency is well suited to frequency band localization.

A second Experiment (D ′ ) is carried out with the same configuration of sources in order to study lower gravity wave modes from the WEof computation at T 2 . A scale based approach similar to the one proposed for the Reference simulation is used. The gravity wave mode studied over the same grid scale is Mode 1 (λ h = 9.5 km and λ v = 500 m). The associated horizontal phase speed is c h,1 =1ms -1 , which is five times the phase speed for Mode 5. Thus, one way to get Mode 1 is to change the analysis time step t. The new value is t ′ = t/4 = 116 s and WEof analysis is applied over a period of 4 h. The observed resulting wavelength is λ h = 8.6 ± 1 km (not shown), which must be compared to the theoretical 9.5 km. The difference can be largely explained by the fact that the PCA is performed between 1.7 and 3.6 h. About 3.6 h are necessary for Mode 1 to cross the domain (13 km), so this mode, which is generated at the domain center, has not crossed the complete domain during the PCA analysis time.

Time localization

Frequency localization is only one side of the problem. Indeed, transient processes can appear in the same frequency band and have clearly different forcing mechanisms. In such a case, the analysis cannot be solved by a single wavelet analysis. It would, however, be unrealistic to consider that any process can be efficiently separated from the rest of the dynamics on the basis of a single PCA. In such cases, WEof offers a complementary approach where the wavelet analysis is used to extract the chosen frequency band, and PCA is used to decorrelate the dynamical processes.

In the following, Experiments E and F are designed to study the time and space localization of internal waves using WEof analysis. In the first experiment (E), we evaluate the ability of PCA to reconstruct transient processes when two transient processes of the same nature are generated within the same frequency band but at different locations. In Experiment F, forcing mechanisms are not of the same nature and the resulting wave pattern is different.

Surface elevation and vertical velocity fields are extracted in the forcing frequency band and are further decomposed in WEofs using PCA analysis. Each of the first three WEofs is then investigated to understand its behavior and verify that the wavelengths and propagation angles found are associated with the modes observed in the Reference experiments.

Waves emanating from two internal forcing sources

Experiment E is based on a grid resolution of 300 m. The analysis parameters of vertical velocities and surface elevations are given in Table 3. The modeling period is 5 days (120 h) but WEof analysis begins after 20 h, to focus on the separation of the processes by the WEof method. Indeed the main purpose is to analyze a time-independent forcing associated with a time localized forcing: the propagation of the lower modes of the former are thus not crucial in the present study.

As stated in Table 1, two forcing sources are applied in the inner ocean at the same period of 2.6 h. Source 1 continuously generated waves at the domain center at a depth of 125 m. The second source is applied at the centre of the top right quarter of the domain. It is localized in time between 57.5 and 66.5 h and the amplitude of the forcing is twice the amplitude of time-independent forcing. Wavelet coefficients are computed in the 2.6-h band and the first three WEofs are shown in Fig. 7a-c. Fig. 7d gives the temporal variation of WEofs.

WEof 1 accounts for 68.8% of the data field variance. As in the previous experiments, several modes are present propagating from Source 1 (Fig. 7a). A mixing of modes higher than 6 that looks like Mode 11 can be observed. Theoretical values for the associated wavelengths are λ h = 868.5 m and λ v = 45.5 m while the observed values are λ h = 859 ± 210 m and λ v =44± 4 m. A focus on the area where Mode 11 seems visually to be present, i.e. along the diagonal between 8 and 9.8 km from the origin and from 50 to 250 m deep, reveals that 31% of the points correspond to Mode 11 propagation (λ h = 868 ± 210 m). Another maximum can be found at the wavelengths corresponding to Modes 6 and 7. Actually the lower mode observed on Fig. 7aisMode6(λ h = 1592 m and λ v = 83 m). Moreover, between 9.3 and 11.9 km over the diagonal section and for depths ranging from 50 to 256 m, 25% of the points show a propagation of the observed Mode 6. The remainder of the meaningful wave propagation occurs at Mode 9 (21% of the area), Mode 5 (18%), Mode 10 (16%) and Mode 4 (12%). The waves propagate with an angle of 86.9 ± 0.5 • , which has to be compared to the theoretical value of 87 • associated with a 2.6-h-period internal wave. A map of the surface elevation anomaly (not shown) for the same time period confirms both the localization of the source and the amplitude of the wave described by WEof 1.

Based on Fig. 7d, the amplitude of WEof 2 (explaining 22.8% of the field variance) is maximal during a period slightly longer than the period when Source 2 is active. The waves propagate at an angle of 87 ± 0.5 • with the vertical (Fig. 7b), in good agreement with the theoretical value of 87 • for an internal wave at 2.6 h. The dominance of Mode 6 with λ h = 1681.7 ± 210 m and λ v = 82.3 ± 8 m after 51 h is evident in Fig. 7b. These results agree well with the theory which gives λ h = 1592 m and λ v = 83 m for the same mode (Table 2). Focusing on the area along the diagonal of the domain between 11.9 and 12.3 km from the origin and 50 and 255 m deep, 56% of points are associated with Mode 6. Here, the statistical approach provides useful information about modal patterns: Mode 5, for instance, appears to explain 30% of the propagation in the previous domain. We can thus conclude that WEof 2 is associated with the transient process localized in time.

The decorrelation of the two waves is not complete however. From Fig. 7b one can sense interactions, which are confirmed by the amplitudes shown in Fig. 7d. A glance at the temporal evolution of WEof 1 amplitude indicates that WEof 1 contains information from the Source 2 waves generated between 57.5 and 66.5 h. To confirm and study these interactions, an additional "Reference" experiment is carried out with only one continuous source of wave generation at the domain center (not shown). A comparison with Fig. 7a reveals obvious differences near the Source 2 generation area, which confirms an interaction between the two gravity waves in this area.

Finally, from Fig. 7d, WEof 3 amplitude can be observed to reach a maximum at the rise and the fall of Source 2 forcing, indicating that this WEof is directly linked with transient processes. It is associated with the appearance and disappearance of the temporally localized process. Moreover, the vertical velocity pattern shows waves generated by Sources 1 and 2 (Fig. 7c).

In conclusion, Experiment E shows the ability of PCA to extract processes generated at the same frequency but at different locations when one of them at least is localized in time. There are limitations on separating two processes at the same frequency, which are clearly associated with the ability of the PCA to decorrelate physical processes. Wavelet analysis in this case provides an interesting tool in frequency space to increase the chances of success of PCA.

Waves generated by internal and external forcing sources

Experiment F is very similar to Experiment E except for the fact that the time localized forcing is applied to surface elevation (Source 2, leading to an anomaly of 2 m). Permanent forcing is thus applied in the middle inner fluid, which leads to a smaller anomaly (Source 1). An interesting result of Experiment F is the ability of the multivariate WEof method to analyze processes when perturbation energy is mainly carried by transient processes of different natures: Mode 0 (the external gravity wave) propagates most of the perturbation energy in Experiment F while, in Experiment E, most of the energy is associated with the lower internal modes. Results for vertical velocities are given in Fig. 8a (WEof 1 reconstruction) and Fig. 8b (WEof 2 reconstruction). The temporal variation of the two WEofs is shown in Fig. 8c.

WEof 1 in Fig. 8a represents waves generated by Source 2 forcing localized in time (also evidenced by surface elevation maps not shown here). The associated wavelengths are found to be associated with Modes 4-9. For Mode 4, we obtain λ h = 2395 ± 210 m and λ v = 119 ± 7 m, which are close to theoretical values (λ h = 2389 m and λ v = 125 m). This mode is the most energetic mode that can be observed near the wave generation area. About 16% of the points located over the diagonal of the domain between 9.7 and 12.3 km from the origin at depths between 21 and 217 m (nearly all the water column to the left of Source 2) have wavelengths close to the theoretical wavelength of Mode 4 after 61.9 h (λ h = 2389 ± 210 m). Other modes are rather difficult to identify. However, Mode 9 is observed close to Source 1 in a region where several interferences can also be found. For this mixing of modes λ h = 1144 ± 210 m and v =53± 6 m (to be compared with the theoretical values for Mode 9: λ h = 1061 m and λ v = 56 m). Between 8.5 and 10.6 km from the origin over the diagonal and for depths from 21 to 217 m, Mode 9 is obtained for 33% of the points. The propagation angle of the rays is 86.7 ± 0.7 • , not far from the theoretical propagation angle of 87 • .

WEof 2 (Fig. 8b) clearly shows waves propagating from Source 1, the time independent forcing. This is confirmed by the surface elevation maps (not shown) where the amplitude is maximum at Source 1. The associated modes range from Modes 9 to 11. For Mode 9, measured wavelengths are λ h = 1002 ± 286 m and λ v =55± 4 m, not far from theoretical values (λ h = 1061 m and λ v = 56 m). Statistics performed over the diagonal between 8.5 and 10.2 km from the origin and for depths of 21 to 217 m show that 65% of the area exhibits propagation with Mode 9 wavelength characteristics (λ h = 1061 ± 210 m). The remainder of the area shows smaller wavelengths (i.e. higher baroclinic modes).

For higher modes, wavelengths are more difficult to extract because of the model resolution. One finds λ h = 858 ± 210 m and λ v =46± 4 m, which are close to the theoretical values for Mode 11: λ h = 868.5 m and λ v = 45.5 m. A focus on the diagonal between 9.3 and 11.5 km and for depths from 21 to 217 m reveals that Mode 11 propagation affects 62% of this part of the domain. These modes propagate along rays making an angle of 86.9 ± 0.5 • with respect to the vertical, in good agreement with the theoretical value.

WEof 1 explains 73% of the field variance while WEof 2 accounts for 17.4% of this variance. The fact that the time localized process appears in WEof 2 in Experiment E, but in WEof 1 in Experiment F is to be related to the nature of the forcing leading to high amplitude surface elevation and vertical velocities. From Fig. 8c, one can clearly see that, during the period of Source 2 activity, the amplitude of WEof 1 (related to waves generated mainly by Source 2) increases while WEof 2 amplitude decreases. Part of the gravity waves generated by Source 1 can thus be found in WEof 1.

Sensitivity to time localization

Experiment F can be further investigated to study the separation of two processes depending on their duration. Let δt 1 and δt 2 be the duration of the forcing for Sources 1 and 2, respectively. Tests on several configurations with various values of the δt 1 /δt 2 ratio show that the reconstruction results are directly linked with this ratio. It plays an important role in both wavelet analysis and PCA. The larger the ratio, the more accurate the analysis, as δt 2 has to be small enough compared to δt 1 to allow the two processes to be well separated. A lower bound of the ratio is difficult to define because the results deteriorate progressively as the ratio decreases. The results presented for Experiments E and F were obtained for δt 1 /δt 2 = 13.3. Further testing shows that an efficient WEof decomposition in wavelet space can be obtained for a ratio δt 1 /δt 2 > 3. For smaller ratios, the results are rather poor, with a mixing of physical processes within the same WEof. In the case of wavelet analysis, an additional loss of information occurs because of the Heisenberg box resolution.

The case of two forcing sources at the same location and at the same frequency is also studied. Two external forcing sources act at the domain center at the same period (2.6 h), but with different durations (δt 1 for the time-independent waves and δt 2 for the time-localized waves generated at the end of modeling). The time-localized forcing is twice as strong as the time-independent one. Results are just mentioned here and the associated figures are not shown. In this case, WEof analysis separates the two internal wave propagations even if they originate from the same location. For a δt 1 /δt 2 ratio of 3, WEof 1 accounts for the most energetic process related to the higher surface elevation (localized in time). However, it contains parts of the information associated with the time-independent waves (WEof 1 amplitudes are non-zero before the time-localized forcing is applied). WEof 2 represents the waves generated by the time-independent forcing source but contains interactions with waves generated by the time-localized forcing source. WEof 3 represents the transitory period when the time-localized forcing is imposed and removed.

Thus, WEof decomposition is an efficient tool for extracting time-dependent physical scales from a complex ocean background.

Discussion

The main purpose of the study is the analysis of transient dynamical processes in the coastal ocean based on the WEof method. We have chosen to illustrate such transients by what is certainly one of the most fundamental adjustment processes in geophysical fluid dynamics at coastal ocean scales: gravity waves. Indeed, any perturbation of the ocean mass field leads to a dynamical adjustment that precludes any possible equilibrium. These dynamical processes have been under study for a long time both in the atmosphere and in the ocean, where analytical models provide the simplest situations. In coastal oceanography, the dynamics of such transient processes, whether they are restricted to gravity waves or not, is part of the kinematical and dynamical fields.

However, coastal modeling cannot and must not be restricted to process-oriented modeling. [START_REF] Auclair | The penetration of the Northern Current over the Gulf of Lions (Mediterranean) as a downscaling problem[END_REF] show, for instance, that coastal dynamics exhibits strong non-linearity over the shelf, leading to non-Gaussian statistics. We therefore revisited the fundamental concept of scale modeling. The numerical model is clearly seen as an inves-tigation tool suited to well-defined time and space scales, and consequently to a particular frequency band. The difference with process-oriented modeling is that non-linear interactions between dynamical features within the studied frequency band are taken into account. In Experiments A-C, the same surface forcing is applied to a flat bottom, linear stratified, hydrostatic ocean. Each experiment is dedicated to the study of a particular mode or a group of modes of the gravity wave generated and leads to the evaluation of the modeling process when compared to theory.

We show in the present paper that WEof analysis is a key tool in scale-oriented modeling. Experiments D-F were indeed associated with the generation of simultaneous transient gravity waves well-separated in frequency bands, at different locations and even over differing time scales. In every case, it was shown that, with certain restrictive assumptions, it was possible to have access to any physical scale of any transient process: each WEof is thus associated with a theoretical mode or group of modes of the propagating gravity wave.

Several restrictive conditions for WEof analysis have also been identified. They are mainly associated with wavelet analysis as they concern the wavelet analysis Cone of Influence (Section 2.1) or the Heisenberg wavelet box. As far as time localized dynamical processes are concerned, it was further shown that their time scale must differ by a factor of at least 3. A clear understanding of the wavelet analysis time-frequency resolution probably remains one of the more fundamental restrictions. The Heisenberg boxes must be recognized as the fundamental atoms. Indeed, the use of a continuous wavelet analysis could easily lead to the wrong conclusion that any frequency discrimination is possible. Another restriction comes from the fact that any analysis for which the fundamental time scale is not much larger than the record length will encounter difficulties in resolving the dynamics.

The WEof technique should be extended to an automatic identification of transient processes [START_REF] Bakshi | Multiscale PCA with application to multivariate statistical processes monitoring[END_REF]. Such monitoring can be based on a systematic computation of scalograms within the modeled frequency band leading to an alert strategy when the energy spectrum overshoots a chosen threshold.

The consequences in terms of modeling are numerous. First, this transient detection provides a dynamical analysis of the modeling field; second, transient processes are usually difficult to identify among stochastic processes and yet their consequences are crucial to the understanding of any dynamical evolution and most energy transfers within the coastal ocean. A selective scale-oriented control of model trajectory can eventually be based on such monitoring. Instant detection of a transient process whose time scale is clearly known provides support for further data assimilation.

Conclusion

WEof analysis, a combination of wavelet and principal component analysis, has been applied to gravity wave, scale-oriented modeling. This multivariate, fully 3D technique appears to be a suitable tool for extracting transient processes within a given frequency band.

To evaluate its performance and understand the restrictive assumption it is based on, WEof has been successively applied to separate gravity waves generated by simultaneous forcing sources acting either at the same location but at different frequencies or at the same frequency but over different time periods and at different locations. In every case, WEofs perform well as long as a certain number of fundamental rules are obeyed. That is, respect of the wavelet cone of influence, recognition of the wavelet Heisenberg boxes as basic analysis atoms and clear localization of the processes.

Appendix A. Generation of gravity waves

Several techniques can be used to generate numerical gravity waves with well-known characteristics in a 3D coastal ocean model. A possible approach consists in applying an oscillating perturbation on surface elevation. In this case, the equations of motion for the 3D model SYMPHONIE become (for a vanishing Coriolis parameter): where u, v, w are the components of the total velocity, η the surface elevation anomaly η forcing is the applied perturbation of a few tens of centimeters, A the horizontal viscosity coefficient fixed from 1.5 to 250 m 2 s -1 and the term K ∂u,v ∂z stands for the Reynolds stress associated with turbulent flows.

∂u
A second method to generate internal waves is to apply a perturbation in the inner water column. This can be associated with the experiment of an oscillating circular cylinder exciting internal waves in a stratified fluid (for details see [START_REF] Sutherland | Visualisation and measurement of internal waves by "Synthetic Schlieren". Part 1: vertically oscillating cylinder[END_REF]. The forcing is applied via the pressure gradient so that the vertical velocity satisfies:

∂ 2 w ∂z∂t = 1 ρ 0 ∂ ∂x 2 + ∂ ∂y 2 P F (21)

Fig. 1 .

 1 Fig. 1. St. Andrew's cross: propagation of an internal wave in a continuously stratified fluid. Phase velocity c ϕ , wave number k and group velocity ៝ c g. Rays of energy propagation are emanating from the forcing area and make an angle θ 0 with the vertical.

  Fig.2shows a diagonal cross section of the vertical velocity (Fig.2a) and the surface elevation (Fig.2b) for Reference Experiment A after 23 h. The propagation being symmetrical, only one-quarter of the modeling domain is shown. The vertical velocity cross section clearly displays a wave propagating nearly horizontally. Fig.2bexhibits an isotropic horizontal propagation. A direct measurement gives λ h = 465.2 ± 35 km for the wavelength and c h = 48.7 ± 3.6ms -1 for its phase speed. The statistical approach confirms this measurement with 80% of the computed length scale in the area in the range 473 ± 35 km. It can further be noted that all the computed wavelengths are between 370 and 550 km (the horizontal model resolution is 70 km over the diagonal). These results are in good agreement with the theoretical results for Mode 0 given in Table2, i.e. λ h = 473.1 km and c h,0 =49.5ms -1 .

Fig. 3 .

 3 Fig. 3. Experiment B: vertical velocity (m s -1 ) for Mode 1 and higher after 8 h over Section [P 1 P 2 ] between 33.5 and 67.5 km from the origin. Values are restricted to [-2e-05, 2e-05].

Fig. 4 .

 4 Fig. 4. Experiment C: vertical velocity (m s -1 ) over Section [P 1 P 2 ], between 33.5 and 42 km from the origin, for Mode 5 and higher after 19 h. Two locations (M 1 ) and (M 2 ) are considered for horizontal and vertical wavelengths λ h and λ v . Values are restricted to [-2.7e-05, 2.7e-05].

Fig. 5 .

 5 Fig. 5. Experiment D with two surface forcing at T 1 =2.65handT 2 = 6. 12 h. Scalogram for the surface elevation at point S D (Fig. 6a). Values are restricted to [2e-03, 1.1e-02].

Fig. 6 .

 6 Fig.6. Experiment D with two forcing at the same location (point S D ) but at different periods (T 1 = 2.65 h and T 2 = 6.12 h). Reconstructed vertical velocity (m s -1 ) over

Fig. 7 .

 7 Fig. 7. Experiment E with two internal forcing at T = 2.65 h but at different locations (S E,1 and S E,2 ). Reconstructed vertical velocity (m s -1 ) over Section [P 1 P 2 ] between

Fig. 8 .

 8 Fig. 8. Experiment F with an external and an internal forcing both at T = 2.65 h. Reconstructed vertical velocity (m s -1 ) over the Section [P 1 P 2 ] between 8.1 and 15.2 km from the origin. (a) WEof 1 (73% of variance explained) standing for the external forcing, after 61.9 h. Values are restricted to [-6e-05, 6e-05]; (b) WEof 2 (17.4% of variance explained) standing for the internal forcing, after 43.6 h. Values are restricted to [-9e-06, 9e-06]; (c) amplitude vs. time for the first two WEofs.

Table 1

 1 Characteristics of the oscillating sources and studied processes in Experiments A-F

		Studied processes	Source S 1 Source S 2	Spatial localization
	Experiments A, B and C Experiment A: Mode 0	Period T 1
		Experiment B: Modes 1-2	Permanent
		Experiment C: Modes 5-10		-
	Experiment D	Modes 4-5-10	Period T 1	Period T 2
			Permanent Permanent
	Experiment E	Modes 6-9-11	Period T 1	Period T 2 = T 1
			Permanent Localized over t 2
	Experiment F	Modes 4-9-11	Period T 1	Period T 2 = T 1
			Permanent Localized over t 2

Table 2

 2 Theoretical characteristics of the waves for several modes at periods T 1 = 2.65 h and T 2 = 6.12 h: horizontal and vertical phase speeds (c h and c v ), horizontal and vertical wavelengths (λ h and λ v )

		Mode 0	Mode 1	Mode 4	Mode 5	Mode 6	Mode 9	Mode 10	Mode 11
	λ v (m) (T 1 , T 2 )	-	500	125	100	83	56	50	45
	c v (m/s) (T 1 )	-	0.0523	0.0131	0.0105	0.0087	0.0058	0.0052	0.0048
	c v (m/s) (T 2 )	-	0.0227	0.0057	0.0045	0.0038	0.0025	0.0023	0.0021
	λ h (km) (T 1 )	473.1	9.554	2.389	1.911	1.592	1.061	0.955	0.868
	λ h (km) (T 2 )	1091	22.039	5.509	4.408	3.673	2.449	2.204	2.003
	c h (m/s) (T 1 , T 2 )	49.52	1	0.25	0.2	0.17	0.11	0.10	0.09

Acknowledgments

The authors would like to thank the Laboratoire d'Aérologie (Toulouse, France) computer team: Serge Prieur, Didier Gazen, Juan Escobar and Laurent Cabanas for their support. No internal wave modeling would have been possible either without the support of Patrick Marsaleix (LA, Toulouse, France) in the use of the SYMPHONIE model. We would also like to thank the reviewers for their meticulous review of the paper.

with P F the forcing pressure taken to be of the form P F (x, y, z) = P 0 exp -(xx 0 ) 2 + (yy 0 ) 2 R 2 cos(ωt) sin 2π zz 0 8 (22) where P 0 = 0.05, R = 1 and x 0 , y 0 and z 0 give the source location.

Appendix B. Determination of wavelength

Horizontal and vertical wavelengths can be determined by two methods. They are first classically and objectively measured directly over the surface elevation and vertical velocity plots as the average distance between two wave crests in the area of interest. The modal number is then deduced from Table 2 as the mode corresponding to the measured wavelength.

A second method based on a statistical approach is proposed to confirm the previous objective determination of the mode number. If we consider the evolution of the vertical velocity for a plane wave:

where W is the wave amplitude, φ the phase, k the total wave number and r the distance from the source. The horizontal and vertical phase speed can be obtained as:

Horizontal and vertical wavelengths are then directly computed from the components of the wave vertical velocity. However, such a determination may seem rather inappropriate as it is based on local properties (the phase and its gradients), while the wavelength is a global property in so far as it is defined for a wave propagating in space (the application of the local frequency and wave number can be found in [START_REF] Bretherton | The general linearised theory of wave propagation[END_REF][START_REF] Phillips | The Dynamics of the Upper Ocean[END_REF][START_REF] Leblond | Waves in the Ocean[END_REF]. The same conclusion could be drawn from the previous measurement of the wavelength, but the measurement is carried out objectively, i.e. in a region, which is visually identified as showing the most salient oscillations. The rationale is to replace this visual objective determination by a statistical approach in the oceanic region where the measurement was made or where some information is crucially needed. Histograms of wave velocities obtained can consequently be computed, the width of each range being associated with an a priori uncertainty.

To compute uncertainties, we consider the maximum and minimum values possible for a given wavelength. Results are then given for the mean value plus or minus its difference with an extreme value. In all cases, uncertainties are at least equal to the grid resolution, the result varying between the mean value plus or minus at least half of the grid resolution (vertical or horizontal).