Modeling the ocean and atmosphere during an extreme bora event in northern Adriatic using one-way and two-way atmosphere–ocean coupling

We have studied the performances of (a) a two-way coupled atmosphere–ocean modeling system and (b) one-way coupled ocean model (forced by the atmosphere model), as compared to the available in situ measurements during and after a strong Adriatic bora wind event in February 2012, which led to extreme air–sea interactions. The simulations span the period between January and March 2012. The models used were ALADIN (Aire Limitée Adaptation dynamique Développement InterNational) (4.4 km resolution) on the atmosphere side and an Adriatic setup of Princeton ocean model (POM) (1/30× 1/30 angular resolution) on the ocean side. The atmosphere–ocean coupling was implemented using the OASIS3-MCT model coupling toolkit. Two-way coupling ocean feedback to the atmosphere is limited to sea surface temperature. We have compared modeled atmosphere–ocean fluxes and sea temperatures from both setups to platform and CTD (conductivity, temperature, and depth) measurements from three locations in the northern Adriatic. We present objective verification of 2 m atmosphere temperature forecasts using mean bias and standard deviation of errors scores from 23 meteorological stations in the eastern part of Italy. We show that turbulent fluxes from both setups differ up to 20 % during the bora but not significantly before and after the event. When compared to observations, two-way coupling ocean temperatures exhibit a 4 times lower root mean square error (RMSE) than those from one-way coupled system. Two-way coupling improves sensible heat fluxes at all stations but does not improve latent heat loss. The spatial average of the two-way coupled atmosphere component is up to 0.3 C colder than the one-way coupled setup, which is an improvement for prognostic lead times up to 20 h. Daily spatial average of the standard deviation of air temperature errors shows 0.15 C improvement in the case of coupled system compared to the uncoupled. Coupled and uncoupled circulations in the northern Adriatic are predominantly wind-driven and show no significant mesoscale differences.


