Investigation of Air-Sea Turbulent Momentum Flux over the Aegean Sea with a Wind-Wave Coupling Model

Near surface turbulent momentum flux estimates are performed over the Aegean Sea, using two different approaches regarding the drag coefficient formulation, a wave boundary layer model (referred here as KCM) and the most commonly used Coupled Ocean–Atmosphere Response Experiment (COARE) algorithm. The KCM model incorporates modifications in the energy-containing wave spectrum to account for the wave conditions of the Aegean Sea, and surface similarity to account for the stratification effects. Airborne turbulence data during an Etesian outbreak over Aegean Sea, Greece are processed to evaluate the simulations. KCM estimates found up to 10% higher than COARE ones, indicating that the wave-induced momentum flux may be insufficiently parameterized in COARE. Turbulent fluxes measured at about 150 m, and reduced to their surface values accounting for the vertical flux divergence, are consistently lower than the estimates. Under unstable atmospheric stratification and low to moderate wind conditions, the residuals between estimates and measurements are less than 40%. On the other hand, under stable stratification and strong winds, the majority of the residuals are more than 40%. This discrepancy is associated with the relatively high measurement level, shallow boundary layer, and the presence of a low level jet.


Introduction
The atmosphere-ocean interface differs significantly from that of the atmosphereland interface due to the presence of surface gravity waves. Several studies [1,2] have demonstrated that Monin-Obukhov (MO) similarity theory, although developed over land, is applicable in the Surface Layer (SL) of the Marine Atmospheric Boundary Layer (MABL). However, the usual application of MO theory should be confined to the upper portion of the SL where turbulent fluxes are assumed to be almost constant with height, above the shallow layer that is directly influenced by the waves, the Wave-influenced Boundary Layer (WBL). In the WBL and outside the viscous sublayer, additional scaling parameters are required for similarity, such as the non-dimensional atmospheric forcing of waves, the angle between wind, and the dominant wave direction [3].
Several studies have demonstrated the wind stress dependence on the sea state (e.g., [4,5]), implying greater stress over young-developing seas compared to older seas which are in 'equilibrium' with the wind flow. Young seas occur with fetch-and/or duration-limited flows. The wave age parameter is usually used as the measure for the sea state evolution. Over the open sea, meteorologists rely on bulk formulas in order to relate the values of wind speed and temperature with their associated fluxes, through transfer coefficients. The most commonly used bulk formula is COARE [6][7][8], which was originally based on the parameterization of [9]. This type of formulation gives a robust relation between observed fluxes and mean In this study, two algorithms are employed over the AS, with different theoretical approaches regarding the drag coefficient formulation. In particular, the well-established bulk algorithm COARE 3.5 and the WBL model proposed by Kudryavtsev, Chapron, and Makin [33], hereinafter named KCM, are examined under various wind conditions and stabilities encountered during the AG. Their main difference is that COARE parameterizes the sea surface drag coefficient in terms of surface roughness, gustiness and atmospheric stability, while KCM is based on momentum conservation within the WBL and relates the drag coefficient directly to waves described statistically in terms of the directional wave spectrum.
The rest of the paper is organized as follows: Section 2 describes the methodology for the air-sea turbulent momentum flux estimation, through the algorithms and models applied in this research. In this section, flux measurements, in-situ wave measurements, the area of interest, the prevailing atmospheric conditions, and the numerical simulations are also described. Section 3 provides analysis of the results according to those observations. In particular, in Section 3.1, the wave simulations are provided in order to investigate the prevailing wave conditions in relation to observations during the AG campaign. In the next section (Section 3.2), the airborne covariance fluxes are discussed in association to their position to the nearest upwind island and the mean flow parameters. Also the two flux algorithms are compared and evaluated against the observations. Finally, Section 4 summarizes the main results and conclusions. Some results of the turbulent sensible heat flux are also demonstrated in the paper, mainly for stability classification purposes, but they were not modeled.

Materials and Methods
The fluxes over the AS during AG are examined by applying the flux-parameterization model KCM and COARE 3.5 algorithm. The default energy-containing wave spectrum of KCM, which is estimated by measured wind, has been replaced by a spectrum calculated by the Simulating Waves Nearshore (SWAN) model [34] in order to better represent the actual energy-containing spectrum under the prevailing conditions during AG. SWAN is driven by the Weather Research and Forecasting (WRF) Model [35]. Also, in a modified (KCM mod ) version of KCM, atmospheric stability effects have been included in the wind speed profile through the use of SL similarity functions that are used in COARE.

