HF Radar Measurements of Surface Waves in the Gulf of Naples (Southeastern Tyrrhenian Sea): Comparison With Hindcast Results at Different Scales

HF radar systems wave measurements are evaluated against numerical simulations in the Gulf of Naples (Southeastern Tyrrhenian Sea). Wave measurements are obtained from three CODAR SeaSonde HF radars installed along the coast of the Gulf of Naples. The numerical models employed are WavewatchIII, implemented on a regional scale with a resolution of about 10 km in longitude and latitude in the whole Mediterranean Sea, and SWAN, implemented with a 200 m resolution in the area of interest. Numerical simulations are also validated against experimental data acquired by a buoy installed offshore the Gulf of Naples. The agreement between HF radar measurements and model hindcasts is evaluated through the estimate of statistical error indices for the main wave characteristics (significant wave height, mean period, and mean direction). The consistency between wave parameters retrieved by HF radars and hindcasted by the models opens the way to future integration of the two systems as well as to the utilization of HF radar wave parameters that could be envisaged for data assimilation in wave models.


INTRODUCTION
Monitoring the physical processes taking place in the nearshore region is a fundamental pre-requisite for a proper integrated management of the coastal zones. In this framework, land-based remote sensing by HF radars presently provides a challenging opportunity for simultaneously measuring surface currents and wave parameters (Paduan and Rosenfeld, 1996;Rubio et al., 2017;Capodici et al., 2019). Crombie (1955) described the physics underlying the acquisition by HF radars of the backscatter echo coming from the rough surface of the sea produced by the resonant first-order Bragg waves. First-order echoes provide measurements of surface currents (Paduan and Graber, 1997;Gelpi and Norris, 2003). Second-order echoes are weaker and noisier than the first-order ones (Gurgel et al., 2006) but information on wave parameters can be derived from this part of the spectrum using methods of integral inversion (Barrick, 1977(Barrick, , 1979Lipa and Nyden, 2005).
An advantage of such instruments is that they can collect continuous real-time wave measurements over wide areas, while traditional wave buoys can record parameters only relative to a fixed position (Lorente et al., 2019). Only over the last years the accuracy of HF radar wave measurements has started to be tested (Lopez et al., 2016), thus the analysis of the spectrum for wave applications is still a developing topic (Gurgel and Schlick, 2005).
In addition to instrumental measurements, several wave models used for both hindcasting and forecasting wave parameters are currently available. Among the most commonly adopted in the scientific literature are: WavewatchIII (WWIII), used as operational wave forecasting system by National Oceanic and Atmospheric Administration (NOAA) (Tolman and Chalikov, 1996); WAM (WAve Model), included in the Integrated Forecast System of European Center for Medium-Range Weather Forecasts (ECMWF) (Komen et al., 1994); and SWAN (Simulating WAves Nearshore), developed by Delft University of Technology (Booij et al., 1999) and aimed at downscaling the wave field up to the shallow waters. Previous investigations revealed how the parametrization of coastal effects such as (but not limited to) wave refraction, damping, shoaling and breaking is not trivial, leading to locally varying wave field and making an absolute comparison between wave models and HF radars difficult .
In a recent work (Saviano et al., 2019), the performance of a HF radar system operating in the Gulf of Naples (Southeastern Tyrrhenian Sea, Western Mediterranean Sea; hereinafter referred to as GoN, see Figure 1) was evaluated in comparison to wave measurements recorded by a fixed buoy. As the buoy was not positioned in the area directly covered by the HF radars, it was possible to evaluate the consistency of the patterns of measurements rather than the specific values. Notwithstanding this limitation, the comparison indicated a general good agreement between the two platforms, in particular with the HF radar installed in proximity of lower edge of the domain (therefore closer to the buoy, Figure 1), scoring a yearly averaged r = 0.73. By contrast, the comparison between the buoy and the two HF antennas located in the innermost sectors of the GoN returned lower numerical consistency, although similar patterns were pointed out (with yearly r values of 0.61 and 0.58). This result was expected, not only owing to the non-negligible physical distance between the buoy and the two innermost sites, but also due to the wider fetch possibly leading to more intense sea states and to bathymetric issues ( Figure 1A). On the basis of this significant variation, additional investigation was needed; as such, the present study provides further investigations by comparing HF radar wave measurements in the GoN with the outputs of two wave models, WWIII and SWAN. WWIII provides information at a regional scale (Mentaschi et al., 2013a(Mentaschi et al., , 2015 while SWAN allows to downscale the wave features over finer resolutions, accounting for the complex morphology of the coastline (Booij et al., 1999). A first validation of the wave models was obtained through the comparison with wave measurements collected by a wave buoy installed at the outer edge of the GoN, which allows a direct comparison between the modeled and the sampled data. HF radar measurements were subsequently compared with the numerical model outcomes over three years, from January 2010 to December 2012. Such a long-term approach is, to the best of the present knowledge, unique, as previous studies covered shorter periods (usually from days to months) (see e.g., Wyatt et al., 2003;Siddons et al., 2009;Waters et al., 2013;Orasi et al., 2018;Lorente et al., 2018;Lopez and Conley, 2019). The paper is structured as follows. Section 2 presents the data used for the study and describes the locations of the wave measurement devices, while the indices employed for validating the models outcomes and evaluating the reliability of HF radars measures are presented in section 3. Sections 4 and 5 show and discuss the results of the aforementioned comparisons. Finally, in section 6 conclusions are summarized and possible future developments are introduced.

Observations: HF Radars and Wave Buoy
Waves were both measured indirectly and directly. In the former case, a SeaSonde HF radar network covering the internal part of the GoN has provided spatial information on the main wave characteristics, integrated along range bins at roughly 1 km resolution (Figure 1). In addition, direct measurements from a directional buoy moored off Capri island (southwestern edge of GoN; Figure 1A) have been used as a benchmark in order to validate both the regional wave hindcast (developed with WWIII) and the downscaled local waves propagation inside the GoN (developed with SWAN).

HF Radar Wave Observations
The GoN is a coastal embayment along the Southeastern Tyrrhenian Sea, featuring key natural ecosystems and archeological sites, together with anthropized areas and brownfileds (e.g., Ferrigno et al., 2018;Armiento et al., 2020;Aucelli et al., 2020;Margiotta et al., 2020). Wave measurements inside the GoN were collected by a 25 MHz SeaSonde HF radar system (manufactured by CODAR Ocean Sensors Ltd), operated by the Department of Science and Technology at the "Parthenope" University of Naples. The network, used for the simultaneous measurement of the surface current (Menna et al., 2007;Uttieri et al., 2011;Cianelli et al., 2012Cianelli et al., , 2015Kalampokis et al., 2016;Ranalli et al., 2018;Bagaglini et al., 2020) and wave fields (Falco et al., 2016;Saviano et al., 2019Saviano et al., , 2020 was set up in 2004 with two monostatic antennas, in Portici (PORT) and Massa Lubrense (SORR), and in 2008 it was integrated with a third site in Castellammare di Stabia (CAST, see Figure 1). Such data have also been used in operational contexts and for physical-biological applications Cianelli et al., 2017). The GoN network is the first and longest-running installation in Europe (Rubio et al., 2017), and is presently operational, although since 2013 the SORR site has been discontinued. The SeaSonde system recorded and averaged single spectrum characteristics every 10 min along 1 km equally spaced circular annuli (range cells, RCs) centered on the antenna using a CODAR proprietary software (Seasonde Radial Suite R7u2) (for specific cylindrical handling of HF radar data see Sciascia et al., 2018). For each RC, information on spatially averaged wave parameters (significant wave height, centroid period, and direction, referred to as H HF s , T HF c , and θ HF , respectively) were then extracted from the second-order spectrum by applying a Pierson-Moskowitz model (Lipa and Nyden, 2005). To estimate the directional ocean wave spectrum, SeaSonde assumes that the ocean wave spectrum is homogeneous over the radar RC (Lipa et al., 2018). For the analysis we used the results of the range cell located between 5 and 6 km from the coast (RC5), considered as the best operational selection, and we filtered the data to remove spikes and spurious data (details can be found in Saviano et al., 2019). In particular, only data relative to waves higher than half a meter were considered. Owing to intrinsic limitations in the physics behind the functioning of HF radars, when sea waves score H s < 0.50 m the second-order echo is too small compared to spectral noise, and cannot reflect the radar transmitted electromagnetic signal efficiently (Wyatt and Green, 2009;Wyatt, 2011;Atan et al., 2016;Lopez et al., 2016). Such condition leads to erroneous readings or data gaps (Basañez et al., 2020). As such, measurements below this lower threshold were discarded from the analysis. In the GoN, sea states with significant wave heights lower than 0.5 m typically occur during the spring-summer periods (Pugliese Carratelli and Sansone, 1987;Benassai et al., 1994a,b;Buonocore et al., 2003).

Buoy Wave Observations
A SEAWATCH Wavescan Buoy manufactured by Fugro OCEANOR (Trondheim, Norway) operated from August 2009 to mid-December 2010, and from mid-March to November 2012. The buoy (PC buoy, Figure 1) was installed at the southwestern limit of the basin of interest. It was equipped with a directional wave sensor and collected information about significant wave height (H BUOY s ), mean period (T BUOY m ), peak period (T BUOY p ), and directions (θ BUOY ), providing data every 30 min; more details can be found in Saviano et al. (2019). Data on H BUOY s , T BUOY m , and θ BUOY were used for the calibration of the regional and local models. The integration time of buoy measurements was 30 min.

Wave Numerical Modeling: From Basin to Local Scales
Information from the large scale hindcast (WWIII) were downscaled inside the GoN (SWAN) in order to obtain a reliable description of wave characteristics influenced by the bathymetry and orography of the GoN.

Regional Model: The Mediterranean Basin and WWIII
The regional hindcast was developed by the Department of Civil, Chemical and Environmental Engineering of the University of Genoa (DICCA,www.dicca.unige.it/meteocean/hindcast), with the WWIII numerical model (Tolman, 1989(Tolman, , 2009). The hindcast was done in the Mediterranean Sea from 01/Jan/1979 to 31/Dec/2018, with an approximate 0.1 • resolution both in longitude and latitude. The wave model was forced by 10 m wind fields obtained from the non-hydrostatic mesoscale model Weather Research and Forecast (WRF-ARW) version 3.3.1 (Skamarock et al., 2008), based on NCEP Climate Forecast System Reanalysis (CFSR), for the period from January 1979 to December 2010 and CFSv2 (Climate Forecast System version 2) for the period from January 2011 to December 2018. The output was hourly recorded at the nodes of the computational grid, providing the spectrum integrated quantities (H REG s , T REG m , θ REG m ), the peak wave features (T REG p , θ REG p ), along with the u REG w and v REG w wind velocity components (Mentaschi et al., 2013a(Mentaschi et al., , 2015Cassola et al., 2016). Though the information provided by the hindcast is discretized (i.e., one dataset for every node of the computational grid), wave features can be continuously defined within the boundaries of the model through a timespatial interpolation of the grid nodes outcomes. Nevertheless, the spatial resolution of the model may not allow to get the waves data at fine scales in the near-shore areas, as it actually happened for the GoN (see Figure 1B). In order to evaluate the wave features at the HF radar RCs, the SWAN model (Booij et al., 1999) was nested within the frame of WWIII. A finer regular grid was developed, with a 0.0025 • lon/lat resolution (almost 200 m at these latitudes), using a bathymetry of the area provided by the Hydrographic Institute of the Italian Navy (see Figure 1C). Then, hindcast data were used to set the wave conditions at the grid boundaries and to define the wind forcing over the whole domain. At this stage, it is worth mentioning that in SWAN sea states are modeled through a JONSWAP spectrum (Hasselmann et al., 1973), whose shape particularly depends on two parameters: the peak enhancement factor and the directional spreading (referred to as γ and θ s , respectively). The former parameter affects the sharpness in frequency of the spectrum, i.e. high values of γ imply that the majority of sea states energy relies in proximity of the peak period; the latter parameter has a similar meaning, but with respect to the waves incoming direction: low (high) values of θ s refer to sea states characterized by incoming directions weakly (highly) spread around the peak one. Therefore, in order to select the most performing SWAN setting, we tested four different sets: the values of γ taken into account are those characteristic of wind-waves dominated and swell dominated sea states (2.2 and 3.3, respectively), while θ s was varied between 10 • and 20 • . The values of the two parameters were varied at a time, resulting in four possible model settings. The modeled wave features due to each of the tested sets were than evaluated against the buoy measurements, in order to retain the most efficient γ /θ s combination (i.e., the one leading to the modeled wave features closer to the buoy data). A summary of the simulations is reported in Table 1.

STATISTICAL DESCRIPTORS FOR COMPARATIVE ANALYSES
As a first step, before comparing HF radar measurements with numerical model reanalysis, a validation of WWIII and SWAN performances was carried out against buoy data. Hence, we first evaluated the hourly hindcast outcomes at the buoy location within the period 25/Aug/2009-23/Dec/2010; at a second time, we compared the observations and the numerical simulations. The statistical parameters used for computing the correlation between the data are: where S i and O i are simulated and observed data, respectively; σ stands for the data standard deviation, whereas the overstanding bar indicates the average. ρ represents the correlation index, spanning in the ±1 range, where +1 (−1) indicates perfect correlation (anti-correlation) for two investigated series, while 0 indicates no-correlation. The statistical significance of the correlation is based on the criteria discussed in Rumsey (2016). NRMSE (Normalized Root Mean Square Error) is a measure of the distances between two datasets, and is one of the most employed indexes for evaluating the goodness of fit between the observed values of a variable and the correspondent values predicted by a model; then, substracting the average component of the error from the NRMSE, the scatter index (SI) is obtained. NRMSE and SI combine information on both the average and the scatter errors between two series, and should provide values as small as possible (ideally, if a perfect match exists, both the indexes should attain values equal to zero). NBIAS is a measure of the average component of the error, which should result in close-to-0 values for series in good agreement. slope measures the best linear interpolate for the resulting scatter of two series; under the hypothesis that two datasets are correlated, it attains positive values that may be lower or higher than 1, depending if the data to be validated are generally smaller or higher than the benchmark, respectively. Then, the index proposed by Hanna and Heinold (1985), referred to as HH, was employed. HH allows to overcome drawbacks that may arise for the simulations negatively biased, i.e. that underestimate the quantities measured (see Mentaschi et al., 2013b, who illustrate that in this case NRMSE and SI cannot be reliably used). Finally, when dealing with circular quantities such as the wave directions, the aforementioned indexes cannot be directly computed. In this case, an analysis should be performed on the differences between the simulated and the observed data, previously normalized in the [−π; π] space (i.e., in the worst case two waves may happen to be antiparallel). The directional indices have been therefore adjusted as follows: where the modulus operator mod −π ,π implies to subtract a 2π angle if θ i > π, on the other hand, if θ i < π a 2π angle is added to the difference.
The same test was performed for the hindcast data with respect to the HF radar data at PORT and SORR for the period 01/Jan/2010-30/Dec/2012. No analysis was available for CAST, as here the radar measurements fall outside the hindcast domain (see Figure 1B). Subsequently we analyzed the results linked to the SWAN model through the same procedure. In this case, the comparison between the HF radar and the model data was also carried out for CAST that is located within the SWAN domain. For the comparison, the HF data were averaged on 1 h. A Taylor diagram (Taylor, 2001) was also employed in order to select the best setting among the tested ones (see section 2.2.2).

Validation of Numerical Models Against Buoy Measurements
Comparisons between WWIII and buoy data are shown in Figure 2. Correlation indexes (ρ) show values close to one, whereas the error indexes are characterized by close-to-zero values, indicating coherence among measurements, especially with respect to the significant wave heights (H s ) and the mean periods (T m ). The scatter plots show that only in few cases the parameters related to the different sources are significantly diverging from the 45 • line (which indicates perfect correlation). Indeed, the color in the scatter plots is scaled according to the frequency of occurrence, i.e. the darker colors refer to conditions that seldom take place either in the buoy and WWIII outcomes, while the brighter colors are related to sea states more frequently encountered in both the numerical model and the buoy measurements. A higher uncertainty characterizes the mean directions, which can be anyway considered consistent as most of the computed θ (about the 70% of data) lie in the [−45 • ; 45 • ] range. As for the local scale analysis (i.e., the wave downscaling in the GoN with the high resolution model SWAN), Table 2 shows the correlations for each trial set of the model, evaluated versus the wave heights of the PC buoy, H s being taken as the leading parameter (data not shown). Results over the whole series are sensibly improved with respect to the regional scale ones. However, the values of the error indices suggest that the resulting wave features are weakly sensitive to the selection of the peakedness factor γ and the directional spreading θ s . This is further confirmed by inspection of the Taylor diagram: as Figure 3 shows, distributions relative to the different simulations almost overlap each other. If an overview of the total statistics does not reveal a significant dependence on the peak enhancement parameter and the spectrum directional spreading, zooming to a particular wave height series shows how different choices for these parameters may actually affect the model outcomes. Generally, θ s equal to 20 • results in wave heights slightly higher than the 10 • ones: the modeled wave heights can be therefore either closer or farther with respect to the buoy ones, depending on the sea state under investigation (i.e., if the sea states are overestimated or underestimated with respect to the sampled ones, see Figure 4). However, these differences are damped in the overall simulations, which embed almost 2 years of hourly sampled sea states.
For this study the following SWAN simulations were performed setting γ and θ s equal to 2.2 and 20, respectively, according to values typically suggested for the Mediterranean   (Figure 6). Finally, it can be noticed how the scatters of the distribution of H s and T is low-bound by the thresholds set to filter the radar data, as explained in section 2. In terms of θ m , the agreement between HF radar and model data is remarkable in PORT, while in SORR the results are more scattered, but mostly within [−50 • ,50 • ]. As for the mean wave directions, it should be mentioned though how it is not trivial to accurately model them, especially for moderate sea states (as shown by the color bar of Figure 2). As such, this could add further uncertainty in the comparison between the modeled and the HF radar sampled θ . Figures 7-9 show the comparison of the H s time series of the HF stations with the simulations. At PORT location, downscaling the wave field over a high resolution grid leads to H s values generally higher than those computed with the regional model, resulting in a general overestimation of the SWAN performance at this location. Furthermore, at CAST and SORR locations in some sea storm events an overestimation of HF acquisitions is recorded. Figures 10-12 show the scatter plots of wave parameters (H s , T, θ ) between the radar network acquisition and the simulations performed with SWAN. In PORT (Figure 10), the correlation ρ remains substantially unchanged, reaching values close to 0.70 for both comparisons. This implies that sea states characterized by similar relative intensities among the respective series occur at the same time. Nevertheless, results are more scattered with respect to the main bisector, in particular for the most intense sea states whose H s are generally higher in the local model. This can be seen in the decreasing/increasing values of slope/NBIAS (from 0.87 to 0.70 and from −0.05 to −0.22 for the regional and local models, respectively). The overall increment of the sea state magnitudes leads to better comparisons for the local model mean periods, with lower values NBIAS (which changes from 0.24 to 0.17) and slope (from 1.22 to 1.13). The wave incoming directions also show a slightly increasing deviation with respect to those of the regional analysis, with the NBIAS index decreasing from −3.70 to −5.26. As for SORR (Figure 11), a direct comparison with the regional counterparts underlines how the local numerical model leads to a higher consistency with the HF radar wave measurements for all the parameters under investigation, with the error indices (NRMSE, NBIAS, SI, and HH) switching to lower values, while ρ increases from FIGURE 3 | Taylor diagram of H s between PC buoy and SWAN for different γ and θ s tested in the model implementation (see Table 1). The panel on the right hand side zooms into the Taylor diagram to better show the overlap of the points.

Local Model
0.65 to 0.71 and slope reaches a value closer to 1 (1.41 vs. 1.52 of the regional analysis). These considerations equally apply to T and θ . Finally, in CAST (Figure 12), results present a reasonable consistency for all the investigated parameters. As for SORR, the modeled H s are slightly underestimated with respect to those sampled by the HF radars, though in this case the scatters are even closer (see, e.g., the positive NBIAS and slope, The first subrange corresponds to the H s values most frequently scored in the GoN, the second one to rough sea and storms, while the latter one to extreme, particularly severe cases. This was verified at yearly scale in Saviano et al. (2019) but also at pluriannual scale in a more recent work . As a general outcome, the statistical descriptors are more robust in the range 0.5 − 2.0 m. PORT and SORR scored the lowest confidence with the model at H s > 4.0 m, a result that could mirror the lower performance of HF radars at values beyond this critical value. For CAST the lowest confidence was found in the 2.0−4.0 m interval; this result is inconsistent with the other sites, and may be biased by a limited sensitivity of the radar due to infrastructural issues (see Supplementary Material).

DISCUSSION
A performance assessment of wave measurements was carried out for the GoN by comparing HF radar measurements and model data over a three-year period. A network of three HF radars was operational in the area of interest, providing hourly data of the main characteristics of the wave field. In a previous work (Saviano et al., 2019), a year-long assessment of HF radar derived wave measurements was qualitatively carried out by FIGURE 4 | Comparison of H s between PC buoy and SWAN simulations during a sea storm. The plot reports model performance using set1 and set2 cases (see Table 1).
comparison with simultaneous recording of a not co-located buoy with low correlation mainly at PORT and CAST locations (0.61 and 0.58, respectively). In the present work, the performance of the HF radars was further investigated by comparing wave measurements with model outputs over a 3-year period. The comparison between data at same location of HF radar and models improve the results obtained in the previous study. The numerical simulations were developed through two different models: (i) WWIII, used with a coarse resolution over the investigated area; and (ii) SWAN, used to downscale the information of WWIII over a finer resolution. Both models were validated through a direct comparison between their outcomes and data sampled by a wave buoy, moored at the edge of the GoN.    Figure 2 proved how the regional model WWIII is able to provide reliable outcomes for H s and T m , while θ m require further deepening. The same consideration can be reported for the SWAN simulations, developed by employing the regional wave information off Capri island and propagating them to the buoy site. In this case, different settings have been tested, varying the parameters most likely affecting the local analysis: the peakedness γ and the directional spreading θ s , which rule the shape of the modeled wave spectra along the frequency and the wave directions, respectively.   The simulations linked to the trial sets showed that the results are not significantly affected by the selection of parameters (see Table 1 and Figure 3), resulting in sound comparisons in all cases. The reason is related with the features investigated: the mean parameters are in fact computed by integrating the 2Dspectrum in frequency and direction. The shape of the spectrum is not affecting its underlying area as much as it may actually do for single points referred to specific characteristics (i.e., the peak  The comparison between WWIII and HF radar measurements yields a statistical agreement at PORT (Figure 5), while at SORR the wave features of the two platforms show more significant deviations (Figure 6), although still being consistent. Previous studies on the comparison of mean wave direction and peak period show that the detected discrepancies are linked to lower sea states, where directional spectra can be contaminated by noise due to spurious features such as those associated with antenna side lobes, or to the sensitivity of these parameters to non-sea signals (i.e., interference or ships) in the radar backscatter Lorente et al., 2018).

Results shown in
Looking at wave direction in SORR (Figure 6, right panel), the majority of θ HF values is higher than θ REG m , i.e. waves sampled by the HF radar are rotated to the East with respect to the model ones. This could be due to the coarse resolution of the WWIII grid, not allowing to efficiently resolve the refraction induced by the southern cape of the GoN for the waves coming from S-SW (see Figure 1B).
When zooming into the local model (SWAN), the comparisons between the measured and the modeled waves lead to an improvement in SORR (Figure 11). In WWIII, the wind growth effect modeled over the wide GoN may have led to over-estimated intensities in the wave magnitude, resulting in more scattered comparisons especially for the H s parameter (see Figure 10, right panel). On the contrary, in SWAN the higher resolution of the computational grid allows to better resolve the refraction of SW/NE waves around the cape of the basin.
The results from CAST are consistent with those recorded for the other two sites. In this case, however, no comparison is possible with hindcast data, as this sub-basin is not covered by WWIII.
The additional comparisons of the Hs retrieved by the radars and by SWAN over three Hs ranges (0.5-2.0, 2.0-4.0, and > 4.0 m, in Supplementary Material) sheds additional light on "windows of confidence" in the evaluation of the agreement between the two platforms. Generally speaking, the correlation between the radars and SWAN is higher at lower H s , but drops as H s increases. This indicates that, on a general basis, the comparison between the two platforms is robust over the ranges of H s most frequently encountered in the GoN, while some flaws emerge in the case of extreme waves. The analysis of this outcome points to different plausible mechanisms leading to such discrepancies. From the radar point of view, a crucial issue is represented by the discontinuity in functioning (Wyatt et al., 1999), requiring advances in background theory (Wyatt, 2011) and in software implementation for wave retrieval (Lopez et al., 2016). The lower bound set by waves with H s = 0.5 m does not seem to impair HF radar efficiency in comparison to wave models. More critical seems to be the detectability of waves above 4.0 m. Although previous studies from the GoN show a consistent pattern between the radars and the buoy in harsh conditions (Saviano et al., 2019), further investigations will be needed to assess if the H s overestimation is explained by the inversion method used, currently based on an ideal Pierson-Moskovitz spectrum. This is sufficiently sound to describe unimodal energy spectra from fully developed seas, though future implementations of the proprietary software will be capable of handling bimodal distributions as well (Lipa et al., 2018). Promising improvements along this line are also provided by the development of modified inversion methods (Lopez and Conley, 2019) and by the use of neural networks (Hardman and Wyatt, 2019). In addition, a re-processing of the HF radar data with new releases of the proprietary software may return more consistent results, as discussed in Basañez et al. (2020).
From the model perspective, it should be considered that both WWIII and SWAN have been validated using the PC buoy, located outside the coverage area of the radars and working at a lower frequency compared to the HF radars. It might be possible that the models performance could be affected by spatial variability, reducing the validity of models in the coastal area. Further validations will be performed when in situ data from coastal buoy will be available for the GoN. The wave field in the shadow zone behind the island of Capri may be more difficult to resolve, because of the complex morphology of the coastline. Some of the deviations for modeled low sea states associated with high sea states sampled by the radars could therefore be explained by looking at their positions along the coastline, with respect to S-SW incoming waves. The problems related to complex morphological features of the GoN were underlined also by Inghilesi et al. (2016), using both WAM and SWAN models. The results also indirectly show coherent values for the wave features modeled by WWIII and SWAN, even if the simulations were produced with different approaches (i.e., time-spatial regression for the regional analysis, physical wave propagation for the local analysis). Similarities in the modeled fields are expected, because the bottom depths corresponding to these RCs are deep enough for the large waves not to be dramatically affected by the bathymetry of the area, i.e. the sea is almost 100-m deep at both PORT and SORR.
In the present work the comparison between HF radars and wave models was carried out exclusively on RC5. This procedure was justified considering previous results carried out in the GoN (Falco et al., 2016;Saviano et al., 2019Saviano et al., , 2020De Leo et al., 2020) showing that RC5 guaranteed the optimal balance between data retrieval (minimum data gaps) and distance from the coast (sufficient echo intensity, deep water conditions). Moreover, this RC can be used as a robust estimator of the wave filed dynamics developing in the farthest RCs (Saviano et al., 2019). Future implementations will include the analysis of any potential variability in the wave field by performing radar-model comparisons on a larger number of RCs. Further investigations are in progress using satellite altimetry data and buoys recently deployed inside the zone covered by the radars to be able to understand and solve the shortcomings in the comparison. The datasets of the different platforms employed in this study (PC buoy, HF radars and model simulations) are coherent with the wave climate of the Tyrrhenian Sea and with the different dynamics taking place in both the internal and external parts of the GoN (Benassai et al., 1994a,b;Buonocore et al., 2003;Falco et al., 2016;Morucci et al., 2016;Saviano et al., 2019).
It is worth mentioning that the performance of HF radars is limited by numerous sources of error (e.g., environmental noise, interpretation methods, basin structure; Laws et al., 2010). Nevertheless, the results presented here suggest that the HF radar wave measurements are reasonable and reliable at least as a first approximation of the sea state. However, it is necessary to implement new algorithms to broaden the range of wave parameters observed by HF radars. Moreover, the possibility of integrating HF radar data into wave models (Caires and Wyatt, 2003;Siddons et al., 2009;Waters et al., 2013) can be a useful tool for optimizing the initial conditions of a model and therefore improving the accuracy of its estimates. Some assimilation schemes with WWIII, SWAN and WAM are proposed in literature (Caires and Wyatt, 2003;Siddons et al., 2009;Waters et al., 2013), and currently represent a promising scientific perspective for the near future.

CONCLUSIONS
A 3-year (2010-2012) characterization of the wave height, period and direction in the GoN was carried out by comparing HF radar derived wave parameters with WWIII and SWAN model predictions. The comparison allowed a characterization of the range of wave variations in the region, adding value to the HF radar wave estimates. Results demonstrate that HF radars might be employed as remote sensing tools to retrieve wave parameters in coastal areas, even though the signal should be further analyzed and possibly validated against wave measures provided by other devices; equally, skill metrics indicate that WWIII and SWAN performances were consistent with a wide range of sea states. These results show that in coastal areas a synergistic observational and modeling approach can provide a complete characterization of the wave conditions suggesting possibilities to expand the wave monitoring networks.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
SS, FD, MU, and GB: methodology and writing review and editing. FD and SS: software, validation, and data curation. GB, MU, and EZ: project administration. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
The Department of Science and Technology (formerly the Department of Environmental Sciences) of the Parthenope University of Naples operates the HF radar system on behalf of the AMRA consortium (formerly CRdC AMRA), a regional competence center for the analysis and monitoring of environmental risks. Our radar remote sites are hosted by the ENEA Centre of Portici, the Villa Angelina Village of High Education and Professional Training and La Villanella resort in Massa Lubrense, and the Fincantieri shipyard in Castellammare di Stabia, whose hospitality is gratefully acknowledged. The authors gratefully acknowledge the Civil Protection Department of the Campania Region for providing buoy data, the Port Authority of Capri Island for general information on the installation of the PC wave buoy, and the CAE company for providing technical information on sensors installed on the PC wave buoy. The authors would also like to thank the Reparto Geofisica Marina e Oceanografia of the Istituto Idrografico della Marina (IIM) for providing us with high resolution bathymetry employed for the numerical simulations in the GoN. Mention of trade names or commercial products does not constitute endorsement or recommandation. The authors thank the editor and the Reviewers for constructive feedback.