Introduction
The Adriatic Sea is a semi-enclosed basin in the north of the central Mediterranean, oriented in the SE-NW direction (Fig. 1).The northern part of the Adriatic consists of a shallow shelf with depths not exceeding 60 m.Further southeast the depths in the middle Adriatic reach 260 m in two Jabuka pits, while the southernmost part, the south Adriatic basin, is the deepest part of the Adriatic with depths over 1200 m.The eastern coast of the Adriatic (Croatia) is mostly rocky, steep, with an abundance of rocky islands, while the western Adriatic coast (Italy) is low and sandy with several lagoons (Artegiani et al., 1997).
The northernmost part of the northern Adriatic is the shallow, semi-enclosed Gulf of Trieste.The main source of freshwater in the Gulf of Trieste is the Soča/Isonzo river, while other important riverine inputs to the northern Adriatic water stem from the Tagliamento, Livenza, Piave, Adige and, most importantly, Po, rivers.Buoyant spreading of their low salinity waters tends to generate a classical coastal current, which, under the influence of Earth's rotation, remains more or less confined to the western Adriatic (Italian) coast where it flows south as a part of Western Adriatic Current.Adriatic river climatological discharges, used ubiquitously in the Adriatic modeling community, were compiled in Raicich (1994), while another, more recent, work investigated circulation, dense water formation and its spreading along the Adriatic basin during the February 2012 bora, using new Croatian river climatologies (based on 2009-2011 river-runoff observations) as hydrological forcings in their circulation model (Janeković et al., 2014).
Dominant wind forcings can be separated into two classes: the cold northeasterly bora wind and the warm southeasterly scirocco.The bora denotes a predominantly winter occurrence of strong, cold and dry continental air flowing over the Dinaric orographic barriers over the Adriatic sea.The generating mechanism of the bora was thought to be of katabatic origin (Kuzmić et al., 2007), but several investigations now show that severe episodes of bora are generated by orographic wave-breaking (Grisogono and Belušić, 2009;Tudor and Ivatek-Šahdan , 2010).The bora seaward airflow is channeled mainly over Trieste, Senj and other gaps in the Di-naric Alps (see Fig. 1 for jet locations) on the eastern coast of the Adriatic (Cushman-Roisin et al., 2001), producing the well known wind jets over the ocean surface, which generate surface offshore water removal in the east and downwelling on the west Adriatic coast (Kourafalou, 2001).Bora wind episodes lead to intense air-sea interactions, which increase the net upward heat fluxes through the ocean surface and induce negative buoyancy flux on the northern Adriatic shelf (Raicich et al., 2013), thus making it one of the most important Mediterranean dense water formation sites.
The Adriatic sea has often been investigated using oneway coupled atmosphere-ocean modeling chains.A series of observations and simulations have been used routinely to analyze the occurrence of the well-known double gyre system, occurring north of the 45 • N parallel during intense bora events (Kuzmić et al., 2007;Zore-Armanda and Gačić, 1987).Measurements (Zore-Armanda and Gačić, 1987) suggested a positive surface current curl to the north of and a negative current curl south of Rovinj.This pattern was reproduced numerically by several modeling groups (Kuzmić et al., 2007;Paklar et al., 2001).It was also reported that the NW cyclonic gyre, the so-called Trieste gyre, has a barotropic current field while the anticyclonic gyre south of Rovinj is baroclinic in nature (Kuzmić et al., 2007).These findings indicate a controlling influence of the shallower topography along the Italian coast.This was further explored to quantify the extent of topographic control of wind-driven circulation during bora and scirocco (Malačič et al., 2012).It was shown that the topographic influence during bora episodes is substantial.
Interannual simulations of the Adriatic deep-water formation indicate that it is the atmospheric forcing on synoptic spatial and temporal scales that determines the characteristics and variability of the Adriatic dense water (Mantziafou and Lascaratos, 2008).Consequently incorporation of an ocean feedback into the atmospheric model, namely, the two-way atmosphere-ocean coupling, should lead to improvements in modeled ocean and atmosphere response during synoptic events and their consequent relaxation on longer timescales.Two-way atmosphere-ocean coupling approach to modeling the Adriatic was to the best of our knowledge first implemented to achieve realistic simulations of the Adriatic baroclinic circulation during bora episodes in the winter and spring 2001 (Pullen et al., 2003(Pullen et al., , 2006) ) and to study marine atmospheric conditions over the Adriatic during February 2003 (Dorman et al., 2006).A coupled atmospherewave-ocean modeling system (without ocean feedback to the atmosphere) over the Adriatic domain was recently used to simulate the formation and relaxation of Adriatic dense water during and after February 2012 bora event (Benetazzo et al., 2014).All these studies indicate that two-way atmosphereocean coupling could be beneficial for modeling the air-sea state during periods of intense atmosphere-ocean interactions.
As already noted, the February 2012 bora was a landmark event (Raicich et al., 2013;Mihanović et al., 2013) due to both its strength and its duration.Consequently, there were several routine (and ad hoc) field campaigns performed in the Gulf of Trieste and in mid-Adriatic immediately before and immediately after the strongest bora episode, and the event itself was subsequently widely reported (Raicich et al., 2013;Mihanović et al., 2013) as well as modeled (Janeković et al., 2014;Benetazzo et al., 2014).Atmospheric forcing caused exceptionally large surface heat fluxes, leading to an almost unprecedented water cooling and dense water formation in the Gulf of Trieste, where the sea temperatures dropped to 4 • C and density anomalies exceeded 30.55 kg m −3 (Raicich et al., 2013;Mihanović et al., 2013).Intensive air-sea interactions and the substantial number of observations make this event a suitable candidate for verifying two-way atmosphere-ocean coupled simulations.
In this paper we present two-way atmosphere-ocean coupling simulations of the February 2012 bora, along with comparisons with an identically set up, one-way coupled system.In Sect. 2 we present the two numerical systems in detail, focusing first on the one-way coupled configuration and then on the two-way coupled configuration.In Sect. 3 the observational data are described, together with dates and locations of the CTD (conductivity, temperature, and depth) casts in early 2012.In Sect. 4 the results of both one-way and twoway coupled numerical systems are presented and discussed, focusing on air-sea fluxes and air and sea temperatures, and we compare these with observations.Concluding remarks are presented in Sect. 5.The version of the model used for the experiments in this paper is that of the currently operational version used at Slovenian weather service.It runs on a 432 × 432 horizontal Lambert conic conformal grid with 4.4 km resolution and 87 vertical levels with the model top at 1 hPa and model integration time step of 180 s.The model domain spans 0.7 • W, 28.6 • E in longitude and 37.4 • N, 55.0 • N in latitude.The physics package used in the model uses Modular, Multi-scale, Microphysics and Transport (3MT) structure (Gerard et al., 2009).Initial conditions for the model are provided by atmospheric analysis with 3-hourly threedimensional variational assimilation (3D-Var) (Fischer et al., 2005;Strajnar et al., 2015) and optimal interpolation for surface and soil variables.Most of the conventional and satellite observations are assimilated, including surface stations (Synop) and ship reports, radiosondes, Automated Meteorological Data Relay (AMDAR), atmospheric motion vectors (AMV), Advanced Microwave Sounding Unit (AMSU)-A and -B, Metop-A and -B polar orbiters and water vapor channels from the geostationary Meteosat Second Generation Spinning Enhanced Visible and Infrared Imager (SE-VIRI).Observations within ±1.5 h from the analysis time are assimilated.Information at the domain edge is obtained from the global model by applying Davies relaxation (Davies, 1976).Lateral boundary conditions (LBCs) are provided by the ECMWF boundary conditions optional project.LBCs are applied with a 1 h period in the assimilation cycle and a 3 h period during model forecasts.Boundary condition information is interpolated linearly for time steps in-between these times.Within one-way coupled ALADIN (from now on referred to as ALADIN1w, while the two-way coupled version will be referred to as ALADIN2w) itself there is no implementation of ocean component; ocean variables values are therefore kept constant during the model run.In the operational (one-way coupled) ALADIN1w in Slovenia, sea surface temperature (SST) is initialized from the most recent host model analysis of the ECMWF model that uses Operational Sea Surface Temperature and Sea Ice Analysis (OS-TIA; Donlon et al., 2012), supplied by the National Environmental Satellite, Data and Information Service (NESDIS) of the American National Ocean and Atmospheric Administration (NOAA).OSTIA combines satellite data with in situ measurements by applying optimal interpolation.In AL-ADIN1w the SST remains constant over each 24 h of the model run.
In both ALADIN1w and ALADIN2w, turbulent fluxes over the ocean surface are parameterized in the usual fashion as follows.Wind stress is calculated as where ρ a is the moist air density, c pa is the specific heat of air at constant pressure, C D is the drag coefficient and V is the surface wind velocity.Sensible heat flux is parameterized as where C H is the sensible heat transfer coefficient, T air is air temperature at the lowest model level (10 m approximately) and T sea is sea surface temperature.Latent heat flux is calculated as where L is latent heat of vaporization, C E is the vapor transfer coefficient, q air is the atmospheric specific humidity at the lowest model level and q sat (T sea ) is the saturated specific humidity at the sea surface temperature.The turbulent transfer coefficients C D , C H and C E are computed using a modified Louis scheme that takes into account different roughness lengths for heat and momentum (Mascart et al., 1995).
Net radiation at the surface, R n , is computed as where j SW is the net shortwave radiation at the ocean surface, = 0.96 is the atmosphere and ocean long-wave radiation emissivity factor, R atm↓ is the atmospheric downward longwave radiation, σ is the Stefan-Boltzmann constant and T sea is the sea surface temperature.
The total heat flux Q T presented in this paper was obtained from (5) Throughout the paper, fluxes are depicted as positive when they represent an energy gain for the ocean and as negative when they represent energy loss from the ocean.

