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

Abstract

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.

1. Introduction

The atmosphere-ocean interface differs significantly from that of the atmosphere-land 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 flow variables, but it does not consider the structure of the WBL, and there is no feedback accounting for the dynamic wind-wave interactions.

Analytical and numerical based studies have been developed in order to add information regarding the WBL processes and their effect on the momentum and energy flux [10,11,12,13]. Within the WBL, apart from the turbulent fluctuations, there is a component coherent with the underlying waves [14,15]. The wave-induced fluctuations attenuate as the vertical distance from the sea-surface increases, and, in case of stationarity and horizontal homogeneity, the boundary layer far above the waves resembles that over a flat surface [14]. The total surface stress in WBL, defined as the sum of wave-induced and shear (turbulent) stress, is considered most of the times positive (i.e., downward transfer of momentum). However, under low wind speed with strong incoming swell, the wave-induced stress can be negative (i.e., upward transfer of momentum), according to [16]. Under such conditions, the turbulent stress may approach zero at some height, which can lead to the formation of a wave-driven wind jet at the top of WBL [17] and impact the wind and turbulence profile in the MABL [18]. WBL height is estimated to be of O (1 m) [19,20,21] depending on the peak of the wave spectrum (dominant wave) for wind-driven sea waves (wind-seas) but it can extend to much greater heights in swell dominated conditions [22,23]. Furthermore, wave-breaking contributes significantly in the dynamics of the WBL; [24] estimated that a substantial portion of the wave-induced stress is supported by breaking waves (up to 50% depending on the wind) due to air-flow separation at the wave breaking fronts. Apart from this, the wave-breaking controls the energy dissipation of waves as well as the heat transfer through white-capping and sea spray.

The coastal regions are usually dominated by young-developing seas due to the limited fetch. On top of this, the changes of bottom-topography near the shore, influences the approaching waves (depth-limited) in terms of wave propagation, height, steepness and breaking point and consequently the WBL processes. In addition, the significant heterogeneity near the shoreline limits the applicability of the assumptions valid over the open sea. Heterogeneity could impact areas further downstream of coastlines, because well-defined layers in the vertical may not exist due to the formation of an Internal Boundary Layer (IBL). Ref. [25] recognized that the stress dependence on wave age may be partially due to the advection of land turbulence offshore, decaying with the downstream distance, and [26] traced the influence of heterogeneity downstream and found residual turbulence advected from upstream, detached from the sea surface under strong wind conditions. Also, [27] argued that the linear formulation, used to extrapolate flux measurements at higher levels to the corresponding surface values, does not apply to coastal zones due to advection and the IBL formation. For example, during the Risø Air-Sea Experiment (RASEX), an elevated stress maximum was observed and maintained by the horizontal advection of land turbulence [28]. According to [29], no significant dependence on the MO stability parameter was found in a study at a coastal site and the universal functions for the non-dimensional gradient wind or temperature did not comply with the classical expressions.

The present study attempts to estimate and evaluate the air-sea turbulent momentum flux over the Aegean Sea (AS), in the eastern part of the Mediterranean, under Etesian atmospheric conditions. The AS region has a very large number of small and large islands and it is characterized by a complex shoreline and bottom topography. Only a few studies regarding turbulence have been performed in this region, to our knowledge. In particular, apart from the Aegean-GAME (AG) aircraft campaign, where the atmospheric turbulent quantities are described and analyzed by [30] and further investigated by [31], there is also another study with micro-meteorological measurements by [32]. However, in the study of [31] no wave-state or topographical wake effects were considered regarding the variation of the fluxes. On the other hand, the study of [32] was carried out at two sites on land, 30 m from shoreline, and their observations were influenced by local characteristics of the airflow due to the proximity to the shore.

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.

2. 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.

2.1. 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 stereo-photograph 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].

2.2. 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 $\mathsf{\Delta }{\theta }_v={\theta }_v\left(sea\right)-{\theta }_v\left(air\right)$ 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”.

2.3. 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).

2.4. 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 wave simulations configuration.

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 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.

3. Results and Discussion

3.1. Wave Conditions

3.1.1. SWAN Sensitivity

