A Multimodal Wave Spectrum – Based Approach for Statistical Downscaling of Local Wave Climate

Characterization of wave climate by bulk wave parameters is insufficient for many coastal studies, including those focused on assessing coastal hazards and long-term wave climate influences on coastal evolution. This issue is particularly relevant for studies using statistical downscaling of atmospheric fields to local wave conditions, which are often multimodal in large ocean basins (e.g., PacificOcean). Swell may be generated in vastly different wave generation regions, yielding complex wave spectra that are inadequately represented by a single set of bulk wave parameters. Furthermore, the relationship between atmospheric systems and local wave conditions is complicated by variations in arrival time ofwave groups fromdifferent parts of the basin.Here, this study addresses these two challenges by improving upon the spatiotemporal definition of the atmospheric predictor used in the statistical downscaling of local wave climate. The improved methodology separates the local wave spectrum into ‘‘wave families,’’ defined by spectral peaks and discrete generation regions, and relates atmospheric conditions in distant regions of the ocean basin to local wave conditions by incorporating travel times computed from effective energy flux across the ocean basin. When applied to locations with multimodal wave spectra, including SouthernCalifornia and Trujillo, Peru, the newmethodology improves the ability of the statistical model to project significant wave height, peak period, and direction for each wave family, retaining more information from the full wave spectrum. This work is the base of statistical downscaling byweather types, which has recently been applied to coastal flooding and morphodynamic applications.


