Fabrice Collard

Fabrice Ardhuin

Bertrand Chapron

Routine monitoring and analysis of ocean swell fields using a spaceborne SAR

Satellite Synthetic Aperture Radar (SAR) observations can provide a global view of ocean swell fields when using a specific "wave mode" sampling. A methodology is presented to routinely derive integral properties of the longer wavelength (swell) portion of the wave spectrum from SAR Level 2 products, and both monitor and predict their evolution across ocean basins. SAR-derived estimates of swell height, and energyweighted peak period and direction, are validated against buoy observations, and the peak directions are used to project the peak periods in one dimension along the corresponding great circle route, both forward and back in time, using the peak period group velocity. The resulting real time dataset of great circle-projected peak periods produces two-dimensional maps that can be used to monitor and predict the spatial extent, and temporal evolution, of individual ocean swell fields as they propagate from their source region to distant coastlines. The methodology is found to be consistent with the dispersive arrival of peak swell periods at a mid-ocean buoy. The simple great circle propagation method cannot project the swell heights in space like the peak periods, because energy evolution along a great circle is a function of the source storm characteristics and the unknown swell dissipation rate. A more general geometric optics model is thus proposed for the far field of the storms. This model is applied here to determine the attenuation over long distances. For one of the largest recorded storms, observations of 15 s period swells are consistent with a constant dissipation rate that corresponds to a 3300 km e-folding scale for the energy. In this case, swell dissipation is a significant term in the wave energy balance at global scales.

Introduction

Storms over the ocean produce long surface gravity waves that propagate as swell out of their generation area. In deep water, the wave phase speed C and period T are proportional. As the phase speed of the dominant waves Cp does not exceed 1.2 times the wind speed at 10 m height U10 [START_REF] Pierson | A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii[END_REF], the longest period waves must be generated by very intense winds. For example, the generation of waves of period T larger than 16 s requires winds with speeds over 18 m s -1 blowing over a distance of the order of 1000 km, to produce a significant energy, or yet stronger winds over a shorter fetch [START_REF] Munk | Directional recording of swell from distant storms[END_REF]. Such a large region of high winds is generally associated with a smaller storm center from which the long swells radiate. Waves further evolve after their generation with an important transfer of energy towards both high and low frequencies, due to nonlinear wave-wave interactions. Away from that core, nonlinear interactions become negligible [START_REF] Hasselmann | Part 3. evaluation on the energy flux and swell-sea interaction for a Neuman spectrum[END_REF] and long period swells have been observed to propagate over very large distances, up to half-way around the globe [START_REF] Munk | Directional recording of swell from distant storms[END_REF], radiating a large amount of momentum and energy across ocean basins. This measurable long-distance propagation is made possible by a limited loss of energy.

Copyright 2009 by the American Geophysical Union. 0148-0227/09/$9.00

