Impact of currents on surface flux computations and their feedback on dynamics at regional scales

A twin numerical experiment was conducted in the seas around the island of Sardinia (Western Mediterranean) to assess the impact, at regional and coastal scales, of the use of relative winds (i.e., taking into account ocean surface currents) in the computation of heat and momentum fluxes through standard (Fairall et al., 2003) bulk formulas. The Regional Ocean Modelling System (ROMS) was implemented at 3 km resolution in order to well resolve mesoscale processes, which are known to have a large influence in the dynamics of the area. Small changes (few percent points) in terms of spatially averaged fluxes correspond to quite large differences of such quantities (about 15 %) in spatial terms and in terms of kinetics (more than 20 %). As a consequence, wind power input P is also reduced by ∼ 14 % on average. Quantitative validation with satellite SST suggests that such a modification of the fluxes improves the model solution especially in the western side of the domain, where mesoscale activity (as suggested by eddy kinetic energy) is stronger. Surface currents change both in their stable and fluctuating part. In particular, the path and intensity of the Algerian Current and of the Western Sardinia Current (WSC) are impacted by the modification in fluxes. Both total and eddy kinetic energies of the surface current field are reduced in the experiment where fluxes took into account the surface currents. The main dynamical correction is observed in the SW area, where the different location and strength of the eddies influence the path and intensity of the WSC. Our results suggest that, even at local scales and in temperate regions, it would be preferable to take into account such a contribution in flux computations. The modification of the original code, substantially cost-less in terms of numerical computation, improves the model response in terms of surface fluxes (SST validated) and it also likely improves the dynamics as suggested by qualitative comparison with satellite data.