Introduction
At a given time, the wave state of the ocean surface is a composite of wind seas and swell.Wind seas are generated by and strongly coupled with local winds, whereas swell is generated remotely and might have propagated over large distances.Though multiple definitions exist, swell can be distinguished from wind seas when the wave phase speed exceeds the overlaying wind speed by 20% (Semedo et al. 2011).Swells and seas are functions of both the intensity and frequency of atmospheric systems (Young et al. 2011).The wave state of the ocean surface represents a multitude of wind seas and swell trains, which each have a particular set of bulk wave statistics (significant wave height H s , peak wave period T p , mean wave direction D m , and directional spreading s).Often, the wave spectrum exhibits multiple wave energy peaks (e.g., California deep-water wave climate; Crosby et al. 2016), with contributions of energy generated locally and energy propagated from distant regions.However, using the full wave spectrum from numerical wave models or wave buoys produces large volumes of data that pose management and analysis challenges.As a result, many studies examine bulk wave parameters of the largest energy peak of a spectrum while virtually eliminating secondary spectral peaks, which leads to incomplete and potentially misleading results.Recent advances allow for simple statistical representation of multimodal wave spectra but have yet to be widely applied (Portilla-Yandún et al. 2015).
Statistical downscaling methods for projecting local wave climate exploit the relationship between wave conditions and the magnitude and frequency of atmospheric systems.Statistical downscaling defines a linear (Wang et al. 2010) or nonlinear (this study) relationship between atmospheric variables, such as sea level pressure or wind speed, and wave parameters.These methods can be compared to dynamical downscaling, where a numerical wave model is forced by spatiotemporally varying winds (Erikson et al. 2015).Though statistical downscaling methods have been applied extensively with success in the Atlantic Ocean and Mediterranean Sea (Wang et al. 2012;Laugel et al. 2014;Camus et al. 2014a;Rueda et al. 2016a), the ability to use these methods in larger ocean basins (e.g., Pacific Ocean) is still being explored (Graham et al. 2013;Camus et al. 2014a;Rueda et al. 2017).Statistical downscaling relies on the quality of the definitions of the predictor (e.g., regional sea level pressure) and predictand (e.g., local wave conditions).One commonly used method for defining these fields is presented in Camus et al. (2014a).In that work, daily regional sea level pressure (SLP) fields averaged over an optimal number of days defined the predictor for daily total bulk wave parameters at a particular location (predictand).The spatial range of the predictor encompassed the primary areas of wave generation for the location, identified using Evaluating the Source and Travel Time of the Wave Energy Reaching a Local Area (ESTELA; Perez et al. 2014).The temporal range of the predictor was the average travel time of wave energy generated inside the wave generation region to the location, also identified using ESTELA.
In smaller ocean basins, such as the North Atlantic Ocean and Mediterranean Sea, the region of wave generation is relatively small.Thus, the travel times of wave energy generated within the basin (i.e., far from the location versus close to the location) differ by only a few days, usually less than 5 days.Additionally, wave spectra are often unimodal, with swell arriving from one discrete generation region.In these cases, the definition of the predictor and predictand following Camus et al. (2014a) is successful in reproducing historical or projecting future wave climates.However, in large ocean basins, such as the Pacific Ocean, the spatial and temporal definitions of the predictor and total bulk parameter definition of the predictand yield less successful statistical downscaling results.The reasons are twofold.First, waves may be generated in multiple discrete generation regions, yielding mixed sea states of local seas and multiple swell trains.Total bulk parameters and a single wave generation region do not adequately represent the wave climate.Second, travel times of wave energy generated in different parts of the region may differ by a few weeks.As a result, sea level pressure fields averaged over a time period do not physically relate to waves arriving on a particular day.In this work, we seek to improve the definitions of the predictor and predictand of Camus et al. (2014a) by 1) introducing ''wave families'' to better model multimodal spectra and multiple wave generation regions and 2) using isochrones of travel time to account for the vastness of wave generation regions.
The improved methodology is presented by defining the predictor and predictand for a location offshore of Southern California.Local wave conditions along Pacific coastlines are influenced by waves generated and propagated over very large distances, often the entire extent of the Pacific basin (Adams et al. 2008).Waves generated by distant storms in the central and western North Pacific and those generated in the South Pacific and Southern Ocean contribute to the bimodal wave spectrum of California (Crosby et al. 2016), making it difficult to form statistical relationships between atmospheric conditions and waves using existing methodologies for the reasons identified above (Espejo et al. 2014).The success of statistical downscaling for use in large ocean basins would allow for a more rapid proliferation of regional projections using numerous global climate models and forcing scenarios, creating an ensemble that would better estimate future wave conditions (Perez et al. 2015).Additionally, recent advances in statistical downscaling allow for simulation and projection of extreme bulk parameters, such as daily maximum significant wave height (Rueda et al. 2016a), coastal flooding (Rueda et al. 2016b), and coastal morphodynamics (Antolinez et al. 2016).
In this work, an improved methodology for defining the optimal predictor for statistical downscaling of multivariate (H s , T p , and D m ), multimodal wave climatology is presented, built on the methodology defined in Camus et al. (2014a).Wind speed, SLP, and fetch are common variables used in predictor definitions for waves (Wang et al. 2010).Though near-surface winds are the drivers of wave growth, SLP and squared gradients of SLP fields contain information of both wind direction and magnitude, with spacing of isobars describing wind speed and orientation of isobars describing wind direction (Espejo et al. 2014).The methodology presented here utilizes the relationship between SLP and wave parameters while improving the ability to represent multimodal wave climates through wave families and incorporating physical intuition into the temporal and spatial definition of the predictor by using isochrones of travel time.The paper is structured as follows: Data and methods are presented in sections 2-5 through an application to Southern California, followed by an additional application to Trujillo, Peru, in section 6.This work concludes with brief comments in section 7.

Methodology and data
Statistical downscaling defines a relationship between a predictor, here daily sea level pressure and the spatial gradients of sea level pressure, and a predictand, here daily multivariate wave parameters for multiple wave families at a particular location.This work seeks to improve the ability to statistically model multimodal wave climate.The methodology (Fig. 1) can be separated into three broad steps: step A, parameterization of spectral data; step B, spatiotemporal definition of the predictor; and, step C, a multiple multivariate linear regression model to assess the relationship between the variables.These steps are subdivided as follows.
A1: Identify wave generation areas using ESTELA (Perez et al. 2014).A2: Partition the wave spectral data into sea and swell components.B1: Define the temporal range of the predictor using isochrones.B2: Identify wave families using wave generation areas and sea and swell partitions.B3: Construct a daily predictor of SLP for each wave family using the isochrones to create a temporal relationship between atmospheric fields and waves.C1: Use multiple multivariate linear regression analysis to assess the skill of the predictor to reproduce the predictand (daily wave conditions for each family).
The methodology is applied to a location in deep water off the coast of Southern California, United States (338N, 1208W), as a demonstrative tool.
The historical SLP data used for statistical downscaling in this work were extracted from the National Centers for Environmental Prediction's Climate Forecast System Reanalysis (CFSR) dataset (Saha et al. 2010).CFSR SLP data are available at 6-hourly resolution from 1979 to 2009 on a 0.58 global grid.A global ocean wave reanalysis provided hourly directional wave spectra for the period 1979-2009 on a 1.58 longitude 3 18 latitude global grid (Perez et al. 2015).To generate this reanalysis, the third-generation spectral wave model WAVEWATCH III (WW3; Tolman 2009) was forced at a global scale by CFSR near-surface winds.Bathymetry and shoreline data were populated with ETOPO1 (Amante and Eakins 2009) and National Geophysical Data Center Global Self-Consistent, Hierarchical, High-Resolution Shoreline (Wessel and Smith 1996).Wave spectra were computed with 158 directional resolution and 32 frequency bins ranging nonlinearly from 0.0372 to 0.714 Hz with a factor of 1.1.