Numerical Models
KCM is a wind-over-wave-coupling model with two components, the water-side and the air-side, which are coupled via the wind forcing in terms of the growth rate (the energy rate transferred to the wind-waves by the wind). The form drag and the skin-friction drag are analyzed separately in order to derive the sea surface drag [33]. The form drag contribution is further split into the non-breaking (regular) and breaking waves. The growth rate functional form is based on [36] empirical parametrization, with modifications on the basis of the rapid distortion theory of [37,38].
The semi-empirical spectrum model (water-side) follows the approach of [12,39]. The total wave-number spectrum is defined as a composition of the energy-containing wave spectrum (long wind waves) [40] and the equilibrium spectrum (short wind waves). The equilibrium spectrum results from the solution of the wave action balance equation at uniform conditions (no surface current and steady wind). However, it has been constrained by [41] to reproduce the statistical properties of the sea surface, based on field stereophotograph measurements of short wind-waves.
Regarding the air-side, the momentum conservation equation integrates the impact of the surface-waves narrow band in terms of a wave-number vector. Regular waves contribute to the form drag through the relation between the surface pressure and the waves, while breaking waves contribute through the action of the pressure drop on the surface slope discontinuity (breaking front). Then, the momentum conservation is integrated over all wave components utilizing a simple turbulent closure scheme valid only for slow waves relative to the wind (no swell). The stability functions employed in KCM mod are identical to the ones used in the bulk-flux algorithm COARE, as mentioned above.
The COARE algorithm follows the standard MO approach and was initially developed based on the international TOGA-COARE field program [42]. Subsequent improvements to the algorithm [7,8] have been done since its initial release [6]. COARE parametrizes the sea surface drag in terms of surface roughness, gustiness, and atmospheric stability. The latest version (i.e., COARE 3.5) includes sea-state-dependent parameterizations regarding the surface roughness. The parameterization of the Charnock coefficient as a piecewise linear function of wind speed [8] was used in the present study. The use of parameterizations based on wave age as almost linear empirical functions of wind speed, which have been found to hold on average over the open ocean, cannot give significant improvement. The stability function, used for unstable conditions, is a combination of Kansas forms that blends smoothly for convective conditions [43]. The modifications for stable conditions are based on [44], with profile data taken over the Arctic ice cap [45].
The SWAN Cycle III wave model (version 41.10), developed by Delft University of Technology, is based on the wave action balance equation, including sources and sinks. It simulates the wave parameters in coastal regions given the bottom topography and driven by the wind speed at 10 m above mean sea level (msl) and the sea surface currents [46,47]. It also calculates the friction velocity employing the drag formulation of [48].

