During late January and early February 2012, a persistent cyclonic circulation associated with an exceptional cold anomaly dominated the Mediterranean region. Among the resulting effects, the northern Adriatic sea basin (NA) experienced a very large energy losses, mostly related to the intense and cold Bora winds blowing from north-east. Sea water temperature along the Italian coast dropped down to 6 • C, while part of the Venice lagoon got frozen. These series of exceptionally cold air outbreak episodes, as well as their effects on the NA circulation and dense water formation, are investigated by means of the Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) Modeling System, where the oceanographic model ROMS, the atmospheric model WRF and the wave model SWAN are coupled via MCT. In this specific application to the NA sea configuration, lasting from January 23 to February 23, 2012, particular emphasis was devoted to the analysis of the atmosphere-ocean-waves interactions. First, we employ the "stand alone" WRF atmospheric model in 4 different modes ("zero mode", i.e. using the skin temperature from the global atmospheric model without updating in the Sea Surface Temperature (SST); "static mode", i.e. retaining the January 23 radiometer SST; "dynamic mode", updating every 6 hours the SST as derived from radiometer data at 0.83 deg resolution; "OML mode", as above, but using a simple Ocean Mixed Layer model available within WRF to predict the temperature evolution). Second, the WRF-ROMS one-way forced case is analyzed, where no feedbacks to the atmosphere are provided from the ocean model ROMS, but momentum and heat fluxes are determined by WRF model. Then, the WRF-ROMS two-way coupled case is implemented (where the atmosphere model exchanges momen-tum and heat, and the ocean model exchanges SST with the Atmospheric model). Finally, the WRF-ROMS-SWAN two-way coupled case for waves-ocean-atmosphere is performed, where common variables are exchanged every 1200 seconds. While WRF-ROMS communicate ad described above, the atmospheric model sends the wind components to the wave model and receives back the wave parameters, from which it derives the surface roughness. On the other hand, the wave model sends to the ocean model the wave parameters, and ROMS provides velocity components, surface height and bathymetric variation towards the wave model. To our knowledge this is the first example of such a fully two-way coupled approach to wave, atmosphere and ocean domains in the NA basin; the work aims therefore at disentangling the different contributions of the interplay between ocean-atmosphere. Although uncoupled models were able to capture the cold air intrusions and the evolution of the heat fluxes, their magnitude resulted overestimated when compared to values derived from observations. The heat fluxes and air temperature at sea level, heavily affected by the SST used to force the atmospheric model (especially in a shallow, semi-enclosed basin such as the NA sea), reconciled better to observations when the atmosphere and the ocean were mutually coupled. More specifically, the atmosphere-ocean coupled model produces air temperatures that are in better agreement with measurements with respect to those of the uncoupled case. Besides, the estimate values of the heat fluxes during the peak of the event improved significantly, going from -1200 W/m2 to about -900 W/m2, very close to those obtained using data collected at the Ocean TOwer "Acqua Alta". We also noted a more spatially coherent distribution of the heat fluxes. The fully coupled case including also wave dynamics provides a further improvement in the estimate of the above mentioned parameters, leading to a better description of extreme events during the simulated period (namely, previous cases failed in reproducing the highest and lowest values of the heat fluxes and temperatures).