Investigation of model capability in capturing vertical hydrodynamic coastal processes: a case study in the North Adriatic Sea

In this work we consider a numerical study of hydrodynamics in the coastal zone using two different models, SHYFEM and MITgcm, to assess their capability to capture the main processes. We focus on the North Adriatic Sea during a strong dense water event that occurred at the beginning of 2012. This serves as an interesting test case to examine both the models 5 strengths and weaknesses, while giving an opportunity to understand how these events affect coastal processes, like upwelling and downwelling, and how they interact with estuarine dynamics. Using the models we examine the impact of setup, surface and lateral boundary treatment, resolution and mixing schemes, as well as assessing the importance of nonhy-drostatic dynamics in coastal processes. Both models are able to capture the dense water 10 event, though each displays biases in different regions. The models show large differences in the reproduction of surface patterns, identifying the choice of suitable bulk formulas as a central point for the correct simulation of the thermohaline structure of the coastal zone. Moreover, the different approaches in treating lateral freshwater sources affect the vertical coastal stratiﬁcation. The results indicate the importance of having high horizontal resolu-15 tion in the coastal zone, speciﬁcally in close proximity to river inputs, in order to reproduce the effect of the complex coastal morphology on the hydrodynamics. A lower resolution off-shore is acceptable for the reproduction of the dense water event, even if speciﬁc vortical structures are missed. Finally, it is found that nonhydrostatic processes are of little importance for the reproduction of dense water formation in the shelf of the North Adriatic Sea.


Introduction
Coastal hydrodynamic processes play an important role in ocean dynamics.Being at the interface between land and sea, they are strongly influenced by the input of freshwater through river discharge, tides, topographic features, as well as human activities and are affected, at the surface, by the winds and by heat and water fluxes.
The hydrodynamics typically observed in coastal areas involve processes interacting on a wide range of spatial and temporal scales, as well as slowly and rapidly varying features (Moum et al., 2008).The scale interaction seen in the coastal zone is driven by a number of local (wind, sources of freshwater) and large-scale (pressure, surface heat and mass fluxes) forcings.Along the coast, the surface wind affects the dynamics of freshwater from river inputs with different wind regimes causing the buoyant flow to narrow or thicken, leading to increased upwelling or downwelling (Magaldi et al., 2010).
Moreover, the interaction with coastal water bodies leads to the identification of "regions of freshwater influence" (ROFI hereafter) and interaction zones in the proximity of lagoons and transitional areas (Garvine, 1995).The presence of complex coastal morphologies, embayments, promontories, and sudden bathymetric changes can interact with coastal currents producing small-scale features (filaments) with specific temporal and spatial variation (Doglioli et al., 2004; Published by Copernicus Publications on behalf of the European Geosciences Union.Burrage et al., 2009).Islands, as well, can enhance smallscale features and are characterized by vertical movements during specific tide and wind conditions (Orlić et al., 2011;Orlić and Pasarić, 2011).In fact, the majority of horizontal structures observed in the coastal zone are characterized by Rossby and Richardson numbers of around 1 (submesoscale), representing areas of frontogenesis where vertical fluxes and buoyancy are enhanced (Thomas, 2008).In such a complex environment, sudden changes in the forcings can trigger strong hydrodynamic events, such as the formation of dense water (DW), wind-driven upwelling, and peak river floods.
Modeling the coastal zone and the specific hydrodynamic processes occurring there is challenging due to the number of spatial scales involved and the complex morphologies (Wolanski et al., 2003).In particular, these processes can produce strong vertical motions, which are difficult to model, requiring high resolution and an accurate representation of the underlying physics, perhaps even requiring the inclusion of nonhydrostatic processes.Improving modeling skills for reproduction of coastal processes is a balance between trying to capture the full range of physical processes involved (turbulence, mixing, non-hydrostatic vertical motion), while at the same time introducing suitable numerical approaches for efficient simulation of the processes.Modeling tools, with appropriate horizontal and vertical discretization, are needed (finite difference -finite volumes -finite elements; structured -unstructured grids).Also the choice in numerical parameterization schemes, particularly concerning vertical mixing, play a central role (Durski et al., 2004).
When modeling vertical processes, one issue to consider is whether nonhydrostatic processes are important for reproducing them.Several studies investigated this issue: Mahadevan (2006) studied the effect on submesoscale processes, stating the difficulty in identifying specific vertical features connected with nonhydrostatic process modeling.A major effect of the choice of resolution in pattern reproduction is stressed.Jiang et al. (2011) found that nonhydrostatic effects do not play a major role in coastal upwelling, but, interestingly, they identify their impact on the horizontal patterns (enhanced meandering).Magaldi and Haine (2015) stressed how the nonhydrostatic processes can affect the energy transfer between scales and they pointed out the need to investigate the possible role of nonhydrostatic processes in quantifying the modulation of scale interaction (on the horizontal) along the coast.
The Adriatic Sea is an example of a water body that is strongly linked to its coastal system, being a semi-enclosed basin with a particular topography, i.e., having a very shallow northern area becoming deeper towards the south, and a large number of freshwater sources (Russo and Artegiani, 1996;Cushman-Roisin et al., 2001).This makes it prone to DW events, when cold north-easterly winter winds induce water sinking in the shallow northern Adriatic (Vested et al., 1998;Vilibić and Supić, 2005;Mihanović et al., 2013;Vilibić and Mihanović, 2013;Durrieu de Madron et al., 2013).These extreme DW events have many complex influences and thus are particularly challenging to understand and model, though their impact on the general circulation has made them an important topic of research (Querin et al., 2013;Janeković et al., 2014;Benetazzo et al., 2014).Here we focus on one particularly strong DW formation event that occurred in the beginning of 2012.The extreme intensity of this DW event motivated many studies, with a large collection of in situ data, providing insights on the hydrodynamic features that occurred.Therefore, this case serves as an interesting test to assess the models, allowing us to compare our results with previous study efforts (Mihanović et al., 2013;Vilibić and Mihanović, 2013;Benetazzo et al., 2014), while having an opportunity to complement the understanding of how these events affect coastal hydrodynamic processes and, in particular, probing into what are the most suitable modeling strategies to reproduce them.
There are still many aspects of coastal dynamics that are not well understood and there are limits to the information garnered from in situ observations and measurement campaigns.Much about the dynamics of these processes must be studied through the use of numerical models.With this in mind, our approach here is to use two very different numerical models, SHYFEM (shallow water hydrodynamic finite element model) and MITgcm (Massachusetts Institute of Technology general circulation model), in order to compare their strengths and weaknesses in representing these processes.In particular we assess their ability to capture the DW event, its formation and propagation, as well as associated coastal upwelling and downwelling.In addition to a comparison between the two different models we also compare two simulations, one imposing hydrostatic balance, the other fully nonhydrostatic, in order to determine what impact nonhydrostatic processes have in regional coastal processes and DW phenomenon.
In Sect. 2 we describe the models used and simulation setup, as well as providing a list of observational data used for comparison with the models.In Sect. 3 we present the results, beginning in Sect.3.1 with a validation of the models against observational data.In Sect.3.2 we take a broad look at how the models represent the coastal dynamics, namely the DW formation and propagation, followed by an analysis of the coastal upwelling (Sect.3.3), and the impact of estuarine dynamics (Sect.3.4).We discuss the results in Sect. 4 and draw our conclusions in Sect. 5.

