Estimating Coastal Winds by Assimilating High-Frequency Radar Spectrum Data in SWAN

Many activities require accurate wind and wave forecasts in the coastal ocean. The assimilation of fixed buoy observations into spectral wave models such as SWAN (Simulating Waves Nearshore) can provide improved estimates of wave forecasts fields. High-frequency (HF) radar observations provide a spatially expansive dataset in the coastal ocean for assimilation into wave models. A forward model for the HF Doppler spectrum based on first- and second-order Bragg scattering was developed to assimilate the HF radar wave observations into SWAN. This model uses the spatially varying wave spectra computed using the SWAN model, forecast currents from the Navy Coastal Ocean Model (NCOM), and system parameters from the HF radar sites to predict time-varying range-Doppler maps. Using an adjoint of the HF radar model, the error between these predictions and the corresponding HF Doppler spectrum observations can be translated into effective wave-spectrum errors for assimilation in the SWAN model for use in correcting the wind forcing in SWAN. The initial testing and validation of this system have been conducted using data from ten HF radar sites along the Southern California Bight during the CASPER-West experiment in October 2017. The improved winds compare positively to independent observation data, demonstrating that this algorithm can be utilized to fill an observational gap in the coastal ocean for winds and waves.


Introduction
Regional coastal modeling of the ocean and the atmosphere have come a long way in the last 20 years but most models still suffer from errors due to parameterization or inaccurate model inputs (i.e., bathymetry, initial or boundary conditions). These model errors can be addressed with the assimilation of local observations. In this study, we examine coastal wind estimation by assimilating observations of Doppler spectra from coastal high-frequency (HF) radar sites into a regional wave model with a variational assimilation approach. The wave model used here is SWAN ( [1,2]) which is a state-of-theart third-generation wave model used predominantly in coastal and inland waters.
HF radar data assimilation studies commonly utilize the surface current measurements traditionally associated with this data source. A few of these studies are [3][4][5]. The use of HF radar Doppler spectra to correct wave models is less widely used. Siddons et al. (2009) [6] examined the use of three different approaches for HF data assimilation into SWAN: (1) Ensemble Kalman Filter, (2) Ensemble Optimal Interpolation and (3) 3D variational scheme. This study makes use of beam-forming HF radars. Waters et al. (2013) [7] implemented an optimal interpolation assimilation approach that makes use of the windsea and swell parts of the spectrum separately. To the authors' knowledge, the 4D-variational approach presented here has never been applied to HF radar data for wave modeling.
The methodology for the variational assimilation framework in the SWAN model was first laid out by Walker (2006) [8] and then used by Veeramony et al. (2010) [9]. The approach used only examined the spatial and spectral advection terms while ignoring source terms. Orzech et al. (2013Orzech et al. ( , 2016 [10,11] expanded the implementation to include all the source terms with a new version of the adjoint of the wave-action-balance equations. Walker and Brunner (2021) [12] applied a stationary buoy assimilation setup with a similar methodology as presented here and showed good model skill when compared to both assimilated and unassimilated buoy observations. Shore-based HF radars are routinely used to map surface currents by measuring the Doppler shifted backscatter from ocean waves. The backscattered energy from the ocean surface is due to Bragg scattering of the transmitted electromagnetic wave by ocean surface waves. Coherent reflections in the backscattered spectrum are generated at ocean wavelengths exactly 1 2 the wavelength of the transmitted wave. Most applications make use of HF radars as a tool for producing high-resolution maps of surface currents. It is also possible to extract wave information from the second-order portion of the spectrum, where reflections are generated from waves at all frequencies not just the Bragg peaks. See [13] for a brief introduction to HFR theory.
The objectives of this work were to develop a system for assimilation of HF radar data for general nearshore domains, to produce improved wave nowcasts/forecasts, and provide wind forcing corrections. Initial implementation of our HF radar assimilation algorithm in SRI's swanX 5DVAR system (based on SWAN 41.20) is complete. It includes a ground-wave HF radar model, based on first-and second-order Bragg scattering. The model uses the wave spectrum output from SWAN and the winds and ocean current field from COAMPS-OS along with the HF system parameters. The adjoint of the HF radar model translates errors in the HF radar spectrum into errors in the wave spectrum. The adjoint of SWAN then produces an effective wind error, consistent with the ST6 wind growth model. The whole system is converged, and wind field updates are made to the COAMPS wind inputs. Here, we present an initial demonstration of the SWAN HF radar assimilation system for CASPER-West experiment. Comparisons between the wind estimates from the SWAN HF assimilation and independent buoy data have shown the algorithm capable of correcting winds inputs for this coastal domain.
The rest of the paper has the following structure: Section 2 describes the SWAN and HFR forward scattering models; Section 3 presents the data assimilation framework including the cost function, adjoint models and implementation. The results of a realworld application of the assimilation are discussed in Section 4 followed by conclusions in Section 5.

The Models
The algorithm used here is based on the open-source SWAN wave spectrum model and makes use of HFR Doppler spectra data. We briefly describe the model, and then discuss the observed HFR Doppler spectra data. Then, we discuss the HFR forward model and the assimilation framework. Then, we examine the wind model and finish with a brief discussion about implementation.

SWAN Wave Model
SWAN is a third-generation wave model for coastal regions [1,14]. The basic equation used is the spectral balance equation for wave action N(x, s, t) = E(x, s, t)/σ, where E is the variance density spectrum, σ is the intrinsic wave radian frequency, x = (x, y) spatial location, s = (σ, θ) spectral location, and t is time. The action balance equation is given by where ∇ = ( ∂ ∂x , ∂ ∂y , ∂ ∂σ , ∂ ∂θ ), and C = (C x , C y , C σ , C θ ) are the wave-energy propagation velocities in physical and spectral space, consisting of wave group velocities plus ocean currents for physical space and the effects of refraction and wave-current interaction for spectral space. These are described in detail in [1], as is the source term on the right-hand side of Equation ( 1). This source term includes models for the effects of wind growth and loss of energy due to white-capping, bottom friction, and depth-induced breaking, as well as energy transfer within the spectrum due to nonlinear wave-wave interactions (resonant triad and quartet interactions). These source-term models can be written collectively as where A, the so-called "linear" wind-growth term, is independent of E, and all other terms are proportional to E, but with a collective coefficient S E that depends at least weakly on E. This set of equations can be solved for the action spectrum for spatial region Ω subject to appropriate initial conditions on Ω and boundary conditions on ∂Ω. For portions of the wave spectrum with propagation velocities that carry energy into Ω, the "incident" wave spectrum E(x, s, t) on ∂Ω must be specified. The spectrum is taken as periodic in θ and the σ boundaries are located far above and below the energy-containing region of the spectrum, yielding E = 0 at these boundaries. In addition to these boundary and initial conditions, complete specification of the mathematical problem requires that the winds and bathymetry, including any tidal offset, be prescribed for Ω.
SWAN version 41.20 (used here) has been updated from 40.85 to include dissipation of swell energy, non-breaking dissipation, optional observation-consistent wind input, and white-capping physics; the complete suite of improvements is known as ST6. There are two options for the implementation of non-breaking dissipation. The first uses the formulation by Young et al. (2013) [15], updated by Zieger et al. (2015) [16], and the second uses the formulation by Ardhuin et al. (2010) [17]. We utilize the latter, with the default values. The details of the ST6 physics can be found in Rogers et al. (2012) [18]. Notably, this formulation changes the standard Komen et al. (1984) [19] wind speed scaling to a similar scaling with a free parameter, i.e., U = S ws u * . Other concurrent changes introduced include changes to the wind drag formula and wave growth term, described in the SWAN documentation [20].

HF Radar Scattering Model
We start with the following expression for the observed energy in the Doppler spectrum σ, as a function of Doppler frequency −∞ < ω < ∞, and range r > 0, from a cross-loop HF radar system as the sum of the first-and second-order Bragg scattering contributions. σ(ω, r) =σ 1 (ω, r) +σ 2 (ω, r) .
The first-order scattering is given bŷ where (r, φ) is the spatial polar coordinate system, centered on the radar, and k 0 is the electromagnetic wavenumber. Here, S(x, k) is the variance spectral density (the wave spectrum, equivalent to E(x,s) at a fixed time), k is the hydrodynamic wave vector, is the wave vector associated with the Bragg waves propagating toward the radar (here, waves at −k b are propagating away from the radar) and direction is idicated by m = ±1 . The intrinsic radian wave frequency ω k = gk tanh kh, where k = |k| and h = h(x) is the water depth, and ω b = ω k evaluated for k = k b . The illumination patterns w n (φ) for the loop antennas are where φ 0 is a reference for the antenna orientation.
At each spatial location, there is a near-surface ocean current u(x) and this will impart a Doppler shift ∆ω = k · u. This shift can be applied as follows, here, ω is Doppler-shifted so that ω =ω + k b · u. Applying to Equation (5) and integrating inω yieldŝ this is the final form of the first-order Bragg-scattering model, where n = 1, 2 defines the antenna. For second-order scattering, the contribution to the Doppler spectrum, equivalent to Equation (3), iŝ where k = −m k b − k, with m = ±1, is an additional term defining propagation direction of the Bragg waves, and Γ is a coupling coefficient where ∆ = 0.011 − 0.012i for seawater, and (k, ϕ) are polar coordinates in the wave spectrum k-plane. It should be noted that the last term in Γ goes to zero for deep water (kh → ∞).
Accounting for the Doppler shift associated with near-surface currents yieldŝ

The Assimilation Framework
Assimilation of HFR data is carried out in a variational framework based on the models described above. This requires definition of a cost function and one or more control variables that are adjusted to minimize the cost function. In this case, the cost function is a measure of the error between the predicted and observed HFR Doppler spectrum, and the control variable is the surface wind field used in the SWAN model. The adjustments to the wind field are guided by the solution of an adjoint version of the SWAN model with input from the adjoint HFR model which, in turn, has as its input the error in the present HFR prediction. The gradient of the cost function with respect to the wind, used to adjust the wind field, is computed from the adjoint wave action spectrum, using an expression derived from the wind-growth term in the SWAN model. In the sections below, we first describe the cost function, followed by the adjoint forms of the models, and then the expression for the gradient. We then describe the sequence of operations for the assimilation algorithm.

Cost Function
The Doppler spectrum at second order has two components that are summed to produce the complete spectrum. Since they are summed, we can treat them separately and then combine the results. First, we must examine our data. The model for the first-or second-order returns produces the noise-free Doppler spectral density in units of normalized radar cross section (square meters per square meter) per radian (i.e., the units are rad −1 ). The data s, on the other hand, are uncalibrated and contain noise, and can be represented as where C n is a fixed (for channel n) calibration constant and e is a fixed noise level. Rearranging yields In practice, e can be estimated from the extreme upper and and lower frequency bands where the signal energy is unlikely to appear; however, this can lead to negative spectral density values in regions where s ≈ e. For this reason, we will use For the systems of interest, the calibration constant C n is unknown, so construction of a cost function based on the error variance such as whereσ is the predicted spectrum, is problematic. We could attempt to estimate C n from the data by requiring the total backscattered energy, or the peak backscatter, to be the same for the model and data, but those approaches have issues of their own. If we define our cost function using the derivative in frequency of the difference in the log of the predicted and observed spectral density, we obtain We will use this form, which is insensitive to the unknown calibration constant, going forward.

Adjoint Models
Here, we describe the adjoint models used in the algorithm. We begin with the adjoint form of the SWAN model, followed by the adjoint forms of the first-and second-order HFR backscatter models.
The adjoint form of the SWAN model, described in Veeramony et al., 2010 [9], Orzech et al., 2016 [11], and Walker and Brunner 2021 [12], is where A(x, s) is the adjoint wave action spectrum, and S D is an 'observation' term, discussed below. The adjoint SWAN model is solved backward in time over the period of interest, subject to homogeneous boundary conditions. The HFR spectrum model has first-and second-order parts; this will result in two observation terms in Equation (15), S D = S D1 + S D2 . For S D1 , we setσ =σ 1 in J to get J 1 ; taking the first variation with respect to E yields Integrating by parts yields where Note that in Equation ( 17), integration by parts produces an additional term evaluated at the upper and lower boundaries of ω, where the spectrum, and hence this term, will vanish. Then, combining Equations ( 17) and ( 18) and re-arranging the order of integration results in where we have combined the r and φ integrations into a single 2D spatial integration in x (i.e., dx/r = dr dφ). The resulting contribution to the adjoint SWAN observation term is a portion of the integrand in this expression multiplying δE, where we have integrated in ω, evaluating the integrand at the shifted Bragg frequencies (ω = ±ω b + k b · u). The gradient term is assigned at the Bragg wave-vector locations in the wave spectral domain k = ±k b and, in the spatial domain, at x = (r, φ). For use in the SWAN assimilation algorithm, S d (k, x) = δS(k, x). We now compute the the second-order contribution to the observation term S D2 . We start by inserting the expression for the second-order Doppler spectrum into the cost function and then take the first variation with respect to E and proceed as above for S D1 . This yields where k = −mk b − k, ω d = k b · u + mω k + m ω k and m, m = ±1 where, again, these are for approaching/receding waves (respective halves of the k plane). In S D1 and S D2 , the error in the Doppler spectrum at ω d is mapped to a correction in the wave spectrum at k and the correction is scaled by the scattering coefficient Γ and, for the second-order contribution, wave energy spectral density at the resonant wavenumber.
The final piece needed for the assimilation algorithm is an expression for the gradient of the error in the Doppler spectrum (the cost function) with respect to the wind that can be used to guide adjustments to the wind. The source function S on the right-hand side of the SWAN model Equation (1) contains the wind-growth source term, S in (detailed in [14]). The desired gradient is obtained by taking the first variation with respect to the vector wind. The first variation of the cost function J with respect to the vector wind U w is where the vector nature of the wind is acknowledged. This equation governs the reduction in J that will result from a change in U w . To reduce J, we set δU w the incremental change in the wind to where α is positive. In the SWAN model, there are many choices for S in , and they all are somewhat complex empirical parameterizations. For this reason, we compute δS in /δU w by perturbing the vector wind components by a small amount, and calculating the change in S in , divided by the component perturbation. The result is integrated over the spectral domain to yield a spatially-resolved time-varying gradient to guide the correction to the vector wind components.

Implementation
The initial guess for the wave spectrum over the entire region is a previously corrected SWAN model background that has corrected the swell due to assimilation of National Data Buoy Center (NDBC) buoy data. The adjoint Equation (15) is solved using the HF radar observations as input, the gradient. The gradient is then used in a conjugate-gradient procedure [21] to calculate the wind field that results in the minimum of the cost function J. In practice, the following sequence of steps is followed:

•
The adjoint solution is calculated using the error in the most recent prediction as input; • Using the adjoint solution, the gradient is determined; • The conjugate-gradient descent algorithm is used to calculate a new estimate of the wind field and U 10 ; • The SWAN model is run with corrected inputs and a new wave-spectrum prediction for the region is generated; • The forward HF radar model is run with the new spectrum as input and a new prediction of the data is generated.
A graphical description of the complete algorithm structure is shown in Figure 1. The algorithm makes use of operation forecast data from the COAMPS model (as a first guess for the winds), models and adjoint models (the SWAN wave spectrum model and the HF radar forward spectrum model and corresponding adjoints), ancillary data such as bathymetry, and HF radar Doppler spectrum data (from CODAR HFR sites). The algorithm outputs are wind and wave products, improved wave spectra and improved estimates of the wind field. Estimates begin with a first guess wind field for the region obtained from operational wind forecasts. Estimates of the wave spectrum and the HF radar Doppler spectrum are calculated. The HF radar Doppler spectrum is compared to that for the data (from the CODAR sites) and the difference is fed back through the adjoint HFR and SWAN models. The gradient of the error in the estimated HFR spectrum with respect to the wind field is calculated from the adjoint wave spectrum. This gradient is used to adjust the wind field using a descent algorithm and the steps are repeated until the wind field converges and the HFR spectrum is a best fit. The algorithm makes use of operation forecast data (as a first guess for winds), models and adjoint models, ancillary data (such as bathymetry), and HF radar Doppler data. The algorithm outputs are wind and wave products, improved wave spectra and improved estimates of the wind field.

Results
The HFR inverse modeling framework developed here was applied for estimation of the 10-m winds in the Southern California Bight during the CASPER-West field experiment in October 2017. The aim here was to utilize HF radar Doppler spectra observations to estimate the winds near the coast. The SWAN model (version 41.20) was run with all the relevant model physics enabled (i.e., wind generation, bottom friction, three-wave and four-wave interactions, depth-induced breaking, dissipation of swell energy (non-breaking dissipation) and optional observation-consistent wind input and white-capping physics). The complete suite of improvements are known as ST6 with all input parameters set to their default values. The details of the ST6 physics can be found in Rogers et al. [18]. The initial guess for each 24-h assimilation window was taken from a previously computed SWAN assimilation using NDBC buoy observations to correct the open boundary conditions. The algorithm product is a wave and wind estimate for the region from a non-stationary SWAN model run for the time of interest that is a best fit to the observational data.

HF Radar Data Description
Here, we utilize the Doppler spectra returns of sea state from 10 SeaSonde type Coastal Ocean Dynamics Applications Radar (CODAR). The locations of HFR sites can be seen in Figure 2. Due to the nature of the HF radar, observations are generated on a polar coordinate grid centered on the antenna locations. The range cells on this polar coordinate grid are spaced at 1.5 km. The Doppler spectra data from ocean waves are received at each HF site and are calculated by spectral processing of the signal within each range bin. The received signal was sampled at approximately 2 Hz and Fourier-transformed with a 512 sample window, resulting in a Doppler spectrum being recorded approximately every 4 min, each spectrum containing 512 Doppler frequency samples from −1 to +1 Hz. The data used for this study include these raw Doppler spectra, which are not published online, for the 10 CODAR antennas. Each site consists of two cross-loop antennas and one omnidirectional antenna. The omnidirectional data was not considered to be useful for this study. For purposes of comparison to model predictions, the data were averaged over 3 range cells and a roughly 1-h time interval. Thus, hourly predictions of the Doppler spectrum were made for 8 range samples from 10-42 km.  Figure 3 shows an example comparison of range-Doppler plots for the Refugio State Beach (RFG1) HF radar site. The 1st order Bragg peaks are visible in both the observed and predicted data. The 1st order peak for the positive Doppler frequencies (waves coming towards the antenna) is stronger than for the negative Doppler frequencies (wave receding from the antenna). This trend is more prominent in the predicted spectra than the observed spectra although it is still present. The second-order peaks are also visible in both the predicted and observed plots although much less obvious than the first-order peaks.

Comparisons to Buoy Wave Data
In this section, we present comparisons between wave parameters from NDBC buoy 46217 and the predicted wave spectra from the same location. Comparison at buoys 46025 and 46221 are similar, so they are not shown in detail. We limit our comparisons to the significant wave height H s , the mean wave period T m , and the mean wave direction θ m . The figures presented here show the general performance of the algorithm at these buoy locations before we move on to the wind comparison. Figure 4 shows a comparison of the estimated wave parameters to the observed data from buoy 46217 in the Anacapa Passage near the Channel Islands. The water depth at this buoy location is 114 m. The top panel of Figure 4 is a time series of the estimated H s , compared to the buoy observations. There is good agreement between the estimated and observed values for the majority of the study period. There are some times where the predicted H s overestimates the observations. This is most likely due to a poor specification of the boundary conditions during several assimilation cycles. The middle panel of Figure 4 shows the time series comparisons for T m . Again, there is good agreement for a large portion of the month with the prediction underestimating the T m for several times. The bottom panel of Figure 4 is the time series comparison for θ m . The observations of θ m show a wide spread of values with estimates generally following the same pattern as the observations. Quantitative error statistics for H s , T m , and θ m for buoy 46217 are presented in Table 1. Table 1. Error statistics for significant wave height (H s ), mean wave period (T m ) and mean wave direction (θ m ) from comparison of the SWAN assimilation results to observation data used in the background. The mean error in quantity X is X , the RMS error is rms X . For comparison, the mean and standard deviation of the observations (X and σ X ) are included on the second line, in parentheses.

Comparisons to Buoy Wind Data
Here, we compare the 10-m winds from the SWAN HF assimilation algorithm against winds from NDBC buoy 46053. The buoy winds have not been incorporated in any way into our assimilation framework, so, they can be treated as independent observations. We restrict our comparisons to the U and V components of the wind. Figure 5 shows the comparisons for buoy 46053 located East of Santa Barbara in a depth of 405 m. This continuous record of 10-m winds was generated using the middle 12 h from each of the 24-hour assimilation cycles. The comparison shows that the SWAN wind estimate solution closely matches the winds from the NDBC buoy for the majority of this study period. The COAMPS first guess solution does not agree with the buoy data as well as the HF SWAN solution. Notably, in the U time series (top panel), there are several times on 3 October and 10 October where there is a significant difference between the first guess and the assimilation estimate where the assimilation estimate agrees with the observations. This points towards errors in the coastal wind field coming from COAMPS. Figure 6 shows a side-by-side field plot of the 10-m winds on 3 October 2017 at 2300 z. The left panel shows the COAMPS first guess winds that contain a strong eastward (positive U) wind feature in the Santa Barbara Channel. The right panel shows the SWAN HF assimilation 10-m winds with much weaker winds in the Santa Barbara Channel relative to the COAMPS first guess.
The observations from buoy 46053 agree with the SWAN HF assimilation winds and show the potential for corrections to errors in coastal wind forecasts.

Discussion and Conclusions
In this work, we have described the development, implementation, and testing of an HF radar assimilation algorithm for the SWAN model. This framework utilizes a cost function to determine the error between the observed and predicted Doppler spectra from an HF radar first-and second-order forward model. Adjoint models of both the HF radar forward model and the SWAN model are then solved. The first-and second-order gradients of the cost functions are calculated from the adjoint spectrum solutions. Lastly, a conjugate-gradient algorithm is used to determine the appropriate Doppler spectrum and wind conditions that minimizes the cost function.
The algorithm was applied and validated using real-world observations from coastal wave buoys and HFR sites during the CASPER-West field experiment in October of 2017. The background wave field was taken from a buoy assimilation experiment where the open-boundaries were corrected. Predicted and observed significant wave height H s , mean wave period T m , and mean wave direction θ m agreed well over the study period with a few short periods of disagreement.
In conclusion, a variational HF radar Doppler spectrum assimilation algorithm for the SWAN model has been developed. Corrections in the near surface winds during the CASPER-West experiment were generated. A comparison of these winds to independent buoy observations shows good agreement and the skill of the algorithm in adjusting errors in coastal wind forecasts. Additional testing of the algorithm in other regions of interest is necessary before this tool can be utilized for operational forecasting. Further expansion of the algorithm to potentially account for bistatic or skywave HF radar configurations is ongoing and should allow for more general utility and relocatablity.  Data Availability Statement: All wave buoy and bathymetry used in this study are publicly available from the U.S. National Oceanic and Atmospheric Administration. The HFR Doppler spectra dataset is a level 0 data product that is not routinely reported, but can be obtained directly from the site operators.