POM: one-way coupled system setup
The ocean model for one-way and two-way coupling is the Princeton ocean model (POM) (Blumberg and Mellor, 1987;Mellor, 1998).Our configuration of POM was set up for the Adriatic domain with 271 longitude cells spanning 12.0-21.0• E and with 209 latitude cells spanning 39.0-45.9375• N. Grid horizontal angular resolution is 1 • /30 × 1 • /30 and its vertical discretization contains 25 σ layers.We will refer to this one-way coupled Adriatic setup of the model as POM1w and the two-way coupled system as POM2w.As noted below, POM2w setup inherits most of its features from the one-way coupled setup.
River outflows, depicted in (Fig. 1), were included in both POM1w and POM2w using climatologies for the rivers Stella, Tagliamento, Livenza, Piave, Sile, Adige, Reno, Marecchia, Vibrata, Ofanto, Bojana, Drini, Shkumbi and Seman from (Raicich, 1994).Climatologies for rivers Mirna, Raša, Crikvenica, HPP (hydroelectric power plant) Senj, Zrmanja, Krka, Jadro, Cetina, Neretva, Ombla and HPP Dubrovnik were taken from (Janeković et al., 2014).Po River discharges were included using daily means of hourly discharge observations from Pontelagoscuro station.Soča/Isonzo river discharge was implemented using hourly operational HFS (hydrological forecasting system) forecasts from the hydrological forecasting system for the Soča/Isonzo River (see subsection 2.4 below).Rižana climatology was created from measurements provided by the Slovenian Environment Agency (ARSO).River temperatures were not defined and at the river mouth river discharges attain the temperature of the surrounding ocean, but with zero salinity.Po River discharge was split into five separate neighboring cells, partly to mimic the Po River delta but also because high Po River discharges (above 2000 m 3 s −1 ), when applied to a single discharge cell, sometimes cause violations of the Courant-Friedrichs-Lewy (CFL) stability criteria.
Lateral boundary conditions at the open boundary are applied for temperature, salinity, sea surface elevation and zonal and meridional velocity components.All these are taken daily from MyOcean Mediterranean Forecasting System (MFS;Tonani et al., 2009) and the model is forced with them for 24 h of model runtime.Every 24 h of model runtime, the simulation run was terminated and ocean models were hot started from their previous-step restart files.Both coupled and uncoupled ALADIN setups are initialized from the same 3-hourly 3D-Var data assimilation cycle using OSTIA SST as the surface boundary condition.Tidal forcing in both ocean models is provided by the Oregon State University Tidal Prediction Software (OTPS), set up on the POM1w/2w grid, but tidal effects are not considered in this paper.
Both setups were initialized on 6 November 2011, 00:00 UTC, with temperature, salinity, sea surface elevation and zonal and meridional velocity components of the My-Ocean MFS (Tonani et al., 2009).
POM1w mode splitting is inherited from POM: the external time step, used for the barotropic mode computation, is set to 1 s; the internal time step, used for the baroclinic mode calculations, is set to 90 s.These time steps and grid resolutions were set in order to satisfy the CFL condition for numerical stability.The vertical turbulence closure scheme is the usual 2.5 level Mellor-Yamada; horizontal diffusion is treated using the Smagorinsky formula.
Atmospheric input in POM1w is implemented through the input of hourly ALADIN1w fields contained in ALADIN1w gridded binary files.POM1w bulk flux computations were disabled and the POM1w model obtains all its fluxes from the ALADIN1w model.

Two-way atmosphere-ocean coupled system setup
Apart from the coupling interface, the two-way coupled setup of POM (referred to as POM2w) is identical to the POM1w setup.The spatial oceanic grid and vertical discretization of the POM2w model domain remain the same as with POM1w, as do the river discharges and turbulence closures.The baroclinic time step in POM2w was adjusted to match the AL-ADIN2w computational time step of 180 s.The exchange of coupling quantities between ALADIN2w and POM2w was enforced at each computational time step, namely, every 180 s of coupled system runtime.As in POM1w, the bulk heat flux computation in POM2w was disabled to be consistent with the POM1w setup and to obtain a balanced heat flux budget in the atmosphere-ocean system.Heat fluxes are computed only by ALADIN2w and are sent to POM2w via a coupling exchange scheme.
The atmosphere-ocean coupling itself is implemented using the OASIS3-MCT model coupling toolkit (Valcke, 2013) with the data flow as shown on Fig. 2. In the Figure, the rounded rectangles represent distinct OASIS objects, which are effectively treated by OASIS as independent models, even though that is only true in the case of POM2w and ALADIN2w, and is not true for pseudo-MFS and pseudo-MERGER, as will be explained immediately below.
The coupling scheme uses four OASIS objects with domains shown on Fig. 3.A more detailed description of all OASIS objects is as follows: 1. ALADIN.Atmospheric model.Receives the SST field from the pseudo-MERGER pseudo-model and sends the computed mean sea-level pressure, air temperature, precipitation, wind speed (u and v directions), humidity, sensible and latent heat fluxes, and shortwave and longwave downward radiation fields to the POM model.

POM. Ocean circulation model.
Receives mean sealevel pressure, air temperature, precipitation, wind speed (u and v directions), relative humidity, sensible and latent heat fluxes, and shortwave and longwave downward radiation fields from the ALADIN2w    N parallel between Italy and Greece, see Fig. 3).( 4) This merged and remapped MFS-POM2w SST field is sent from pseudo-MERGER to ALADIN2w, completing one coupling time step.
The ALADIN analyses used to initialize 24 h two-way coupled runs are taken from an independent data assimilation cycle, computed in advance.The same analyses are used in one-way and two-way coupled setups.This means that the initial conditions for ALADIN2w do not contain oceanographic information from the previous POM2w run; only OSTIA analyses are used for SST in the ALADIN assimilation cycle.The full two-way ocean-atmosphere assimilation system would consist of a fully coupled ocean-atmosphere model for integration during the data assimilation cycle of ALADIN as well: such a strategy was successfully implemented in (Pullen et al., 2007) and its implementation is left for future work.

The hydrological forecasting system on the Soča/Isonzo river basin
To the best of our knowledge, Soča/Isonzo river-runoff measurements at the river mouth are not available for the February 2012 bora event or for the preconditioning phase.We therefore decided to employ a HFS on the Soča/Isonzo river basin to remedy this lack of measurements.HFS has been operative at the Slovenian Environment Agency since early 2012 and provides Soča/Isonzo river hourly discharges for ocean circulation models in the numerical experiments presented in this paper.The most important input to the HFS system is that of the meteorological forecasts for precipitation and air temperature at ground level from the ALADIN and ECMWF atmospheric models, which are updated four and two times a day, respectively.The HFS is based on DHI MIKE11 engines -the hydrological model for rainfallrunoff simulations and the 1-D hydrodynamic (HD) model for river routing simulations.The rainfall-runoff model (Nielsen and Hansen, 1973)  The rainfall-runoff model was calibrated on the measured precipitation and discharge data in the period between 1998 and 2007.The precipitation data set contained daily and hourly values from approximately 115 meteorological stations from the entire river basin.The discharge data set was prepared for 19 Slovenian hydrological stations only as the discharge values for the Italian part of the basin (approximately one third of the basin area) did not exist.The Italian subcatchments were therefore assigned parameter values similar to the parameters of the calibrated Slovenian subcatchments that have common topographic and hydrological characteristics.The model calibration was performed on low and high flow conditions as well as on long-term water balance on the downstream-lying hydrological stations Soča/Isonzo Solkan and Vipava Miren in southwestern mainland Slovenia (not shown on the map).Model parameters were verified on the Slovenian hydrological stations with the 2008-2010 data set.