Methods
In this study we use two different three-dimensional (3-D) hydrodynamic models, SHYFEM and MITgcm.Both are designed for oceanographic studies and both have been previously applied in the open sea and in the coastal area of the Ocean Sci., 12, 51-69, 2016 www.ocean-sci.net/12/51/2016/Adriatic Sea (Bellafiore and Umgiesser, 2010;Querin et al., 2013).

SHYFEM
The SHYFEM model (Bellafiore and Umgiesser, 2010;Umgiesser et al., 2004)  Bottom stress is applied using a constant bottom friction coefficient (0.0025).We use a TVD (total variation diminishing) scheme for both the horizontal and vertical advection in the transport and diffusion equation for scalars, with constant diffusivity (0.2 m 2 s −1 ).Horizontal advection of momentum is discretized by an upwind scheme and horizontal eddy viscosity is computed by the Smagorinsky's formulation.For the computation of the vertical viscosities and diffusivities, a k-turbulence scheme is used, adapted from the GOTM (General Ocean Turbulence Model) model described in Burchard and Peterson (1999).
On the surface, a constant value for the wind drag coefficient is used (0.0025).To reproduce the surface heat fluxes, shortwave radiation from the atmospheric model is imposed, whereas the long-wave radiation is computed according to the Clark et al. (1974) formula.Bulk formulas are computed considering the sea surface temperature, the winds at 10 m height, the dry air temperature and the air pressure at 2 m, and the relative humidity as inputs.The latent heat flux and the sensible heat flux are computed according to the Kondo (1975) bulk formula.Cloud cover is taken from the atmospheric model.