Parameterization of spectral data a. Wave generation areas
Wave climate along the California coast is a function of locally generated seas and swell generated in the North and South Pacific (Wingfield and Storlazzi 2007).To minimize the study area, a method developed by Perez et al. (2014) is used to identify the areas of wave generation contributing to local wave conditions at the particular location of interest.The method is referred to as ESTELA and uses geographic criteria and two-dimensional wave spectra to map areas of wave generation for a target point.
ESTELA reduces computational expenses by eliminating global grid cells that are blocked from the target point by landmasses, assuming that swell waves propagate along great circle paths (Snodgrass et al. 1966).A full discussion of the methodology and limitations can be found in Perez et al. (2014).Most importantly, ESTELA provides an estimate of the average effective energy flux toward the target point and the average travel time.
The ESTELA effective wave energy flux map for the location offshore of Southern California reveals two discrete wave energy generation regions: the North Pacific and an area in the western South Pacific near Australia and New Zealand (Fig. 2).These wave generation regions will contribute to the definition of the wave families.The travel times of wave energy reveal the time scales of wave propagation in the Pacific basin.Maximum travel time for swell generated in the North Pacific is 18 days, while the maximum travel time for swell generated in the South Pacific is 21 days (Fig. 2).Large differences in energy travel times across the basin (e.g., North Pacific swell arriving in 3 days generated in the eastern North Pacific and North Pacific swell arriving in 12 days generated in the western North Pacific) inspired improvements to the temporal definition of the predictor of Camus et al. (2014a).
b. Wave spectral partitions Camus et al. (2014a) defined the statistical predictand as daily bulk wave parameters (e.g., H s ), but the information of the full, directional wave spectrum is critical for areas affected by multimodal wave climates.To efficiently represent multimodal wave conditions, energy of the full directional wave spectra is split into three partitions.Spectral partitions were defined using an algorithm adapted from terrestrial watershed delineation (Hanson and Jensen 2004).Local seas (zeroth partition) were identified as energy traveling in directions consistent with concurrent wind directions with a wave speed less than 1.5 times that of the wind speed, and two swell partitions were identified by condensing the energy surrounding spectral peaks.The first and second partitions compose wave energy from the dominant and secondary swell bands, respectively.Bulk wave parameters for each of these three spectral partitions were saved hourly on a global scale, such that wave conditions at a location are defined by (H s0 , T p0 , D m0 , H s1 , T p1 , D m1 , H s2 , T p2 , D m2 ).It is important to note that this partitioning algorithm does not use a cutoff frequency to differentiate between seas and swell but instead characterizes seas and swell based on the ratio of wave and wind speeds.Portilla-Yandún et al. (2015) argue that partitioning based on spectral peaks is more physically accurate than use of a cutoff frequency.However, this results in energy associated with potentially improper partitions (swell with wave periods ,7 s).Improvements to the partitioning are not within the scope of this work.