Observations
The oceanographic buoy Vida is a coastal observation platform, operated by the National Institute of Biology (NIB).It is located in the southern part of the Gulf of Trieste at 13.55505 • E, 45.5488 • N, as depicted in the inset of Fig. 1 (marked with a red cross).Data from the buoy are multifaceted (air temperature, air humidity, currents, waves, sea temperature, salinity, dissolved oxygen, chlorophyll concentration, etc.) and are transferred in real-time to the NIB server via an Ethernet link and are publicly available (http://www.nib.si/mbp/en/buoy/).A Nortek AWAC acoustic Doppler cur-rent profiler is situated on the sea bed beneath the buoy to monitor profiles of current (at 1 m intervals along the water column) and the sea floor temperature at a depth of 22.5 m.The sea temperature at 2.5 m beneath the surface is measured using a SeaBird 16plus Seacat salinity and temperature sensor with a sampling period of 300 s.Temperature measurements used for comparisons in this paper were downsampled in time to 1 h temporal resolution using hourly instant values.
An acoustic wind gauge (Gill instrument, WindMaster Pro Ultrasonic Anemometer) is also installed on the oceanographic buoy at approximately 5 m above sea level, and is used to measure wind speed and direction with a sampling frequency of 10 Hz.Wind speed measurements presented in this paper were downsampled in time to a 30 min temporal resolution using half-hourly instant values.
In 2012 the NIB made regular monthly CTD measurements in the Slovenian part of the Gulf of Trieste.These campaigns took place in the vicinity of buoy Vida (13.55505 • E, 45.5488 • N) on 26 January, 16 February, 27 February and on 12 March, providing temperature, salinity and density anomaly vertical profiles at buoy Vida before and after the bora episode.
The Paloma platform, located in the center of the Gulf of Trieste (45.618 • N, 13.565 • E) is operated jointly by the Regional Environmental Protection Agency (ARPA FVG) and National Research Council, Institute of Marine Sciences (CNR-ISMAR) (for details see Raicich et al., 2013).The observed quantities from Paloma used in this paper were seawater temperature at 3 and 15 m depths, air temperatures, relative humidity, mean hourly wind velocity and hourly mean pressure.
The Acqua Alta platform is located in the Gulf of Venice (45.314 • N, 12.519 • E) and is operated by the CNR-ISMAR.The following observational quantities from Acqua Alta were used to estimate in situ bulk fluxes during January and February 2012: air temperature, air pressure, wind speed and seawater temperature at 2 m depth.
To compute the bulk fluxes ALADIN uses an implicit scheme, which takes into account several prognostic quantities (like wind gustiness and friction velocity) that were not available within the existing observational data.This flux scheme therefore could not be used on observational data in the same manner as implemented in ALADIN itself.To compute the fluxes from measurements, we therefore used a COARE 2.0 algorithm (Fairall et al., 1996) that needs no additional data apart from those available.

Results and discussion
In the February 2012 bora episode, there was a strong bora with wind speeds above 20 m s −1 blowing from 5 to 9 February 2012, on 9 February the wind speed fell to about 7 m s −1 for about 10 h, and rose again to reach wind speeds of 15-20 m s −1 , peaking between 10 and 14 February.Between 14 and 19 February the bora subsided and wind speeds gener-

Air-sea fluxes
As noted above, both coupled and uncoupled system use the following ALADIN fluxes at the air-sea interface: sensible heat flux Q H , latent heat flux Q E and net radiation R n (composed from solar and infrared contributions) at the ocean surface.The only difference between coupled and uncoupled flux calculation in ALADIN stems from the sea surface temperature boundary condition in ALADIN: ALADIN1w calculates its fluxes using sea surface temperature from the OS-TIA product while ALADIN2w computes these fluxes using POM2w sea surface temperature.OSTIA SST in the northern Adriatic was consistently overestimated throughout the event (see green lines in all the panels of Fig. 11), which led to an overestimation of fluxes in the one-way coupled setup, as shown below.
Figure 5 depicts latent, sensible and total fluxes, averaged in time over the period of the bora event between 28 January-16 February.Latent heat flux (top row in Fig. 5) reaches a time-averaged value of approximately −200 W m −2 in the Trieste jet, about −450 W m −2 in the Senj jet, −400 W m −2 in the Karlobag jet and −230 W m −2 in the Sukošan jet (see Fig. 1 for the locations of these jets).The coupled system exhibits somewhat lower latent heat losses, leading to positive values of Q E2w − Q E1w (top row, right-hand panel in Fig. 5).This happens especially in the Gulf of Trieste, where the coupled system appears to lose about 40 W m −2 (amounting to approximately 20 %) less heat due to evaporation.In the Trieste jet wake, along the latitudes of Istria, the coupled system loses about 20 W m −2 (about 10 %) less than the uncoupled one.
In the Senj jet there are notable differences between the coupled and uncoupled setups in the bay of Kvarner and between the island of Cres and the mainland (see Fig. 1 for location).The coupled setup exhibits up to 75 W m −2 (amounting to approximately 20 %) less total heat loss in those regions.
In the coupled-uncoupled heat flux differences (see righthand column of Fig. 5) a bipolar zonal pattern occurs, reaching all the way to Italy at the latitudes 44.5-44.8] from ALADIN2w.First row right: difference in latent heat flux between ALADIN2w and ALADIN1w.Second row: the same as first row, but for sensible heat flux Q H . Third row: the same as first row, but for total heat flux Q T .(Note that the fluxes are generally negative during the event, implying ocean heat loss.Positive difference Q 2w − Q 1w thus means that the uncoupled system is losing more heat, while negative difference implies that the coupled system is losing more heat.)ing just below the tip of the Istrian peninsula: there is a northern zonal belt where the uncoupled system loses more heat than the coupled (the blue region starting right at the tip of Istria) and a southern zonal belt where the opposite is observed (the red region starting at Lošinj island, see Fig. 1 for location).This arises from the fact that the direction of the Senj jet in the two systems differs slightly.When comparing the two jets this leads to regions where the central flow of Senj jet form one model overlaps with a weaker lateral flow from the same jet from the other model.The same pattern is encountered when comparing sensible heat losses.
A significant difference between the latent fluxes is noticeable along the Italian coast.There is a coastal belt where an uncoupled system loses up to 50 W m −2 more heat than the coupled system.This happens due to OSTIA SST overestimation in the northern Adriatic during the bora outbreak as will be discussed in more detail in Sect.4.3.
Sensible heat losses (second row of pronounced sensible losses is in general approximately 30 nautical miles shorter than the zonal extension of evaporative losses.These described patterns dominate the time-averaged total heat flux Q T , which reflects them to a large extent as seen from the bottom row of Fig. 5.The time-averaged total flux can be multiplied with the respective time window (28 January-16 February) to obtain quantitative estimates of the differences of total energy lost in the two systems (not shown in the paper).Over the entire period of the bora event, the Gulf of Trieste, the Gulf of Kvarner and the western Adriatic (Italian) coastal belt lost approximately 100 MJm −2 less energy per unit area in the coupled system than in the uncoupled system.On the other hand, the zonal belt extending from the tip of Istria to the Italian coast lost about 20 MJm −2 more energy per unit area in the coupled than in the uncoupled system.The northern Adriatic west of Istria lost on average 30-50 MJ m −2 more energy per unit area in the uncoupled system than in the coupled system.
To estimate flux differences between the two setups in the time domain, we compared the flux time series from both setups at three observation platform locations in the northern Adriatic, the Vida buoy and the Paloma and Acqua Alta platforms (see the inset of Fig. 1 for their locations).We further compared the modeled fluxes from the two setups with bulk fluxes computed from observations at all three platforms.These results are shown in Figs.6-8.
At the Vida and Paloma platforms, the following observational parameters were used to compute the heat fluxes: sea temperature at 2 m depth at Vida or 3 m at Paloma (note that  the water columns at platform locations exhibit negligible stratification during strong bora events), air pressure, wind speed, relative humidity and air temperature.At Acqua Alta station, measurements were used for all parameters except relative humidity.The latter, for flux computations at Acqua Alta location, was taken from the ALADIN1w model.Error estimates for fluxes, computed from observations at the platforms, were determined in the following way.Uncertainties in observed air temperatures, sea temperatures, relative humidity, atmospheric pressure and wind speed were adopted from Raicich et al. (2013) (note that the very same data sets are used at Vida and Paloma as in Raicich et al., 2013).The uncertainties of air temperature and SST were thus set to be 0.3 and 0.1 • C, respectively.Uncertainty in relative humidity was set to 3 % while that of sea level air pressure was set to 0.3 hPa.A relative error of 12 % was set for wind speed.Sensible and latent heat flux error estimates were then computed using the Coupled Ocean-Atmosphere Response Experiment (COARE) algorithm on an ensemble of observed air temperature, sea temperature, relative humidity, atmosphere pressure and wind speed measurements that were perturbed by the values of their respective uncertainties stated above.When comparing sensible and latent heat fluxes from AL-ADIN setups to those from COARE 2.0, it has to be remembered that flux comparisons are generally highly sensitive to parameterization, but Louis and COARE schemes agree reasonably well (Eleuterio, 1998).Such was also the outcome of our comparisons, as shown in Figs.6-8.On the other hand, the fluxes in this paper differ from those presented in (Raicich et al., 2013) in which the Kondo scheme was used for transfer coefficients, by as much as 30 % during the peak periods of bora  in this paper to Fig. 6 in Raicich et al., 2013).As demonstrated below, however, ocean cooling in the northern Adriatic is well reproduced in POM2w (and indeed even overestimated in POM1w) using our values of fluxes.
As judged from observations at Vida, Paloma and Acqua Alta, it appears that air-sea coupling has a positive impact on the sensible heat flux Q H (see middle left panes in Figs.6-8): at all three stations the sensible heat fluxes, as computed from observations, are better reproduced in ALADIN2w.The latent heat fluxes, computed from observations, are on the other hand somewhat better reproduced by ALADIN1w.The root mean square differences (RMSD) between ALADIN1w(2w) setups and fluxes, computed from observations, are shown in Table 1.
It should be noted that the fluxes in these images are depicted as negative when the ocean is losing heat.A positive total heat flux difference Q T2w − Q T1w > 0 between a coupled and an uncoupled setup therefore implies that the uncoupled system is losing more heat than the coupled one.This was the case at the Paloma and Acqua Alta stations during the bora event.Total heat flux difference during the event reached 150 W m −2 at Paloma (bottom right panel in Fig. 7) and 100 W m −2 at Acqua Alta (bottom right panel in Fig. 8), amounting to as much as 20 % of the total flux during the bora outbreak.Before and after the bora event the air-sea coupling does not appear to have created a significant impact on the fluxes.
At the Vida location, fluxes from the two setups do not differ significantly, even during the bora event (bottom right panel in Fig. 6).This is most probably due to its location, which, unlike Paloma and Acqua Alta, is very close (within one ALADIN or POM grid point) to the land, and the impact of the ocean feedback is consequently less pronounced in the atmospheric model.