MITgcm
The MITgcm solves both the hydrostatic and nonhydrostatic Navier-Stokes equations under the Boussinesq approximation for an incompressible fluid with a spatial finite-volume discretization on a curvilinear computational grid.The model formulation, which includes implicit free surface and partial step topography, is described in detail by Marshall et al. (1997a, b).The model domain, that covers the entire Adriatic Sea, is discretized by a non-uniform curvilinear orthogo-nal grid of 432 × 1296 points.The model has 100 vertical z levels with a thickness of 1 m in the upper 23 m gradually increasing to a maximum of 17 m for the remaining 64 levels.The bathymetry used by MITgcm is provided by the National Group of Operational Oceanography (GNOO; http://gnoo.bo.ingv.it/bathymetry/).As in Sanchez-Garrido et al. (2011) and Sannino et al. (2014), an implicit linear formulation of the free surface is used.The model uses constant horizontal eddy coefficients for momentum (viscosity: 10 m 2 s −1 ), temperature, and salinity (diffusivity: 2 m 2 s −1 ).Vertical eddy viscosity and diffusivity coefficients are computed in the MITgcm using the turbulence closure model developed by Bougeault and Lacarrere (1989) for the atmosphere and adapted for the oceanic case by Gaspar et al. (1990).
The river runoff is considered explicitly and modeled as a lateral open-boundary condition.As in Querin et al. (2006), the rivers are included by introducing small channels in the bathymetry that simulate the river bed close to the coast.Velocity is imposed at the upstream end of each channel, with the prescribed discharge rate being obtained by multiplying the velocity by the cross sectional area of the channel.
No flux conditions for either momentum or tracers and no slip conditions for momentum are imposed at the solid boundaries.Bottom drag is expressed as a quadratic function of the mean flow in the bottom layer: the (dimensionless) quadratic drag coefficient is set equal to 0.002.
The net transport through the southern open boundary is corrected during run-time at each time step to balance the effects of river discharge and of the evaporation minus precipitation budget on the surface level.This solution prevents any unrealistic drift in the sea surface elevation.Tides are imposed as a barotropic velocity at the southern boundary.
At the surface, the wind drag coefficient is computed following the default MITgcm formulation: where U 10 is the wind speed at 10 m.The treatment of surface heat forcing is done with the same bulk formula used in SHYFEM, except for the sensible and latent heat fluxes, where the formulation proposed in Large andPond (1981, 1982) is used.

Simulation setup
Two numerical experiments were carried out.In the second experiment the nonhydrostatic version of the MITgcm model is run, again over the same time period, to assess the importance of nonhydrostatic processes.
As the two models have different grids their resolutions are considerably different.In Fig. 1 we show maps of the difference in resolution of the two models, with red and blue indicating where MITgcm is more or less resolved than SHYFEM, respectively.As can be seen overall the MITgcm has higher resolution.Only in coastal regions do the models have comparable resolution, with the blue regions in the right panel indicating where SHYFEM is more resolved.

Observational data
In order to validate the model simulations, a number of observational data sets are used.Figure 2 shows their location.
i. CTD (conductivity, temperature, and depth) transects of temperature and salinity are provided from a cruise with the R/V DallaPorta, along the Senigallia transect (Fig. 2), where temperature and salinity profiles were acquired with a SeaBird Electronics SBE 911-plus CTD, on the 27 March 2012.These sets of data are part of a larger data set, collected in the bimonthly monitoring activity along that transect.
ii.Sea surface temperature (SST) from satellite data obtained using Moderate Resolution Imaging Spectroradiometer (MODIS).MODIS is a key instrument aboard the NASA Terra and Aqua satellites, which acquires measurements in 36 spectral bands.It can provide a wide range of atmospheric, land, and oceanic products: specifically for the ocean, MODIS SST is retrieved from radiometric measurements at 11 and 4 µm wavelengths with 1 km of spatial resolution.We selected the MODIS-Aqua SST for the Adriatic Sea, acquired during daytime on the 26 January, 5 and 16 February 2012, and available on the OceanColor web page of the Goddard Space Flight Center of NASA (http://oceancolor.gsfc.nasa.gov/).Since the SST products were highly affected by clouds, they have been adequately cloud masked, with the MODIS atmosphere products (http:// modis-atmos.gsfc.nasa.gov).To prevent the loss of river plume information in the area close to the coast, a quality mask was not applied to the SST, because the effects on the SST fields were negligible as tested by Reinart and Reinhold (2008).Finally, SST fields were remapped to geographical lat-long coordinates.
iii.In the Gulf of Trieste, time series of surface (2.5 m) temperature and salinity and bottom (22.5 m) temperature from the Vida buoy (location 45

Results
The backdrop for our simulations is the extreme DW outbreak that occurred during the winter of 2012.In the Adriatic Sea during the period January-February, there was an unprecedented generation of DW with record breaking density anomalies of above 5 kg m −3 relative to a value of 1025 kg m −3 (Mihanović et al., 2013).The event took place after a particularly warm and dry year, resulting in a reduction of coastal freshwater supply, in the backdrop of an already long-term trend in increasing salinity.The event was then triggered by an extended period of cold weather with strong Bora winds that lasted for about 3 weeks in the coastal eastern Adriatic region, between the 25 January and the 14 February 2012 (Mihanović et al., 2013).In what follows we will show how the two different models capture various aspects of this phenomenon, beginning first with a comparison with the measurements available for this period and then showing how the models reproduce a number of the specific processes affecting the coastal zone, namely the timing of the DW outbreak, its formation and evolution, coastal upwelling and riverine processes.We also look at a comparison of the hydrostatic and nonhydrostatic simulations.

Model validation
Here we provide an assessment of how well the two hydrostatic models do in reproducing the hydrodynamics of the coastal zone, as seen in the observations.In Fig. 3 the time series of surface and bottom temperature as well as surface salinity from the Vida buoy is shown.After the strong cold Bora wind event, the surface waters lose heat leading to DW sinking.The monotonic decrease of temperature continues until the beginning of February when the sudden loss of surface heat produces an abrupt drop in temperature from about 10 to 6 • C (2-3 February), before climbing up to 8 • C a couple of days later (Fig. 3a).On the bottom, temperatures reach even lower values (5 • C around the 15 February), suggesting an injection of cold water down from the surface from the areas in the vicinity during the Bora event.In the lead up to the event, up to the end of January, both models reproduce well the surface and bottom temperatures.Also they both capture the onset of the event, registering the starting moment of the cold water sinking at the beginning of February.However both models overestimate the minimum temperature values reached during the event.In the case of SHYFEM, the surface values are well represented; however, bottom temperatures are overestimated, with their values being close to those of the surface, indicating that the model mixes the water column too quickly and does not show any significant unstable stratification due to the DW sinking.MITgcm reproduces the surface temperature before the event, with just a small underestimation in the first simulated month.However, bottom temperature is closer to the observations during the event, and the greater difference in the temperature values between surface and bottom, indicates that MITgcm has a more unstable stratification and less mixing than SHYFEM during the DW event.
From the statistical analysis of the whole temperature time series, which is shown in Table 1, it is evident that the two models well reproduce measurements, with biases always lower than 0.2 • C. Correlation is higher than 0.96 for SHYFEM surface and bottom temperature data, while slightly lower for MITgcm.The two models show higher errors in reproducing the time series variability, as expressed by root mean square error (RMSE) values of around 1 • C.
The Vida buoy also shows a general increase of surface salinity during the first months of 2012 (Fig. 3b), probably connected with the low discharge of freshwater characterizing the whole north Adriatic Sea (NAS) in that period, in particular the closest river Isonzo (Mihanović et al., 2013).Both models overestimate the surface salinity during the whole period, as shown in Table 1 (bias around 0.3 psu for both models), with SHYFEM showing a lower variability compared to MITgcm results.Correlation for surface salinity is higher for  SHYFEM (0.84) than for MITgcm (0.73), but high enough for both models to state that the reproduction capability of the haline temporal evolution is matched (Table 1).
The AA Platform surface (−2 m) and bottom (−12 m) data for density anomaly, temperature, and salinity (Fig. 4) have a similar trend as that seen in the Vida buoy time series, with the density anomaly peak reached at the end of the first week of February 2012.The measurements reveal the stable stratification, just before the DW formation event (density anomaly difference, between surface and bottom, around 1 kg m −3 , even if the water column is thermally unstably stratified), the passing of the well-mixed DW (until the 22 March), and the subsequent re-stratification.SHYFEM does match the stable density stratification before the event, just for a few days, even if it has a slightly more homogeneous water column, compared with measurements (Fig. 4a).The DW signal is registered by SHYFEM perfectly matching the density anomaly values in the month of February.The major discrepancy is in the reproduction of the surface density anomaly after the event when a mass of lighter water is measured on the surface.MITgcm, as well, matches the general trend, with density anomalies closer to the measured ones, before the event, compared with SHYFEM, but overestimating the values during the event, in February.
Temperature trends are well matched by both models: Table 2 shows correlation values, for SHYFEM, of 0.88 and 0.97, for surface and bottom temperature, respectively.Also MITgcm shows high correlation values for temperature, even if slightly lower than SHYFEM (0.77 and 0.86 on surface  and bottom, respectively).SHYFEM better reproduces the temperature variability before the event and has a better match with observations in the post-event period, compared with MITgcm.Due to the lack of measurements of bottom temperature in the period just before the event, it is not possible to state whether the unstable thermal stratification, reproduced by SHYFEM or the well-mixed thermal structure, simulated by MITgcm at the AA platform, represents the real process.The salinity time series indicates that the density anomaly discrepancy in the models is due to the freshwa-ter dynamics, specifically the lack of direct measurements to impose as input for the studied period, as was the case at the Vida buoy.In fact higher salinity biases are registered by both models on the surface (1.18 and 1.07 psu for SHYFEM and MITgcm, respectively) while a better match is seen at the bottom (Table 2).Clearly river inputs that provide an incorrect amount of freshwater discharge, can directly affect the salinity variation (seen in Table 2, with very low correlation values for the two models).The salinity mismatch also affects the surface density anomaly (Table 2).
In Fig. 5 we show comparisons of model temperature and salinity with three CTD profiles from the Senigallia transect, moving away from the coast the profiles are indicated with dash-dot, dashed and solid lines respectively (location indicated in Fig. 2) for the 27 March 2012.The DW signal, produced at the beginning of February in the northern end of the basin, flowed along the Italian shelf and can be detected at the bottom around 20 km offshore, from CTD profiles (dashed profiles).SHYFEM reproduces the salinity profiles with a general underestimation of 0.5 psu.Generally MITgcm overestimates coastal salinity, while the more offshore profiles show similar differences from measurements for both models.Coastal haline stratification could be missed as a consequence of the Po River plume mismatch, which is discussed below.SHYFEM shows quite a clear underestimation of surface temperature, particularly offshore, with a bias of about −1 • C. However it matches better the data along the entire water column.MITgcm overestimates the surface temperature by 2 • C while it underestimates the bottom values by 1 • C.
In Fig. 6 we show maps of the SST from MODIS satellite observations and model minus satellite differences for both models, for three different times: before, during, and just after the DW period (26 January, 5 and 16 February).The comparison shows generally common behavior for both models, with differences in small-scale features.The 26 January satellite SST reveals a bulk of cold water, as is typical of the winter season, flowing out from the Po River that produces a clearly identifiable strip of coastal cold waters along the Italian littoral, just south of the river mouths.The cold discharge from the northern rivers is detected and the whole coast is characterized by a SST lower than 6 • C. SHYFEM and MITgcm overestimates the surface temperature of these waters, by around 2 • C. MITgcm tends to slightly overesti-  mate (between 0.5 and 1 • C) the whole area of NAS, except for a small area in the Gulf of Trieste, where also SHYFEM displays a slight underestimation (−0.5 • C).If results are considered in a basin-wide perspective, SHYFEM tends to underestimate SST in the deepest areas offshore, located in the center of the Adriatic Sea, while the biggest errors for MITgcm, are detected along the Italian littoral in the middle Adriatic Sea, with a strip of coastal waters underestimated by more than 2 • C (Fig. 6 top).Better performances of both models can be seen in the comparison with satellite images for the two dates during and just after the DW formation event.Still there is an overestimation of temperatures in the narrow strip in the proximity of rivers but the overall bias is reduced, in the range [1, +1] • C. A major discrepancy is seen in the reproduction of a cold structure just offshore of the Croatia littoral.Both models overestimate SST there and this suggests that a specific process is not reproduced that can be linked with atmospheric forcing as well as lateral freshwater sources.In fact it is possible that the amount of cold water injected into the system from the Po River, in the period preceding the DW event, and not reproduced by the models, enters into the general circulation of the basin and also affects the coastal area in front of Croatia (Fig. 6, center).The comparison with the satellite image from the 16 February shows the zone of highest bias (positive for SHYFEM and negative for MITgcm) just south of the Po River delta, crossing longitudinally the basin and in the transition zone between colder and warmer waters.It seems that in this frontal zone, delimiting the area with DW where vertical mixing would occur, the two models behave differently.Another major discrepancy between the two models can be seen in the narrow strip along the northernmost littoral, where SHYFEM overestimates SST by around 1 • C and MITgcm underestimates it by the same quantity.Generally SHYFEM has a bias in the range [−0.5, +0.5] • C in the offshore area of the NAS, while MITgcm tends to overestimates SST by 1 • C there (Fig. 6, bottom).To correctly interpret the outcomes from the model-satellite comparison, we should highlight that there are several factors, which might affect the performance of SST satellite-derived results.Satellite-derived SST is the skin layer temperature and it provides information on only a few microns of the sea surface.SST measured by buoys or derived by models are generally collected at depths from 0.5 to 5 m below the sea surface.These SSTs are called bulk SSTs.Therefore, the skin SST can be significantly different from the bulk SST.Referring to Donlon et al. (2002), surface thermal stratification can induce differences of some degrees between the skin and the bulk temperatures.In the western Adriatic Sea shelf, where the majority of river discharges occurs, the buoyancy flux due to river runoff at the sea surface causes a significant increase in the difference between the skin and the bulk temperatures.In addition, spatial variations in the near-coast surface winds might induce different levels of heating in different areas and generate spatial gradients in SST (Otero et al., 2009).It has to be stressed that also water turbidity due to river runoff can affect the SST: a modeling implementation in the Black Sea, done by Kara et al. (2004) and Kara et al. (2005) demonstrated that high turbidity affects the depth corresponding to solar radiation extinction and consequently the calculation of SST.Kara et al. (2005) demonstrated that using a clear-water constant attenuation depth assumption (as done also in the modeling work here proposed), as opposed to turbid water type values in the modeling implementation, produced monthly SST biases as large as 2 • C in the winter period in the Black Sea.Not being possible to apply different values of depth corresponding to solar radiation extinction, based on the presence of sediments (dynamics not simulated in the models), we had to take into account a possible bias in simulating SST close to river inputs of 2 • C.

Dense water formation and propagation
In Fig. 7 we show, for SHYFEM and MITgcm (both in the hydrostatic and nonhydrostatic implementation) time series of depth profiles of the average density anomaly, temperature, salinity, and root mean square (RMS) vorticity averaged over the NAS (in the shelf region above latitude 44 • N, for depths lower than 40 m -area shown in Fig. 2).The DW formation event is marked by a strong increase in the density anomaly at the beginning of February, with values reaching +5 kg m −3 for SHYFEM, and slightly less for MITgcm.The SHYFEM values are in agreement with those measured in Mihanović et al. (2013).SHYFEM describes the sudden formation of DW in the NAS and its sinking/mixing over the whole water column (Fig. 7a).Also the subsequent increase of the density anomaly at larger depths, during and after the event, is detected by SHYFEM.From Fig. 7, the moment of DW formation is clearly identified, for SHYFEM, after the first week of February, lasting for 1 week and then the progressive decrease of density anomaly marks its flow southward just out of the NAS.The temperature profile time series (Fig. 7b) for SHYFEM identifies the cold waters produced at the beginning of the event.Interestingly, the bulk of cold water changes its characteristics and temperature while sinking, even after the event.This stresses the fact that DW characteristics are evolving, being influenced by the mixing taking place with the surrounding warmer waters.SHYFEM simulates an increase in surface salinity during the DW event, suggesting evaporative processes due to the effect of the cold Bora wind.Therefore, during the DW formation, SHYFEM has a more homogeneous haline environment, highly thermally unstably stratified, that leads to DW sinking.MITgcm, like SHYFEM also registers the rapid increase in density anomaly at the beginning of February, even if the rate of increase and the peak reached by MITgcm is lower than SHYFEM.Similar temperature variations, with an unstable stratification characteristic of the DW formation, are seen by both models but it is less pronounced in SHYFEM.MITgcm has a lighter water environment at the beginning of the simulation, compared with SHYFEM, probably due to the presence of less saline waters on the surface (Fig. 7c).This can be responsible for the lower density anomaly simulated during the DW event by MITgcm.Another major difference between the models is in the evolution of the bottom salinity: MITgcm shows an increase just after the event, with an higher stable haline stratification.
From Fig. 7 we can also compare the hydrostatic and nonhydrostatic MITgcm simulations.For temperature, the nonhydrostatic run stratifies thermally, just after the beginning of the simulation, slightly more than the hydrostatic run.Also slightly colder waters are present near the bottom, during and after the DW event, in the nonhydrostatic MITgcm.This results in a relative increase of density anomaly close to the bottom in February for the nonhydrostatic run, though these small differences do not lead to any significant change in the vertical dynamics between the two runs (Fig. 7d).
The DW formation corresponds with a strong increase in vorticity detected by both models (Fig. 7d), though SHYFEM has a much stronger and more prolonged vorticity intensification relative to the MITgcm runs.

Circulation and vertical dynamics
In Figs. 8 and 9, we show maps of surface vorticity (surface currents overlaid) and net vertical velocity over the water column (wind vectors overlaid), respectively, from SHYFEM and MITgcm for the 26 January, 5 and 16 February, and 27 March 2012.Additionally, maps of the difference between MITgcm hydrostatic and nonhydrostatic runs are shown.The 26 January is characterized by spatially variable wind over the NAS.SHYFEM has a narrow band of positive vorticity along the Italian littoral.Just off the Po River delta, a band of negative vorticity is seen, probably due to the advection of freshwaters out of the main branch (Fig. 8).A cyclonic circulation in front of the Venice Lagoon is detected by SHYFEM, linked to upwelling there (Fig. 9), directly induced by the wind.MITgcm has the same patterns seen by SHYFEM but their values are about 1 order of magnitude lower, due to the much less energetic currents (0.1 m s −1 vs. 0.2 m s −1 in SHYFEM).The direct effect of wind forcing is seen in the surface vorticity map for the 5 February, during the DW formation event, both in SHYFEM and in MITgcm.Clearly there is a strong enhancement of coastal upwelling along the eastern coastline during the DW outbreak, as a result of the strong Bora winds driving Ekman suction.In this case the vertical dynamics in SHYFEM is due to the Ekman transport (coastal upwelling in the Gulf of Trieste coastal area) as well as to the surface cooling by the Bora wind.Strong negative vertical velocities indicating sinking are seen in the center of the Gulf of Trieste and in the whole NAS coastal zone.Interestingly, as stated by Mihanović et al. (2013), other sources of DW are seen along the coast of Croatia and in specific areas in the archipelago in front of it.MITgcm has a general cyclonic circulation in the NAS, bordered by littoral negative vorticity in the northern end of the basin.As with SHYFEM, the area of DW sinking is seen but with a lower magnitude of vertical velocity, even if a number of small-scale features are reproduced, showing higher horizontal variability of vertical processes.On one hand, the higher resolution of MITgcm over the NAS, allows for the reproduction of more smallscale vortical structures, identifying a wider spatial range of processes, compared with SHYFEM.On the other hand, the larger structures reproduced by SHYFEM (NAS gyre during DW event) seem more energetic, with lower dissipation along the vortices boundaries and lower large-to-small-scale energy turbulent cascade, increasing the net vertical transport.The stronger horizontal surface dynamics, registered in SHYFEM, can lead to energetic vortical structures that enhance the larger-scale vertical dynamics connected to them.
Just after the DW formation, on the 16 February, Bora wind starts to be weakened but is still present.SHYFEM shows a general surface circulation, in the whole NAS, in the direction east-west, directly following the wind curl.No specific downwelling is seen by SHYFEM, while the coastal area of the Gulf of Trieste and the Croatia littoral show upwelling, probably due to local effects of wind stress along the basin border.MITgcm still has negative velocities in the NAS and surface currents seem mainly directed along the north-south axis, with meandering and small-scale patterns.
The low wind on the 27 March produces surface currents in the NAS with different behaviors in the two models.SHYFEM shows weak but well-defined geostrophic circulation, going from east-west along the coastline in the northern end of the basin.Coastal vertical movements are not enhanced, except for a slightly positive vertical velocity in the offshore areas on the NAS.MITgcm shows a counter current, in the anticyclonic direction, with almost zero net vertical velocity.
Overall from Fig. 8, we see that the two models produce different surface current patterns.SHYFEM has more energetic coastal currents flowing southward, enhancing the freshwater transport out of NAS.MITgcm, throughout the DW event, has weaker surface dynamics, increasing the residence time of freshwaters in the NAS.
Figures 8 and 9 also provide insight into the role played by nonhydrostatic processes.From the difference plots it is clear that nonhydrostatic processes have no impact on the dynamics in the shallowest coastal area of the NAS.Only in the deeper basin further south do differences between the hydrostatic and nonhydrostatic simulations appear, in particular along the slopes of the sills of the south Adriatic.
To clarify how the surface forcings affect the two model simulations, Fig. 10 shows the total heat flux, in terms of gain/loss, for the four dates presented above.As a general picture, for the dates before and after the DW event, a lower gain of heat by MITgcm is seen, compared with SHYFEM, while during the 5 February the heat loss is more diffused for SHYFEM than for MITgcm.MITgcm has more local areas of heat loss, directly connected with the Bora wind jets, but also shows specific areas with small heat gain.In any case the values of heat loss seen by the two models correspond with the ones also simulated by Janeković et al. (2014).
In order to look closer at the coastal upwelling taking place, we examine the vertical velocity focusing on the coastal area of the Gulf of Trieste within 18 m depth (Fig. 11).SHYFEM and MITgcm time series of vertical velocity profiles, averaged over this sub area are shown (Fig. 11a).Here we present only the hydrostatic simulation because no significant differences are seen in the nonhydrostatic run for this area.A strong signal of positive (upward)  velocity is detected by both models during the DW event, though it is much stronger in SHYFEM than MITgcm.This upwelling is the result of the net Ekman suction induced by the Bora wind while there is a general DW sinking in the rest of the Gulf of Trieste.This is evident in Fig. 11b, where we show the comparison between the daily maximum Ekman induced vertical velocity (estimated from the two dif-ferent wind-stress formulations used in the two models) and the daily maximum net vertical velocity.Both models show an Ekman induced upwelling during the DW event, though the magnitude of the vertical velocity in SHYFEM is more comparable to the Ekman values, whereas MITgcm shows a lower net vertical velocity in the same period.The small differences in the Ekman velocities computed by the two mod- els are connected with the different formulation of wind drag coefficients, though overall they are very similar.

Riverine processes
Among the different coastal processes interacting during the DW event, the riverine inputs have an important role to play that must be taken into account.Different models can behave differently in reproducing the river plumes shape, in terms of both horizontal spreading and vertical mixing.Figure 12 infers both of these aspects, showing, for the four dates considered above, the depths where the highest haline gradient (freshwater above saline water) occurs, within the isoline of 37 psu that is chosen as a limit to border the ROFI environments.This choice mimics the approach proposed in Falcieri et al. (2014) that identifies the plume limit at 36 psu.We chose a slightly higher value, in order to include the bulk zone of the plume and the relative mixed area in its proximity.It has to be mentioned that, due to the low discharge characterizing the simulated period that enhances the DW formation, the ROFI is limited to a narrow coastal strip, except for the 27 March, when wind is weak and the discharges of rivers are relatively higher than in the preceding period.As for the previous images, the nonhydrostatic run is not shown due to the negligible differences in the NAS, compared with the hydrostatic run.
There is always a wider extension of surface freshwater for the MITgcm run than for SHYFEM.This is particularly evident along the northern littoral of the basin and can be seen throughout the whole period (Fig. 12).Focusing on the major river in the area, the Po River, it seems that the two models mixes differently along the water column, with a higher freshwater stratification for MITgcm than for SHYFEM, during periods of low wind and higher discharge (27 March).During the DW event, the effect of surface wind stress is high and leads to a more confined strip of freshwater in both models, though more so in SHYFEM, where the simulated stronger coastal current is enhancing the freshwater flow southward along the littoral (Fig. 12, 5 February).The major differences between the two models are seen in the 16 February, when there is a much larger surface spreading of freshwater in MITgcm.

Discussion
The reproduction of the majority of coastal processes, such as the DW formation and its spreading southward, coastal upwelling, and estuarine processes, requires taking into account a number of issues from a modeling point of view.SHYFEM and MITgcm are two very different models in terms of their numerical approach, parameterization, and their treatment of boundaries and forcings.
The two models demonstrated major differences in reproducing the correct amount of water with density anomaly higher than 5 kg m −3 (Fig. 7).The density anomaly produced by SHYFEM during the DW event is higher and closer to the measurements.The energy balances of the two models are different, as can be deduced by the total heat maps shown in Fig. 10.Even before considering how the dynamics acts on the water masses, it is important that the correct energy is injected into the system that will then be transferred into the vertical dynamics.The sinking processes would lead to higher vertical velocities and initiate stronger mixing with the surrounding waters due to the higher thermohaline gradient.
The validation section revealed the importance of having the correct setup in order to reproduce the predominant drivers of this phenomenon, i.e., the mechanical action of wind, acting on the sea surface, and the thermal flux due to the sudden cooling of the air-sea interface.Thus, the availability of correct and adequately resolved data set to force the models and a suitable treatment of surface boundary stress and heat-mass fluxes are required.The two models are forced with the same atmospheric data, but still have different features on the surface: Fig. 7 shows how SHYFEM is able to capture values of density anomaly comparable with the ones found by Mihanović et al. (2013).This implies that the surface forcings are realistic enough for the investigation of these phenomena and the differences between the models' results can be linked partially to the treatment of these forcings.
The two models apply different formulations in treating the wind stress, inducing slightly different Ekman transports.However, despite these differences, the Ekman velocities computed by the two models have the same timing, particularly during the DW event (Fig. 11b).On the local coastal scale, it seems that, even if the formulation is different, the dynamics connected with wind stress is similarly reproduced.Therefore, differences can be ascribed more to other factors, namely the bulk formulas used to compute the heat and mass surface fluxes.
In fact the choice of suitable bulk formulas that take into account the specific processes connected with the air-sea interaction and the heat and mass transfer through the interface strongly influence the capability to reproduce the DW formation.The comparison both with the Vida buoy and the Acqua Alta platform data shows different behaviors by the two models, concerning the SST.The different choice adopted by SHYFEM and MITgcm in dealing with the parameterization of the sensible and latent heat fluxes give rise to the different results of the two models, particularly after the DW event, approaching springtime when there is an increase in the heat gained by the sea at the surface.MITgcm has a warmer SST on the whole NAS area (Figs. 3, 4 and 6).The slightly lower heat flux through the surface computed by MITgcm (Fig. 10) can result in a lower injection of energy and the lower dynamics seen also in Figs. 8 and 9. Figure 7 also reveals that the surface salinity in MITgcm during the DW formation is lower compared with SHYFEM, which corresponds with the lower density anomaly registered by the model.Different parameterizations of latent heat adopted by the two models can play a role in this, in particular for the computation of evaporation.
Moreover, the differences seen in the salinity are affected by the lateral boundaries (i.e., river inputs), both in terms of forcing availability and boundary treatment.Due to the lack of measured discharge data for a number of rivers in the simulated area and the use of river water temperature from 2007 at the Tagliamento River, applied on all the lateral freshwater sources, the models show discrepancies in matching the surface salinity temporal variability (Figs. 3 and 4) and the spatial surface thermal pattern close to river mouths (Fig. 6).What is known is that winter 2012, unlike the climatology, is colder than the average.Therefore, the limit of both models, in providing the real temperature of freshwater discharges in the NAS, and the use of climatological discharge time se-ries is responsible for the mismatch found.Even if the use of the climatology is reasonable, at least for what concerns the salinity, for a general characterization of the dynamics as the model-measure biases are around 1 psu, this choice can deeply affect the reproduction of suddenly varying or extreme events like the DW formation in winter 2012.
Moreover, even though the two models are forced with the same data sets at the lateral boundaries, their surface biases have different signs, suggesting that there is also a substantial difference in the momentum laterally injected into the system and in the vertical mixing simulated (Figs. 6 and 12).The two models apply different approaches in dealing with lateral boundaries, with differences in the reproduction of river mouth morphology, and in the momentum applied for the freshwater discharge.The difference in resolution, just along the coastline, that shows a higher-resolved river channel shaping in SHYFEM, can affect the river discharge inflow as a consequence of the geometry of the input points.The advection induced in the transition zone between the narrow, well-defined river channel and the open sea, as reproduced by SHYFEM, could be more pronounced and could lead to different plume shaping.MITgcm includes, as well, river channels but the lower resolution, just in the proximity of estuaries and deltas, can affect the dynamics.Moreover, MITgcm imposes velocity values at the upstream end of each channel, and discharges are computed multiplying the values by the cross sectional area.SHYFEM directly imposes the measured discharges and momentum is injected into the system as a consequence of the water level gradients between the boundary nodes and the surrounding nodes.It is not possible to state if one of the two approaches is more suitable, but as it has a possible effect, even if minor compared with the geometric effect due to the different resolutions, it is worth mentioning.
The different temperature and salinity fields simulated by the two models close to the river mouths, mainly due to the freshwater sources, provide different baroclinic gradients, affecting the coastal thermohaline circulation.Additionally, the horizontal schemes used in the models, for horizontal advection and diffusion of scalars (i.e., temperature and salinity), would then lead to different baroclinic currents and to the higher amount of freshwater in the sub area of NAS for MITgcm, compared with SHYFEM (Fig. 7).The model biases computed at the Senigallia CTD transect (Fig. 5) are another example of this situation: coastal haline stratification could be missed as a consequence of the Po River plume mismatch, which is evident in the thermal bias shown in Fig. 6 that lead to the different riverine water spreading.
The model differences in the more offshore area of the NAS can also be influenced by the different resolution of the two models: on one hand it seems that higher resolution is needed along the coast, to reproduce the complexity of the coastal morphologies; on the other hand, a higher resolution offshore, as is the case for the MITgcm grid, leads to a number of small-scale vortical structures, generally missed by SHYFEM (Fig. 8).It seems that SHYFEM is horizontally more diffusive than MITgcm, but less dissipative (from an energetic point of view).Therefore, less horizontal fronts are seen by SHYFEM but they are steeper.An open point, that cannot be easily discriminated from the obtained results, is the relative importance of horizontal advection and mixing parameterization in the reproduction of these processes.In order to add some information on the relative effect of horizontal viscosity parameterization, a preliminary test simulation with MITgcm was carried out, though not shown here.The variable horizontal viscosity parameterization, proposed by Leith (1968) and already implemented in Sannino et al. (2014) for the Gibraltar Strait, was checked.Enhanced vorticity patterns, compared with the hydrostatic MITgcm run discussed in the paper were seen, suggesting the horizontal viscosity plays a role in correctly reproducing horizontal features.Consequently, effects on the vertical transports were detected.Further investigation of this topic is left for future work.
Analyzing the models' capability to reproduce the vertical hydrodynamic field structure, as a general comment, it seems that SHYFEM outputs are characterized by higher vertical mixing compared with MITgcm, that can be ascribed to both the higher thermal energy gradient between water masses simulated in the process and the vertical mixing connected with the turbulence closure schemes used.
Finally, Figs. 8 and 9 suggest the neglectable importance of nonhydrostatic processes in producing the dense water formation and coastal upwelling in the NAS.These processes are governed by the mass balance, in terms of horizontal transfer of water due to wind vs. vertical suction of water to compensate, which are already reproduced with the hydrostatic approximation.

Conclusions
The coastal zone of the NAS is characterized by a number of hydrodynamic processes that interact and evolve on different spatial and temporal scales.The present work demonstrates the complexity of modeling these specific processes and identifies a number of issues needed for choosing the most suitable modeling strategy for this typology of study.The main findings are listed below.
-The two models use different bulk formulas for surface latent and sensible heat, accordingly to their state-ofthe-art setups already tested in the area.This leads to different heat transfer at the surface, giving rise to an overall different energy balance.Lower convective dynamics over the water column is reproduced in the case of MITgcm relative to SHYFEM.Therefore, the choice of suitable bulk formulas, specifically in the coastal zone, is a central point for modeling implementations.
-There are differences in the small-scale hydrodynamic structure in the offshore area of NAS that are connected with higher resolution over the whole domain in MITgcm.However, these fine-scale features in MITgcm have little impact on the overall reproduction of the dense water formation; therefore, the presented implementation identifies the spatially variable minimum resolution adequate to reproduce the investigated processes, that span from less than 500 m in the nearshore area up to 1-2 km offshore.Higher resolutions do not add information on the main investigated dynamics.
-A highly resolved coastal zone, with the possibility to reproduce the complex morphology connected with lateral freshwater inputs, can provide the correct momentum injection into the system and affects the capability to reproduce buoyant processes in the coastal area.
-Nonhydrostatic processes have little impact on the coastal features seen on the shelf of the NAS, suggesting that the hydrostatic models are adequate for simulating DW formation in the shallow areas of the basin.
There are a number of outstanding issues that are not tackled in this present work.It is important to point out that in this paper we used different wind stress formulation and bulk formulas for the two models, which result in uncertainties in the model comparison.Further studies, perhaps using direct forcing with heat and mass fluxes provided by meteorological models, would be helpful in reducing these uncertainties.Other open questions not considered in this work are the effects that different horizontal advection, mixing, and turbulence closure schemes have on coastal hydrodynamic processes, such as the dense water event considered here.Such a study would require using a single model with different implementations of these schemes to precisely characterize and attribute their impact on the coastal dynamics.Also, here we found that nonhydrostatic processes have little impact in the shallow coastal shelf of the NAS, though there were differences from the hydrostatic case seen in the deeper part of the basin.Exploring the entire Adriatic basin may reveal if the nonhydrostatic dynamics plays any part in the wider propagation of dense water through the basin, particularly in the southern Adriatic pit.Despite these outstanding questions, this work provides some clarity on the chosen setups that were already used for implementations in the Adriatic Sea, giving suitable suggestions for improvements.These modeling implementations, mainly devoted to process investigations, can be used to guide choices made for possible future operational products.

Figure 1 .
Figure 1.Maps of the difference in grid resolution (in km) between MITgcm and SHYFEM (SHY minus MIT) in the entire Adriatic (top left), and a close up of the north-eastern coastal area, the Gulf of Trieste (bottom left).Red indicates where MITgcm is higher resolved than SHYFEM while blue indicates the reverse.The two grids, SHYFEM (top right) and MITgcm (bottom right) are shown for the area of the Gulf of Trieste.

Figure 2 .
Figure 2. Map showing location of data sources used for validating the model: CTD data (red dots), Acqua Alta CNR Platform (yellow dot), and Vida buoy (green dot).Purple dots show the location of river inputs.

Figure 3 .
Figure 3.Time series of (a) surface (left panel) and bottom (right panel) temperature and (b) surface salinity for SHYFEM (blue), MITgcm (red) and the Vida buoy observations (black).

Figure 4 .
Figure 4. Timeseries of surface (solid) and bottom (dashed) (a) density anomaly, (b) temperature, and (c) salinity for SHYFEM (blue), MITgcm (red), and the CNR Platform observations (black).Note there is a gap in the data for the bottom observations between the 8 and 25 January.

Figure 5 .
Figure 5. Senigallia transect: profiles of salinity (left) and temperature (right) for three CTD profiles along the transect for the 27 March 2012 (big red dots shown in Fig.2).The inset shows comparison between observations (black), SHYFEM (blue) and MITgcm (red) for the innermost, shallower profile (dash-dot), one in the center of the transect (dashed), and the outermost, deeper profile (solid).

Figure 6 .
Figure 6.MODIS SST images (left column) and bias maps showing the difference between model SST from MODIS satellite observations for SHYFEM (center column) and MITgcm (right column), for times before, during, and just after the dense water event, namely 26 January (top row), 5 February (center row), and 16 February 2012 (bottom row).Units are [ • C].

Figure 7 .
Figure 7. Time series of depth profiles of the (a) density anomaly, (b) temperature, (c) salinity, and (d) RMS of vorticity averaged over the north Adriatic area, for SHYFEM and MITgcm, both in the hydrostatic and nonhydrostatic implementation.

Figure 8 .Figure 9 .
Figure 8. Maps of surface vorticity with surface current overlaid for SHYFEM (left panel) and MITgcm (central panel), in the hydrostatic implementation, and differences in vorticity between MITgcm hydrostatic and nonhydrostatic implementations (HY-NH, right panel), for the dates indicated.

Figure 10 .
Figure 10.Maps of surface net heat flux for SHYFEM (top) and MITgcm (bottom) for the dates indicated.Units are [W m −2 ].
Figure 11.(a) Time series of depth profiles of vertical velocity averaged over the coastal area of the Gulf of Trieste, with depth lower than 18 m, and (b) time series of maximum daily Ekman wind curl and daily maximum of net vertical velocity for SHYFEM and MITgcm.

Figure 12 .
Figure12.Maps of the depth at which the highest salinity gradient (freshwater above, saline waters below) occurs, within the 37 psu isoline that identifies ROFI in the NAS, for the dates indicated.
Institute of Marine Science -National Research Council) within the last 15 years in the area in front of the Venice Lagoon.Water levels are set at the mean sea level as an initial condition and are then adjusted to the computed values.3-D velocity values are initially set to zero.The main open boundary is located at the Otranto Strait.3-D temperature and salinity and tidal water level time series force the open-boundary section.At the lateral open boundaries, corresponding to river inflows, discharge time series are imposed.

Table 1 .
Statistical analysis of simulated water temperature and salinity time series computed at the Vida buoy.Analyses provided are the difference between mean of observations and simulations (Bias), the root mean square error (RMSE), and the correlation.

Table 2 .
Statistical analysis of simulated water temperature and salinity time series computed at the Acqua Alta CNR Platform.Analyses provided are the difference between mean of observations and simulations (Bias), the root mean square error (RMSE), and the correlation.