Introduction
The assessment of the fluxes at the air/sea interface is an issue of crucial relevance for many topics in geophysics.A correct parametrization of such exchanges is relevant for climatic studies, climate change, weather and ocean forecasting and their applications in marine and maritime sciences.Wind stress, which is the medium of the momentum flux between atmosphere and ocean, is one of the main drivers of the ocean circulation for a wide range of spatial and temporal scales.The wind stress (τ ) in ocean models, when not directly provided by atmospheric forecasts, is usually computed through the so-called bulk formulas as described by Fairall et al. (1996), where τ is equal to the square of the wind speed at 10 m times the air density multiplied by a dimensionless drag coefficient (usually function of wind speed).Fairall et al. (2003), updating their previous work, suggest the use of relative wind vectors to compute the wind stress, i.e., to take into account ocean currents subtracting them from the absolute wind vectors.The contribution of the ocean currents in the computation of the wind stress has been neglected in ocean modelling for many years.The fastest ocean current is 1-2 order of magnitude smaller than the stronger wind: for this reason the surface current contribution was often neglected in applying bulk formulas, even if an estimation of surface currents is often easily available from the ocean model itself.Considering that the computation of the wind stress contains a squared velocity term, it can be easily un-derstood that the relative contribute of ocean currents is also squared, becoming really relevant for low-wind conditions.Further, as the drag coefficient is also a function of the wind speed, the inclusion of the surface currents also affects the drag term, further increasing the impact of such a component.Heat fluxes may also be impacted by including the surface currents, in this case with a linear relation with wind-current velocities.By taking into account the surface current component, bulk formulas for momentum, sensible and latent heat fluxes can be written as: (1) where ρ a is the air density, u a and u s are the vector velocities for air and sea surface respectively, t a − t s is the difference in temperature between air (at 10 m) and sea surface, q a − q s is the difference in humidity, C pa and L e are the specific heat of air and the latent heat of water evaporation respectively, while C d , C s and C l are the coefficient for momentum, sensible heat and latent heat transfer respectively.Some recent papers provided evidence of a moderate but actual impact of such a modification on fluxes at global/oceanic scales.Kara et al. (2007) showed that the inclusion of ocean currents and dominant waves into the drag computation leads to a daily reduction of the drag of about 10 % at daily scale and for the entire globe, with large variability between mid-latitude (smaller impact) and tropics.Another model study (Dawe and Thompson, 2006) found that, for the North Pacific, heat fluxes and wind stress changed about 1-2 % at basin average, while localized changes (in the tropics) reached up to a 10 % reduction of both momentum flux and surface currents.In that study the wind power input to ocean surface is reduced by 27 %, quite in good accordance with previous findings of Duhaut and Straub (2006).In the Gulf Stream region, this reduction of the wind work was estimated to be around 17 % (Zhai and Greatbatch, 2007).Deng et al. (2009) also assessed the effect of coupling currents with winds.They found a 10 % change in surface currents when considering surface currents velocities in the bulk formulas, quite in agreement with other authors.All authors found that in the tropics such changes are more relevant than for mid-latitudes.This is a valid generalization for large scales, while an insight of what happens at midlatitude and at regional and coastal scales was not provided yet.To address this issue we focused our attention in the seas around Sardinia (Western Mediterranean sea), which is a highly variable and dynamic area interested by several (sub-)mesoscale structures of different origin (Fuda et al., 2000;Puillat et al., 2002;Ribotti et al., 2004;Testor et al., 2005;Olita et al., 2013) and also affected by strong wind events characterized both by seasonality and high-frequency peaks.
The area is quite heterogeneous in terms of circulation and dynamical characteristics.The Sardinian Sea, i.e., the continental shelf and slope area west of Sardinia, is part of the Algero-Provençal Basin.From the basin scale circulation perspective, the Sardinian sea is located in between the Algerian Basin to the south, dominated by the inflow of Atlantic water from Gibraltar advected eastward by the Algerian Current, and the Provençal Basin to the north characterized by the path of the Northern Current (Millot et al., 1999) moving south-westward along the Italian and French continental shelf and where a surface cyclonic gyre drives the northern sub-basin circulation Lévy et al. (1998).The southern branch of this cyclonic gyre contributes to the formation of the North Balearic front (Fuda et al., 2000;Testor and Gascard, 2003;Olita et al., 2014) which represents the separation between the Atlantic water reservoir of the Algerian Basin and the saltier and denser waters of the Provençal basin (e.g.Olita et al., 2014).In a recent paper (Olita et al., 2013) we suggested, through the analysis of the outputs of a 3-D assimilative model, that the upwelling occurring along the SW Sardinian coast was pre-conditioned by the presence of a quasipermanent southward current (Western Sardinian Current -WSC) whose origin was in part due to the approaching of anticyclonic eddies to the western Sardinia shelf.This was also supported by the findings of Pinardi et al. (2013) where the same current (they called Southerly Sardinia Current -SSC) is described as permanent at low-frequency scales (decadal) bordering a northern branch of the Atlantic water flow in the Western Mediterranean.In the southern part of the model domain the Sardinian Channel connects Tyrrhenian and Algerian sub-basins.Here the Algerian Current (e.g.Millot et al., 1999) transports Atlantic water towards and across the Sicily Channel.North of Sardinia the Bonifacio Strait (∼ 15 km wide) separates Sardinia and Corsica and connects, with its narrow passage, Algero-Provençal and Tyrrhenian basins.Winds crossing the strait contribute to the generation of a wind-driven quasi-stable cyclonic gyre (Perilli et al., 1995;Millot et al., 1999) in the northern Tyrrhenian sea, east of Sardinia, that represents the most energetic mesoscale structure of the northern Tyrrhenian sea (Iacono et al., 2013).All these different characteristics make this domain a good test case to study the impact of the inclusion of surface currents on the surface fluxes and their feedback on circulation at regional and coastal scales.
The aim of the present work is to study the impact of the surface currents in the computation of the surface momentum and heat fluxes, through the bulk formulas (Fairall et al., 2003), in turn driving the surface dynamics and surface temperature.The latter can be modified both through changes in surface heat fluxes and as a consequence of changes in vertical and horizontal motions.To evaluate such an impact, we performed a twin experiment with the Regional Ocean Modelling System (ROMS).ROMS was implemented in the central area of the Western Mediterranean sea, around Sardinia Island, at 3 km resolution.The chosen resolution al- lowed to respect the suggested 1 : 3 ratio (e.g.Debreu and Blayo, 2008) between child and parent grid resolution, as well as to well resolve mesoscale processes (considering that the smallest Rossby radius of deformation for this area is of the order of about 10 km).Details on the model implementation and experimental setup are provided in Sect.2, together with information on the data and analyses performed.Two experiments were conducted, reproducing the circulation of the year 2012, with and without the contribution of surface currents in the computation of the momentum and heat fluxes.In Sect. 3 we validate the model vs. satellite SST and compare the outcomes of the two setups under different points of view.Finally, concluding remarks are drawn in Sect. 4.