Atmosphere temperatures
The temperature at the lowest model level (approximately 10 m) fields over the Adriatic part of the computational domain was determined for various forecast ranges on 10 February 2012 (Fig. 9).The differences between near-surface temperatures obtained from ALADIN1w and ALADIN2w, shown in the right column of Fig. 9, increase with time, which is an expected consequence of the fact that both AL-ADIN1w and ALADIN2w are initialized from the same initial conditions at the beginning (00:00 UTC) of each 24 h forecast run.At a forecast range of 1 h (top row of Fig. 9), the differences are limited mostly to the shallow coastal regions (especially along the Italian Adriatic coast) where the differences in SST between OSTIA and POW2w are largest.At this time the near-surface air temperatures over land in the two cases are almost identical.Later on during the forecast these differences increase, as shown in the right-hand column of Fig. 9).The air masses from the two setups are advected along with the dominant wind, i.e., towards the southwest, and the pattern of their respective differences resembles the wind field (e.g., pronounced jet structures perpendicular to the eastern Adriatic coast in the right column of Fig. 9).These temperature differences are also seen across the Italian Peninsula and over the Tyrrhenian Sea, where the differences in near-surface temperature exhibit a characteristic dipole-like pattern, which is related to displacement of the precipitation features around the center of a low-pressure area (shown as a "+" character in the middle and bottom panels of the right column of Fig. 9).This difference, arising as a non-linear response of the atmospheric models to the differences between OSTIA and POM2w, is shifted eastward into the southern Adriatic at a forecast time of 23 h, in accordance with the trajectory of the pressure low.
This progressive impact during the forecast is further confirmed by objective verification of 2 m temperature forecasts.The mean bias and standard deviation of errors scores for 2 m temperatures, computed for 23 downwind meteorological stations (see right panel of Fig. 10 for their locations) in the eastern part of Italy for all days in February 2012 are shown in Fig. 10.The two-way coupled experiment is considerably and systematically colder than the one-way coupled setup, which is an improvement for prognostic lead times up to 18 h.The difference in standard deviation of errors is rather small after 24 h of forecast, but still shows a slight improvement in case of ALADIN2w compared to AL-ADIN1w.
The above suggests that the improvements due to two-way air-sea coupling are not due to more detailed ALADIN2w model dynamics but mostly due to improved forecast of the temperature of the air advected across Adriatic.The impact on other meteorological variables and/or other regions is much smaller.This is somewhat expected, since the AL-ADIN atmospheric component communicates asymmetrically with the ocean model in the coupled system -only SST is fed back from POM2w to ALADIN2w.The improvement of the temperature error variance and the impact of ocean on other meteorological fields might be greater during the warm part of the year, when small-scale atmospheric processes are more important and their initialization could be a direct consequence of ocean dynamics.