Measurements
The AG airborne data, from 31 August to 7 September 2011 are used to evaluate the simulations of this research. These observations, at a 32 Hz sampling frequency, were obtained with the UK Facility for Airborne Atmospheric Measurements BAe-146 research aircraft and cover the whole Aegean inside and above the Planetary Boundary Layer (PBL) (up to 4.5 km). In this study, we used the wind components measured with a radome probe, air temperature measured with a Rosemount sensor, and humidity (specific humidity) measured with a Lyman-Alpha hygrometer. Sea-surface skin temperature (SST) measurements were made with an infrared Heimann sensor for turbulent fluxes estimation. There were also available atmospheric measurements made by many slow response sensors. More details are provided in [30]. Three selected flights are used: (a) b640 (0559-0954 UTC), (b) b641 (1113-1535 UTC), both performed on 4 September 2011, and (c) b643 (0854-1408 UTC), performed on 7 September 2011 [30]. The atmospheric turbulent fluxes (777 samples flow level straight segments) have been calculated using the eddy-correlation method with a horizontal averaging length of 3 km, in order to include most of the energy-containing eddies [49] and to avoid spatial non-homogeneity effects. In most cases, the measuring height (~150 m), which was imposed by safety regulations, is above SL, given that the average height of MABL has been estimated from 200 to 400 m for stable conditions and from 500 to 900 m for unstable conditions [30]. Therefore, the measured turbulent momentum flux is expected to deviate from the corresponding surface value due to vertical flux divergence. In order to account for the reduction of the fluxes with height and to be comparable with the simulated surface fluxes, a linear height-correction is applied. The correction is based on the average behavior of stable and unstable cases between the surface and the zero-flux height [50], as presented by [49] and applied in [31]. Under unstable conditions, the zero-flux height was set to 700 m (average height of the top of well mixed MABL capped by a temperature inversion), while under stable conditions this was set to 250 m, which is the average height of the low level wind jet (LLJ) (near zero shear and low turbulence) that was usually observed in sounding legs under these conditions [30]. A higher order interpolation method would require more information than only an average MABL height, which was not available. The random error of this correction can be high, but on average this correction will bring the measured fluxes closer to the actual near-surface values. The potential virtual temperature difference between the sea surface and air ∆θ v = θ v (sea) − θ v (air) has been used to define the atmospheric stability.
Wave observations from the buoys network, operated by the Hellenic Centre for Marine Research (http://www.oceansites.org, http://www.poseidon.hcmr.gr accessed on 14 September 2021) are also available at 3-h intervals at 6 different locations (Figure 1b). The buoy data used in this paper were the significant wave height, first moment mean wave period, peak wave period, and mean wave direction. The water depth is less than 300 m in all six examined locations, except for buoy near Crete where the depth is 1073 m. Buoys' distance from the nearest shore ranges from 3.6 km (Lesvos) to 35.8 km (Crete), and only the Crete buoy location can be classified as "open-sea".

Meteorological Conditions
The AS is a semi-enclosed basin characterized by a complex shoreline with numerous islands of various sizes. The AS is mostly dominated by wind-seas [51] with intense sea surface temperature spatial variability [52]. The spatiotemporal variations of the wind and wave fields are strongly influenced by this complex topography.
The prevailing meteorological conditions during the AG were representative for an Etesian winds outbreak, which is a persistent northerly wind over the AS during a warm period, as described in detail by [30]. The measured wind speed at the observational level often exceeded 15 m/s with gusts over 20 m/s. Moving from west to east of the AS, the wind direction turned from north to north-western directions. The air temperature varied from 20 to 25 • C, while the SST from 17 to 27 • C. For the given period, the AS can be spatially divided in three distinct regions based on the predominant conditions: the northeastern area, with persistent strong winds (~16 m/s) under stable conditions due to lower SST (~21 • C), the central AS (Cyclades), with near neutral conditions, the strongest winds (up to 20 m/s), and the largest disparity in terms of wind speed and direction due to topographical effects, and the south and southeastern area, with persistent moderate winds (~12 m/s) and unstable conditions due to higher SST (up to 27 • C).

Numerical Applications
Area-averaged wind speed, air temperature, SST, and relative humidity from the AG flights (b640, b641, and b643) are used as input to both KCM versions and COARE 3.5, in order to simulate the surface fluxes. At the same time, SWAN non-stationary simulations of high temporal and spatial resolution (see below), driven by the surface wind fields from the WRF model (at 10 m above msl) without surface current information, are used to calculate the long-wave spectrum. The long-wave spectrum is then projected along the flight path by applying spatial bi-linear interpolation. For the WRF model set up, we follow [31]. The considered PBL parameterization scheme is the Mellor-Yamada-Nakanishi-Niino [53] combined with the Community Land Model version 4 land-surface parameterization scheme. This land-surface scheme is selected because it follows the COARE 3.5 formulation over water surfaces [8]. The wave simulations use the configuration for the physical process, as summarized in Table 1. SWAN has been applied from 4 to 7 September 2011 and the domain covers the whole AS ( Figure 1). Bathymetry data with a spatial resolution of 30 arc-seconds from the General Bathymetric Chart of the Oceans, GEBCO (information on the data set is available from the GEBCO project web pages, [54]). A curvilinear computational grid has been implemented with 995 × 895 meshes in "longitudinal-direction" (mean length 765 m) and "latitudinal-direction" (mean length 920 m), respectively. The computational spectral domain consists of 50 logarithmically distributed discrete frequencies in the range [0.01-5 Hz], and 72 equally distributed directions in the range [0-360 • ]. Stationary computations have been applied in order to provide initial conditions for the non-stationary runs. Spatial interpolation of the wind field is performed at every SWAN grid point, and the adopted wind-input time-step is set to 1 h.

Physical Process Formulation
Wave growth [11,55] White-capping [11] Quadruplet wave-wave interaction Fully explicit computation of the nonlinear transfer with Discrete Interaction Approximation per iteration [56] Bottom friction [57] Triad wave-wave interactions Lumped Tried Approximation [58] Depth-induced wave breaking Breaker index scales with both the bottom slope and the dimensionless depth [59] Numerical scheme Scheme for non-stationary runs [60]; Second Order Upwind scheme for stationary runs also known as a BDF scheme [61]

SWAN Sensitivity
The performance of the wave simulations is examined through the comparison with buoys and in particular the significant wave height (swh), the mean wave direction (mwdir), and the first moment mean wave period (T m01 ). The most energetic (peak) wave period T p has also been examined for the accuracy of the simulated spectral distribution. Furthermore, the criterion of [40] is applied to determine whether the peak wave train is a wind-sea. Initially, the performance of SWAN is examined with respect to computational time-step and space resolution as well as nesting approach.
Given that SWAN performance declines when wave energy travels many grid cells per time step, the time step should be chosen according to the peak wave frequency [62]. Therefore, 3 simulations with time steps of 60, 30, and 10 s are performed, all with the same initial wave conditions. The spatial or temporal differences between 10 s and 30 s are found negligible. In case of 60 s, however, the differences are significant in the leeward side of the islands (e.g., swh differ up to 0.5 m downwind of Crete), probably due to the additional resolving small scale waves that are involved in growing seas. Thus, the 30 s is considered to be an optimal choice in terms of accuracy and computational time.
As the AS is a region with complex shoreline and seabed topography, the spatial resolution is of critical importance, especially when island-induced wave-energy-blocking is present. Therefore, two runs are performed, a fine one with 765 × 920 m and a coarser one with 1.5 × 1.8 km spatial resolutions. In terms of integral wave parameters, no significant dependence on grid resolution is found, as a resolution improvement can be masked by possible errors due to the wind spatial interpolation [63]. However, the higher spatial resolution is chosen, as directional differences (up to 20 • ) are noticed regarding the peak wave frequency.
Finally, the SWAN model has been applied with and without nesting. In the case of nesting, the area of interest (innermost domain, Figure 1b) is nested in a parent domain comprising the Mediterranean Sea (Figure 1a). The analysis demonstrated that there are no significant spatial or temporal differences, except at the outer boundaries of the innermost domain, due to the absence of boundary forcing in the case without nesting. This is justified by the dominant northerly flow (Etesian flow) and the topography of the AS (semi-enclosed basin), hence no energy-containing wave-trains were traveling from south to north. Therefore, the non-nesting approach has been chosen due to the reduced computational time requirement.

Simulated Sea State during Etesian Conditions
The wave conditions during the AG campaign calculated by SWAN were determined by the Etesians. As displayed in Figure 2, both days exhibit similar wave condition; low swh values (<1 m) are found at the northern and northwestern part of the Aegean due to low winds and limited fetch, within the Cyclades complex due to wave blocking and spreading, and at the deceleration zone just north of Crete [30]. High swh values (>1.75 m) are found downstream of Chios and just north of the Cyclades region, downstream of the Kythira and Crete straits (southwest of Crete), and at the southeast of Crete, and associated with flow acceleration due to channeling effects between Crete and Karpathos islands [64].  Estimates of the full range saturation spectrum have been calculated by the composition of the SWAN and the equilibrium s, calculated by KCM mod using the airborne measurements. Figure 3, displays two indicative simulated spectrums from the two flight segments at the east-central (west of Chios) and southeast AS (near Karpathos). As displayed in Figure 3, the sea at the east-central Aegean is under-developed due to short fetch and low WRF surface winds, and the dominant waves (peak wave frequency) are propagating from north directions. At the southeast AS, a more developed wave spectrum with two peaks is found, as the prevailing surface wind moving southward is changing from north to west-northwest directions. In addition, the secondary peak located in the gravity-capillary range (~370 rad/m or~1.7 cm wavelength) is much greater in the east-central AS due to stronger winds at the observational level compared to the southeast AS.

Wave Simulation Comparison with In-Situ Measurements
A quantitative evaluation of SWAN and WRF models in relation to buoy data is summarized in Tables 2 and 3 for the period of 4 days. The statistical metrics applied for scalar (a, R, and snrmse) and directional quantities (mae and nrmse θ ), are described in the Appendix A.
WRF wind speed has also been extrapolated to the height of 3 m using similarity profiles, so that it can been compared with the buoy data at the same level. The WRF model appears to overestimate the wind speed at most locations. The largest overestimation is found at Skyros, while the model performed better at Mykonos and Lesvos with less than 0.3 m/s differences (Table 2). Regarding the directional wind pattern, the model, on average, reproduce well the Etesian outbreak. In terms of mean wind speed the model performed better at Lesvos, while at Saronikos located at the sheltered zone downstream of Attica region, it exhibits the highest differences, due to a mismatch in wind direction.  SWAN model tends to overestimate swh and particularly peak values, at all locations. This is partly due to the narrower frequency range of buoy measurements compared to the frequency range of SWAN, and the wind forcing uncertainty in terms of duration, strength and direction. The highest swh mean values are found at Mykonos (Obs = 1.50 m, Sims = 1.52 m) where the highest wind velocity was recorded. The best correlation between simulated and observed swh values is found at Saronikos, Lesvos, and Crete (0.79, 0.75, and 0.75, respectively). This is expected, as the wind forcing was most accurate at Lesvos, while there is large fetch for waves to develop for north-northwest directions at Crete. The lowest model performance at Skyros (a = 0.81, R = 0.47, and snrmse = 0.340) is related to the wind speed overestimation.
SWAN estimated mwdir closely follows the wind forcing, due to the absence of swell (unimodal seas). Mykonos has the best performance (mae = 10 and nrmse θ = 0.034) but also the least observed variation (20 • ) in the northwest direction. The Saronikos location displays the worst performance (mae = 35 and nrmse θ = 0.151), which is directly related to directional wind forcing. Regarding T m , the closest match is displayed at the open-sea location of Crete (a = 0.94, R = 0.73, and snrmse = 0.115), followed by Mykonos (a = 0.90, R = 0.66, and snrmse = 0.131), and Santorini (a = 0.94, R = 0.64, and snrmse = 0.089). At Saronikos, the wind direction discrepancy favors the development of young waves and the dissipation of the older ones (higher and lower frequencies, respectively), that finally determine the value of T m01 [65].
The comparison demonstrates that the model overestimates T p everywhere except at Saronikos. The highest agreement between observed and modelled T p values is found at the Mykonos location (a = 0.98, R = 0.75, and snrmse = 0.095), where the highest wind speed is observed, and Lesvos (a = 0.95, R = 0.73, and snrmse = 0.102). At Saronikos, SWAN has the poorest performance (a = 1.10, R = 0.22, and snrmse = 0.606), as well as for the mean wave period. Finally, in all buoy stations, wind-sea waves were the dominant wave type.

Island Effects to the Turbulent Fluxes
Atmospheric turbulence calculated from airborne measurements generally displays significant scatter relative to that from other sources (towers), mainly caused by the spatial average over areas with strong heterogeneity (land-sea or SST fronts), which do not allow for longer averaging to reduce random errors [66]. In order to exclude the possible effects of the islands, we used the prevailing wind direction at each data point to find the nearest upwind shoreline, and we kept for analysis the data where the nearest shoreline was within the quadrant (i.e., within ±45 deg) of the prevailing wind direction. Then, the fluxes of momentum-friction velocity at the height of measurements (u h * ) and sensible heat (SHF) were distributed in bins (of 5 km) of distance; their averaged values are presented as a function of upwind shore distance in Figure 4. From Figure 4a, it is apparent that the measurements position has a significant role and the momentum flux at the measurement height is generally increased in the lee side of islands due to wind eddies shed by the islands. At a distance less than 5 km, the mean momentum fluxes and its standard deviation is much greater relative to the measurements taken further offshore. In particular, the highest mean values (0.35 ± 0.28 m/s and 0.32 ± 0.13 m/s) are found at~5 km for positive and negative ∆θ v respectively. For distances from 5 to 20 km, the fluxes gradually decrease as the offshore distance increases, though the reduction rate with distance is greater for the negative (stable), compared to the positive ∆θ v (unstable) samples.
The island effect on the airflow is particularly evident over the Cyclades region, the area with the highest density in number of islands in AS and strong winds (not shown), while over the southeast AS, under unstable conditions, the impact of very few islands is weak (not shown). Under the prevailing Etesian regime, there is significant variability of atmospheric conditions over the AS. Therefore, in some areas, the variation of fluxes with the offshore distance can also be attributed at least partly to the change of atmospheric stability intensity. Substantial differences between stable and unstable conditions for offshore flows are also reported by [26]. Under stable conditions in particular, the momentum flux was found to decrease and to approach equilibrium at~10 km offshore. A fraction of the turbulence decrease with distance could also be explained by the fact that the surface stress over the young seas is higher compared to older ones [5]. However, [25] found no correlation of the momentum flux with the underlying waves in the first few kilometers downwind of land, and they argued that the flux decrease with distance is mainly due to the decreasing influence of the upstream land.
A similar pattern to u h * is shown for SHF, as both fluxes involve the vertical velocity variance. At distances less than 5 km, the mean SHF and its standard deviation are −45 ± 55 Wm −2 and −17 ± 27 Wm −2 for positive and negative ∆θ v respectively (Figure 4b). Furthermore, the positive ∆θ v samples appear to be considerably affected by the island presence, compared to the negative ∆θ v samples; nevertheless, these cases are very few to draw further conclusions.
It is also known, either from simulations or observations, that oceanic wake eddies with a kilometer's size can be induced by atmospheric or oceanic flows downwind in the lee side of islands [67][68][69]. These eddies, if they are cyclonic, lead to locally cold SST (upwelling) and, if anticyclonic, to warm SST (down-welling). The effect of island wakes on SST can be detected in high resolution SST satellite images. Further, wind induced wake eddies from islands with a significant mountain may lead to a warm SST patch in the lee side of the island. Persistent wind stress curl associated with the air flow distortion may also lead to local changes of SST [68][69][70]. In the case of AS, with many islands close to each other, the situation is complicated, but some cold SST eddies have been observed close to the islands in satellite images during the AG (not shown). In such cold SST eddies downwind of islands, the SHF become negative, which is indicated by the observations in Figure 2b.
We conclude that the direct influence on the fluxes and the decay of the advected upstream residual turbulence can be concealed by the conditions met over the sea, the observational height, the size of the island, and the relative position of the measurements. More specific, the data collected within the distance of 10 km are more likely to have been affected by the nearest upwind island and therefore have been excluded from the final dataset, resulting in 668 samples that will be used hereafter.

The Effect of Mean Flow Parameters on Fluxes Variability
In this section, both momentum and sensible heat fluxes are examined in relation to the mean flow variables (i.e., wind speed and ∆θ v ). Unlike the fast coupling of momentum between the atmosphere and the sea surface, the interaction between the air and sea surface temperature is rather slow [71], especially under stable conditions, which can be forced by warmer air advection [72,73]. According to [74], the difference between air temperature and SST could reflect the atmospheric stability regime qualitatively (based on its sign), but not necessarily quantitatively the atmospheric stratification. This is due to the fact that the interconnection between the air and sea is mainly through molecular thermal conduction, which is much slower than the turbulent mixing. Thus, the temperature gradient within the air which includes short-term variations of atmospheric stratification, if available, is preferable to the gradient based on SST. In our data, concurrent measurements of air temperature at two different heights were not available, and, thus, we used the air-sea virtual potential temperature difference ∆θ v for atmospheric stability classification.
For the AG dataset, the measured wind speed decreases as the airflow decelerates, moving southward from colder (north-east and east-central AS) to warmer (southeast AS) waters, and the conditions become more unstable. Figure 5a illustrates the bin median friction velocity u h * as a function of U (in bins of 2 m/s) for two ∆θ v groups. Under strong shear, like the conditions met during the AG (weakly stratified and near-neutral), the wind velocity is expected to be the deciding factor for the momentum flux, yet, u h * (U) depends significantly on stratification (Figure 5a). u h * (U) increases almost linearly up to 15 m/s (which implies an almost constant drag coefficient), while, for higher winds (U > 15 m/s), it decreases for both groups (Figure 5a) and, as expected, friction velocity at the same wind speed is greater for unstable compared to stable conditions (as implied by stability effects in MO similarity). It must be noted that the relatively large measurement height (150 m) does not allow the use of MO similarity, as in COARE, to estimate neutral drag coefficients. It is also worth mentioning that, on average, samples collected under quite stable conditions (∆θ v < −1 • C) show very small dependence on wind speed and u h * retain very low values (<0.15 m/s) (not shown). In contrast, ref. [71] found a strong dependence of the momentum flux with the wind speed under stable conditions; however, their data was taken at~6 m above msl. At very high wind speeds (as in hurricane conditions) it has been observed that the neutral drag coefficient values ceases to increase with wind speed [75,76]. In case of AG, however, the higher winds are mainly related to shallow boundary layers and more stable conditions and, thus, the measurement level (150 m above msl) is well above the surface layer and close to the MABL top. In these cases, a linear with height model of flux divergence (and the corresponding correction) may not hold well and turbulence could decay faster with height. Also, a near-logarithmic wind profile assumed from similarity theory in SL does not hold in the presence of LLJ, and, thus, wind speed at the measurement height is higher than expected from similarity profile.  The SHF as a function of U∆θ nearly passes through the origin (i.e., near zero flux for vanishing ∆θ, which is actually the method used for SST bias calibration). As indicated in Figure 5b, the bin SHF scatter increases as |U∆θ| increases, this appears from the more irregular behavior of the SHF collected under strong winds (>15 m/s), compared to lower winds conditions. Additionally, for U∆θ < −30 Km/s, the heat fluxes become almost zero at the measurement height, and they are related to very low turbulence primarily in the northeastern part of AS [30].

Momentum Fluxes Comparison
In this section, the simulated friction velocities u * from both algorithms, KCM and COARE, are examined. Estimates of the surface momentum fluxes have been performed from the airborne measurements. The frequency distribution of the symmetrically normalized relative residual (rr) between the two algorithms is displayed in Figure 6 (narrow distribution around zero indicates better agreement). The data are grouped with respect to the bulk atmospheric stability. The data consist of 431 samples under near-neutral to unstable conditions and 237 samples under near-neutral to stable conditions. The results presented in Figure 6 indicate that stability effects are significant and have to be taken into account even in the marine environment, where the turbulent momentum flux is generally small and near-neutral conditions mostly exist. In general, KCM estimates higher values (more than 5% and up to 30%) than COARE under stable conditions. Under unstable conditions, the differences are small and around zero value (|rr| < 10%) (Figure 6a). However, the KCM mod simulated fluxes are systematically higher (up to 10%) than COARE under both unstable and stable atmospheric conditions (Figure 6b). This is probably due to the wave-induced momentum fluxes which are calculated more accurately in KCM, but they are implicitly parameterized in COARE. It should be noted, however, that the wave effect appears to be less significant than the stability effect for the specific dataset.
Furthermore, the frequency distribution (%) of the symmetrically normalized relative residuals between the simulated friction velocities u * from the KCM mod , with the observed height-corrected (i.e., reduced to surface level using linear flux divergence) ones u cor * , relative to either ∆θ v or U, presented in Figure 7. For the near neutral to stable regime (∆θ v < 0 • C), the distribution displays two peaks (non-symmetric bimodal distribution) in the range of 60 to 80% and −20 to −40%. As displayed in Figure 7a, the majority of the simulated u * (~65%) are higher than the height-corrected measurements (rr > 0). These data have been collected under more stable conditions and near the MABL top and, thus, weaker turbulence is observed. The rest of the simulated u * (~35%), are lower than the measurements (rr < 0). These data are reported downwind or near islands, and the turbulence observed is greater than expected. For the near-neutral to unstable regime (∆θ v > 0 • C) the model performs better; however, the airborne fluxes are systematically lower (~85% cumulative frequency of occurrence, Figure 7a) than the model estimations (right skewed distribution), yet the residual distribution peak is found within the limits of 0 to 20%. The moderate to high wind group (U < 15 m/s) displays good performance in terms of simulation efficiency, and the simulated u * differ less than 40% (|rr| < 40%) compared to the observations, with a cumulative frequency of occurrence near 50% (Figure 7b). For the group with the stronger winds (U > 15 m/s), simulations significantly overestimate the measurements and the distribution exhibits a displaced residual peak value found in the interval of 60% to 80% [77], demonstrating that COARE performed better under strong winds than under weak winds. In our case, the stronger winds are related to the presence of LLJ and the measured wind speed is higher than predicted from the similarity profile. Figures 8 and 9 display a time series of the simulated surface and the height-corrected observed momentum fluxes, along with the basic meteorological variables, for two selected flight paths (displayed in Figure 2): in the east-central Aegean (west of Chios), with near-neutral to stable conditions and high winds, and in the southeast Aegean (near Karpathos), with unstable conditions (due to higher SST) and moderate winds. As expected, the simulated fluxes closely follow the wind speed variation in both flights, while the measured ones (u cor * ) do not display such a clear dependence on U. Under stable stratification, the simulated fluxes are substantially greater than u cor * (by a factor of~2, Figure 8). The inhibition of vertical mixing and turbulence under stable conditions reduces the height of the MABL top; under these conditions, the wind aloft may be frictionally decoupled from the surface [78], stimulating the formation of a LLJ. Unfortunately, the smaller effect of other factors (like wave effects on the structure of the lower part of the MABL) cannot be distinguished and analyzed with the present data, due to significant flux divergence in such cases. On the contrary, under unstable stratification, the simulated fluxes are closer to the u cor * values that clearly depend on U (Figure 9). The increased turbulent mixing deepens the MABL under unstable conditions, thereby increasing the interconnection between the surface and the wind aloft above the SL and in the mixed layer.   Figure 8a and (b) as in Figure 8b, but for the flight segment at the southeast AS, on 9 September 2011 (as displayed in Figure 2b).

Conclusions
In this study, surface momentum fluxes are estimated using KCM and COARE algorithms and evaluated against airborne measurements over the AS. In addition to the modification of the energy-containing wave spectrum, the effects of stratification were included in the KCM mod model. The KCM spectral input is provided by the SWAN model, which is driven by the WRF model.
During the Etesian outbreak, spatially consistent wind-induced wave conditions are sustained throughout the AS. SWAN displays good agreement with buoy observations at most locations regarding the mean wave parameters. A higher agreement is found for the significant wave height, as it is less sensitive to the spectral distribution, than for the mean or the peak wave periods.
The observed momentum fluxes (at~150 above msl) are found to be significantly affected by the islands (up to 10 km downwind from islands), especially under stable conditions. For near-neutral to weakly stable and unstable conditions, an increasing tendency of the momentum fluxes as a function of wind speed is observed (up to U = 15 m/s). The averaged sensible heat fluxes (upward or downward) are found to be in the near-neutral regime (|SHF| < 10 Wm −2 ).
The two algorithms (KCM and COARE) display better agreement on the unstable compared to the stable regime. The stability effects are of greater significance than the wave effects for the AG dataset (likely due to the measurement height) and COARE estimates are lower than KCM mod (up to 10%), likely due to the unresolved wave-induced fluxes in COARE. Under unstable conditions, KCM mod momentum fluxes are consistently higher than the measured fluxes, but the normalized residual peak value is close to zero. For stable conditions, the majority of the friction velocity estimates are much higher than the measurements (mostly found on northeast and east-central AS). Notably, the aircraft measurements used in this study are above SL, and the linear flux divergence assumption, which was used to reduce flux values to surface level, may not be valid, especially under stable atmospheric conditions. The agreement between KCM mod and the airborne momentum fluxes is better in cases of lower winds (<15 m/s) compared to stronger ones (>15 m/s). However, most of these wind speed values (>15 m/s) are higher than anticipated from the similarity profile due to the frequent occurrence of a LLJ.
In conclusion, it seems that both the COARE and KCM models performed reasonably in the semi-closed basin of AS under complex atmospheric and wave conditions. However, the wind-over-wave-coupling KCM model should be further tested with field measurements closer to sea surface, acquired under conditions where the impact of waves overweighs that of stratification in the estimation of momentum fluxes.  Data Availability Statement: Data will be made available on request, for collaboration purposes.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
For quantitative evaluation of the degree of accuracy of the models' results, the statistical indicators employed in this study are as follows: For scalar quantities: The slope of the regression line (a) through the origin, which indicates the presence of biases (a < 1 model overestimates, a > 1 model underestimates).
The Correlation Coefficient (R), which is a measure of the degree of linear dependence between predicted P i and observed O i values, and an indicator of the scatter component of the error: The symmetrically normalized root mean square error (snrmse) introduced by [79] combines information for the average and scatter components of the error. According to [80], widespread indicators, such as the scatter index and the root mean square error, identify negative bias simulations as better performing relative to unbiased simulations: For directional quantities: The mean absolute error (mae) of the convex angle between observed and predicted directions: The normalized root mean square error (nrmse θ ) of the convex angle between observed and predicted directions: Criterion for the most energetic wave (peak) train identified as wind-sea: 1.2·U 10m · cos(θ d ) > ·c p and θ d < 45 o (A5) c p is the phase speed of the surface waves at the spectral peak. θ d is the angle between the reference wind and peak wave train direction.