The performance of the wave simulations is examined through the comparison with buoys and in particular the significant wave height $\left(\mathrm{swh}\right)$, the mean wave direction $\left(\mathrm{mwdir}\right)$, and the first moment mean wave period $\left({\mathrm{T}}_{\mathrm{m}01}\right)$. The most energetic (peak) wave period $\left({\mathrm{T}}_{\mathrm{p}}\right)$ 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., $\mathrm{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.

3.1.2. 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 $\mathrm{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 $\mathrm{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.

3.1.3. Wave Simulation Comparison with In-Situ Measurements

A quantitative evaluation of SWAN and WRF models in relation to buoy data is summarized in

Table 2. Mean observed and simulated values for atmospheric and wave variables from 4 to 7 September 2011 at buoy locations.

	Depth at Buoy Location (m)	Distance from Nearest Shore (km)	Wind 3 m above msl (m/s)		Prevailing Wind Direction		$swh$(m)		$T_{m01}$(s)		$T_{p }$(s)		$mwdir$
			Obs	Sims	Obs	Sims	Obs	Sims	Obs	Sims	Obs	Sims	Obs	Sims
Skyros	83	10.2	4.8	6.1	NNE	NE	0.68	0.82	3.3	3.7	4.0	5.1	NNE	NNE
Lesvos	123	3.7	7.2	6.9	N	N	0.74	0.82	3.5	3.5	4.5	4.8	NNW	N
Saronikos	216	9.6	4.5	5.1	E	SE	0.43	0.45	2.9	2.5	3.7	3.2	ENE	NNE
Mykonos	99	5.7	9.5	9.3	NNW	NNW	1.50	1.52	4.2	4.6	5.6	5.8	NNW	N
Santorini	303	9.4	5.9	7.0	NW	NNW	0.81	0.86	3.4	3.6	4.3	5.1	NNW	N
Crete	1073	35.8	5.7	6.3	WNW	NW	0.89	0.94	3.7	3.9	4.7	4.9	NNW	NNW

and

Table 3. SWAN statistical evaluation against buoy measurements from 4 to 7 September 2011.

	$swh$			$T_{m01}$			$T_p$			$mwdir$
	$a$	$R$	$snrmse$	$a$	$R$	$snrmse$	$a$	$R$	$snrmse$	$mae$ (°deg)	$nrms{e}_{\theta }$
Skyros	0.81	0.47	0.340	0.89	0.55	0.146	0.80	0.51	0.267	18	0.058
Lesvos	0.88	0.75	0.224	0.99	0.56	0.095	0.95	0.73	0.102	16	0.050
Saronikos	0.91	0.79	0.254	1.10	0.15	0.241	1.10	0.22	0.606	35	0.151
Mykonos	0.99	0.70	0.172	0.90	0.66	0.131	0.98	0.75	0.095	10	0.034
Santorini	0.93	0.58	0.220	0.94	0.64	0.089	0.85	0.44	0.200	21	0.070
Crete	0.96	0.75	0.252	0.94	0.73	0.115	0.95	0.33	0.236	18	0.058

Statistical indicators: slope of the regression line $\left(\mathrm{a}\right),$ correlation coefficient $\left(\mathrm{R}\right)$, symmetrically normalized root mean square error (snrmse), mean absolute error $\left(\mathrm{mae}\right)$, and directional normalized root mean square error (${\mathrm{nrmse}}_{\mathsf{\theta }}$).

for the period of 4 days. The statistical metrics applied for scalar ($\mathrm{a}$, $\mathrm{R}$, and $\mathrm{snrmse}$) and directional quantities ($\mathrm{mae}$ and ${\mathrm{nrmse}}_{\mathsf{\theta }}$), 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 $\mathrm{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 $\mathrm{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 $\mathrm{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$(\mathrm{a}$ = 0.81, $\mathrm{R}$ = 0.47, and $\mathrm{snrmse}$ = 0.340$)$ is related to the wind speed overestimation.

SWAN estimated $\mathrm{mwdir}$ closely follows the wind forcing, due to the absence of swell (unimodal seas). Mykonos has the best performance ($\mathrm{mae}$ = 10 and ${\mathrm{nrmse}}_{\mathsf{\theta }}$ = 0.034) but also the least observed variation (20$°$) in the northwest direction. The Saronikos location displays the worst performance ($\mathrm{mae}$ = 35 and ${\mathrm{nrmse}}_{\mathsf{\theta }}$ = 0.151), which is directly related to directional wind forcing. Regarding ${\mathrm{T}}_{\mathrm{m}}$, the closest match is displayed at the open-sea location of Crete ($\mathrm{a}$ = 0.94, $\mathrm{R}$ = 0.73, and $\mathrm{snrmse}$ = 0.115), followed by Mykonos ($\mathrm{a}$ = 0.90, $\mathrm{R}$ = 0.66, and $\mathrm{snrmse}$ = 0.131), and Santorini ($\mathrm{a}$ = 0.94, $\mathrm{R}$ = 0.64, and $\mathrm{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 ${\mathrm{T}}_{\mathrm{m}01}$ [65].

The comparison demonstrates that the model overestimates ${\mathrm{T}}_{\mathrm{p}}$everywhere except at Saronikos. The highest agreement between observed and modelled ${\mathrm{T}}_{\mathrm{p}}$values is found at the Mykonos location ($\mathrm{a}$ = 0.98, $\mathrm{R}$ = 0.75, and $\mathrm{snrmse}$ = 0.095), where the highest wind speed is observed, and Lesvos ($\mathrm{a}$ = 0.95, $\mathrm{R}$ = 0.73, and $\mathrm{snrmse}$ = 0.102). At Saronikos, SWAN has the poorest performance ($\mathrm{a}$ = 1.10, $\mathrm{R}$ = 0.22, and$\mathrm{snrmse}$ = 0.606), as well as for the mean wave period. Finally, in all buoy stations, wind-sea waves were the dominant wave type.

3.2. Observed and Simulated Fluxes

3.2.1. 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 $\mathsf{\Delta }{\theta }_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 $\mathsf{\Delta }{\theta }_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⁻² and −17 ± 27 Wm⁻² for positive and negative $\mathsf{\Delta }{\theta }_v$ respectively (Figure 4b). Furthermore, the positive $\mathsf{\Delta }{\theta }_v$ samples appear to be considerably affected by the island presence, compared to the negative $\mathsf{\Delta }{\theta }_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 $\mathrm{SST}$ (upwelling) and, if anticyclonic, to warm $\mathrm{SST}$ (down-welling). The effect of island wakes on $\mathrm{SST}$can be detected in high resolution $\mathrm{SST}$ satellite images. Further, wind induced wake eddies from islands with a significant mountain may lead to a warm $\mathrm{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 $\mathrm{SST}$ [68,69,70]. In the case of AS, with many islands close to each other, the situation is complicated, but some cold $\mathrm{SST}$ eddies have been observed close to the islands in satellite images during the AG (not shown). In such cold $\mathrm{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.

3.2.2. 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 $\mathsf{\Delta }{\theta }_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 $\mathrm{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 $\mathsf{\Delta }{\theta }_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 $\mathsf{\Delta }{\theta }_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\left(U\right)$ depends significantly on stratification (Figure 5a). $u_{*}^h\left(U\right)$ 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 ($\mathsf{\Delta }{\theta }_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.

Figure 5b, illustrates the bin median $SHF$ as a function of $U\mathsf{\Delta }\theta$, relevant to the bulk formulation. The $SHF$ as a function of $U\mathsf{\Delta }\theta$ nearly passes through the origin (i.e., near zero flux for vanishing $\mathsf{\Delta }\theta$, which is actually the method used for SST bias calibration). As indicated in Figure 5b, the bin $SHF$ scatter increases as $\left|U\mathsf{\Delta }\theta \right|$ 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\mathsf{\Delta }\theta$ < −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].

3.2.3. 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 $\left(rr\right)$ 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 $\mathsf{\Delta }{\theta }_v$ or $U$, presented in Figure 7. For the near neutral to stable regime ($\mathsf{\Delta }{\theta }_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 ($\mathsf{\Delta }{\theta }_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.

Figure 8 and Figure 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 $\mathrm{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.

4. 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⁻²).

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.