The wave field at any time t, latitude φ and longitude λ, is described by its spectral densities G as a function of frequency f and direction θ. In the limit of geometrical optics, the spectral density is radiated at the group speed Cg in the direction of wave propagation, and can be expressed as a function of G at any previous time t0. Allowing for a spatial decay at a rate µ, the spectral energy balance is [e.g. [START_REF] Munk | Directional recording of swell from distant storms[END_REF] G(t, φ, λ, f, θ) = G (t0, φ0, λ0, f, θ0) exp t t 0 -µCgdt .

(1)

In deep water without current, the initial position φ0, λ0 and direction θ0 are given by following the great circle that goes through the point of coordinates φ, λ with a direction -θ over a distance X = (t -t0)Cg = (t -t0)g/(4πf ). This corresponds to a spherical distance α = X/R along the great circle, where R is the Earth radius. Equation (1) can be used to invert µ from wave measurements. For swell periods shorter than 13 s, [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF] have measured an e-folding scale Le = 1/µ = 5000 km (this number corresponds to a 0.1 dB/degree attenuation in their analysis). For larger periods, [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF] could only conclude that Le is larger, possibly infinite. In the past 40 years, little progress has been made on these conclusions [WISE Group, 2007]. Yet this question if of high practical importance, either for wave forecasting [e.g. [START_REF] Rascle | A global wave parameter database for geophysical applications. part 1: wave-current-turbulence interaction parameters for the open ocean based on traditional parameterizations[END_REF] or other geophysical investigations regarding air-sea fluxes or microseismic noise [e.g. [START_REF] Grachev | Upward momentum transfer in the marine boundary layer[END_REF][START_REF] Kedar | The origin of deep ocean microseisms in the north Atlantic ocean[END_REF].

A theoretical upper bound for Le is given by the viscous theory [START_REF] Dore | Some effects of the air-water interface on gravity waves[END_REF], see also Appendix A for a simple derivation]. According to theory, the largest shears are found right above the water surface, and the air viscosity dominates the dissipation of swells, giving, in deep water,

Le,max = ρwg 2 4ρaσ 3 √ 2νaσ , ( 2
)
where νa is the air viscosity, σ = 2π/T . For T = 13 s this gives Le,max = 45000 km, which means that over a realistic propagation distance of 10000 km the energy of 13 s swells is only reduced by 25%. Swells are thus expected to be very consistent over distances that are only limited by the size of ocean basins. The analysis of swells at this global scale should provide insights into their dynamics, including propagation and dissipation, but also into the structure of the generating areas, in a way similar to the use of the cosmic microwave background for the analysis of the early universe.

The present paper provides two important intermediate steps toward that goal. First, we demonstrate in section 2 how sparse data from a single space-borne synthetic aperture radar can be combined dynamically to provide a consistent picture of swell fields. This internal consistency reveals the quality of the SAR-derived dataset which we further verify quantitatively with buoy data. In section 3, we discuss and derive the asymptotic far-field swell energy evolution. Numerical investigations are performed to check the validity of the asymptotic solutions. This result provides a tool to interpret measured swell heights in terms of propagation and dissipation. This method is illustrated with one example that corresponds to a strong swell dissipation. Conclusions follow in section 4.

Space-time consistency of space-borne swell observations

Investigations by [START_REF] Holt | Tracking storm-generated waves in the northeast pacific ocean with ERS-1 synthetic aperture radar imagery and buoys[END_REF] and [START_REF] Heimbach | Development and application of satellite retrievals of ocean wave spectra[END_REF] have shown that space-borne synthetic aperture radar (SAR) data can be used to image the same swell field over 3 to 10 days as it propagates along the ocean surface. These preliminary studies have shown that the combination of SAR data at different places and times yields a position of the generating storm, and predictions for the arrival time of swells with different periods and directions. [START_REF] Heimbach | Development and application of satellite retrievals of ocean wave spectra[END_REF] have further pointed to shortcomings in the wind provided by an atmospheric circulation model for a given Southern Ocean storm, based on systematic biases in wind-forced wave model results compared to SAR observations. Unfortunately, the systematic analysis of such data has been very limited, and generally confined to data assimilation in wave forecasting models [e.g. [START_REF] Hasselmann | An optimal interpolation scheme for the assimilation of spectral wave data[END_REF][START_REF] Breivik | Assimilation of ERS SAR wave spectra in an operational wave model[END_REF][START_REF] Aouf | On the combined assimilation of RA-2 altimeter and ASAR wave data for the improvement of wave forcasting[END_REF]]. This narrow use of SAR data is due to three essential difficulties.

First, a SAR image is not a picture of the ocean surface and the relationship between the spectrum of the SAR image and that of the ocean surface elevation is nonlinear and fairly complex [e.g. [START_REF] Krogstad | A simple derivation of Hasselmann's nonlinear ocean-synthetic aperture radar transform[END_REF]. Sophisticated methods have been developed in order to estimate the surface elevation spectrum [e.g. [START_REF] Hasselmann | An improved algorithm for the retrieval of ocean wave spectra from synthetic aperture radar image spectra[END_REF][START_REF] Schulz-Stellenfleth | A parametric scheme for the retrieval of two-dimensional ocean wave spectra from synthetic aperture radar look cross spectra[END_REF]. These methods had to be implemented by the user of the data, and generally required some a priori first guess of the wave field provided by a numerical model. For longer wave systems, the imaging mechanisms are essentially quasi-linear, making possible a simpler methodology used by the European Space Agency (ESA) to generate a level 2 (L2) product. The method is fully described by [START_REF] Chapron | Wave and wind retrieval from SAR images of the ocean[END_REF]. It uses no outside wave information, and builds on the use of complex SAR data developed by Engen and [START_REF] Engen | Sar-ocean wave inversion using image cross spectra[END_REF] to remove the 180 • directional ambiguity in wave propagation direction. The quality of the L2 data has been repeatedly analyzed [e.g. [START_REF] Johnsen | Envisat ASAR wave mode products -quality assessment and algorithm upgrade[END_REF][START_REF] Collard | Extraction of coastal ocean wave fields from SAR images[END_REF]. Because long ocean swells have large wavelengths and smaller steepnesses, the L2 products corresponding to this spectral range have higher relative quality, confirming that the imaging mechanism is well described under the quasi-linear assumption.

All SAR data used here are such L2 products, provided by ESA and obtained with the L2 processor version operational since November 2007, and described by [START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF]. The data for times before that date were reprocessed with this same processor. In previous realtime data, frequent low wavenumber artefacts were caused by insufficient filtering of non-wave signatures in the radar images. This filtering is necessary to remove the contributions of atmospheric patterns or other surface phenomena like ships, slicks, sea ice, or islands, with spectral signatures that can overlap the swell spectra. The L2 product contains directional wave spectra with a resolution of 10 • in directions and an exponential discretization in wavenumbers spanning wavelengths of 30 to 800 m with 24 exponentially spaced wavenumbers, corresponding to wave periods with a 7% increment from one to the next.

The second practical problem is that the data obtained from an orbiting platform are sparse and with a sampling that makes a direct analysis difficult. Hereafter, we show that the space-time consistency of the swell field can be used to fill in the gaps in the observations and produce continuous observations of swell periods and directions in space and time.

Third, and last, for a simple use of SAR data, some parameters that are not affected by the variable resolution is SAR scenes [START_REF] Kerbaol | Analysis of ERS-1/2 synthetic aperture radar wave mode imagettes[END_REF]. [START_REF] Schulz-Stellenfleth | An empirical approach for the retrieval of integral ocean wave parameters from synthetic aperture radar data[END_REF] have proposed to produce parameters representing the entire sea state by extending the resolved spectrum with an empirical windsa contribution. Here we take the opposite approach and restrict the resolved part of the spectrum by using a spectral partitioning (see Appendix B for details) to retrieve the swell significant wave height Hss, defined as four times the square root of the energy of one swell system, and the peak period Tp and peak directions θp. Thus one SAR typically produced one or two swell parameters for distinct swell systems.

For Hss, only very limited validations have been performed [START_REF] Collard | Coastal ocean wave retrieval from ASAR complex images[END_REF]. We thus perform a thorough analysis of SAR-derived Hss comparing to co-located buoy data (see Appendix B for details).

The bias on Hss derived from ESA Level 2 products is found to be primarily a function of the swell height and wind speed, increasing with height and decreasing with wind speed. Variations in standard deviation are dominated by the swell height and peak period, with the most accurate estimations for intermediate periods of 14 to 17 s.For wind speeds in the range 3 to 8 m s -1 , Hss has a bias of 0.24 m and the standard deviation of the errors is 0.29 m.

We thus corrected the values of Hss by subtracting a bias model given by b h = 0.11 + 0.1Hss -0.1 max{0, U10SAR -7}

(3)

where Hss is in meters and the wind speed U10SAR is in m s -1 . From now on, all the reported values of Hss will be corrected using this expression. After correction, the standard deviation of Hss estimates is reduced to less than σ h = 0.10 + min {0.25Hss, 0.8} (4)

where σ h and Hss are in meters.

The quality of Tp and θp has already been carefully studied and are routinely monitored. The root mean square (r.m.s.) error on Tp is less than 10% of the measured value for Tp > 12 s, and the r.m.s. error on θp is 22 • for these same swells [START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF], with little bias. Because few directional in situ measurements are available, we demonstrate here an original semi-quantitative dynamical validation of these two parameters.

Virtual wave observers

Given these SAR-derived estimates of θp and Tp, linear dispersion relationship and the principle of geometrical optics can then be exploited to predict arrival times and locations of the swell.

virtual buoy position

t 1 (1,1) t 2 (1,1) t 0 (1) interrogation window t 0 (2) t 1 (2,1) t 2 (2,1) t 1 (2,2) t 2 (2,2) t 0 (3) t 0 (4) Figure 1
. Schematic definition of a virtual wave observer. Any SAR observation i is available at a time t0(i) on the black dots. All swell partitions (i, j) (here indicated by the arrows) are propagated and may cross the interrogation window from time t1(i, j) to t2(i, j).

In order to obtain swell conditions at the location of a "virtual wave observer", we define an interrogation window covering 2 by 2 degrees in latitude and longitude. According to the SAR-derived peak parameters, swell partitions from the entire ocean basin are propagated, both forward and backward, along great circles in space and time. These theoretical trajectories are followed with a constant group speed gTp/(4π), starting off in direction θp from the observation point. From any observation time t0, these great circles may cut through the interrogation window from times t1 to t2 (thick solid lines in figure 1, with different colors for different partitions). As the maximum value of |t2 -t0| and |t1 -t0| is increased from 6 hours to 6 days, the timeevolution of the peak frequencies and peak directions at the virtual observer gradually reveals similar ridges to the one observed in real buoy measurements (figure 2).

In fig. 2b-d, each horizontal colored segment corresponds to one swell partition that crosses the spatial window between times t1 and t2. Some segments are very short, corresponding to trajectories that barely cut one corner of the window.

Clearly the SAR detects the direction of the most energetic part of the wave spectrum measured by the buoys (fig. 2a). At frequencies above 0.1Hz, the virtual observer patterns appear rather noisy. Shorter scales are not so correlated. These shorter components are often observed as part of the wind sea for which the propagation with a single group speed and direction is not a good approximation. Also, propagated high frequency swells, such as the 0.12 Hz waves coming from direction 200 on July 10, do not show up in the real buoy record. This is possibly the result of a relatively high dissipation rate for these swells. For frequencies below 0.08 Hz, the virtual observer shows ridge-like structures similar to those observed in situ due to the dispersive arrivals of swells from remote storms [e.g. [START_REF] Munk | Directional recording of swell from distant storms[END_REF][START_REF] Gjevik | long period swell wave events on the Norwegian shelf[END_REF]. Even the faintest events are well detected, such as the 0.06 Hz arrival on 23 July, even though that 15 s swell of 0.5 m is dwarfed by a another 0.8 m and 12 s swell, and a 2 m and 8 s wind sea.

Swell detection with the virtual observer reaches its limits when the swell height is very low, such as on 3 July with a 20 s 0.3 m swell well detected by the buoy (the green-orange ridge at 0.05 Hz). Because we use a single trajectory emanating from one observed swell partition the relatively small interrogation window can easily be missed after 10000 km of propagation.

Our discrete propagation technique suffers from a randomized version of the "garden sprinkler effect" that, if not corrected for, can create unrealistic flower-like patterns in the far field of storms in numerical wave models that use a discretized spectrum [e.g. [START_REF] Tolman | Alleviating the garden sprinkler effect in wind wave models[END_REF]. Our choice of a single group speed and direction, because a narrow swell spectrum is not resolved by the SAR, produces a discrete wave field (the dots in figure 3). With the present processing this is smoothed by the finite size of our interrogation window (figure 1). An extension of the present technique could use neighboring group speeds and directions to take into account the frequency and directional spread of the sea state, which would allow the use of a smaller window. Just like the estimation of propagated wave heights, discussed below, the estimation of a spectra width that cannot be resolved by the SAR may use some further information on the structure of the generating storm.

Other errors in the present technique can also be attributed to the SAR processing. In particular, a maximum value is defined for the transfer function used to obtain the wave spectrum from the SAR image [START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF]. Although this is designed to prevent the amplification of measurement noise, long swells such as this 20 s event have very small slopes, and it is likely that they are underestimated in the wave spectrum due to this threshold in the processing.

When propagated for 6 days, without any other information than the peak frequency and direction at the time of observations, the waves are remarkably consistent with the latest local observations. For the southern swells arriving at Christmas Island between July 16 and 21 (figure 2), the difference in arrival times given by the virtual oberver and real buoy is typically less than 12 h. This is less than 10% of the maximum time between the SAR observation and the virtual observer record. This implies that the accuracy of the peak period estimate for each SAR partition must also be less than 10%, consistent with previous validation studies [START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF][START_REF] Johnsen | Envisat ASAR wave mode products -quality assessment and algorithm upgrade[END_REF]. The consistency of the arrival directions along the ridges also suggests that the root mean square (RMS) error in peak direction estimates must be close to 20 • , comparable to the 22 • RMS difference between mean wave directions obtained from SAR wave mode and a numerical wave model for waves with periods longer than 12 s [START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF].

Although it cannot replace the spectral resolving power of a buoy, the performance of the virtual observer is therefore comparable or better to that of human observers in terms of peak period and direction [START_REF] Munk | Refraction of ocean waves: a process linking underwater topography to beach erosion[END_REF]. The really missing bit is a wave height estimate along the swell propagation path. We will show that this may be obtained by estimating the source storm characteristics and the dissipation rate of swell energy.

Storm source identification and "fireworks"

Along the estimated trajectories, virtual observations can further be produced in a similar fashion. The animation of these propagated swells confirms the very well organized nature of storm swells crossing large ocean basins.

X -5

intense storms that occur over a period of several days (see figure 3 and auxiliary animation). In large ocean basins where swells are likely to be imaged several times by the same satellite, these fireworks can be used to estimate the time of arrival of swells from any given storm [e.g. [START_REF] Holt | Tracking storm-generated waves in the northeast pacific ocean with ERS-1 synthetic aperture radar imagery and buoys[END_REF]. For this reason, these animations have been produced routinely every day since August 2007, for the Pacific, Indian, and Atlantic oceans (see http://www.boosttechnologies.com/esa/images/, e.g. nrt_pac.gif for the Pacific Ocean).

Using backward trajectories, the location and date of swell sources can further be defined as the spatial and temporal center of the convergence area and time of the trajectories. These positions have been verified to correspond to high wind conditions observed by scatterometers and reproduced by ECMWF wind analyses. We consider these storms to be the source of all the swell partitions that produce trajectories that pass within 12 hours and 2000 km of their center. This processing, similar to the one performed by [START_REF] Heimbach | Development and application of satellite retrievals of ocean wave spectra[END_REF], provides a global view of swell fields in both space and time, extending the coverage of similar techniques based on buoy data [START_REF] Hanson | Automated analysis of ocean surface directional wave spectra[END_REF]. In figure 4, a swell covers one Earth quadrant away from the storm, with a large detection gap that extends from the Southern Pacific to California. This blank area is the long shadow cast by French Polynesia where wave energy is dissipated in the surf [e.g. [START_REF] Munk | Directional recording of swell from distant storms[END_REF]. Observations were restricted to swell partitions with periods close to 17 s, but the full dataset typically covers swells with periods of 12 to 18 s, as shown in figure 3.

The apparent self-consistency of both the virtual buoy plot (figure 2.d) and the fireworks animations, are the result of the large auto-correlation length of the swell fields, which was expected from the in situ measurements of [START_REF] Darbyshire | The generation of waves by wind[END_REF]; [START_REF] Munk | Directional recording of swell from distant storms[END_REF]; [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF]. Yet, these plots could not exhibit such patterns without a good accuracy of the SAR-derived peak periods and directions, used in the propagation methodology.

Far-field swell energy

All these illustrations of forward-backward ray tracing indicate the potential to use a simple Geometrical Optics (GO) strategy. The next goal is then to determine the strength of the far-field radiated swell energy. This requires the definition of a swell source, and an estimation of the swell energy dissipation scale. For this we define the time t0 as an initial condition after which there is no significant wave generation or non-linear evolution, for frequencies less than fmax. Namely, at t0 all the wave components with smaller frequencies have already been generated, so that the radiation of these waves is essentially fossil and fully governed by geometrical optics. The possible effects of diffraction and scattering are discussed by [START_REF] Munk | Directional recording of swell from distant storms[END_REF], and, together with dissipation, will cause deviations from the G.O. model outlined below. We therefore make no restrictive hypothesis on what happens before t0, and thus the motion of the generating storm has no direct effect on our results, but it obviously modifies the spatial distribution of the energy at t0, which will be relevant.

In reality, swells evolve over the course of their propagation as the result of their interactions with the local winds, mutual wave-wave interactions, interactions with other wave systems, including the local wind sea. Swells are also expected to evolve according to interactions with other oceanic motions that affect the upper ocean, namely surface currents, internal waves [e.g. [START_REF] Kudryavtsev | The coupling of wind and internal waves[END_REF] and turbulence. Depth and island scattering effects must also carefully be taken into account [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF]WISE Group, 2007]. Compared to these different mechanisms, frequency dispersion and angular spreading effects are certainly the first leading order phenomena to take into account for the major part of the decrease in the height of the swell systems. Indeed, as the long swell systems will be characterized by relatively small steepness parameters, nonlinear mutual wave-wave interactions do not appear to be important in scattering surface wave energy more than a few storm radius distances outside an active generating area Hasselmann [1963]. Furthermore, the level of turbulence in the ocean does not appear to significantly affect the waves Fabrikant and Raevsky [1994]; [START_REF] Ardhuin | On the interaction of surface waves and upper ocean turbulence[END_REF], and the conversion of surface wave energy into internal gravity wave energy by wave-wave interactions does not seem to be a leading order sink term for the energy balance of surface gravity waves. Finally, for the very long trans-ocean fast propagating swell components, surface current bending effects, proportional to the ratio between vertical current vorticity and the group velocity, may also be considered as residual effects. Away from island obstructions, the ratio between the angular width along the great-circle observatory points at very large distances from the generating area, and the mean spread in the generating area is approximately proportional to 1/ sin α with α the spherical distance from the storm. center The change in spectral density F defined by

F (t, φ, λ, f ) = 2π 0 G(t, φ, λ, f, θ)dθ (5)
follows the spatial expansion (close to the source) and contraction (as waves approach the antipodes for α > π/2) of the energy front. This transversal dispersion is associated to a narrowing of the directional spectrum G, for α < π/2, and a broadening for larger distances. This approximation applies to large distances, and relatively small source regions. Closer to the source, the approximation does not hold. Swell amplitudes radiating from large extended sources will decrease more slowly than swell amplitudes emanating from compact sources.

Moreover, we can represent swell waves as a linear superposition of harmonic waves in narrow spectral band.

E(t 0 ,x) "Source" at t=t 0 x E(t 1 ,x) E(t 2 ,x) F(t 0 ,x 0 ,f) f f=f 0 F(t 1 ,x 1 ,f) F(t 2 ,x 2 ,f) f=f 0 f=f 0 f f x 0 x 1 =x 0 +(t 1 -t 0 ) C g (f 0 ) x 2 =x 0 +(t 2 -t 0 ) C g (f 0 ) x P x Q Figure 5
. Dispersion of linear waves in one dimension. At any given time the spectrum is given by a propagation of the spectra at t = t0. Taking x1 = x0 + (t1 -t0)Cg(fp) the spectral density at fp is the same as for time t0, but the spectrum is narrower which gives a smaller elevation variance, E(t1, x1) < E(t0, x0).

Quite naturally, through the method of stationary phase, the group velocity is defined and the slowly-varying wave envelope is found to decay. This decay is inversely proportional to the square-root of the distance (figure 5). Accordingly, far away from the generating sources, and in the absence of dissipation, the swell energy

Es(t, φ, λ) = ∞ 0 F (t, φ, λ, f )df. ( 6
)
should decrease asymptotically like 1/(α sin α) when following a wave group (see Appendix C for a detailed proof).

3.1. The [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF] method

Using measurements with a limited or no directional resolution, [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF] assumed that wave propagation was completely blocked by waters shallower than 60 fathoms (approximately 110 m), and that diffraction could be neglected. For example, in figure 6, the island would be represented by the 60 fathom depth contour. These authors then estimated a loss of swell energy from the deviation of the ratio of directionally-integrated spectra, e.g. F (x2, f0)/F (x1, f0), compared to what is expected from the lateral dispersion effect, taking into account islands. Currents, shallow water areas and diffraction effects around islands are neglected here. These effects are discussed by [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF].

"Source" at t=t 0 x G(t 0 ,x 0 ,f 0 ,θ)

θ θ=θ p G(t 1 ,x 1 ,f 0 ,θ) G(t 2 ,x 2 ,f 0 ,θ) θ=θ p θ θ Island x 2 =x 0 +(t 2 -t 0 ) C g (f 0 ) x 1 =x 0 +(t 1 -t 0 ) C g (f 0 ) "Source" at t=t 1 "Source" at t=t 2 x 0 A B θ p θ 2 θ 1 θ 3 θ 4 θ 5 θ 6 θ 5 θ 6 θ 4 θ 3 C D Figure 6
. Dispersion of linear waves in two dimension, represented here on a flat surface for simplicity. The variable y has been dropped as the spectra shown here are all at y = 0. We call "source" the region of the ocean where waves with frequencies smaller than fmin can be found. As time goes by, the source expands in space due to both frequency dispersion (like in 1D), and geometrical dispersion. The wave energy with frequency f0 that will be observed at point x1 (respectively x2) at time t1 (respectively t2) is, at time t0, along the thin arc circle AB (respectively the thick dotted arc circle CD). Due to the small island between x1 and x2, the energy that would have been recorded at x2, if the island were not present, is actually dissipated on the shore of the island. As a result the local energy density E(x2) is reduced. At frequency f0, contributions to E(x2) only come from angles θ5 to θ6. The directional spectra (bottom) are thus affected by the blocking effect of islands, and the directional narrowing as one goes further from x0 (on the Earth this narrowing reverses after 10000 km of propagation, due to the sphericity).

Rigorously, their method is inexact because the recording stations x1 and x2 measure wave groups that had neither exactly same propagation directions nor the same position when they were near x0. Yet, because the wave field in the neighborhood of x0 is the superposition of many independent wave trains, one can assume that the spectral density F is a smooth function of the direction. Then we may say that over the intervals θ3 to θ4 or θ5 to θ6, G does not change so much, i.e. in figure 6, G(t0, φD, λD, f, θ6) ≃ G(t0, φD, λD, f, θ5). On the sphere the application of eqs.

(1) and (5) yields

F (t, φ, λ, f0) = 1 sin α D C G(t0, φ, λ, f0, θ0) ds R 1 + O 1 sin α , (7)
where the integral is performed over the line segment joining C to D. The error relative to the asymptote is the sum of two terms. One is proportional to the spatial gradient of F , due to the change from the arc circle to the segment, and the other corresponds to the relative variation of G with θ over the range θ5 to θ6, which is small in the far field, provided that the directional spectrum is smooth enough. Under these two smoothness assumptions, and for large sin(α), the ratios of spectral densities F at x1 and x2, as used by [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF], can be compared to the conservative equations ( 7) and used to diagnose the dissipation of swell energy. This does not apply if x1 is in the vicinity of the storm or its antipode, where the observed arrival direction span a large range.

In practice, this method can be very sensitive to the correct estimation of the island shadowing. For that reason, the measurement route chosen by [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF] was far from ideal. Because they needed land to install most recording stations for the wave measurements, they used an almost north-south great circle that extends from the south of New-Zealand (Cape Palisader) to Alaska (Yakutat), a route peppered with islands in its southern part, and partially blocked by the Hawaiian chain in its northern part. Also, storms typically refuse to line up with any measurement array. Their pre-defined great circle, although designed to follow a typical Southern winter swell propagation path, always deviated by some extent form the actual track followed by the most energetic swells they recorded. For the Indian ocean storms, this difficulty was reduced by the relatively narrow range of angles that allows propagation from the Indian to the Pacific ocean.

A method using global swell heights

Now using an instrument with a global coverage, we can carefully avoid both problems by choosing propagation paths far away from the smallest island, and by exploiting only observations well aligned with the storms. However, due to the limited spectral resolution inherent to the SAR wave mode image size and processing, we cannot use the spectral distribution G or F of the energy, and can only use the energy Es integrated over a swell partition. We thus need for Es the equivalent of eq. ( 7) for F . For simplicity, and without loss of generality, we take the source storm centered at time t0 = 0 on the pole S defined by a colatitude ϕ = φπ/2 = 0, so that the distance from the storm center is r = Rϕ. We consider the swell energy Es observed at a position defined by the spherical distance α and we take the reference meridian to be in the direction of the observation point (figure 7). We will later assume that the source area is relatively small with a size R∆α, where ∆α is the maximum value taken by ϕ (figure 7). In all the following derivations, we have chosen a fixed frequency f0 and we follow a wave group of that frequency. The time of observation t is thus related to α by t = Rα/Cg(f0), so that the variable t will be omitted.

In appendix C, we prove and verify that, in the absence of dissipation the swell energy Es decreases like 1/(α sin α) for large values of α, and that typical storms should produce swells within 20% of this asymptote at distances larger than 4000 km from the storm center. A much larger observed deviation should thus reflect a gain or loss of energy by the propagating swells. The expected departures from the asymptotic evolution should be compared to those due to swell dissipation or generation. Even with perfect SAR observations, this is the intrinsic limit of the present method. A 20% error in energy conserving conditions may be misinterpreted as a dissipation or generation with an e-folding scale of the order of 20000 km, which gives a 20% energy change as waves propagate from 4000 to 8000 km away from the storm source.

Illustration

To illustrate the method described above, we analyse of one of the most powerful swell field recorded over the past 4 years by ENVISAT's ASAR. The swell case illustrated in figure 4 is not well suited due to the islands in the southnorth swell tracks and the poor sampling of ENVISAT for the tracks going north-east from the storm, we have thus chosen another source, found on February 12 2007 at 18:00 UTC, and located at 168 E and 38 N. This swell was generated by a storm moving Eastward with maximum Westerly winds of 26 m s -1 at the indicated date, and subsiding to less that 22 m s -1 and veering to South-Westerly 12 h later (according to ECMWF 0.5 • resolution analyses). The storm motion in the direction of wave propagation certainly helped to amplify the local windsea, with a maximum significant wave height Hs of 14.1 m at 00:00 UTC on February 13 (according to a numerical wave model configuration that is otherwise verified to produce root mean square errors less than 9% for Hs > 8 m).

Using SAR-measured wave periods and directions at different times and locations, we follow great circle trajectories backwards at the theoretical group velocity. The location and date of the swell source is defined as the spatial and temporal center of the convergence area and time of the trajectories. We chose a central peak period, here 15 s, and track the swells forward in space and time, starting from the source center at an angle θ0, following ideal geodesic paths in search of SAR observations. Along each track, SAR data are selected if they are acquired within 3 hours and 100 km from the theoretical time and position. Great circle tracks are traced from the source in all directions, except for angular sectors with islands.

In order to obtain enough SAR data, we repeat this operation for regularly spaced values of θ0 with a step of 2 • .

In our case, when varying θ0 from 74 to 90 • (counted clockwise from North), this procedure produced 58 SAR measurements with one swell partition that had a peak wavelength and direction within 50 m and 20 • of the expected theoretical value.

If no energy is lost by the wave field, Es decreases asymptotically as 1/[α sin(α)] away from the source. Among the 58 swell observations, we further removed all the data within 4000 km of the source center, to make sure that the remaining data are in the far field of the storm, and data with a significant swell height Hss less than 0.5 m, after bias correction based on the error model. This makes sure that the signal to noise ratio in the image is large enough so that the wave height estimation is accurate enough.

We then have 35 observations for which we assume that Es is only a function of α, and we define the dissipation rate

µ = - 1 EsR d (Esα sin α) dα . ( 8
)
Positive values of µ correspond to losses of wave energy (figure 8). We then fit an analytical function Hss(α) to the data, defined by a constant µ and Hss(α = π/5), i.e. the swell height at 4000 km from the source. Here the couple Hss(α = π/5) = 4.4 m and µ = 3.7 × 10 -7 m -1 gives the least square difference between the decay with constant µ and the observed swell decays. Further, the uncertainty of that dissipation rate may be estimated from the known uncertainty of the SAR measurement of Hs, given by eq. ( 3)-( 4). A more simple error model, with larger errors based on the Hs12 analysis by [START_REF] Johnsen | Envisat ASAR wave mode products -quality assessment and algorithm upgrade[END_REF], does not significantly alter this analysis. Using that error model and neglecting other sources of error in the present analysis, we perturbed the observed swell heights independently to produce a 400 ensemble of synthetic data sets. Taking the 16% and 84% levels in the estimation of µ, that would correspond to one standard deviation for a Gaussian distribution, we find that 3.1 × 10 -7 < µ < 4.0 × 10 -7 at the 68% confidence level. This is the first ever estimation of the uncertainty on an observed swell dissipation rate. These values of µ are more than twice larger than reported by [START_REF] Snodgrass | Propagation of ocean swell across the Pacific[END_REF] for smaller amplitude swells. The formidable height of 4.4 m at a distance of 4000 km was observed by the SAR for all outgoing directions from at least 74 to 106 • . This same swell was also recorded by buoys in the North-East of Hawaii (NDBC buoy 51001), also with a peak period of 15 s, and a height of 3.4 m on 16 February 2007 at 0:00 UTC. That buoy is located 3300 km away from the center and in a direction close to 112 • . Looking in the North-East quadrant of the storm, one also finds a trace of the swell at buoy 46005, off the Washington coast (4900 km in direction 59 • ). There the swell was also observed with a 15 s period and a maximum height of 3.2 m on February 17 at 17:00 UTC, similar to the SAR observations for the same distance. For directions closer to Northbound, either the generation was weaker or the Alaskan islands sheltered the coastal buoys. For example, the same swell was also recorded by NDBC buoy number 46075, off Shumagin Island, Alaska, at a distance of 3000 km from the source, in the direction 42 • . At that buoy, the peak period was 15.0 s with a maximum swell height of 1.3 m on 15 February 2008 at 18:30 UTC.

Thus the power radiated by the storm is of the order of is 0.5 TW at 4000 km from the storm center, spread over a 50 • angular sector. This power is about 16% of the estimated 3.2 TW annual mean flux that reaches the world's coastlines [START_REF] Rascle | A global wave parameter database for geophysical applications. part 1: wave-current-turbulence interaction parameters for the open ocean based on traditional parameterizations[END_REF]. However, the observed dissipation rate corresponds to an e-folding scale of 3300 km for the energy. Taking an average propagation distance of 8000 km, only 160 GW would make it to the shore. If the same dissipation rate prevailed closer to the source, then the power radiated at 1000 km form the storm center was 1.4 TW.

We thus expect that the far field dissipation of swells, in spite of the small steepness of these swells, plays a significant role in the air-sea energy balance. This effect probably explains the systematic positive bias for predicted wave heights in wave models that neither account for swell dissipation nor assimilate wave measurements [see e.g. [START_REF] Rascle | A global wave parameter database for geophysical applications. part 1: wave-current-turbulence interaction parameters for the open ocean based on traditional parameterizations[END_REF].

Conclusions

Taking advantages of the satellite observations of unprecedented coverage and quality, investigations can repeat and complement the pioneering analysis of swell evolution performed in the 1960s. Severe storms can generate relatively broad spectra of large surface waves. But rapidly, the redistribution of energy, through linear dispersion and nonlinear interaction mechanisms, becomes very effective. The initial wind waves become swells outrunning the wind, leading to the apparition long-crested systems. The propagation properties of these surface gravity waves have been found to closely follow principles of geometrical optics. The consistent patterns of swell fronts dispersing over thousands of kilometers was shown to be useful to provide time series at "virtual wave observing stations", filling gaps in space and time in between the orbit cycles of observation. When compared to buoy measurements, the present results give an explicit dynamical validation of the SAR-derived spectral parameters. As the speed of waves in deep water is proportional to their period or wavelength, information carried by the SAR-derived period and direction distributions, observed at a fairly large distance from the generating area, pertains to the wind conditions existent up to 15 days before.

We also discussed how the swell energy should, in the absence of dissipation, decay in the far field of the storm like 1/(α sin α) where α is the spherical distance between the storm center and observation point. Exploiting that property allowed us to estimate a dissipation rate µ of swell energy with unprecedented accuracy, establishing that swell dissipation can be a significant term in the global wave energy budget.

The proposed methodolgy performed here requires data far enough from the source, typically more than 4000 km, in order to approach this simple asymptotic behavior. At the same time, the swell amplitude should be large enough to be accurately measured by the SAR. Some knowledge of the spectral shape and its spatial distribution inside the storm may be useful to provide better estimates of µ for low dissipation cases, or closer to the storm centers. These further analyses will likely benefit from the joint use of data from altimeters, SARs, and other sources of spectral wave information.

A more systematic analysis and interpretation of this dissipation will be reported elsewhere [Ardhuin et al., 2009a], with applications to wave forecasting models [START_REF] Ardhuin | Spectral wave dissipation based on observations: a global validation[END_REF][Ardhuin et al., , 2009b]]. The parameterization of the dissipation rate could also be used to produce a data-based forecasting system, extending our virtual buoy technique to the estimation of swell heights, with a forward propagation of observations. Going in the opposite direction, toward the storm source, it is possible that the analysis of swell fields could provide a "new" way of looking into the poorly observed structure of severe storms. Because the usual remote sensing techniques for estimating wind fields either do not work for very high winds or are not well validated [e.g. [START_REF] Quilfen | Altimeter dual-frequency observations of surface winds, waves, and rain rate in tropical cyclone isabel[END_REF][START_REF] Quilfen | The potential of quikscat and windsat observations for the estimation of sea surface wind vector under severe weather conditions[END_REF], the use of far-field swell information may provide an interesting complement to the local wind speeds and wave heights. The is inverse problem has already been formulated by [START_REF] Munk | Directional recording of swell from distant storms[END_REF] who already proposed an elegant heterodyning technique to push the spatial resolution for the estimation of storm location from swell data, while [START_REF] Heimbach | Development and application of satellite retrievals of ocean wave spectra[END_REF] have proposed to use wave models to correct wind field errors. The quality of the SAR-derived swell parameters that are coming out of today's ENVISAT and tomorrow's Sentinel-1, together with a good understanding of the swell energy budget, including its dissipation revealed here, may finally enable this vision.

Appendix A: Viscous theory for air-sea interaction

For the sake of simplicity we will consider here the case of monochromatic waves propagating in the x direction only, and we will neglect the curvature of the surface. For the small steepness swells considered here that latter approximation is well founded and a more complete analysis is given by [START_REF] Kudryavtsev | Impact of swell on the marine atmospheric boundary layer[END_REF]. For deep water waves, the free stream velocity above the waves, just outside of the boundary layer is u+(x, t) = -σa cos(kxσt), where a is the swell amplitude and σ = 2π/T is the radian frequency. The sub-surface velocity is u-(x, t) = σa cos(kx-σt) (figure 9). Due to the oscillations that propagate at the phase velocity C, the horizontal advection of any quantity X by the flow velocity u, given by u∂X/∂x, can be neglected compared to its rate of change in time ∂X/∂t since the former is a factor u/C smaller than the latter, which is typically less that 0.1 for the swells considered here. Defining u(x, z, t) = u(x, z, t)u-(x, t), where the brackets denote an average over flow realizations for a given wave phase. The horizontal momentum equation is thus approximated by,

∂ u ∂t = - 1 ρa ∂p ∂x - ∂u- ∂t + G (A1)
where G represents the divergence of the vertical viscous and turbulent fluxes of horizontal momentum,

G = ν ∂ 2 u ∂z 2 + ∂ u ′ w ′ ∂z . (A2)
Because the boundary layer thickness δ is small compared to the wavelength, the pressure gradient in the boundary layer is given by the pressure gradient above the boundary layer, in balance with the horizontal acceleration. This is another way to write Bernoulli's equation (e.g. [START_REF] Mei | Applied dynamics of ocean surface waves[END_REF]), -∂p/∂x/ρa = -σ 2 a sin(kxσt) = ∂u+/∂t. (A3) This yields

∂ u ∂t = 2 ∂u+ ∂t + G (A4)
with the boundary condition for z ≫ δ, u goes to 2u+(x, t).

The equation for the horizontal momentum is thus exactly identical to the one for the oscillatory boundary layer over a fixed bottom with wave of the same period but with an amplitude twice as large. In the viscous case, G = ν∂ 2 u/∂z 2 and one recovers, after some straightforward algebra, the known viscous result, i.e., for

z > ζ, u(x, z, t) = 2σa e -z + cos (kx -σt -z+) -cos (kx -σt) + O(ρa/ρw) (A5)
where z+ = (zζ)/ 2ν/σ, with the surface elevation ζ(x, t) = a cos(kxσt).

The dissipation rate of energy is given by the mean work of the viscous stresses,

β = Cgµv = ρaνu∂u/∂z ρwga 2 /2 . (A6) Equation (A5) gives, ∂ u ∂z = -2σa e -z + 2ν/σ [cos (kx -σt -z+) -sin (kx -σt -z+)] (A7)
This gives the low frequency asymptote to the viscous decay coefficient [START_REF] Dore | Some effects of the air-water interface on gravity waves[END_REF],

µv = -2 σ 2 gCg ρa ρw √ 2νσρa/ρw. (A8)
This result was previously obtained using a Lagrangian approach without all the above simplifying assumptions [START_REF] Weber | Effect of the air on the drift velocity of water waves[END_REF], their equation 7.1 in which our µ correspond to their 2α]. The full viscous result is obtained by also considering the water viscosity νw, which gives the O(ρa/ρw) correction for the motion in the air, and the classical dissipation term with a decay µvw = -4k 2 νw/Cg, which dominates for the short gravity waves. The total disspation rate is simply µ = µv + µvw. For a clean surface with νw = 3 × 10 -6 m 2 s -1 , ν = 1.4 × 10 -5 m 2 s -1 , and ρa/ρw = 0.0013, the two terms are equal for waves with a period T = 1.3 s and a wavelength of 2.6 m. The air viscosity dominates for all waves longer than this, which is typically the range covered by spectral wave models for sea state forecasting. For a comparison with fixed bottom boundary layers, the Reynolds number based on the orbital motion should be redefined with a doubled velocity and a doubled displacement, i.e. Re= 4u orb a orb /ν. For monochromatic waves a orb = a and u orb = aσ = 2πa/T . For random waves, investigations of the ocean bottom boundary layer suggest that the boundary layer properties are roughly equivalent to that of a monochromatic boundary layer defined by significant properties [START_REF] Traykovski | Geometry, migration, and evolution of wave orbital ripples at LEO-15[END_REF].

Although the wind was neglected here, it should influence the shear stresses when its vertical shear is of the order of the wave-induced shear. Taking a boundary layer thickness δ and wind friction velocity u⋆, and assuming a logarithmic wind profile, this should occur when u⋆/(κδ) exceeds 2u orb /δ, where κ is von Kármán's constant. This corresponds to, roughly, u⋆ > u orb . For swells with T < 15 s and Hss > 2 m (i.e. u orb,s > 0.4 m s -1 ), and winds less than 7 m s -1 (i.e. u⋆ < 0.2 m s -1 ), the wind effect on fe may be small and the previous analysis is likely valid. In general, however, the nonlinear interaction of the wave motion and wind should be considered, which requires an extension of existing theories for the distortion of the airflow to finite swell amplitudes.

Finally, the above result for the air viscosity is easily generated to any water depth D by dividing the free stream velocity u+(x, t) by tanh(kD), so that the dissipation rate is a factor 1/ tanh 2 (kD) larger in intermediate water depth.

Appendix B: Quantitative validation of H ss A classical analysis of SAR estimation errors is provided by a direct comparison of swell parameters, estimated from level 2 products, with buoy measurements at nearly the same place and time [START_REF] Holt | Tracking storm-generated waves in the northeast pacific ocean with ERS-1 synthetic aperture radar imagery and buoys[END_REF][START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF].

Previous validations were presented for the total wave height Hs [START_REF] Collard | Extraction of coastal ocean wave fields from SAR images[END_REF] or a truncated wave height Hs12 defined by chopping the spectrum at a fixed frequency cut-off of 1/12 Hz. For that parameter, [START_REF] Johnsen | ASAR wave mode processing -validation of reprocessing upgrade[END_REF] found a root mean square (RMS) difference of 0.5 m, when comparing SAR against buoy data, including a bias of 0.2 m. In the present study, we use Hss values obtained from both SAR and buoy spectra.

For each wave spectrum observed in the world ocean, swell partitions are extracted providing estimations of Hss, Tp, and θp. In practice, the L2 spectra are first smoothed over 3 direction bins (30 • sectors) and 3 wavenumber bins, in order to remove multiple peaks that actually correspond to the same swell system. The peaks are then detected and the energy associated to each peak is obtained by the usual inverted water-catchment procedure [START_REF] Gerling | Partitioning sequences and arrays of directional ocean wave spectra into component wave systems[END_REF]. The swell peak period is defined as the energy-weighted average around ±22% of the frequency with the maximum energy. Likewise the peak direction θp is defined as the energyweighter direction within 30 • of the peak direction. Table 1. Error statistics for swell partitions heights and peak periods derived from SAR wave mode data (after bias correction) against buoy-derived data, for subsets A, B and C of the co-located database. Subset A contains 2399 observations. Subset B is restricted to U 10SAR ≤ 8m/s and contains 1936 observations. Subset C is further restricted to SAR-buoy distances less than 100 km, and contains only 460 observations. RMSE stands for root mean square difference, while the NRMSE is the RMSE normalized by the root mean square observed value. The scatter index (S.I.) is equivalent to the NRMSE computed after bias removal. Finally, r is Pearson's linear correlation coefficient. A preliminary validation of Hss was performed by [START_REF] Collard | Coastal ocean wave retrieval from ASAR complex images[END_REF], using L2 processing applied to 4 by 4 km tiles from narrow swath images exactly located at buoy positions. That study found a 0.37 m r.m.s. error. This smaller error was obtained in spite of a 4 times smaller image area that should on the contrary produce larger errors due to statistical uncertainties. This suggests that a significant part of the "errors" in SAR validation studies are due to the distance between SAR and buoy observations. The swell height validation has been repeated here, using all buoy data from 2004 to 2008, located within 200 km and 1 hour of the SAR observation. These co-located data are made publically available as part of the XCOL project on the CERSAT ftp server, managed by Ifremer. Because we wished to avoid differences due to coastal sheltering and shallow water effects, we restricted our choice of buoys to distances from the coast and the 100 m depth contour larger than 100 km. As a result, most selected buoys are not directional, and we use partitions in frequency only. For the present validation, differently from other section of this paper, we thus define a swell partition as the region between two minima of the frequency spectrum. The corresponding energy Es gives the swell height Hss = 4 √ Es. The buoy swell height is then defined from the energy contained within the frequency band of the SAR partition. The peak period is then estimated as the period where the buoy spectrum is maximum. The database includes 15628 swell partitions observed by the SAR, with matched buoy swell partitions.

Hss

Many of these observations correspond to relatively short swells, for which the waves are poorly imaged. We have thus defined a subset of the database by imposing the following conditions. First the image normalized variance, linked to the contrast intensity and homogeneity, should be in the range 1.05 to 1.5, which limits the dataset to 6651 observations. This removes SAR data with non-wave features (slicks, ships ...) that would otherwise contaminate the wave spectra. Second, both the SAR and buoy peak periods are restricted to the 12 to 18 s range, which reduces the dataset to 4136 observations, and removes most of the problems related to the azimuthal cut-off. Third and last, the SARderived wind speed U10SAR is limited to range from 3 to 9 m s -1 in order to remove low winds with poorly contrasted SAR images and high winds which may still cause some important azimuthal cut-off and contamination of swell spectra by wind sea spectra. This gives subset A, with 2399 observations. The resulting heights are compared to buoy measurements in figure 10.

When the maximum wind is reduced to 8 m s -1 , giving subset B, the differences between SAR and buoy data is reduced, with further reductions when the maximum distance between SAR and buoy data is reduced from 200 to 100 km to give subset C (table 1).

Appendix C: Derivation and verification of the asymptotic swell energy without dissipation C1. Derivation

The swell energy Es is an integral of the local spectrum G over both frequencies f and arrival directions θ,

Es(α) = 2π 0 ∞ 0 G t, φ ′ , λ ′ , f, θ ′ dfdθ ′ . (C1)
Using eq. ( 1) this local integral, can also be written as an integral over the entire source area Ω. The spherical distance between any point P (ϕ, λ) in the source region and the observation point O(α, 0) is α ′ . The observed frequency that is due to this source point is

f = gt/(4πRα ′ ) = f0 α α ′ . (C2)
We may replace f by α ′ in eq. (C1),

Es(α) = f0 α 2π 0 α 2 G (t0, φ, λ, f, θ) α ′2 dθ ′ dα ′ . (C3)
For a circular uniform storm of radius r with isotropic spectra, as used in figure 11, Es(α) is given by the integral over α ′ weighted by the directional width of the spectrum ∆θ ′ . That width is given by the spherical trigonometry relationship

∆θ ′ = 2θ ′ max = 2 arccos cos ϕmax -cos α cos α ′ sin α sin α ′ , (C4)
with ϕmax = r/R. For general spectral distribution, we may transform the integration variables (α ′ , θ ′ ) which are the colatitude and longitude coordinates on the sphere with a pole at the observation point, to coordinates (ϕ, λ) with a pole in the center of the swell field at t0. The transformation Jacobian is simply given by the equality of the elementary area on the unit radius sphere dA = sin ϕdφdλ = |cos φ| dφdλ = sin α ′ dα ′ dθ ′ . We thus have

Es(α) = f0 α sin α Ω α 2 sin α α ′2 sin α ′ G (t0, φ, λ, f, θ) |cos φ| dφdλ (C5) = f0 α sin α Ω G (t0, φ, λ, f0, θ0) dA × 1 + O ∆α α + Ω [G (t0, φ, λ, f, θ) -G (t0, φ, λ, f0, θ)] dA + Ω [G (t0, φ, λ, f0, θ) -G (t0, φ, λ, f0, θ0)] dA . (C6)
where θ is the direction of the great circle at the generation point that goes through that point and the observation point. θ is thus a function of φ, λ, α and θ0. ∆α is maximum value of |α ′ -α|, i.e. the radius of the source region divided by the Earth radius.

For continuous spectra, the second integral on the right hand side of eq. (C1) goes to zero as α goes to infinity (which, on the sphere is limited by π) since the part of the source spectra that contribute to the observations shrink to a smaller and smaller neighborhood around f0 and θ0. The observed frequencies f are limited by

|f -f0| f0 ≤ ∆α α (C7)
This is enough to guarantee that this second integral also contributes a deviation ε2 from the asymptote, limited by

ε2 ≤ A ∆α α f0 max ∂G ∂f (C8)
where the maximum is taken over all the contributing components.

Similarly, the outgoing directions θ received at the observation point are also limited to a narrow window as α increases, giving another deviation term ε3. Using the sine formula in the triangle OP N , sin θ/ sin(πα) = sin λ/ sin α ′ . Thus, in the far field of the storm and its antipode, θ is close to θ0 = λ. Thus sin θsin θ0 = sin λ (sin αsin α ′ ) / sin α ′ , which is less than ∆α/ sin α ′ , and therefore |θ -θ0| is less than arcsin (∆α/ sin α ′ ). If one does not get too close to the storm or its antipode (say, ∆α < α < π -2∆α) then we can give an upper bound of 1/ sin α ′ as a function of α and obtain ε3 ≤ A arcsin 2 ∆α sin α max ∂G ∂θ .

(C9)

The deviation from the asymptote due to ε3 is thus of the order of ∆α/sin α and may increase close to the antipode, if waves from a wide range of directions can reach that point. In practice this does not happen since continents and island chains block most of the arrival directions at the antipode, leaving only a small window of possible arrival directions [e.g. [START_REF] Munk | Directional recording of swell from distant storms[END_REF]. Directional wave spectra in active generation areas are generally relatively broad with ∂G/∂θ/G typically less than 2 for directions within 30 • of the main wave direction. On the contrary, typical storm spectra can give f ∂G/∂f /G as large as 10 for frequencies within 30% of the peak frequency. We thus expect the deviation ε2 to dominate over ε3.

C2. Verification

For a useful comparison with observations, the asymptotic swell height evolution should be approached on a scale smaller than the ocean basin scale. This is easily tested for storms with spatially uniform spectra over a radius r, by evaluating the integral in eq. (C3). We chose a center frequency f0 and consider the swell energy at a distance Rα + ∆x and a time t(α) such that Rα = gt/(4πf0). ∆x is thus an error relative to the theoretical position of the wave group of frequency f0. The relative difference of Es(α) and its asymptotic evolution 1/(α sin α) depends only on the distance from source (1000 km) Relative deviation from 1/[α sin(α)] (%) γ=1.7, f 0 /f p = 1.0 γ=1.7, f 0 /f p = 1.2 γ=1.0 ,f 0 /f p = 1.0 Source size: 1000 km . Convergence of the swell energy Es integral (C1) towards the asymptote 1/(α sin α), as a function of spectrum width for a storm diameter of 1000 km. The result is independent of the choice of fp. In practice the calculations were made for fp = 0.07 Hz (Tp = 13 s). The three lines for each case correspond to position errors ∆x of -200, 0 and 200 km relative to the great circle trajectory. For all cases considered here the deviation is less than 20% beyond 4000 km from the storm center. spectral shape in the storm, the relative frequency f0/fp where fp is the peak frequency, the distance α, the storm size α/(r/R), and the position error ∆x. Results are shown in figure 11 for isotropic directional spectra, in which case ε3 = 0. Although an isotropic spectrum is not realistic at all, it allows for simple calculations. As discussed below, the effect of directional spreading is expected to be less important than the shape of the frequency spectrum.

As indicated by eq. (C8), the contribution of both the spectral shape and f0/fp comes through the maximum relative variation of F in the frequency interval that contribute to Es. Here we take the spectra in the storm to have a JONSWAP shape [START_REF] Hasselmann | Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project[END_REF], that we adjust by varying the peak enhancement factor γ. A relatively broad [START_REF] Pierson | A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii[END_REF] spectrum is obtained with γ = 1. If we chose f0 = fp, this spectrum will give smaller deviations of Es from the asymptote (solid lines in figure 11), than narrower spectra with larger values of γ. [START_REF] Young | Directional spectra of hurricane wind waves[END_REF] showed that wave spectra in Hurricanes generally fall in between these two categories, with a typical value γJ = 1.7. Larger deviations from the asymptote are obtained for f0 < fp, since the forward face of the spectrum is very steep, while smaller errors are obtained for f0 > fp due to the more gentle decrease of F towards the high frequencies.

Similarly, large deviations are produced if the observations are made in a direction far from the peak generation direction in cases when the directional spectrum is narrow. For observation directions 30 • from the peak direction, and spectra with a cos 4 directional distribution, the deviations are still dominated by the dispersive term ε2, as expected.

The other important factor is the distance α relative to the storm size ϕmax = r/R. A faster convergence is obtained for smaller storms. If observations do not correspond exactly to a theoretical propagation at a group speed Cg but are within a distance ∆x of the theoretical position, the values of Es will also be affected in a way that depends on the spectral shape. In the absence of energy gains or losses, and for realistic storm sizes and spectral shapes, the deviation of observations from the asymptote should be less than 20% for x > 4000 km when ∆x < 200 km is enforced.

Figure 3 .Figure 4 .

 34 Figure3. Snapshot of the "fireworks" animation given in the auxiliary material, for July 17 2004 at 0:00 UTC. Each of the 1071 colored dots represent one observed swell partition, within 6 days of its observation, displaced along a great circle with the group speed corresponding to the detected peak period in the direction of the detected peak direction. Only swells with tracks that passes within 1000 km of the storm center (red disk) have been retained.

Figure 7 .

 7 Figure7. Geometry of the "fossil' swell field distribution at time t0 = 0 (shaded area) and observation conditions. Any point P of colatitude ϕ and longitude λ inside of the storm, generates waves that are observed at point O. At time t the observed waves that come from P have a well defined frequency given by eq. (C2), function of the spherical distance α ′ between P and O, and a well defined direction θ at P , relative to the North, which gives a direction πθ ′ at O. In the triangle OP S the angles λ, θ ′ and πθ are related to the distances α ′ ,ϕ and α by the usual spherical trigonometry relationships, e.g. eq. (C4).

Figure 8 .

 8 Figure 8. (a) Location of SAR observations with a 15 s peak period swell system corresponding to the 12 February source, with outgoing directions of 74 to 90 • . The same swell was also observed at all buoys from 46075 off Western Alaska, to 51001 in Hawaii. The dash-dotted line represent great circles leaving the storm source with directions 42, 59, 74, 90 and 106 • . (b) Observed swell wave height as a function of distance. The solid lines represent theoretical decays using no dissipation (blue), or a the fitted linear dissipation (green), for swells observed in February 2007. Circled dots are the observations used in the fitting procedure. Error bars show one standard deviation of the expected error on each SAR measurement.

Figure 9 .

 9 Figure9. Boundary layer over waves in the absence of wind. Because of the larger inertia of the water compared to the air, most of the adjustment from the sub-surface velocity to the free stream velocity in the air occurs on the air-side of the surface.

Figure 10 .

 10 Figure10. ASAR-derived swell partition heights versus buoy swell partition heights after bias correction using eq. (3), for subset A. The solid line joins the median values from SAR observations in each 0.1 m class of buoymeasured height.

  Figure11. Convergence of the swell energy Es integral (C1) towards the asymptote 1/(α sin α), as a function of spectrum width for a storm diameter of 1000 km. The result is independent of the choice of fp. In practice the calculations were made for fp = 0.07 Hz (Tp = 13 s). The three lines for each case correspond to position errors ∆x of -200, 0 and 200 km relative to the great circle trajectory. For all cases considered here the deviation is less than 20% beyond 4000 km from the storm center.

Acknowledgments. SAR data was provided by the European Space Agency (ESA) and buoy data was kindly provided by the U.S. National Data Buoy Center, and Marine Environmental Data Service of Canada. ESA funded the XCOL co-located database and the initial work on the virtual buoy concept. The swell decay analysis was funded by the French Navy as part of the EPEL program. This work is a contribution to the ANR-funded project HEXECO and DGA-funded project ECORS.


Francesco Pappadà
email: francesco.pappada@psemail.eu.

The Dynamics of Tax Compliance

Keywords: E26, E62, H26 tax compliance, volatility, fiscal policy

The literature on tax compliance has focused on its level but little is known about its dynamics. This paper shows that fluctuations in tax compliance are driven by changes in the state of the economy and the response of tax compliance to them. Tax compliance is markedly volatile and there are large differences in such volatility across countries. A large fraction of these differences (about 70%) is explained by different responses of tax compliance to tax changes and output fluctuations.

Introduction

The literature on tax evasion has clearly established the determinants of tax compliance. In a recent literature review, [START_REF] Slemrod | Tax Compliance and Enforcement[END_REF] summarizes three main factors: i) the level of tax base and rates, ii) the remittance rules, namely the entity which reports tax liability, and iii) the enforcement rules. These factors determine the degree of tax evasion, which is quite heterogeneous across countries.1 While the level of tax compliance has received much attention, far less attention has been paid to its dynamics. However, the dynamics of tax compliance should be relevant for policymakers. Fluctuations in tax compliance may indeed translate in volatile tax revenues, affecting the ability of governments to finance public expenditures, and constraining the conduct of fiscal policy.

This paper provides novel empirical evidence on the volatility of tax compliance and its determinants. The main results are twofold. First, I show that tax compliance is very volatile, and there is large cross-country heterogeneity in such volatility. Some countries experience large fluctuations in tax compliance, whereas other have a much more stable tax compliance. Interestingly, the countries where the tax compliance is more (resp. less) volatile are not necessarily the ones where the informal sector is large (resp. small).

Second, I study the determinants behind the cross-country differences in the volatility of tax compliance. There are two potential forces behind the fluctuations in tax compliance: i) changes in the state of the economy (e.g. tax reforms or output fluctuations) and ii) the elasticity of tax compliance to these changes in the state of the economy. I decompose the contribution of these two factors to the volatility of tax compliance and show that the second factor is predominant. Tax compliance tends to be more volatile in countries where it is more sensitive to the state of the economy. More precisely, the country-specific elasticities of tax compliance to tax policy and output fluctuations explain about 70% of the overall cross-country differences in the volatility of tax compliance, whereas tax policy and output fluctuations explain about 20% of it.

The dynamics of tax compliance

A measure of tax compliance I construct a measure of tax compliance for 35 countries by comparing (i) actual Value-Added Taxes (VAT) revenues, as reported by government sources, and (ii) consumption across good categories subject to different VAT rates, as reported in household consumption surveys. I rely on VAT compliance as measure of tax compliance for two reasons. First, it allows the construction of a credible measure of tax compliance across countries and over time. Second, VAT is more frequently adjusted than other tax instruments. 2My measure of tax compliance compares tax receipts to expected receipts as predicted by tax rates and actual expenditures. Letting tr itc , τ itc and c itc denote VAT revenues, VAT rate and consumption of good i in year t and country c, the measure of VAT compliance is defined as:

γ tc = i tr itc i (τ itc c itc ) .
The gap between tax revenues and expected tax revenues, as captured by the distance between γ tc and 1, reflects imperfect tax enforcement from tax authorities. Such measure cannot however shed light on the nature of tax leakages, whether they come from informal exemptions, corruption of tax authorities or non-cooperative tax evasion from agents. The measure accounts for possible changes in consumption patterns c itc as a response to differential tax rates across goods. Fluctuations in tax compliance can only arise from changes in tax compliance within good categories.

I use distinct data sources for total tax revenues, and for reported consumption of 48 disaggregated good categories between 1995 and 2013 in 35 countries.3 I observe total VAT receipts, tr tc = i tr itc , in national accounts. I use annual household expenditure surveys to create actual consumption, c itc , in each sub-category of good. 4 The information in household surveys comes from the purchaser side thereby alleviating under-reporting of undeclared transactions. I also extract from the European Commission documentation and national sources the different tax rates and I reference good categories that are subject to these rates for each country/year. This allows to account for adjustments in the composition of exempted categories. 5

The volatility of tax compliance

This section provides descriptive statistics on the dynamics of tax compliance. Interestingly, economies that are close in terms of economic development and tax enforcement display stark differences in terms of volatility of tax compliance (see for instance Ireland vs. United Kingdom, or Colombia vs. Chile). Moreover, these differences do not reflect differences in the size of the informal sector. In Figure 2, I compare the volatility of tax compliance to a more standard measure of imperfect tax enforcement, i.e., the size of the informal sector as computed in [START_REF] Schneider | The shadow economy: An international survey[END_REF]. While there is a positive correlation between the two measures, countries with similar incidence of the informal sector markedly differ along fluctuations in tax compliance. What could explain these differences?

ance and the corrections to clean for i) within-year VAT reforms, ii) VAT reforms for sub-categories of goods within a 2-digit classification, and iii) reforms that modify the tax environment without modifying tax rates. 6 In the sample, tax compliance is 0.86 on average, which is arguably quite high-maybe due to the over-representation of rich economies or the built-in incentives to report VAT compared to other taxes.

Explaining the cross-country differences in the volatility of tax compliance

There are two potential determinants of the fluctuations in tax compliance: the changes in the state of the economy Z c -for instance tax rates and output -and the elasticity of tax compliance to such changes ε Zc . This implies that the volatility of tax compliance in country c may be written as var(γ c ) = ε Zc • var(Z c ).

In this section, I decompose the contribution of these two factors to the volatility of tax compliance and then analyze to what extent they explain the differences across countries. I first explore the role of the volatility of tax rates and output. Panel (a) and (b) of Figure 3 show that countries with larger tax changes and output fluctuations tend to display a larger volatility of tax compliance. However, there are substantial differences in the volatility of tax compliance for given fluctuations in output and tax rates. I then quantify the elasticity of tax compliance to tax rates by estimating the following baseline specification in difference:

∆ ln γ tc = ε τ ∆ ln τ tc + ε y ∆ ln y tc + βX tc + δ t + µ c + e tc , (1)
where I consider a country c in year t and denote ∆ ln γ tc as the annual percentage change in VAT compliance; ∆ ln τ tc is the main measure of fiscal policy, i.e., the annual percentage change in the standard VAT rate; and ∆ ln y tc is the annual percentage change in output. The vector X includes the following set of time- varying controls: annual changes in the sectoral decomposition of economic activity across 1-digit sectors and in trade (the ratio of exports and imports over GDP). These controls partly clean for confounding factors co-moving with tax compliance and fiscal policy. Finally, the inclusion of country fixed effects, µ c , conditions the analysis on country-specific trends and δ t are year fixed effects capturing aggregate dynamics in tax compliance. Notes: Robust standard errors are reported between parentheses. The unit of observation is a country in a given year. All specifications include year-and country-fixed effects; we sequentially add the following controls: annual changes in the sectoral decomposition of economic activity across 1-digit sectors in columns 2-3; and annual changes in trade (the ratio of exports and imports over GDP) in column 3. The dependent variable is the annual percentage change in VAT compliance. Tax rate is the annual percentage change in the standard VAT rate; GDP growth is the annual percentage change in output (see specification 1).

Table 1 uncovers an important, yet overlooked, characteristic of tax compliance: it fluctuates markedly with fiscal policy and with the business cycle. I find that the average elasticity of tax compliance to the tax rate is about -0.36, and, as shown in columns 2 and 3, the correlation is robust to the addition of controls (sectoral composition in column 2 and trade in column 3). In parallel, an increase of one percent in output is associated with an increase of 20 percent in tax compliance.

The heterogeneous dynamics of tax compliance across countries may be explained by different fluctuations in the state of the economy and different responses of tax compliance to them. What is the contribution of each factor? In order to answer this question, I introduce the following decomposition:

var( γc ) = var(ε τ c • τ ) + var(ε yc • ȳ) + var(ε τ • τc ) + var(ε y • ỹc ) + u, (2)
where I define xc = var(∆ ln x tc ) and x = avg( xc ) for each variable x tc ; ε τ c (resp. ε yc ) is the country-specific elasticity of tax compliance to changes in tax rates (resp. output fluctuations) and u is the residual. As reported in Table 2, the country-specific responses of tax compliance to tax policy and output fluctuations account for the largest fraction (around 70%) of the differences in the volatility of tax compliance across countries, whereas the differences in output fluctuations and tax policy account for a smaller fraction (around 20%). Notes: The table reports the share of each component of the cross-country differences in the volatility of tax compliance.

Final remarks

This paper shows that the response of tax compliance to the state of the economy explains a large fraction (around 70%) of the differences in the volatility of tax compliance across countries. The determinants of the elasticity of tax compliance to the state of the economy remain unclear, and further research should investigate them. Among the possible explanations, there could be the sectoral composition of tax reforms or the implementation of these reforms over the business cycle. The dynamics of tax compliance have important implications on the conduct of fiscal policy. As shown in Pappadà and Zylberberg (2021), the fluctuations in tax compliance strongly distort the relationship between fiscal policy and default risk as credit markets internalize the cost of future distortions in tax compliance.

ONLINE APPENDIX

We provide here robustness and complements to the analysis presented in Section 2.

The dynamics of tax compliance In Table A1, we replicate the estimation of specification (1) and add-as further controls-the interactions between Tax rate/GDP growth and: a dummy for being classified as an industrial country in [START_REF] Vegh | How is tax policy conducted over the business cycle?[END_REF] (column 1); a dummy for having above-median tax compliance (column 2); and a dummy for having an above-median informal sector [START_REF] Schneider | The shadow economy: An international survey[END_REF], see column 3). The rationale for doing so is that the heterogeneous response to fiscal policy and output may reflect various dimensions of heterogeneity across countries. I am generally quite agnostic as to why tax compliance is more or less responsive across countries of our sample. This exercise however shows that the volatility of tax compliance strongly predicts the size of the tax compliance response to fiscal policy, on top of other, important dimensions of heterogeneity, for instance, the level of development or the average size of the informal sector.