Spatiotemporal definition of the predictor
Statistical downscaling efforts link a multivariate predictor to a multivariate predictand through a function that assumes stationarity of patterns and relationships in time.Here, atmospheric conditions were related to deep-water wave parameters.The predictor field was composed of daily standardized SLP anomalies and squared gradients of SLP anomalies (SLPG) at 28 spatial resolution (Wang et al. 2004;Camus et al. 2014a;Perez et al. 2015;Rueda et al. 2016a).Though CFSR SLP fields are available at 0.58 spatial resolution, comparable results and higher computational efficiency without a loss  2014a) assumed that SLP patterns of the last n days contribute to swell waves that reach the location.According to their methodology, the predictor field for waves on day i was the mean SLP and SLPG over the last n days.Stated differently, waves arriving at the location on day i were statistically linked to the mean atmospheric conditions over days i through (i 2 n).However, as stated above, wave energy arriving on day i could have been generated by atmospheric conditions over the Western North Pacific on day i 2 12 and/or by atmospheric conditions over the eastern North Pacific on day i 2 3 (Fig. 2).Following Camus et al. (2014a) and using n 5 5 for the case of Southern California has the dual effects of averaging out atmospheric conditions associated with storms (cyclogenesis over 3-4 days) and potentially missing the link between atmospheric conditions and wave energy separated in time by several days to weeks.
We attempt to improve upon the temporal definition to account for the vastness of large ocean basins by building a physically meaningful and intuitive predictor.Here, we assemble the predictor for waves arriving at the location on day i as atmospheric conditions over the area within the 1-day isochrone, atmospheric conditions from day i 2 1 over the area between the 1-and 2-day isochrones, atmospheric conditions from day i 2 2 over the area between the 2-and 3-day isochrones, and so on (Fig. 2).For the predictor field P i , P i 5 f . . ., SLP i2t11,V t , SLPG i2t11,V t , . . . .g, and t 5 f1, . . ., Tg, where V t is the area between isochrones t 2 1 and t, and T is the maximum number of isochrones (e.g., 18 days for the North Pacific and 21 days for the South Pacific).
Travel times computed with ESTELA are the average travel times for waves generated in a particular region.As such, this temporal delineation of the predictor is subject to uncertainty.However, in ocean basins as large as the Pacific, our improvement of the predictor is substantial.

b. Wave families
For locations experiencing multimodal wave climate, it is important to split both the wave climate and predictor into components defined by wave families.In this work, complexity is added to unimodal bulk wave parameters by accounting for multimodal spectra through redistribution of wave energy from spectral partitions (H s0 , T p0 , D m0 , H s1 , T p1 , D m1 , H s2 , T p2 , D m2 ) into wave families.Based on a priori knowledge of the wave climate at the location (Fig. 3), wave energy is redistributed into local seas (SEA), swell generated in the Northern Hemisphere (NH), and swell generated in the Southern Hemisphere (SH): where S( f, u) is the full directional spectrum, with f as the frequency and u as the direction.Wave information for each family (SEA, NH, and SH) was restricted to three commonly used descriptive statistics: H s , T p , and D m (e.g., H NH s ).Wave energy in the zeroth partition was local seas.Wave energy in the first and second partitions was split between Northern and Southern Hemisphere swell using D m (Fig. 3) and wave generation regions defined with ESTELA (Fig. 2).FIG. 6.As in Fig. 5, but for SH.
We associated the wave generation area in the North Pacific with NH swell and the area in the South Pacific with SH swell (Fig. 2).Mean wave directions approaching from angles between 2408 and 3608 were considered approaching from the Northern Hemisphere.Mean wave directions approaching from angles between 1408 and 2408 were considered approaching from the Southern Hemisphere.Though wave family parameters are available hourly, daily means are calculated for this statistical downscaling method.
Substantial improvements could be made to more accurately define seas and different swells.We suspect inconsistencies exist due to both the partitioning algorithm and the lack of T p as a criterion for identifying either NH or SH.

c. Construction of the predictor
A predictor was defined for each wave family (SEA, NH, and SH).The generation region for SEA encompassed the area over which wave energy reaches the location in 1 day (Fig. 2).The 2-and 3-day generation regions were also tested, with comparable results.To identify generation regions for swell wave families, a threshold was placed on the ESTELA effective wave energy maps, limiting generation regions to areas where energy flux .2kW m 21 8 21 (Fig. 2).The generation region for NH encompassed the area defined in the ESTELA map in the North Pacific (Fig. 2).The generation region for SH encompassed the area defined in the ESTELA map in the South Pacific (Fig. 2).This division coincided nicely with the previous choice to split NH and SH swell at an incident wave direction of 2408 (Fig. 2).SLP and SLPG over these identified areas were compiled using isochrones to incorporate the improvements to the temporal definition of the predictor.