Ocean temperatures
Figure 11 depicts the ocean temperature time series, while Fig. 12 shows model-observation temperature differences and their respective biases and root mean square errors (RM-SEs) at Vida, Paloma and Acqua Alta stations.Ocean temperatures are reproduced significantly better in the POM2w, with RMSEs up to 4 times smaller than those in POM1w, as shown in Fig. 12.These statistics are also gathered in Table 2.
The latter serves as an illustration of the significance of the atmosphere boundary condition over the sea surface, especially in setups where heat fluxes from the atmospheric model are the only ones used in the ocean model.The OSTIA sea surface temperature boundary condition in the northern Adriatic was consistently warmer than observations by more than 2 • C in the Gulf of Trieste throughout the bora event (see the green line in the first three rows in Fig. 11).On the other hand, OSTIA seawater temperatures at Acqua Alta station reproduces the sea surface temperature fairly well (see the green line in the bottom row in Fig. 11).OSTIA overestimation of SST led to a modeled cumulative ocean heat loss and   cooling in the Gulf of Trieste and along the Venetian coast during the event (see Fig. 12) in the POM1w.The uncoupled system at Paloma location was losing about 130 W m −2 more heat per unit area than the coupled setup (see bottom right panel in Fig. 7).At Acqua Alta, however, the amounts of these losses are lower, of the order of 70 W m −2 (see bottom right panel in Fig. 8), which could explain the absence of negative ocean temperature drift during the bora outbreak at this location (compare the red lines in all the panels of Fig. 12).
A possible way to mitigate such drifts in uncoupled systems could be the introduction of flux corrections of the form where (∂Q T /∂T ) is the temperature derivative of the total flux from Eq. ( 5), T CLIM is climatological sea surface temperature field and T OSTIA is the atmospheric model sea surface temperature boundary condition used to compute bulk fluxes.Such flux corrections would effectively nudge the uncoupled ocean model temperatures towards its climatology and at least in principle prevent excessive cooling.Of course in this case climatology needs to be computed in advance and is itself burdened with errors.Our work shows that twoway atmosphere-ocean coupling presents an alternative, and more consistent solution to temperature drifts caused by airsea flux misestimation.
Sea surface temperatures at the end of the bora event on 14 February 2012, 00:00 UTC, from POM1w and POM2w are shown in Fig. 13.When compared to POM2w (bottom panel of Fig. 13), POM1w exhibits 1-2 • C cooler temperatures in the Gulf of Trieste and along the coastal belt of the Venetian Gulf (top panel of Fig. 13).Comparisons with ocean temperature CTD measurements at Vida location also reflect this discrepancy between POM2w and POM1w, as shown in Table 3. POM1w exhibits a persistent cold bias, which increased substantially after the bora event.This POM1w behavior is proposed to be stemming from the overestimation of sea surface temperatures in the atmospheric model OSTIA boundary conditions (green lines in Fig. 11), which leads to excessive heat losses and excessive cooling in the Gulf of Trieste and along the narrow Italian coastal belt.The models and observations shown in Fig. 11    and that POM1w sea temperatures at the Paloma station were consistently the coolest of all three stations.
A double gyre is related to the wind stress curl in the northern Adriatic, arising from pronounced Trieste and Senj jets to the north and south of Istria and moderate to low wind speeds in between.It was often reported (i.e., Kuzmić et al., 2007) that this wind distribution leads to formation of a cyclonic gyre in the northwestern (Venetian) part of the northern Adriatic, and of an anticyclonic gyre to the south of Rovinj, without any notable offshore surface currents along the west coast of Istria.
Circulation from POM2w at 1 and 12 m depths is depicted in Fig. 14.A cyclonic gyre is positioned at roughly (45.15 • N, 12.9 • E), while its anticyclonic counterpart is much smaller and located closer to Rovinj at about (−0.1 • N, 0.5 • E) to the southeast with respect to the center of the cyclonic gyre.In the northern arm of the cyclonic gyre along the Italian coast and along the Venetian coast surface current velocities attain up to 0.5 m s −1 .At 12 m depth, the cyclonic gyre current velocities are less than 0.25 m s −1 .Currents at 12 m depth further indicate a compensating inflow of water into the Gulf of Trieste along its southern coastline.
The northern Adriatic surface circulation patterns from POM1w are very similar and they are not shown explicitly.This similarity in the POM1w and POM2w circulations is due to the fact that circulation in the northern Adriatic during this bora outbreak was completely wind dominated.This can be shown by estimating the Ekman depth during the bora outbreak.Wind speeds of bora in the Gulf of Trieste during the period 9-14 February 2012 were measured to be well above 10 m s −1 for most of the time; see Fig. 4. Following Malačič and Petelin (2009) the Ekman depth in the Gulf of Trieste can be expressed in terms of wind speed and they estimated its value with the following: where ρ air = 1.29 kg m −3 and ρ sea = 1029 kg m −3 are air and seawater densities, C D = 2.6 × 10 −3 is the drag coefficient,  14).This circulation feature is well known and is also confirmed by high-frequency radar measurements in the Gulf, see Cosoli et al. (2013).Due to land-locked nature of the Gulf, however, this Ekman layer outflow is modified by the compensating inflow, which needs to take place at larger depths along the southern part of the Gulf to replenish it with water.These cool water masses from the Gulf of Trieste and Venetian coast then recirculate in the northern Adriatic due to bora induced cyclonic-anticyclonic circulation (see Fig. 14) and re-enter the Gulf of Trieste some 10 days later.It appears therefore (see Figs. 13 and 14) that higher POM2w temperatures at Vida location are at least partly due to bora induced recirculation, which leads to the intrusion of warmer water along the southern part the of the Gulf of Trieste.Warm water intrusion to the Gulf is clearly visible in POM2w but absent from the uncoupled setup.Judging from POM2w ocean temperatures (bottom panel of Fig. 13), Vida is located well within the warm water tongue while Paloma is on its lateral boundary (bottom panel of Fig. 13).
In this paper we present one-way and two-way coupled modeling of ocean and atmosphere during an extreme bora event in the northern Adriatic in February 2012.Ocean feedback to the atmosphere in the two-way coupled system was limited to the sea surface temperature.The sea surface temperature boundary condition for the atmospheric model in the oneway coupled system was taken from the OSTIA SST product, provided by NOAA.The sea surface temperature boundary condition for the atmospheric model in the two-way coupled system was obtained from the ocean component of the twoway coupled system.Both systems used fluxes computed in the ALADIN atmosphere component.OSTIA SST was overestimated throughout the bora outbreak, leading to overestimation of air-sea flux and excessive ocean cooling in the one-way coupled system.
We used in situ observations from three marine observation platforms (Vida, Paloma and Acqua Alta) in the northern Adriatic to compute in situ fluxes and to compare them to modeled fluxes.When compared to observations sensible heat fluxes were better modeled by the two-way coupled system, while latent heat fluxes were somewhat better represented in the one-way coupled system.An improvement of sensible heat fluxes in the two-way coupled system was not unexpected, since SST was the parameter that was communicated from the ocean to the atmosphere component of the two-way setup.At Paloma and Acqua Alta sensible and latent fluxes differ significantly during the bora event while their differences are less pronounced before and after the bora.At Vida location flux differences between the two setups were not significant throughout the simulations.This is probably connected to its location, which is within 1 grid point from the land in both ocean and atmospheric model.These results indicate that the two-way coupling impact will be largest during periods of intense air-sea interactions and less pronounced otherwise.
We performed an objective verification of near-surface air temperatures from both setups, using measurements from 23 downwind stations in eastern Italy during each calendar day in February 2012.The two-way coupled system shows a lower temperature bias for prognostic lead times up to 18 h.Standard deviation of modeled near-surface temperature errors is also somewhat lower in the atmosphere component of the two-way coupled system.However, the atmosphere in the two-way coupled system appears to benefit mostly from a better representation of air temperatures over the ocean, which are later advected over mainland Italy during the bora, improving the skill score.Ocean temperature feedback in the two-way coupled setup did not induce any new dynamics in the atmospheric model during this extreme winter case.We expect this to happen during the summer when small-scale processes in the atmosphere are more important.We leave this for further investigations.
Using in situ platform and CTD observations we show that the two-way coupled system captures the observed ocean temperature response during the bora event significantly better than the one-way coupled system.This is due to OS-TIA SST overestimation during the bora, leading to excessive heat losses from the ocean, as noted above.These heat losses induced overcooling of water in POM1w and this cold bias persisted in the one-way coupled ocean model long after the cessation of bora.We have demonstrated that two-way atmosphere-ocean coupling offers an efficient and consistent way to address this issue.
Currently the initial conditions of the two ALADIN 24 h integration runs are provided by an independent ALADIN assimilation cycle with 3-hourly 3D-Var, which uses OS-TIA SST analysis for the ocean temperature boundary condition.A possible future improvement of the two-way coupled system would therefore be the inclusion of a two-way coupling in the ALADIN assimilation cycle also.On the other hand, there would be a lack of SST information from measurements, which now enter the two-way coupled system through the ALADIN assimilation cycle (apart from information reaching the area of interest from the lateral boundary conditions).Consequently, there would be obvious danger of a significant model drift in longer timescales, especially over several seasons.To avoid this and to ensure optimal use of SST measurements, a comprehensive data assimilation system for the ocean component would be needed.

