Gulf of Naples Advanced Model (GNAM): A Multiannual Comparison with Coastal HF Radar Data and Hydrological Measurements in a Coastal Tyrrhenian Basin

Abstract

High-resolution modelling systems have increasingly become an essential requirement to investigate ocean dynamics over a wide range of spatial and temporal scales, and to integrate the punctual ocean observations. When applied in coastal areas, they also have the potential to provide a detailed representation of transport and exchange processes at the sub-basin scale. This paper presents a validation exercise between the surface fields generated by the regional ocean modeling system (ROMS), developed for the Tyrrhenian Sea and downscaled for the Gulf of Naples (GNAM Gulf of Naples advanced model), and a 4 year-long (2009–2012) record of high-frequency radar (HFR) data. The comparison between hourly and seasonal model results and HFR surface fields is focused on the Gulf of Naples (GoN), where an observational network of three HFR sites has been operational since 2004, and on a specific subdomain characterized by the presence of the Sarno river, a long-term ecological research station (LTER-MC) and one important canyon area. An evaluation on a transect delimiting inshore–offshore zones in the GoN is also presented. The GNAM model was also compared with in situ hydrological parameters of temperatures and salinities retrieved at the LTER-MC fixed monitoring station. According to the skill metrics, basic circulation features are accurately reproduced by the circulation model, despite some model drawbacks in terms of increment of energy content in the surface current field occurring during specific seasonal events. The results allow us to identify potential model errors and to suggest useful improvements, the outcome also confirms the unique capability of HF radar systems to provide fine-scale measurements for the validation of numerical models and to counterbalance the lack of high-resolution measurements in coastal areas.

1. Introduction

The ocean circulation models can be used as tools to support maritime policy and assist high-stakes decision-making related to marine safety, water quality, and mitigation of both natural disasters and anthropogenic hazards in the vulnerable coastal environment [1]. Applications of ocean circulation models in coastal areas require a high-resolution grid to respond to the demand of estimates of oceanic variables at a fine scale, so numerous challenges and issues must be addressed in the initial and boundary conditions, sub-grid scale parameterizations, topography, and forcing [1,2]. The improvement of regional circulation models in the last decade due to the increasing necessity to represent fine scale processes and to provide operational oceanography products, allows today reproducing the variability of the main currents as well as small-scale features of the coastal dynamics such as filaments, coastal eddies, or river plumes covering a wide range of spatio-temporal scales, which subsequently have a variety of applications (e.g., [3]).

This improvement in ocean processes representation is due to the synergies between models and ocean observations [4,5] accomplished via comparison [6,7] or data assimilation [8,9,10] schemes.

In this context, invaluable instruments for a synoptic, basin-scale observation of oceanographic processes are in the HFR system capable of generating fine-resolution fields of the surface currents over broad coastal areas, making them particularly appropriate for the validation of numerical models [11,12,13]. In fact, in different areas of the Mediterranean Sea, many HFR systems have been used with this purpose, such as in the Ebro Delta area [14,15,16,17,18], the Ibiza Channel [1,15,17,19], the Northern Adriatic [20], and the Strait of Gibraltar [15,21]. A synergistic approach based on the integration of numerical models and observational networks can provide its usefulness to characterize the highly dynamic coastal circulation and the related complex interactions.

In this work, we focus on the validation of GNAM (Gulf of Naples Advanced Model), the regional ocean modeling system (ROMS) developed for the Tyrrhenian Sea and downscaled for the Gulf of Naples (GoN). To this aim we developed a multi-annual comparison and analysis of the surface fields generated by GNAM with a record of HFR data for a 4 year-long period (2009–2012).

2. Gulf of Naples

The GoN is a marginal sub-basin of the Tyrrhenian Sea, characterized by an intricate bottom topography (including two deep canyons) and orography (Vesuvius volcano, hills in the urban area of Naples city, Lattari mountains in the south) (Figure 1). The area is subject to intense anthropic pressure, but at the same time it is particularly rich in terms of biodiversity [22,23] with the presence of habitats with a high value of natural capital stock [24]. The surface circulation of the GoN is mainly driven by the local wind field, with the creation of specific patterns affecting both the physical and biological processes developing in the area [25]. During winter months, the predominant and most intense winds come from NNE and NE directions, which determine off-shore transport and sea storms [26,27,28].