Reverse causality

The estimate ε γ may be contaminated by omitted variation and reverse causality. Fluctuations in a few sectors with higher tax evasion may induce the government to adjust the tax rate. A government observing a negative trend in tax compliance may decide to increase the tax rate to mitigate a revenue shortfall. In Table A2, we shed some light on the latter effect by regressing the change in the standard VAT rate on the lagged change in tax compliance. We find that the lagged tax compliance cannot predict future changes in fiscal policy.

  Figure 1 plots the within-country standard deviation of tax compliance between 1995 and 2013 in my sample of 35 high-and middle-income countries. 6 This Figure shows that there is a large heterogeneity in the volatility of tax compliance across countries.

Figure 1 .

 1 Figure 1. Volatility of tax compliance.

Figure 2 .

 2 Figure 2. Volatility of tax compliance and the size of the informal sector.

Figure 3 .

 3 Figure 3. Volatility of tax compliance, tax rate and output.

Table 1 .

 1 Elasticity of tax compliance.

	Tax compliance	(1)	(2)	(3)
	Tax rate	-.357	-.361	-.358
		(.066)	(.066)	(.067)
	GDP growth	.201	.205	.213
		(.135)	(.139)	(.139)
	Controlling for sectoral composition	No	Yes	Yes
	Controlling for trade	No	No	Yes
	Observations	527	527	527