Figure 2 .
Figure 2. Two-way atmosphere-ocean coupling scheme for one coupling time step.Rounded rectangles denote distinct OASIS models, effectively treated by OASIS as independent.The arrows denote exchanged coupling quantities and direction of the transfer.Time instant of any specific coupling exchanges within one coupling time step grows from left (earliest) to right (latest).

Figure 3 .
Figure 3. Domains of all OASIS models used in the atmosphereocean coupling scheme.White rectangle: atmospheric model AL-ADIN domain.Orange rectangle: ocean circulation model POM2w domain.Black rectangle: remapping pseudo model pseudo-MERGER domain.Blue rectangle: domain defined by the My-Ocean MFS NetCDF files, read by the pseudo-MFS pseudo model.

Figure 4 .
Figure 4. Hourly wind speeds measured at oceanographic platforms Vida, Paloma and Acqua Alta during January and February 2012.

Figure 5 .
Figure 5. Mean fluxes in the Adriatic during 28 January-16 February.First row left: mean latent heat flux Q E [Wm −2 ] from AL-ADIN1w.First row middle: mean latent heat flux Q E [Wm −2] from ALADIN2w.First row right: difference in latent heat flux between ALADIN2w and ALADIN1w.Second row: the same as first row, but for sensible heat flux Q H . Third row: the same as first row, but for total heat flux Q T .(Note that the fluxes are generally negative during the event, implying ocean heat loss.Positive difference Q 2w − Q 1w thus means that the uncoupled system is losing more heat, while negative difference implies that the coupled system is losing more heat.)