Numerical model and experiments
The numerical model is an implementation of the Regional Ocean Modelling System (ROMS Shchepetkin andMcWilliams, 2003, 2005) in its official release from Rutgers (svn revision 705).Such a release of the code does not include an option to switch on/off the surface currents in the surface flux computations, so the original model code was modified in this sense.ROMS is a free surface, hydrostatic, primitive equation, finite difference model widely used by the scientific community for many kinds of applications: large-scale circulation studies (e.g.Haidvogel et al., 2000), ecological modelling (Dinniman et al., 2003), coastal studies (e.g.Wilkin et al., 2005;Iermano et al., 2012), sea-ice modelling and others.The model was implemented in the seas around Sardinia (Fig. 1) in a rectangular grid of 3 km resolution on the horizontal plane and 30 s terrain following levels.The equation distributing vertical levels allows a robust description of surface and subsurface layers where most of the dynamical processes occur.Intermediate and deep layers are discretized with larger (on the vertical) meshes.
Bathymetry was derived from the General Bathymetric Chart of the Oceans (GEBCO), a global 30-arc-second database, and smoothed with a Shapiro filter to remove wavelengths of the order of the grid scale.This is in order to minimize the pressure gradient force error (PGFE) often caused by too steep bathymetric gradients.Stiffness parameters (rx0 and rx1, respectively 0.27 and 5.97) are well within the thresholds suggested by developers (ROMS user forum at https://www.myroms.org/forum/).Initial and boundary conditions were provided by the 1/16 • model of the Mediterranean Sea MFS-1671 (Tonani et al., 2009) retrieved through the My-Ocean (www.myocean.eu)data portal.Daily analyses of 3-D fields of velocities, temperature, salinity and elevation (2-D) have been used for model nesting and initialization.At the boundary the model uses Flather conditions (Flather, 1976) for the barotropic velocities while baroclinic velocities and 3-D tracers (T and S) are clamped to the values prescribed by the outer model.Even though ROMS allows the choice of more advanced BC than a simple clamped solution (e.g.Penven et al., 2006;Mason et al., 2010;Marchesiello et al., 2001), we opted for this solution for the present application as radiative conditions are not recommended by developers when using bulk formulas for flux computations.Short sensitivity studies to different BC confirmed that some issue exists when we use radiative conditions, generating spurious upwelling/downwelling features at boundaries.This suggested to opt for a simple but robust BC as the clamped one.At the free surface the Chapman (Chapman, 1985) boundary condition was imposed.A third-order upstream horizontal advection of 3-D momentum (Shchepetkin and McWilliams, 1998) and the k − turbulence closure scheme (Warner et al., 2005) were used in the present implementation.At surface, which is the focus of the present work, we used the 1/8 • 6-hourly ECMWF ERAinterim analyses fields.10 m air temperature, U and V wind momentum components, air pressure, solar short-wave radiation, air humidity and precipitation were used to compute freshwater, momentum and heat fluxes by using the abovecited bulk formulas.

Experiments
Two experiments were performed: the Bulk Fluxes (BF) experiment did not include surface currents (so in Eqs.(1-3) the u s term was neglected), while in Bulk Fluxes with Currents (BFC) the flux computations interactively took into account the surface currents reproduced by the model.This modification is quite straightforward when done in ROMS source code, just by taking care of the fact that currents and winds run on staggered grids, and therefore they have to be "interpolated" before performing subtraction.As a simple proxy for this, we averaged the i and i + 1 velocity points for each wind i point and then we subtracted such quantities.The simulations were integrated for 1 year, to simulate the 2012 year with boundaries and surface forced by the above described analyses fields.Daily averaged fields are then saved in the output files.The issue related to the data assimilation deserves a little discussion.Considering that the model boundaries are provided by an assimilative model we think that for such a domain the information contained in the boundaries would propagate to the nested model without a substantial loss of information.On the contrary for larger off-line nested domains it was shown in literature (Vandenbulcke et al., 2006;Olita et al., 2012) that assimilation would be needed as information coming from boundaries quickly dissipates.Further, and maybe more important, in the present work we aim to observe changes generated by different parametrizations of surface physics: for this reason we should not hide this signal with any statistical correction.

Validation data and metrics
Model performances at surface were evaluated by using satellite SST fields.Satellite SST used is the MyOcean sea surface temperature operational product for the Mediterranean Sea which is available at http://www.myocean.eu.The daily gap-free maps (Level 4, Optimally Interpolated) at 1/16 • of resolution have been used.Three basic metrics, namely Bias, Root Mean Square Error (RMSE) and Anomaly Correlation Coefficient (ACC) provide a good overview of the quality of the model in reproducing the observed SST.While RMSE and BIAS describe the model error in reproducing observed value of the considered variable, ACC measures the ability of the model in reproducing anomalies of the SST signature detected by the satellite, partly overlooking their absolute value.
The three metrics are formulated as follows: where mod and obs are respectively modelled and observed values of the variable and the overbar stands for a long-term temporal average.In the present paper this long temporal average is the AVHRR monthly climatology .The removal of the climatology allows to filter off the seasonal signal that otherwise would hide the response of this metric to the synoptic features.ACC is a dimensionless number ranging from −1 (worst) to +1 (best).
To further evaluate the quality of the two simulations, with particular reference to the dynamical features produced by the two experiments, we compared model outputs with synoptic observations of the sea surface observed in single swaths (Level 2) satellite images.Ocean colour and SST collected by MODIS sensors (on board of TERRA and AQUA satellites) were used for this purpose.Both typology of products (optical and infrared derived respectively) can provide useful information on surface and subsurface structures.An example of such data comparison is presented in Sect. 3 which tries to emphasize differences between the two model setups and similarities with observed features.

Flow decomposition, kinetics and wind work
In order to investigate the impact of the fluxes on the simulated dynamics, we separated the stable and the fluctuating   part of the velocity field as already described, for example, in Olita et al. (2013).The time-averaged term u = u + u represents the stable part of the flow, while u is its fluctuating part.The fluctuating components can be used to describe both eddy kinetic energy (EKE = 1/2(u 2 + v 2 )) and the Reynolds stress covariance term (RS = u v ), i.e., the eddy momentum flux.Reynolds stress covariance shows where the turbulent part of the flow interacts with the mean flow, accelerating or deflecting it from its mean direction (Greatbatch et al., 2010).It is likely that changes in surface parametrization of surface fluxes would influence both the stable and the fluctuating part of the flow, but in different measure.This suggested us that the two should be investigated separately.Wind stress work P , which is defined as the product of wind stress τ by the surface ocean currents u s was computed in order to assess the differences in terms of wind power input to the ocean between the two model experiments.At oceanic scales, including tropical areas, a reduction of about 20-30 % was recorded when the contribution of surface currents is considered (Duhaut and Straub, 2006;Hughes and Wilson, 2008).

SST validation and intercomparison
Figure 2 shows the three time series for SST RMSE, BIAS and ACC metrics computed vs. the SST data.
Such spatially averaged metrics do not show dramatic differences between the two setups, drawing an almost identical trend.However, BFC (i.e., with currents) setup shows, as the simulation progresses, slightly better performances than BF in terms of BIAS and RMSE metrics (see Fig. 2), while almost no differences are recorded for ACC.RMSE improvements seem to be mainly due to the reduction of the BIAS, which is kept closer to zero by the inclusion of currents term on flux computations.Of course, considering the complexity and the diversity in dynamics of the study area, some spatial variability of RMSEs values would be expected.Actually the comparison of RMSE maps (Fig. 3) for the two experiments shows quite a large variability over space.Although the general distribution of the errors is quite similar, there are many spots where differences are of great magnitude.In general terms, the largest improvement for BFC solution is shown on the western side of the domain (the windy side), with an evident reduction of RMSE for the Algero-Provencal basin.Locally, major improvements for BFC setup are located along the western Sardinia coast where the offshore extent of the patch with large RMSE is smaller than in BF (check figure at about 40 • N, 8 • E).Central and southern Tyrrhenian sea also show noticeable improvements, while in northern Thyrrenian the cyclonic gyre east of Bonifacio strait seems to be better described by the BC setup.

Impact on surface fluxes
All the spatially integrated surface fluxes (momentum, sensible, latent and net heat, see Fig. 4) show an impact of the order of few percentage points (∼ 2 %) by averaging time series values, but with a distinct high-frequency behaviour.
The small differences in terms of time series underneath quite large differences in space because of the very nature of the fluxes and the way they are computed (i.e., interactively during the model integration and with a feedback with ocean currents for BFC experiment).In this regard significant information is provided by the time integrated wind stress difference map shown in Fig. 5.Such spatial differences, for wind stress, reach a low of −8 × 10 −3 N m −2 in the proximity of the southern boundary of the domain where the highly unstable Algerian Current flows.Another low is at the turning point of the Western Sardinia Current in the SW corner of Sardinia (−6 × 10 −3 N m −2 ).Negative patches are quite dominant, as expected.Positive patches are less present, reaching a maximum of ∼ 2×10 −3 N m −2 and almost entirely located along the eastern Tyrrhenian coast.In percentage terms these spatial differences range between −15 % and +20 % on an annual basis, while they are obviously larger considering a daily basis.The values of net heat flux difference (right panel of Fig. 5) are highly patchy and correlated with areas of improved model performances in terms of SST (as shown in Fig. 3).Near the western Sardinia coast (40 • N, 8 • E) the model shows the largest correction in terms of heat fluxes (negative blue patch).

Impact on the mean and turbulent surface circulation
As expected, wind stress changes generated significant modifications in circulation and kinetics.Figure 6 shows time series of total kinetic and eddy kinetic energy for the two experiments.It is evident that the introduction of the currents on stress computation (BFC) led to a large (spatially averaged) reduction of the kinetic energies at surface.Such a reduction, about 23 %, shows a long period maximum between days 150 and 250, i.e., during summer.This is probably because winds reach their minimum intensity, then increase the impact of the correction done by surface currents on momentum flux.The largest part of such a difference between total kinetic energies at surface (about 65 %) is actually due to the turbulent part of the flow.Eddy kinetic energy shows (right panel of Fig. 6) the same trend as the  total one, with a maximum difference between the two experiments also during summer.Time-averaged maps of the above quantities provide an insight of the distribution of such differences.
The mean flow (Fig. 7) reveals some relevant change in terms of averaged path of the Western Sardinia Current which shows, surprisingly, a stronger signature in the BFC configuration.Some important change is evident in terms of mesoscale circulation: stable eddies footprints appearing in the SW side of the domain (with probable influence of boundaries) and west of Bonifacio strait, disappear in the averaged circulation field in favour of streams or a meandering feature.On the other side, EKE maps (Fig. 8) reveal that a large part of dynamical differences can be ascribed to the fluctuating part of the circulation as already argued from the time series.Eke values between the two model solutions shows an averaged reduction of about −23 % for BFC.The largest differences between the two EKE estimates are ascribed in the area of strongest mesoscale activity (Algerian Eddies area).Here a qualitative comparison of modelled fields vs. Level 2 single swath MODIS SST for 29 June 2012 (as an example) elucidates the differences between the two solutions in terms of mesoscale dynamics (Fig. 10).
In this figure the BFC solution better matches, in size and location, the large anticyclonic Eddy centred at about 38 • N, 8 • E, whose signature is also evident in the satellite image.BF solution on this case draws a less clear signature of the eddy and of the front on the east side of the eddy that, in the satellite image, seems also responsible for structuring the path of the WSC south of Sardinia.This is likely a recurrent correction in the new model solution, as this area seems highly impacted by the flux correction both in terms of EKE and heat fluxes.Reynolds stress covariance maps (Fig. 9) also provide useful information on the impacts of flux modifications: largest differences are identified in the area of highest mesoscale activity, i.e., in the Algerian Anticyclonic eddies area.Alternating negative and positive patches appear, in the BFC solution, close to the western Sardinia coast: this feature could help to explain the intensification of the WSC that appears for BFC.In this area BFC provide a solution, in terms of Reynolds stress covariance, close to the one we previously   found (Olita et al., 2013) through an interannual experiment performed with a numerical assimilative model.

Conclusions
In the present work the impact of the surface currents in surface flux calculations at regional/coastal scales was assessed.To do this we performed 1-year long simulation with a new implementation of ROMS in the seas around Sardinia (Western Mediterranean Sea) by using two different setups, with and without the contribution of the currents in the computation of surface fluxes through bulk formulas.
We found, according to bibliography that was mainly related to oceanic and basin scales, that domain-averaged momentum and net heat flux change by some few percentage points while more consistent differences are found for surface kinetic energies (BFC records a 21 % reduction on total surface kinetic energy in respect to BF).Differences can be observed both in the mean and fluctuating part of the flow.The latter showed major changes in the SW area of the domain, where mesoscale eddies are dominant.Inclusion of surface currents determines relevant changes not only in dynamics but also in the prognosed surface temperature by means of the surface heat fluxes.Validation with satellite SST reveals that the solution is generally improved, even if only slightly in spatially averaged terms.Shelf-slope area west of central Sardinia largely benefits by the correction, while some areas shows questionable results, as for example the cyclonic area east of Bonifacio strait.
Central and southern Tyrrhenian also show improvements in the BFC solution.
While quantitative metrics for SST reveal that net heat fluxes and resulting SST are improved, it is purely speculative to ascertain (i.e., not quantitatively validated) if, at these scales, the use of relative winds brings quality to the simulated dynamics or not.Comparison of synoptic satellite infrared and optic observations with modelled results did not solve the issue even if it provided some interesting hint in favour of the BFC solution.
Wind stress work, the product of wind stress and ocean surface currents, provides an insight of the wind power input to the ocean.Such an input is reduced for about 14 % as basin average.A difference map between the two estimates is shown in Fig. 11.The map shows larger differences on the left side, which is more windy and dynamic than the Tyrrhenian sea.The most interesting feature is the localized increase of P in coincidence with the WSC (slightly shifted westward in BFC in respect to BF) justifying the increased signature of the WSC we detected in the averaged flow (cfr.Fig. 7).
The present study provides evidence that the contribution of surface currents should not be neglected in the computa-tion of fluxes even at regional/coastal scales and in temperate regions.This is especially true and important for areas highly populated by (sub-)mesoscale features, which, in turn, are responsible for the modulation of relevant physical-biological processes at sea as the triggering of primary production and the biomass redistribution and export.

Figure 1 .
Figure 1.Left: study area with toponyms and main circulation features as known from literature.Right: model domain and bathymetry.The bathymetry used is the GEBCO at 30 of resolution.

Figure 2 .
Figure 2. Top to bottom: BIAS, RMSE and ACC for BF (blue) and BFC (red dashed) experiments.Units for BIAS and RMSE are •C, while ACC is dimensionless.

Figure 3 .
Figure 3. Map of SST RMSE (whole period) for BF (left) and BFC experiments.Units are • C.

Figure 4 .
Figure 4. Top to bottom: wind stress, sensible, latent and net heat flux differences between the two experiments (BFC -BF).Negative sign indicates lower values for BFC in respect to BF.

Figure 5 .
Figure 5. Difference map (BFC-BF) of the time-averaged wind stress (left) and net heat fluxes (right).Blue values indicate a BFC stress/heat lower than BF.Units are N m −2 and W m −2 respectively.

Figure 6 .
Figure 6.Total (left) and turbulent kinetic energy at surface.Red curve is for BF and green for BFC experiment.

Figure 7 .
Figure 7. Mean flow for BF (left) and BFC experiments.Units are m s −1 .

Figure 11 .
Figure 11.Wind stress work difference (BFC -BF).Units are W m −2 .Blue negative patches indicate where the wind power input is reduced by the feedback of currents on momentum fluxes.