Multiple multivariate linear regression
At this point in the methodology, we have defined a predictor (P NH , P SH , and P SEA ) and predictand D SEA m ) for each wave family.To define a functional relationship between each wave family's multivariate predictor and multivariate predictand, principle component (PC) analysis is first performed on the predictor to reduce the number of dimensions of the problem.We preserve the minimum number of PCs that explain 95% of variance of each wave family.The first four empirical orthogonal functions (EOFs) of the predictors for each wave family are displayed in Fig. 4. The first EOF of the NH predictor describes variation within the strength of the Aleutian low, while the remaining EOFs describe dipolar and tripolar variation over the North Pacific (Trenberth and Paolino 1981).These EOFs are likely correlated to large-scale atmospheric patterns or teleconnections, such as El Niño-Southern Oscillation, Pacific decadal oscillation, or the Pacific-North American pattern, though those correlations were not investigated in this study.
For a statistical downscaling of local wave climate, the predictor in this new PC space can be classified using a K-means algorithm to nonlinearly relate atmospheric and wave conditions for SEA, NH, or SH using a weather-type approach (Camus et al. 2014b).This work only focuses on improvements to the definition of the predictor, but applications of a weather-type approach can be found in Camus et al. (2014b), Perez et al. (2015), Rueda et al. (2016a), andRueda et al. (2017).Here, to assess the ability of the predictor to define daily wave conditions, a multiple, multivariate, linear regression model is applied to the daily PC and wave time series over a calibration period .The regression model is assessed over a validation period (2000-10) using the correlation coefficient R, root-mean-square errors (RMSE), and bias: where the subscripts mod and obs refer to parameters x modeled with the linear regression and ''observed'' parameters from the wave reanalysis, respectively; s is the standard deviation of x, and N is the number of data points in the time series.We test improvements to the predictor in two stages: 1) inclusion of isochrones to account for energy travel time in large ocean basins and 2) splitting of energy into wave families to more completely represent directional wave spectra.First, we assess the predictor with the improvement of the isochrones for total bulk parameters.Second, we incorporate the wave families and isochrones for the fully improved predictor.A summary of these statistics can be found in Figs.5-7 and Table 1.The improved predictors are successful for reproducing H s (R 5 0.86, 0.81, and 0.87 for NH, SH, and SEA, respectively) and reasonably successful for T p and D m (R # 0.7); T p is the most difficult parameter to reproduce, likely due to the limited skill of the partitioning algorithm to split local seas and swells in complex sea states.Particularly, the variance is reduced in modeled T p , meaning the regression is unable to reproduce TABLE 1. Correlation coefficient R, RMSE, and bias for wave parameters significant wave height H s , peak wave period T p , and mean wave direction D m for Southern California and Trujillo, Peru, using the methodology defined in Camus et al. (2014a and the two improvements of this work: 1) inclusion of isochrones and 2) definition of wave families.In the case of C2014, error statistics are calculated for H s .In the case of the isochrones improvement, error statistics are calculated for aggregated H s , T p , and D m .In the case of the isochrone and wave family improvement, error statistics are calculated for H s , T p , and D m for each wave family (NH, SH, and SEA).1).

Additional application
To further test the new model, the methodology is applied to a location in deep water offshore of Trujillo, Peru (88S, 79.58W).Trujillo experiences bimodal wave conditions, with wave energy contributions from both the North and South Pacific, making it an ideal location for the application of this new methodology (Fig. 8).Similar to the application for Southern California, wave energy was split into three distinct wave families.SEA was identified as the first partition.Northern Hemisphere swell and Southern Hemisphere swell were distinguished by D m .Mean wave directions approaching from angles greater than 2708 were considered approaching from the Northern Hemisphere (Fig. 8).Mean wave directions approaching from angles less than 2708 were considered approaching from the Southern Hemisphere (Fig. 8).A predictor was defined for each wave family using ESTELA and the methodology described above.
While the methodology is unable to produce correlation coefficients as high as for the Southern California application, it is successful by RMSE and bias metrics (Table 1).Similar to the previous application, the predictors most successfully reproduce H s (R 5 0.84, 0.83, and 0.84 for NH, SH, and SEA, respectively).RMSE and bias are small for H s , T p , and D m .We compare these results again to those employed in Camus et al. (2014a) because this location was included in their original study.Camus et al. (2014a) defined the predictor as the full extent of the Pacific basin using n 5 12 days.As in the previous application, the improved predictor defined here increases R and decreases RMSE and bias for each wave family (Table 1).

Conclusions
An improved method for defining the optimal predictor for statistical downscaling of the multimodal wave climate was presented in the context of two applications where this method yields better results than previous work.In areas where wave spectra exhibit multiple modes due to seas and swells approaching from different generation regions, using a predictor that spatially encompasses only the generation region of the dominant swell mode or describing the wave climate as bulk wave parameters (H s , T p , and D m ) of the dominant peak of the spectrum is insufficient.Additionally, in large ocean basins where simultaneously arriving waves may have been generated several days apart, averaging atmospheric conditions across multiday time scales leads to errors in the timing of relationships between the predictor and predictand.By redistributing energy from spectral partitions into wave families and accounting for the average travel time of waves generated in different parts of the basin through isochrones, a unique spatiotemporal predictor can be defined to successfully reproduce local multimodal wave conditions.
The improved predictor is tested for locations offshore of Southern California and Trujillo, Peru.Both locations experience multimodal wave spectra, with energy contributions from local seas and swells generated in both the Northern and Southern Hemispheres.For these applications, the methodology is successful, increasing the ability of the atmospheric conditions to reproduce daily multivariate wave parameters compared to previous work.
Coastal scientists and engineers are moving toward representing wave conditions using the full wave spectrum, as parameterizing the dominant peak of the wave spectrum leads to loss of sea and secondary swell energy.The proposed method is the base of statistical downscaling of local wave climate by weather types, which can be used to project historical and future wave climatologies, simulate realizations of time series of wave parameters, project extremes in wave heights or wave run-up, or project historical and future coastal response.Applications of this methodology are as varied as climate change, coastal flooding hazards, and shoreline change.

FIG. 1 .
FIG. 1. Flowchart of the general methodology to define the optimal predictor for statistical downscaling of multimodal wave climate.Improvements to the Camus et al. (2014a) definition of the predictor are highlighted.
FIG. 2. ESTELA effective wave energy flux (color bar) for the location (large red dot) offshore of Southern California.The gray shading over the ocean, outlined generally in white dashed lines, denotes a threshold of 2 kW m 21 8 21 imposed to define important wave generation regions.The isochrones of travel time are gray and black lines, with every third isochrone labeled in days.The red dotted line emanating from the location shows the demarcation of 2408 for North and South Pacific swell families.The shaded area above the red dotted line is the spatial definition of the North Pacific swell family predictor (NH).The shaded area below the red dotted line is the spatial definition of the South Pacific swell family predictor (SH).The shaded area inside of the 1-day isochrone (black dashed line) is the SEA.

FIG. 3 .
FIG.3.This method of defining the optimal predictor requires a priori knowledge of the wave climate at the location.The mean directional wave spectrum from 1979 to 2009 for a location offshore of Southern California exhibits bimodality.Evidence for both North and South Pacific swell wave families is present in the secondary peak of the directional wave spectrum.The red line indicates the division of wave energy at 2408 between the North and South Pacific swell wave families.

FIG. 5 .
FIG. 5. Multiple multivariate regressions for the location offshore of Southern California of NH significant wave height H s , peak wave period T p , and mean wave direction Dir over a validation period (2000-09), where reanalysis parameters are gray and estimated parameters are black.Scatterplots and error statistics are shown for each parameter.

FIG. 8 .
FIG. 8. ESTELA effective wave energy flux (color bar) for the location (large red dot) offshore of Trujillo, Peru.The gray shading over the ocean, outlined generally in white dashed lines, denotes a threshold of 2 kW m 21 8 21 , imposed to define important wave generation regions.The isochrones of travel time are gray and black lines, with every third isochrone labeled in days.The red dotted line emanating from the location shows the demarcation of 2708 for North and South Pacific swell families.The shaded area above the red dotted line is the spatial definition of NH.The shaded area below the red dotted line is the spatial definition of SH.The shaded area inside of the 1-day isochrone (black dashed line) is the SEA.