Figure 6 .
Figure 6.Daily averages of bulk fluxes [W m −2 ] at Vida location from ALADIN1w (red line) and ALADIN2w (blue line).Sensible and latent heat fluxes at Vida location were computed using COARE 2.0 algorithm from measurements at Vida: daily averages of these fluxes are depicted with black squares with vertical error bars (see text for details on error estimation).Top left: Q E -latent heat flux.Middle left: Q H -sensible heat flux.Bottom left: Q T -total heat flux.Top right: Q B -long-wave downward radiation.Middle right: Q S -Solar shortwave radiation.Bottom right: Q T2w −Q T1w -total heat flux difference between ALADIN2w and ALADIN1w.Fluxes are depicted as negative when the ocean is losing heat.

Figure 7 .
Figure 7. Same as Fig. 6 but for Paloma location.

Figure 8 .
Figure 8. Same as Fig. 6 but for Acqua Alta location.

Figure 9 .
Figure 9. ALADIN 10 m temperatures [ • C] in ALADIN1w (left column) and ALADIN2w (middle column) on 10 February 2012 at: 01:00 UTC (top row), 13:00 UTC (middle row) and 23:00 UTC (bottom row).Right column depicts temperature differences [ • C] in both setups.Black + signs in middle and bottom panels of the right column indicate the approximate location of the center of the low pressure area at 13:00 UTC (middle panel) and 23:00 UTC (bottom panel).

Figure 10 .
Figure 10.Left panel: bias (top) and standard deviation of errors (bottom) scores for ALADIN1w (red curve) and ALADIN2w (blue curve), depending on forecast length for temperature at 2 m in February 2012.The green curve, quantified by the right-hand axis, depicts the number of cases (available number of measurements at all stations at each forecast hour combined) used in the statistics.Right panel: red circles mark the geographical locations of all 23 meteorological stations used for verifications.

Figure 11 .
Figure 11.Modeled and observed temperatures at Vida, Paloma and Acqua Alta locations.Observed in situ temperatures are plotted with a black line, OSTIA SST boundary condition for ALADIN1w with a green line, POM1w modeled temperature with red line and POM2w modeled temperature with a blue line.First row: temperatures at Vida location.Second row: temperatures for Paloma location at 3 m depth.Third row: temperatures at Paloma at 15 m depth.Fourth row: temperatures for Acqua Alta locations at 2 m depth.

Figure 12 .
Figure 12.Differences, biases and RMSEs of modeled and observed temperatures at Vida, Paloma and Acqua Alta locations.Top left: sea temperature difference between observations and POM1w (red line) and POM2w (blue line) at Vida location.Top right: the same for Acqua Alta location.Bottom left: the same for Paloma location at 3 m depth.Bottom right: the same for Paloma location at 15 m depth.Respective biases and RMSEs are included below each panel.

Figure 13 .
Figure 13.Sea surface temperature [ • C] in POM1w (top panel) and POM2w (bottom panel) at the end of the bora event on 14 February 2012, 00:00 UTC.Stations locations are drawn in pink color as a square (Acqua Alta), circle (Paloma) and triangle (Vida).

2 Modeling systems 2.1 The ALADIN atmospheric model: a one-way coupled setup
a two time level semi-Lagrangian semi-implicit advection scheme with several choices of physical packages.Its dynamical core is shared with the Integrated Forecasting System of the European Centre for Medium-range Weather Forecasts (ECMWF) and of AROME (Application de la Recherche à l'Opérationnel à Méso-Echelle;Seity et al.,  2011).

www.ocean-sci.net/12/71/2016/ Ocean Sci., 12, 71-86, 2016
POM2w domain and this was achieved using MFS NetCDF files.The time step for exchanging fields is the same for all models -180 s.At the first time step, the models are initialized independently.From the second time step on, the fol-POM2w domain we opted for the use of MFS SST field instead of OSTIA SST field.This was done because POM2w is forced by the MFS and merging POM2w SST with OSTIA SST might introduce substantial SST discontinuities at POM2w domain boundary.To further reduce sharp SST gradients on the pseudo-MFS/POM2w boundary, a linear interpolation between MFS and POM2w SST fields is performed within pseudo-MERGER during the merging and remapping process.The interpolation of SST takes place in all cells located up to five grid cells north of the POM2w open boundary (39 of pseudo-MFS and pseudo-MERGER pseudo models into the OASIS scheme was necessary because the AL-ADIN2w domain extends well beyond the POM2w domain (see Fig.3).We therefore needed to provide an SST estimate to ALADIN2w in regions outside the www.ocean-sci.net/12/71/2016/Ocean Sci., 12, 71-86, 2016 is a deterministic, conceptual, lumped model describing the land phase of the hydrological cycle in a simplified form.The parameters of the physical and semi-empirical model formulations for the snow, surface, root-zone and groundwater storages constitute the average values for each of the sub-catchments.The hydrological Soča/Isonzo river basin model is defined in 30 subcatchments with an average area of 115 km 2 .In the model snow module the sub-catchments are divided into 100 m high elevation zones giving the model setup a semi-distributed character.The river network of the Soča/Isonzo river basin HD model is represented by 17 branches having a calculation point density of less than 2000 m in river length.The upstream and lateral boundary conditions of the HD model are the runoff hydrographs simulated by the model while the downstream boundary condition is the constant water level of the Adriatic Sea.

Table 1 .
Root mean square differences (RMSD) between modeled ALADIN fluxes (modified Louis scheme) and fluxes from in situ observations (COARE 2.0).

Table 2 .
Biases and root mean square errors (RMSE) [ • C] between modeled POM1w, POM2w ocean temperatures and in situ observations at Vida, Paloma and Acqua Alta observational platforms during February 2012.

Table 3 .
Biases and root mean square errors (RMSE) [ • C] between modeled POM1w, POM2w ocean temperatures and CTD observations at Vida location on four different dates before and after the 2012 bora outbreak.