Table 2 .

 2 Decomposition of cross-country differences in the volatility of tax compliance.

	Response to tax policy	var(ε τ c • τ )	.34
	Response to output	var(ε yc • ȳ)	.39
	Fluctuations in tax policy	var(ε τ • τc )	.17
	Fluctuations in output	var(ε y • ỹc )	.04
	Residual	u	.06

See[START_REF] Schneider | The shadow economy: An international survey[END_REF] for cross-country estimates of the informal sector.

The corporate tax rate and the income tax rate are typically adjusted once every 8-12 years; adjustments are then non-negligible. By contrast, the effective VAT rate is more often modified, reflecting frequent changes in the classification of goods across the different brackets. Furthermore, the systematic analysis of fiscal consolidations classified by[START_REF] Alesina | Is it the" How" or the" When" that Matters in Fiscal Adjustments?[END_REF] shows that almost half of the tax surplus is generated through reforms of indirect taxation.

I rely onOECD and Eurostat and their harmonized

COICOP (Classification of IndividualConsumption by Purpose) sub-categories of goods.4 These household surveys are standardized across countries, and they follow similar methodology (i.e., sampling and questionnaire). The aggregate consumption constructed from these surveys strongly correlates with total output, but there exists (standard) measurement error (see[START_REF] Aguiar | Has consumption inequality mirrored income inequality?[END_REF][START_REF] Kolsrud | The value of registry data for consumption analysis: an application to health shocks[END_REF], for a correction method and a comparison with registry data). This error will translate into a measurement error in the denominator of the expression for γ tc .5 See[START_REF] Pappadà | Austerity and tax compliance[END_REF] for additional details about the measure of VAT compli-