In spring and autumn, the principal wind directions are NE and SW, the surface currents associated are characterized by the presence of recirculation structures, with both cyclonic and anticyclonic gyres at basin and sub-basin scales [26,29,30].

Finally, in summer, the setup and reinforcement of the Azore anticyclone, and in the latest years of the African anticyclone [31] determine the onset of a stable, moderate breeze system with winds generally weaker than the rest of the year [32].

This general wind regime is associated with numerous low-pressure systems passing over the basin, with frequent stormy and windy events principally during winter and autumn [33,34].

3. Materials and Methods

In this work, we use a free-run hindcast simulation spanning the period 2009–2012 of the ROMS–GNAM model and HFR data that provides high-resolution synoptic observations of sea surface currents, with a limited geographical coverage but high spatiotemporal resolution, similar to that of the high-resolution model. The HFR hourly derived surface current maps are provided at 1 km spatial resolution.

3.1. Numerical Model

Numerical model simulations were performed using the ROMS, a 3D free-surface hydrostatic primitive-equation finite-difference model, widely used by the scientific community for a diverse range of applications [35,36,37,38,39,40,41]. On the vertical, the primitive equations are discretized over variable topography using stretched terrain-following coordinates, with the stretched coordinates allowing increased resolution in areas of interest. In the horizontal, the primitive equations are evaluated using boundary-fitted, orthogonal curvilinear coordinates on a staggered Arakawa C-grid. Both the vertical and the horizontal stencils utilize centered, second-order finite differences. Advection is evaluated using a third-order, upstream biased scheme.

A one-way nesting has been used in order to obtain high resolution for the Campania coastal area. Simulations are first performed on a grid covering the whole Tyrrhenian Sea, with a 2 km resolution (Figure 2a), using as initial and boundary conditions the data of the Mediterranean Forecasting System, that have a resolution of approximately 6–7 km and are available with a daily frequency [42]. The results of this first simulation are then used as initial and boundary conditions for a finer grid, covering the Campania coast with a 500 m resolution (Figure 2b). All simulations on the finer grid are made using 30 vertical levels, with a resolution at the surface going from less than 10 cm to almost 2 m, depending on the local depth, as typical of terrain-following coordinates.

Boundary conditions on open boundaries are exactly as proposed by [43], with an adaptive algorithm where inward and outward information fluxes are treated separately. Because of the essentially hyperbolic nature of the incompressible, hydrostatic primitive equations, external data are required only for inward boundary fluxes. The outward fluxes are treated with an algorithm for two-dimensional radiation, with an adaptive nudging that adequately incorporates the external information minimizing over- and under-specification problems.

The turbulence model adopted in the vertical direction is the K-profile parameterization (KPP) introduced by [44], that is a first order scheme that has been shown to perform well in open ocean settings. The scheme matches separate parameterizations for vertical mixing of the surface boundary layer and the ocean interior. A formulation based on boundary layer similarity theory is applied in the water column above a calculated boundary layer depth. This parameterization is then matched at the interior with schemes to account for local shear, internal wave, and double diffusive mixing effects.

In the horizontal direction a constant value is given to the diffusion coefficient, as usually done in this type of application, with the value being chosen for its capability of dumping spurious values without introducing an excessive smearing of the solution. In all simulations a value of 10 m² s⁻¹ has been used for velocity, temperature, and salinity.

Surface forcing has been calculated from atmospheric values using the bulk parameterization of [45], based on the coupled ocean–atmosphere response experiment (COARE) algorithm. This algorithm computes bulk air–sea fluxes, including latent and sensible heat flux, net heat flux, and other associated fluxes including evaporation, evaporation minus precipitation, sensible heat flux due to rain, buoyancy flux, and wind stresses, and requires input atmospheric parameters at the ocean surface such as wind velocity, temperature and pressure, cloud cover, rain rates, and surface net solar radiation.

Atmosphere parameters have been obtained from “ERA-Interim” data [46] of the European Center for Medium-Range Weather Forecasts (ECMWF). ERA-Interim is a climate reanalysis dataset that uses a fixed version of a numerical weather prediction system to produce reanalysed data, where observations are blended or assimilated with a previous forecast to obtain the best fit to both. River inflows have been modeled as point sources, with climatological daily flow rate data obtained from a variety of sources, in particular from Italian state agencies in charge of monitoring the environment.

