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

Given these SAR-derived estimates of $\theta_{p}$ and $T_{p}$, linear dispersion relationship and the principle of geometrical optics can then be exploited to predict arrival times and locations of the swell.

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 $gT_{p}/(4\pi)$, starting off in direction $\theta_{p}$ from the observation point. From any observation time $t_{0}$, these great circles may cut through the interrogation window from times $t_{1}$ to $t_{2}$ (thick solid lines in figure 1, with different colors for different partitions). As the maximum value of $|t_{2}-t_{0}|$ and $|t_{1}-t_{0}|$ is increased from 6 hours to 6 days, the time-evolution 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 $t_{1}$ and $t_{2}$. 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. Munk et al., 1963; Gjevik et al., 1988). 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. Tolman, 2002). 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 (Johnsen and Collard, 2004). 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 (Johnsen and Collard, 2004; Johnsen et al., 2006). 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 (Johnsen and Collard, 2004).

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 (Munk and Traylor, 1947). 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.

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.

From the relatively sparse and track-based initial satellite observation sampling, the swell persistency can then be used to capture ”fireworks” patterns exploding from the few 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. Holt et al., 1998). For this reason, these animations have been produced routinely every day since August 2007, for the Pacific, Indian, and Atlantic oceans (see http://www.boost-technologies.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 Heimbach and Hasselmann (2000), provides a global view of swell fields in both space and time, extending the coverage of similar techniques based on buoy data (Hanson and Phillips, 2001). 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. Munk et al., 1963). 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 Darbyshire (1958); Munk et al. (1963); Snodgrass et al. (1966). 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.

Using measurements with a limited or no directional resolution, Snodgrass et al. (1966) 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(x_{2},f_{0})/F(x_{1},f_{0})$, 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 Snodgrass et al. (1966).

Rigorously, their method is inexact because the recording stations $x_{1}$ and $x_{2}$ measure wave groups that had neither exactly same propagation directions nor the same position when they were near $x_{0}$. Yet, because the wave field in the neighborhood of $x_{0}$ 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 $\theta_{3}$ to $\theta_{4}$ or $\theta_{5}$ to $\theta_{6}$, $G$ does not change so much, i.e. in figure 6, $G(t_{0},\phi_{D},\lambda_{D},f,\theta_{6})\simeq G(t_{0},\phi_{D},\lambda_{D},f,\theta_{5})$. On the sphere the application of eqs. (1) and (5) yields

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 $\theta$ over the range $\theta_{5}$ to $\theta_{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(\alpha)$, the ratios of spectral densities $F$ at $x_{1}$ and $x_{2}$, as used by Snodgrass et al. (1966), can be compared to the conservative equations (7) and used to diagnose the dissipation of swell energy. This does not apply if $x_{1}$ 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 Snodgrass et al. (1966) 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.

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 $E_{s}$ integrated over a swell partition. We thus need for $E_{s}$ the equivalent of eq. (7) for $F$.

For simplicity, and without loss of generality, we take the source storm centered at time $t_{0}=0$ on the pole $S$ defined by a colatitude $\varphi=\phi-\pi/2=0$, so that the distance from the storm center is $r=R\varphi$. We consider the swell energy $E_{s}$ observed at a position defined by the spherical distance $\alpha$ 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\Delta_{\alpha}$, where $\Delta_{\alpha}$ is the maximum value taken by $\varphi$ (figure 7). In all the following derivations, we have chosen a fixed frequency $f_{0}$ and we follow a wave group of that frequency. The time of observation $t$ is thus related to $\alpha$ by $t=R\alpha/Cg(f_{0})$, so that the variable $t$ will be omitted.

In appendix C, we prove and verify that, in the absence of dissipation the swell energy $E_{s}$ decreases like $1/(\alpha\sin\alpha)$ for large values of $\alpha$, 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.

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 south-north 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 $H_{s}$ 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 $H_{s}>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 $\theta_{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 $\theta_{0}$ with a step of $2^{\circ}$. In our case, when varying $\theta_{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, $E_{s}$ decreases asymptotically as $1/[\alpha\sin(\alpha)]$ 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 $H_{ss}$ 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 $E_{s}$ is only a function of $\alpha$, and we define the dissipation rate

Positive values of $\mu$ correspond to losses of wave energy (figure 8). We then fit an analytical function $H_{ss}(\alpha)$ to the data, defined by a constant $\mu$ and $H_{ss}(\alpha=\pi/5)$, i.e. the swell height at 4000 km from the source. Here the couple $H_{ss}(\alpha=\pi/5)=4.4$ m and $\mu=3.7\times 10^{-7}$ m^-1 gives the least square difference between the decay with constant $\mu$ and the observed swell decays. Further, the uncertainty of that dissipation rate may be estimated from the known uncertainty of the SAR measurement of $H_{s}$, given by eq. (3)-(4). A more simple error model, with larger errors based on the $H_{s12}$ analysis by Johnsen et al. (2006), 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 $\mu$, that would correspond to one standard deviation for a Gaussian distribution, we find that $3.1\times 10^{-7}<\mu<4.0\times 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 $\mu$ are more than twice larger than reported by Snodgrass et al. (1966) 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^{\circ}$. 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 (Rascle et al., 2008). 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. Rascle et al., 2008).