Notes: Robust standard errors are reported between parentheses. The unit of observation is a country in a given year. All specifications include year-, country-fixed effects and the following controls: annual changes in the sectoral decomposition of economic activity across 1-digit sectors and annual changes in trade (the ratio of exports and imports over GDP). The dependent variable is the annual percentage change in VAT compliance. Tax rate is the annual percentage change in the standard VAT rate; GDP growth is the annual percentage change in output (see specification 1). We augment the baseline specification (1) by adding interactions between Tax rate/GDP growth and a dummy for being a high-response country (see Section ?? for the definition of the High-response dummy).

We add further interactions between Tax rate/GDP growth and: a dummy for being classified as an industrial country in [START_REF] Vegh | How is tax policy conducted over the business cycle?[END_REF] † (column 1); a dummy for having above-median tax compliance (column 2); a dummy for having an above-median informal sector [START_REF] Schneider | The shadow economy: An international survey[END_REF], see column 3). †: In [START_REF] Vegh | How is tax policy conducted over the business cycle?[END_REF]'s classification, industrial countries are Australia, Austria, Belgium, Canada, Denmark, Finland, France, Germany, Greece, Ireland, Italy, Japan, Luxembourg, the Netherlands, New Zealand, Norway, Portugal, Spain, Sweden, Switzerland, United Kingdom; non-industrial countries are Bulgaria, Chile, Colombia, Cyprus, Czech Republic, Estonia, Hungary, Korea, Malta, Poland, Romania, Russia, South Africa. Four countries of our sample are not classified by [START_REF] Vegh | How is tax policy conducted over the business cycle?[END_REF] and six countries are not covered by [START_REF] Schneider | The shadow economy: An international survey[END_REF], hence the lower number of observations in columns 1 and 3 compared to our baseline. Notes: Robust standard errors are reported between parentheses. The unit of observation is a country in a given year. All specifications include year-and country-fixed effects; we add the following controls: annual changes in the sectoral decomposition of economic activity across 1-digit sectors in column 2-3; and then annual changes in trade (the ratio of exports and imports over GDP) in column 3. The dependent variable is the annual percentage change in the standard VAT rate. Tax compliance (lag) is the annual percentage change in tax compliance in t -1; GDP growth (lag) is the annual percentage change in output in t -1.