3.2. High-Frequency Radar

The HFR provides high-resolution synoptic observations of sea surface currents, with a limited geographical coverage but high spatiotemporal resolution, similar to that of the high-resolution model.

HFRs have proved to be excellent tools to provide measurements of coastal currents (e.g., [47,48,49,50,51]). They provide current fields on a spatial regular grid that are also well suited for validating coastal circulation models (e.g., [52,53,54]) and for being assimilated into them (e.g., [55,56]).

The high spatial and temporal resolution surface current data utilized in this study have been provided by such a system. In 2004, the first core of a network of HFR coastal radars was installed in the GoN, permitting real-time synoptic monitoring of the surface current field at the basin scale. This system has been operated by Università degli Studi di Napoli “Parthenope”. HFR transmits an electromagnetic wave towards a target and measures the time taken by the reflected echo to return to the antenna. A coastal radar system emits electromagnetic waves in the 3–30 MHz band and uses gravity waves propagating over the sea surface as a target. When a Bragg scattering coherent resonance occurs [57], the backscattered signal shows a clear peak in the signal spectrum [58]. If a surface current field is present beneath the gravity waves, the peaks in the backscattered signal will be shifted due to the Doppler effect [59], the inversion methods provide information on the direction and intensity of the surface current field. As each HF radar antenna only measures the radial component of surface velocity, at least two transceiving antennas are required [60,61,62]. The system installed in the GoN (Figure 1) is a SeaSonde type manufactured by CODAR Ocean Sensors (Mountain View, CA, USA). It works in the 25 MHz band and measures surface currents relative to the first 1 m of the water column [63]. The temporal resolution of the system is 1 h, while the range is approximately 35 km from the coast. In this configuration, the spatial resolution was 1 km. The HFr network worked continuously in the years 2009–2012, except for a failure at the CAST site from July 2011 to November 2011.

3.3. LTER-Mare Chiara

The Long-Term Ecological Research Program Mare Chiara (LTER-MC) is one of the longest and most detailed plankton monitoring activities in the Mediterranean Sea, with regular sampling since January 1984. The sampling site LTER-MC (40°48.5′ N, 14°15′ E) is located over a depth of approximately 73 m, two nautical miles from the coastline in the GoN, a relatively deep embayment (average depth ~170 m) along the coast of the southeastern Tyrrhenian Sea (Figure 1). Following a series of oceanographic campaigns in the inner GoN in summer 1983 [64], the site was selected in an area that receives municipal inputs from one of the most densely populated Mediterranean coastal regions. At the same time, the GoN is also influenced by the dynamics of the offshore oligotrophic Tyrrhenian Sea waters and occasionally by water masses from the adjacent Gulf of Gaeta [10].

Sampling at LTER-MC took place every fortnight until 1991 and weekly from 1995 to date, with a major interruption from August 1991 to February 1995. In this study, we used the temperature and salinity data from January 2009 (cast MC845 on the 7 January 2009) to December 2012 (cast MC1037 on the 27 December 2012), computed then with the Gibbs Sea Water (GSW) Oceanographic Toolbox to calculate the conservative temperature and the absolute salinity [65].

3.4. Statistical Methods

A preliminary investigation on the spatial coverage of the HFR data in the years (2009–2012) used for the analyses was carried out. In the comparison between HFR and GNAM, HFR data with at least 70% of coverage during the year were used. Figure 2c shows the percentage of coverage of the year 2012.

The statistical metrics used in the present study to compare two data sets x = {x₁, x₂, x₃ … x_N} and y = {y₁, y₂, y₃ … y_N} (where 1, 2, … N represent the time steps) include the mean ($\underset{\_}{x}$), model mean error (bias), scalar correlation ($\rho$), and root mean squared error (RMSE), defined, respectively, as:

$$\underset{\_}{x}=\frac{1}{N} {\displaystyle \sum }_{i=1}^N{x}_i$$

(1)

$$\rho =\frac{1}{\left(N-1\right)}{\displaystyle \sum }_{i=1}^N\left(\frac{x_i-\underset{\_}{x}}{{\sigma }_x}\right)\left(\frac{y_i-\underset{\_}{y}}{{\sigma }_y}\right)$$

(2)

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(u_i^{\left(x\right)}-u_i^{\left(y\right)}\right)}^2}{n} }$$

(3)

where ${\sigma }_x$ and ${\sigma }_y$ are the standard deviation of x and y, respectively.

$\rho$ gives the overall measure of correlation, ranging from 0 to 1.

Energetics analysis has been often used to investigate the temporal variability of ocean currents and eddies (e.g., [66,67,68,69,70]). The eddy state is commonly defined as the deviation from the time-mean state [71]. For a given time-mean state, the flow velocity is split into the time-mean (mean flow) and time-varying (eddy) parts. Accordingly, the total kinetic energy (KE) at a given time, is decomposed into the mean flow and eddy parts as well as a residual term that carries information of both flows. As we do not have currents data in the water column to compare to the model, but only the observations on the surface of the currents from the HFR data, we calculate an indicator of the kinetic energy from the time average velocity u and v [72]. In our study we define this quantity as:

$$KE=\underset{\_}{E_k}=\underset{\_}{\left(u^2+v^2\right)}$$

(4)

To characterize the rotational state of the surface currents at the scale of Gulf of Naples and catch the cyclonic and anticyclonic features of the coastal circulation, we calculate an index of relative vorticity [73] between boxes located symmetrically at approximately 10 km north, south, east, and west from the center of the basin (see Section 4.2). Daily velocities are averaged in each box, to calculate the vorticity index as:

$$VOR=\frac{\left(\underset{\_}{v_{East}}-\underset{\_}{v_{West}}\right)}{\Delta X}-\frac{\left(\underset{\_}{u_{North}}-\underset{\_}{u_{South}}\right)}{\Delta Y}$$

(5)

where ΔX and ΔY are ≃20 km. Direct (indirect) trigonometric directions define the positive (negative) values, indicating a cyclonic (anti-cyclonic) circulation.

4. Results

This section will present the comparisons between the GNAM model and the HFR data. The analyses were initially carried out on the entire GoN. Subsequently, areas of particular interest were identified (three close to the coast and one offshore); finally, an NS transect was defined in the GoN to obtain information on exchanges between the inner and outer part of the basin.

4.1. Time-Averaged Surface Circulation

The performance of the GNAM model is assessed in terms of surface currents in the coastal areas monitored with HFR. The analyses were carried out on the entire GoN using the points where the HFR coverage was greater than the 70%.

Table 1. In table are reported the $\rho$ of GNAM and HFR current component (u, v) in the GoN (with 70% of coverage in HFR data).

Year	2009		2010		2011		2012
Months	u	v	u	v	u	v	u	v
Jan	0.82	0.71	0.97	0.92	0.51	0.56	0.89	0.79
Feb	0.62	0.84	0.68	0.66	0.7	0.83	0.65	0.76
Mar	0.81	0.91	0.87	0.79	0.88	0.84	0.81	0.57
Apr	0.83	0.91	0.76	0.81	0.77	0.25	0.78	0.73
May	0.58	0.43	0.87	0.72	0.78	0.69	0.65	0.62
Jun	0.38	0.72	0.76	0.47	0.60	0.51	0.22	0.57
Jul	0.67	0.69	0.74	0.50	0.17	0.68	0.48	0.47
Aug	0.68	0.14	0.51	0.39	0.76	0.45	0.67	0.27
Sep	0.81	0.49	0.64	0.53	0.57	0.20	0.32	0.39
Oct	0.41	-0.06	0.58	0.60	0.57	0.008	0.57	0.85
Nov	0.31	0.34	0.76	0.78	0.52	0.47	0.73	0.84
Dec	0.71	0.9	0.78	0.63	0.88	0.65	0.77	0.89
Mean	0.64	0.75	0.64	0.63	0.58	0.66	0.51	0.65

shows the monthly correlation coefficient ($\rho$) of the mean zonal (E–W) velocity component (u) and the mean meridional (N–S) velocity component (v) obtained from HFR and GNAM comparison over the year 2009–2012.

This analysis showed a good agreement between GNAM and HFR, with values of $\rho$ in the winter season between 0.51 and 0.97 for the u component and 0.56 and 0.92 for the v component.

Generally, low $\rho$ were found in summer and autumn, mainly for the v component, with a cc minimum value of 0.008 in October 2011, related to a failure of the HF radar (see Material and Methods).

The analysis over quarters divided according to the following scheme (JFM–AMJ–JAS-OND) were conducted according to the GoN climatology. We report the component comparisons of the velocity u and v, the kinetic energy, and the correlation coefficient between GNAM and HFR for the autumn quarters (OND) of all years (Figure 3a–d), as they represent the best response together with the winter period (JFM). All quarters are shown in the Supplementary Materials (Figures S1–S3).

As can be observed in Figure 3a,b, the HFR and GNAM seasonally averaged current patterns exhibit a close resemblance in terms of mean component (u and v) but differ due to a model overestimation of the current intensity, especially for the u component.

The data–model comparison shows a good agreement in terms of current variability, as reflected by the KE maps. A linear structure in the S–N direction can be observed during the years, with cores reaching a similar peak of energy.

The correlation index shows anticorrelation areas that vary during the year but interest mainly the Penisola Sorrentina for the u component and the offshore area for the v component. The v component is the most impacted, generally on summer and post-summer periods (August and October–November 2009, August 2010, April and September-October 2011, and August 2012). The bathymetry used in the model extends to 20 m deep (see Material and Methods), thus it may limit the modelization of the coastal margins in such areas where relatively abrupt topographical gradients exist. Moreover, radar coverage is not optimal in zones very close to the coast, where comparisons to the model are therefore more complex.

During the investigated period, a general underestimation in the kinetic energy plots was found (Figure 3c), the best results on KE comparison are shown in the quarter (OND), while the quarters (AMJ–JAS) show lower agreement.

4.2. Hovmöller Diagrams and Selected Areas

The analyses carried out throughout the GoN were repeated on areas of particular interest and investigations along an N–S transect were implemented to describe in detail the circulation and exchanges between the internal and external part of the basin (Figure 4). The selected coastal areas correspond to the area of the fixed hydrological monitoring station (LTER-Mare Chiara) in the northern part of the GoN. An area in front of the HF PORT station, and the area of Castellammare di Stabia, which is in front of the HF CAST station but also in correspondence with the mouth of the river Sarno, represent the most input of fresh waters in the GoN. Finally, an area corresponding to the deep area of the GoN, the Dohrn canyon, was investigated.

The transect goes from Capo Posillipo to Capri, crossing the entire area under investigation (red line in Figure 4). The Hovmöller diagrams were calculated over quarters for the entire study period, the OND quarters for the years of investigation are shown (Figure 5, Figure 6, Figure 7 and Figure 8). We performed this analysis to document the exchanges between the inner and external parts of the GoN. The results describe well the two opposite patterns visible in the northern and southern portions of the transect, as well as the changes in the orientation of the circulation, according to [6]. In general, the comparison of the spatial distributions of the u component and the v component between the model and HFR along the latitudinal section appears to be complex, but show a closer correspondence during the autumn and winter quarters, though there are apparent interannual variabilities presenting large differences in the model performance, showing repeated and abrupt changes of sign, as well as intermingling situations of meridional coherence of the zonal flux, finally the overall pattern was strongly time-dependent.

We can observe that the kinetic energy of the model is lower than the observations, and that the spatial–temporal differences between u and v are noticeable during some monthly periods (October 2009, December 2010, October and December 2011, and October 2012). The November period shows some better agreement, with the bimodal structures in u (distributed from north to south), and the alternance and intensification of v across the section. All the other quarters are available in the Supplementary Material (Figures S4–S6).

The time series of the velocity components and of the kinetic energy are shown for these four areas (averaged from daily to monthly) in Figure 9.

Generally, the coastal area (Figure 9a–c) shows a lower agreement in the v component and higher values of kinetic energy in HFR retrievals; instead, in the area of Canyon Dohrn (Figure 9d), the comparison shows agreement in the velocity components and kinetic energy.

The systematic underestimation of the velocity modulus in the areas closer to the shore could be due to the complete lack of waves modeling and therefore the absence of wave–current interactions. Indeed, waves induced effects particularly important in the near-shore areas, such as the transfer of kinetic energy to the water column due to breaking waves, which are completely missing. The coupling of ROMS with a wave model appears to be an important step for a future version of the model.

Another possibility is that the atmospheric data used to force the model may be the cause of the problem: the ECMWF data used for the evaluation of surface stresses may lack the necessary resolution to describe accurately the near shore areas, where variations in land topography could be strongly affecting the predicted winds, leading to an underestimation of the wind speeds. The ECMWF-interim data used in this work have a 1/8 degree resolution; higher resolution atmospheric data are currently sought to improve the model.

The possibility that the velocity underestimation may be due to an excessive damping caused by the horizontal diffusion has also been examined. The horizontal diffusion coefficient has in fact been set to an arbitrary value well above its theoretical value in order to provide damping of the numerical noise. Numerical tests, not presented here, have shown that this does not appear to be the cause of the velocity underestimation, as halving the diffusion coefficient does not provide any significant increase in the velocity module.

It should also be noticed that the value used in this work, 10 m² s⁻¹, is smaller than those used in similar applications that have appeared in the literature. For example, [74] modeled the circulation in the GoN with a resolution very similar to the one used here, and adopted a coefficient of 80 m² s⁻¹ after performing a series of sensitivity experiments.

4.3. Vorticity

Vorticity was calculated in the four locations identifying the north, south, east, and west areas around the center of the GoN (see Figure 4). The analyzes were carried out using a moving average on windows of 10, 30, and 90 days over the primary daily estimates, to identify monthly and three-monthly variability. The results are shown on Figure 10. In general, the observations (HFR) and the model (GNAM) both show a tendency toward negative values/anti-cyclonic circulation. The trimestrial smoothing (90 days, bottom) shows a good agreement, except during the spring period of 2010 when the model was more cyclonic and during the winters 2010 and 2011 when the model shows circulation schemes more anticyclonic with respect to HFR observations. The monthly smoothing (30 days, center) allows us to identify when opposed events occurred, with the tendency of the model to sometimes respond opposedly during some months, compared to the observations (e.g., November 2009, June 2010, or August 2012).

4.4. Hydrological Comparison with LTER-Mare Chiara

The GNAM model was also compared with in situ hydrological parameters of temperatures and salinity retrieved at the LTER-MC fixed monitoring station.

The temperature and salinity time series are reported both on the surface and on the bottom for the years of acquisition 2009–2012, Figure 11. Mean values in the water column are presented too, to represent the contribution of the vertical layers to the average (the Hovmöller diagrams of the vertical profiles of temperature and salinity, for both model and observations, can be consulted in the Supplementary Figure S7). To allow comparisons between the in situ hydrology (weekly sampling) and the GNAM model output (daily outputs), we interpolate the observations with a regular daily time step, and we smooth the model’s output with a running mean of 7 days.

Seasonality of the surface temperature seems to be well reproduced with the monthly peaks from spring to autumn that generally phase well between the model and the observations. When considering the averages in the water column, they indicate a good agreement between the model and the observations during the winter months (from December to March). Then, in general between April and October, the model shows a mean underestimate of around 1 °C to the in situ temperature. On the contrary, the deepest layer presents an overestimation of around 1 °C, and up to 2 °C. Looking at the water column bottom (60 m deep), the variability range of the model is more limited and shows a constant increase in temperature from March to October–November, which is absent from observations where the summer-to-autumn cycle is flatter. While surface stratified layers would reproduce better, the variability driven by vertical and horizontal mixing, ultimately controlling the circulation features budgeting the heat content and freshwater trade-offs, the bottom layer in the model appears to accumulate more/export less the thermal content of the water column. This excessive downward propagation of the surface signal during the restratification could be due to the model vertical mixing that is too high. As a consequence, we should expect an underestimate of the surface currents. This is not the case, as discussed above, which implies that at least partially the bias in the model surface currents are due to the wind forcing and, specifically, to the lack of representation of the coastal orographies in the atmospheric model grid.

In terms of salinity, surface layer frames relatively well the maximum values of salt (≃38 $\mathrm{g}·\mathrm{k}{\mathrm{g}}^{-1}$) during the autumn and winter periods, even if the variability lacks to reproduce the local disruptions, probably due to the interannual variability of the freshwater runoffs that are not taken into account in the model (climatological data were used here). The decrease of surface salinity is nevertheless visible during the spring periods, even the minimum values (≃37 $\mathrm{g}·\mathrm{k}{\mathrm{g}}^{-1}$) seen in the in situ observations are not reached. About the variability in the bottom layers, the inter-annual variations seem to generally follow the same trends as the observations. Intra-annual variations are not well reproduced, and in the same way as for temperature, this points out that the model underestimates the horizontal import/export of salty water parcels and their local mixing with the surrounding water mass.

5. Discussion and Conclusions

Multiple platforms collecting measurements of different variables at different scales provide various insights into the model performance. We have evaluated the high-resolution GNAM model using multi-platform observations, HFR, and a fixed monitoring station. The comparison of surface currents retrieved by the HF radar network and obtained by GNAM model on four years (2009–2012) of acquisition was presented. For the same period, a comparison of physical variables of temperature and salinity acquired by a fixed monitoring station located in the GoN was carried out.

The first comparison of surface current components showed a good agreement on the overall period, the monthly correlation coefficient resulted in larger than 0.5 in winter (January–March) for both components. With a maximum value in January 2010 of 0.973 for the u component and 0.917 for the v component, respectively. In all years of comparison, a good agreement was found in the winter months and a lower agreement in those close to the summer, mainly for the meridional component (v). This agreement was found also in the quarter analysis of the years, confirming the good representation of the model. The results show the reliability of the ROMS model to solve the surface dynamics and reproducing the seasonal patterns typical of the study area [6,10,25,26].

The results of the Hovmöller diagrams reproduce the structures and changes in direction of surface currents observed in a previous study [6], showing more complexity in the summer period. The analysis over quarters reproduces the two opposite patterns in the northern and southern portions of the transect, as well as the changes in the orientation of the circulation.

Even though the coastal ocean forecasting system presented here does not currently employ the assimilation of ocean variables, ROMS contains the tangent linear and adjoint formulations necessary for advanced assimilation techniques such as 3DVAR and 4DVAR data assimilation methods. Further improvement of the forecast may be achieved by assimilating data from the acoustic Doppler current profiler (ADCP) and HFR measurements, and the coastal model presented here lays the foundation for such improvements [75].

We show in our study that the surface temperature reproduces a seasonal and inter-annual cycle in good agreement with the observations, while the time series at the bottom of the water column at the LTER-MC coastal station shows a warming cycle that anticipates the observations, every year, and could be due to the mixing scheme implementation that could be modified with stronger turbulent diffusion rates. In general, both vertical and horizontal mixing are difficult to model properly in their amplitudes and dynamics, which could render complex the right representation of the heat budget in such a coastal area. This should require more direct in situ observations of the local turbulence to allow a finer tuning of the model.

In the recent joint studies about the marine ecosystem of the GoN [76], it has been shown that the hydrology at the LTER-MC coastal station was marked by a specific cycle in salinity via the establishment of the surface freshwater layer in spring. This fresh feature appears to be relatively well reproduced in the model with the minimum of surface salinity reached in April–May. Given that the model does not take in account the fresh runoffs but only the main rivers flow, the agreement in the variability for both model and observations suggest a realistic representation of the local horizontal advection of the rivers’ freshwater discharge (e.g., Sarno and Volturno rivers). The authors emphasized the observed delay between the maximum of rain (March) and the minimum of salinity (May–June), which could be driven by the mesoscale activity mitigating the inshore and offshore exchanges, as it is suggested by the vorticity index showing variations at the monthly time scale.

This index represents a very informative proxy to infer the local advection of the various fresh and salty inputs in the GoN by mesoscale features, as it was shown that the GoN area was prompt to be influenced by such structures [77]. A positive/cyclonic mode could favor the northward–westward export and distribution of the fresh runoffs from the Sarno river (located in the east corner of the GoN) along the coastline. On the contrary, a negative/anti-cyclonic mode should lead to another pattern of circulation, by redistributing the coastal runoffs more southward and towards the gulf’s center. Our study highlights this inner feature of the GoN, that should be driven by surface wind forcing. The surface response to such forcing could be then determined with the help of the GNAM model, that should allow to disentangle the respective influence of large-scale, mesoscale, and sub-mesoscale processes in setting the GoN surface and subsurface hydrology, through the trade offs with the open Tyrrhenian Sea [6,10].

During the post-summer and autumn season, the local hydrological cycle was characterized by a salty water layer at mid-depth in September–November, favorably hosting diffusive processes on the vertical through the establishment of a salt-fingering regime [23,63]. This represents a possibly significant input of salty flow, whose importance for biological communities is still to be assessed. The GNAM model is consequently of interest here, given the good agreement between the seasonal comparisons at this particular period, and the possibility of tracking back virtual particles as it includes implemented Lagrangian tools. More generally, the model paves the way for studies dedicated to the connectivity between the various areas of the GoN.