Hydrographic response to winter shamal winds in the Northwestern Arabian Gulf

1 Introduction

Occupying a shallow, semi-enclosed basin, the Arabian Gulf (hereafter referred to as the Gulf) stretches from the northwest to the southeast along the eastern Arabian Peninsula (Figure 1), connecting to the Gulf of Oman and the Arabian Sea via the Strait of Hormuz at its southern outlet. It spans an area of approximately 250,000 km² and holds a water volume close to 8,630 km³ (Emery, 1956; Elshorbagy et al., 2006; Riegl and Purkis, 2012; Vaughan et al., 2019). Water depths average roughly 35 m across most of the basin but deepen markedly toward the Strait of Hormuz, reaching close to 100 m (Chao et al., 1992; Hereher, 2020; Mushtaq et al., 2025; Subeesh et al., 2025). Riverine input is modest and largely restricted to the northern basin, where the Euphrates, Tigris, and Karun together discharge a mean flow near 703 m³s^–1 (Al-Ansari et al., 2019; Janabi, 2010). Located between 24°N and 30°N within the subtropical high-pressure belt, the Gulf experiences predominantly clear skies and arid conditions, resulting in high net evaporation rates estimated at approximately 1.3 myr^–1 (Sheppard et al., 2010; Al-Senafi et al., 2024). This combination of shallow depth, restricted exchange, and strong evaporative forcing renders the water column highly sensitive to episodic atmospheric events.

The prevailing wind regimes over the Gulf are dominated by northwesterly winds (Al-Senafi and Anis, 2020; Aboobacker et al., 2024; Alrushaid and Al-Senafi, 2026), locally known as Shamal (Arabic for “north”), and southeasterly winds, known as Kaus (Figure 2). Shamal winds can intensify considerably during synoptic events, reaching speeds of up to 23 ms^–1. In the present study, following the criteria of Al-Senafi and Anis (2015) and Alrushaid and Al-Senafi (2026), a Shamal event is defined as winds from the west-northwest to north (287°< direction< 360°) with hourly mean speeds of at least 9.85 ms^–1 sustained for three or more consecutive hours in a day. These events occur primarily during the summer (June–September) and winter (December–March) seasons, accounting for approximately 85% of all recorded Shamal occurrences (Al-Senafi and Anis, 2015; Aboobacker et al., 2024). Summer Shamals are typically driven by a quasi-stationary monsoonal low-pressure system centered over northwestern India and southern Saudi Arabia (Rao et al., 2003), while winter Shamals are generally associated with mid-latitude frontal systems originating over the Eastern Mediterranean and propagating southeastward into the Gulf region (Rao et al., 2001; Thoppil and Hogan, 2010). Winter Shamal events typically persist for 2–5 days, depending on the speed and evolution of the associated frontal system (Ali, 1992; Walters, 1990; Al-Senafi and Anis, 2015).

In shallow marine environments such as the northwestern Gulf, where depths are generally less than 37 m, hydrodynamic processes and turbulent mixing are largely regulated by surface forcing mechanisms, including buoyancy fluxes, heat exchange, tides, and wind-driven momentum transfer (Al-Senafi et al., 2020; Li et al., 2020a). Against this background, episodic atmospheric disturbances, notably Shamal wind events, play a pivotal role in modulating vertical mixing, initiating dense water formation, resuspending sediments, and redistributing nutrients and other biogeochemical constituents throughout the water column. In these relatively shallow settings, the cold, dry, and intense Shamal winds are expected to induce convectively driven mixing, comparable to the deep convective overturning associated with Mistral events in the northwestern Mediterranean (Schott et al., 1996; Vreugdenhil and Gayen, 2021) and Bora events in the Adriatic Sea (Lee et al., 2005). In the Gulf, such convection is likely to penetrate to the seabed, promoting the resuspension of benthic nutrients, which may subsequently enhance phytoplankton productivity. In addition, aeolian dust transport associated with Shamal events has been shown to further modulate phytoplankton dynamics (Nezlin et al., 2010), highlighting an additional pathway through which atmospheric processes influence marine primary production in this region.

Characterizing this response requires distinguishing the two turbulent regions of the water column. Vertical mixing in the surface boundary layer (SBL) is primarily driven by air–sea interactions, including wind stress (e.g., Tsagareli et al., 2010), buoyancy and heat fluxes (e.g., Anis and Moum, 1994), and surface wave activity (e.g., Paskyabi et al., 2012), often resulting in a homogeneous or near homogeneous density layer. The SBL depth, defined by the onset of stratification (Anis and Singhal, 2006), responds to vertical mixing intensity (Sutherland et al., 2014), varying diurnally by tens of meters (e.g., Schneider and Müller, 1990) and seasonally by up to 100 m as shown by Sprintall and Roemmich (1999). This layer supports most biological activity and is central to air–sea gas exchange. Tidal and storm–driven currents further enhance mixing through shear production in both the SBL and deeper layers. The bottom boundary layer (BBL), adjacent to the seabed, typically spans 4–50 m (Trowbridge and Lentz, 1991; Gloor et al., 2000) and exhibits elevated turbulence levels driven by bottom friction. Under constant stress conditions, the velocity profile often follows a logarithmic distribution (Doron et al., 2001). Intense mixing within the BBL promotes the resuspension and transport of benthic materials, including nutrients, sediments, and toxins (Dewey et al., 1987).

Despite its environmental and economic significance, the dynamics of turbulence and vertical mixing in the Gulf remain poorly understood, and direct in situ observations of these processes are scarce (Fieseler et al., 2023). Comprehensive observational and numerical investigations of the Gulf by Emery (1956); Brewer and Dyrssen (1985); Chao et al. (1992); Reynolds (1993); Swift and Bower (2003); Thoppil and Hogan (2010), and Bruciaferri et al. (2022) have advanced understanding of its circulation and water-mass structure, yet have paid limited attention to the role of Shamal-driven mixing. In particular, integrated observations linking atmospheric forcing, hydrographic structure, and bottom boundary layer turbulence during Shamal events, together with their evaluation against numerical turbulence closure models, remain to our knowledge largely absent for the northwestern Gulf.

The present study addresses this gap by integrating observational data and numerical modeling to examine turbulence dynamics and vertical mixing in the northwestern Gulf during extreme Shamal wind events. The observational campaign aimed to characterize the intensity and temporal variability of turbulent mixing and to relate these patterns to external forcing and current velocities within the SBL and BBL. The modeling effort evaluates the performance of a one-dimensional turbulence model in reproducing Reynolds stresses, eddy diffusivities, turbulent kinetic energy, and dissipation rates within the BBL under the combined influence of tidal flows and Shamal winds, offering predictive insights into present and future extreme forcing scenarios.

The remainder of this manuscript is organized as follows. Section 2 describes the field setup, the turbulence parameter estimation, and the setup and evaluation of the one-dimensional turbulence closure models. Section 3 presents and interprets the results for the two winter Shamal events, covering the atmospheric forcing, the hydrographic response, and the turbulence within the bottom boundary layer. Finally, Section 4 synthesizes the principal outcomes and their implications for Shamal-induced mixing in the Gulf.

2 Methods

2.1 Observational details

The field campaign was carried out near Qarooh Island, Kuwait, approximately 40 km northeast of Al Khairan and 120 km south of the Euphrates–Tigris delta (Figure 1), between mid-January and mid-April 2013. In the study area, the coastline and bathymetry are aligned northwest–southeast (315°–135°), and water depths reach a maximum of 37 m. Measurements were obtained from a surface meteorological station and an oceanographic mooring station (Figure 3).

2.1.1 Atmospheric

The meteorological station was installed on a pier at Qarooh Island (28°48.973’N, 48°46.457’E), with sensors mounted approximately 7 m above sea level (ASL). The measured variables comprised wind speed and direction, air temperature and relative humidity, incoming solar and infrared radiation, barometric pressure, precipitation, and sea surface temperature (SST). A summary of the instruments and their specifications is provided in Supplementary Table S1.

Sensor readings were recorded every 2 minutes, averaged over 10-minute intervals, and transmitted to a remote server via Iridium satellite (Data Garrison system). The radiometers were cleaned daily to maintain measurement accuracy. Instrument malfunctions during the study period nonetheless caused the loss of wind speed, wind direction, air temperature, and relative humidity records, with the air temperature and relative humidity data missing from mid-January through early March. To fill these gaps, 10-minute measurements of wind speed and direction (Vaisala WS425 Ultrasonic), air temperature (Vaisala HMP45A), and relative humidity (Vaisala HMP45D) were obtained from the Kuwait Meteorological Office (KMO) Qarooh Station, located on the same pier adjacent to the original station. All replacement data were standardized to a reference height of 10 m ASL.

The Coupled Ocean–Atmosphere Response Experiment (COARE 3.0) bulk flux algorithm (Fairall et al., 1996, 2003; Drennan et al., 2003) was employed to calculate turbulent air–sea interface fluxes of latent heat (JqL), sensible heat (JqS), momentum (wind stress, τ), net shortwave radiation (Jqsw), and net longwave radiation (Jqlw). These fluxes were summed to obtain the net surface heat flux (Jq) given by Equation 1

(1) Jq0=Jqsw+Jqlw+JqL+JqS

Following the convention of Anis and Singhal (2006) and Li et al. (2020b), all heat fluxes are defined as positive upward (out of the ocean). A positive net surface heat flux (Jq0>0) therefore denotes ocean heat loss (surface cooling), whereas a negative value denotes heat gain (surface warming). In the text, the magnitudes of the individual flux components are reported with their sense (gain or loss) stated explicitly; the figures follow the signed convention (positive upward). The surface buoyancy flux (Jb0) was then computed from Jq0 using Equation 2:

(2) Jb0=gαTJq0ρoCp

Here, g is the gravitational acceleration, α_T represents the thermal expansion coefficient of seawater, ρ_o denotes the near-surface seawater density, and C_p is the specific heat capacity of seawater computed as described in Fofonoff and Jr (1983). A positive buoyancy flux (Jb0>0) indicates net cooling at the ocean surface, promoting vertical mixing and resulting in a deeper mixed layer. In contrast, a negative buoyancy flux (Jb0<0) signifies net surface heating, enhancing water column stability and thereby suppressing vertical mixing (Dorrestein, 1979).

To quantify the relative influence of wind-induced shear and buoyancy-driven mixing within the SBL, the Monin–Obukhov length (L) was calculated using Equation 3:

(3) L=−u3κJb0, with u=(τρo)1/2

where u∗ is the friction velocity at the water surface, and κ = 0.41 is the von Kármán constant. Under convective conditions (Jb0>0), L takes negative values, and its magnitude provides a vertical scale at which the contributions of buoyancy and shear production of turbulence are approximately balanced (Thorpe, 2005).

2.1.2 Hydrography

Hydrographic measurements were obtained from a mooring deployed approximately 2 nautical miles (∼3.7 km) northeast of Qarooh Island, Kuwait (28° 50.938’N, 048° 47.534’E; Figure 1). The mooring configuration consisted of two primary sections, a vertical mooring array and a bottom-mounted pod, separated vertically by 21 m (Figure 3). Instrumentation on the vertical mooring and bottom pod, including sensor types and sampling intervals, is detailed in Supplementary Table S2. All hydrographic measurements were block-averaged over approximately 10–12 minute intervals and subsequently interpolated onto a common 10 minute time grid to provide temporally homogeneous observational datasets for the analysis and figures.

Instrument malfunction in the three fluorometers restricted fluorescence data collection to the period between Jan 18 and 30. Conductivity measurements from the six CTDs suffered from biofouling despite anti-fouling measures such as copper films and meshes. Consequently, constant salinity profiles were assumed throughout this study, following initial measurements prior to biofouling and observations from Reynolds (1993), who indicated stable salinity profiles off the Kuwaiti coast during Feb and Mar.

Current velocities from the high-resolution ADCP (ADCP-HR) were initially measured in beam coordinates (instrument relative u, v, w) and subsequently converted to the Earth coordinate system (North-South, East-West, Up-Down) using the manufacturer’s provided tilt, roll, and heading corrections (Lohrmann et al., 1990). Acoustic backscatter (echo intensity) recorded by the ADCP is reported as a qualitative proxy for suspended particulate matter; no independent calibration (e.g., water-sample or optical turbidity measurements) was available to convert backscatter to an absolute concentration.

Temperature profiles from 25 sensors (one surface-float sensor, 23 sensors distributed along the mooring, and one bottom-pod sensor) were used to estimate the mixed-layer depth. This followed the approach of Thomson and Fine (2003), employing the split-and-merge segmentation algorithm originally developed by Pavlidis and Horowitz (1974). The algorithm identifies segment boundaries in the temperature profiles through piecewise polynomial fitting and defines the mixed layer as the uppermost segment. Because this estimate was derived solely from temperature profiles, it represents a temperature-defined mixed-layer depth rather than an independently determined density mixed-layer depth. Temperature-based estimates have been widely applied to observational profiles, including in the saline, semi-enclosed Mediterranean Sea and Red Sea (de Boyer Montégut et al., 2004; D’Ortenzio et al., 2005; Abdulla et al., 2018). Initial uncontaminated measurements from the CTDs indicated that temperature was the principal contributor to density stratification at the study site.

2.2 Estimating turbulence parameters

The stability of a stratified shear flow is typically quantified using the Richardson number (Ri), defined as the ratio of the buoyancy frequency squared to the shear frequency squared:

(4) Ri=N2S2

In Equation 4 the term S2=(∂u∂z)2+(∂v∂z)2 represents the squared vertical shear of the horizontal velocity, and N2=−gρ0∂ρ0∂z is the squared buoyancy frequency. The classical threshold value Ric = 0.25 marks the critical limit below which turbulence is likely to be sustained (Miles, 1961). However, other observational studies suggest that Ric can range between 0.2 and 1 depending on local conditions (Galperin et al., 2007; Zaron and Moum, 2009). In general, flows with Ri< Ric are prone to turbulence, while those with Ri > Ric tend toward stability.

The ability to quantify turbulence is essential for identifying the dominant physical mechanisms, such as wind stress, surface buoyancy fluxes, waves, or current shear that drive mixing, and the vertical flux of constituents such as nutrients and sediments. Four key turbulence parameters are estimated in this study: the turbulent kinetic energy (TKE), the TKE dissipation rate (ϵ), BBL Reynolds stress (τ_b), and the vertical eddy diffusivity (K_Z). These parameters are later used to assess and compare the performance of the numerical model simulations (see Section 2.3).

Time series of vertical velocity fluctuations (w′), measured by a bottom-mounted high-resolution ADCP (ADCP-HR), were used to compute the power spectral density (PSD), from which both ϵ and TKE were derived. The vertical velocity component is preferred for estimating low ϵ values due to its lower noise floor (1.1×10–⁶ m² s^–2 rad^–1 m) compared to that of the horizontal component (2.6×10–⁶ m² s^–2 rad^–1 m) and its reduced sensitivity to contamination from surface waves (Stapleton and Huntley, 1995; Gordon et al., 1999). The most common approach for estimating ϵ involves fitting the PSD in the inertial subrange to Kolmogorov’s –5/3 power law (Kolmogorov, 1941; Davidson, 2015), as expressed in Equation 5:

(5) Eww(K)=CKolϵ2/3K−5/3

Here, Eww(K) is the wavenumber spectrum of the vertical velocity fluctuations, and C_Kol = 1.5 is the Kolmogorov constant for the transverse spectrum (Pope, 2000).

In shallow environments, orbital velocities associated with surface waves may contaminate the velocity spectrum and lead to overestimation of the turbulence dissipation rate (Lumley and Terray, 1983). To reduce this influence, the dissipation rate was estimated from the frequency spectrum of the vertical velocity fluctuations using the modified inertial-subrange formulation of Trowbridge and Elgar (2001):

(6) Sww(f)=1255CKolϵ2/3(V2π)2/3f−5/3,

In Equation 6 S_ww is where S_ww(f) is the frequency spectrum of the vertical velocity fluctuations, f is frequency, and V=(u¯2+v¯2)1/2 is the burst-mean horizontal current speed.

The upward-looking ADCP-HR sampled 15 cells with a vertical resolution of 75 mm. Each burst contained 512 samples collected at 8 Hz, corresponding to a duration of 64 s, and bursts were recorded every 6 min. Following the quality-control procedure of Rusello (2009) and Al-Senafi and Anis (2020), bursts failing the 70% beam-correlation criterion were excluded before calculating the power spectral density of the vertical velocity fluctuations. During the five-day SE1 period, 1,128 of the 1,200 recorded bursts (94%) satisfied this criterion, with 72 bursts (6%) rejected. During the six-day SE3 period, 115 of the 1,440 recorded bursts (8%) were rejected, leaving 1,325 bursts (92%). Overall, 187 of the 2,640 recorded bursts (7%) were rejected, leaving 2,453 bursts (93%) for analysis. Because spectra were calculated for 15 cells in each retained burst, the analysis comprised 36,795 cell–burst spectra.

The inertial-subrange fitting interval was selected separately for each burst and cell to exclude the surface-wave contribution at lower frequencies and the instrumental-noise floor at higher frequencies. Across the accepted spectra, the lower fitting limits ranged from 0.0313 to 1.6718 Hz, whereas the upper fitting limits ranged from 0.0625 to 2.5626 Hz. The spectral parameters were estimated using robust regression (Huber, 1964; Huber and Ronchetti, 2009).

Frequency was converted to streamwise angular wavenumber using Taylor’s frozen-turbulence hypothesis (Equation 7),

(7) K=2πfV,

where K is expressed in radm^–1. Because V varied among bursts, the corresponding wavenumber limits also varied.

TKE was obtained by integrating the vertical-velocity wavenumber spectrum over the identified inertial subrange after subtracting the instrumental-noise contribution, as given in Equation 8:

(8) k=32[∫KlKhEww(K)dK−∫KlKhNnoise(K)dK],

where K_l and K_h are the lower and upper limits of the inertial subrange, respectively, and N_noise(K) is the instrumental-noise spectrum. The integral was evaluated over the inertial subrange only, so the resulting quantity represents the turbulent kinetic energy associated with the resolved inertial-subrange scales rather than the full-spectrum total; it therefore constitutes a lower bound on the total TKE. The factor 3/2 follows from the assumption of local isotropy within the inertial subrange (Kolmogorov, 1941), under which k=32w′2¯; consistent with this assumption, the transverse form of the Kolmogorov spectrum (C_Kol = 1.5) was used. Local isotropy is expected at the inertial-subrange scales over which the spectrum was integrated, although the larger energy-containing scales near the seabed may be anisotropic.

The third parameter is the BBL Reynolds stress (τ_b), which is a critical metric for sediment resuspension and transport. Once the bottom shear stress exceeds a critical threshold, sediments and nutrients on the seabed can be entrained into the water column (Kim et al., 2000). Following Lohrmann et al. (1990), the Reynolds stress is estimated using Equation 9:

(9) τb=ρbu′w′¯

where ρ_b is the near-bed density.

The fourth parameter, the vertical eddy diffusivity (K_z), plays a vital role in the vertical transport of heat, nutrients, and other properties. Such transport affects primary productivity and can trigger harmful algal blooms or fish kills (MacIntyre, 1993; Glibert et al., 2002). Following Osborn (1980), the vertical eddy diffusivity is computed using Equation 10:

(10) Kz=ΓϵN2

The mixing efficiency was set to Γ = 0.2 following Oakey (1985). The Osborn relation was applied only under statically stable conditions (N² > 0); K_z was not calculated when N² ≤ 0 because the formulation is not applicable under convectively unstable conditions. No additional lower positive threshold for N² was imposed. Consequently, relatively large values of K_z may occur when elevated TKE dissipation coincides with a small but positive N².

Sensitivity to the assumed mixing efficiency was evaluated using Γ = 0.05, 0.10, and 0.20, consistent with the range reported in previous studies (Oakey, 1985; Yamazaki and Osborn, 1993). Because K_z varies linearly with Γ, estimates obtained using Γ = 0.05 and 0.10 were, respectively, 25% and 50% of those obtained using Γ = 0.20. The assumed mixing efficiency therefore affects the absolute magnitude of K_z but not its normalized temporal variability or the timing of enhanced mixing.

Finally, the vertical heat flux in the BBL can be estimated using the computed Kz (Equation 11):

(11) Jqbbl=ρbCpKz∂T∂z

This estimate (Jqbbl) is used to assess how far Shamal-induced surface cooling penetrates into the BBL. Following the sign convention adopted above, Jqbbl is positive upward, so a negative value denotes a downward heat flux carrying surface cooling toward the seabed.

2.3 Numerical model setup and evaluation

Capturing turbulence in aquatic environments remains a substantial challenge due to the episodic nature of turbulent motions and the necessity for high-resolution observational techniques. When constrained by appropriate observational datasets, numerical models offer a practical framework for simulating and evaluating turbulent processes. In the present analysis, a one-dimensional (1-D) implementation of the General Ocean Turbulence Model (GOTM) (Burchard and Bolding, 2001; Burchard, 2002) was employed to replicate vertical turbulent mixing and assess model outputs against field measurements. The model was driven by site-specific hydrographic data (temperature, salinity, velocity profiles) and by surface momentum and heat fluxes derived from meteorological observations recorded at Qarooh Island.

The simulation focused on evaluating two well-established second-moment turbulence closure schemes embedded in GOTM: the k–ϵ model (Launder and Spalding, 1972) and the k–kl model (Mellor and Yamada, 1982). These schemes solve prognostic equations for turbulent kinetic energy (k), coupled with either its dissipation rate (ϵ) or the turbulence length scale (l), to represent vertical mixing processes. Because the one-dimensional model cannot resolve horizontal advection and other three-dimensional processes, temperature and horizontal velocity were relaxed toward the observed profiles (Burchard and Petersen, 1999; Burchard et al., 2002, 2008; Anis and Singhal, 2006). The observed profiles were interpolated onto the model grid and updated hourly over the observed portion of the water column. Temperature was relaxed over approximately 4.0–26.5 m (30 observation levels) and the horizontal velocity components over approximately 5.3–25.8 m (43 ADCP levels); both were interpolated onto the 270 uniformly spaced model layers (0.1 m vertical resolution). Following Burchard et al. (2002), the relaxation (nudging) tendency was specified as:

(12) ∂X∂t|R=Xobs−XmodτR

where in Equation 12, X represents temperature or either horizontal velocity component, X_obs and X_mod are the observed and modeled values, respectively, and τ_R = 1 hour is the relaxation timescale. Temperature and velocity were therefore constrained model variables. Turbulent kinetic energy, its dissipation rate, Reynolds stress, and vertical eddy diffusivity were not prescribed or assimilated and were used to evaluate the turbulence closure schemes.

The two-equation turbulence models are derived from the Navier–Stokes equations using the Boussinesq approximation (Boussinesq, 1903), and they incorporate the Reynolds (1895) decomposition to isolate mean and fluctuating components of the flow (e.g., Kundu and Cohen, 2010). The computation of k is described by Equation 13:

(13) ∂tk=Dk+P+B−ϵ

where D_krepresents the combined effects of viscous and turbulent transport, P shear production, and B buoyancy production. The k–ϵ model computes a separate equation for ϵ, using a similar structure to the one used for k (Equation 14):

(14) ∂tϵ=Dϵ+(ϵk)(cϵ1P+cϵ3B−cϵ2ϵ)

with c_ϵ₁ = 1.44 and c_ϵ₂ = 1.92 being empirical coefficients, while c_ϵ₃ is computed dynamically depending on stability conditions. In contrast, the k–kl model relates ε to a prognostic equation for l (Equation 15):

(15) l=(cμ0)3k3/2ϵ

which is related to ε by Equation 16:

(16) ∂t(kl)=l2(cl1P+cl3B−(1+cl2(l/lz)2)ϵ)

Here, cμ0=0.55, c_l₁ = 1.8, and c_l₂ = 1.33 are model coefficients, c_l₃ is computed dynamically, and l_z is a wall-proximity length scale (Burchard, 1999). Both closure schemes yield estimates of the eddy viscosity (KM) and eddy diffusivity (KH), which are essential for closing the second-moment turbulence equations (Burchard et al., 2008):

(17) KM=cμk2ϵ

(18) KH=cμ′k2ϵ

The coefficients c_μ in Equation 17 and cμ′ Equations 18 are stability-dependent functions derived from the non-equilibrium turbulence formulation of Canuto et al. (2001). These functions extend earlier parameterizations by Mellor and Yamada (1974) and Kantha and Clayson (1994), incorporating a revised treatment of pressure–velocity and pressure–temperature correlations, as well as increasing the critical Ri from 0.23 to 0.85. These stability functions depend on both the shear and buoyancy numbers, as well as the ratio of k² to ϵ², as described by Canuto et al. (2001), and are defined by Equations 19, 20:

(19) cμ=0.107+0.01741αN−0.00012αS1+0.26αN+0.029αS+0.0087αN2+0.005αNαS−0.000034αS2

(20) c'μ=0.1120+0.004519αN−0.00088αS1+0.26αN+0.029αS+0.0087αN2+0.005αNαS−0.000034αS2

where the non-dimensional parameters αS and αN represent the shear and buoyancy influences on turbulence, respectively (Equations 21, 22):

(21) αS=k2ϵ2S2

(22) αN=k2ϵ2N2

The numerical configuration, forcing, initial and boundary conditions, relaxation settings, and turbulence- closure parameters used in the simulations are summarized in Supplementary Table S3.

To assess the skill of the turbulence models in reproducing observed parameters (k, ϵ, τ_b, and K_z), four statistical metrics were employed. The mean bias error (MBE) and root mean square error (RMSE) quantify the average and squared differences, respectively, between modeled (x) and observed (y) values (Equations 23, 24):

(23) MBE=1n∑i=1n(xi−yi)

(24) RMSE=1n∑i=1n(xi−yi)2

The Nash–Sutcliffe model efficiency coefficient (Ns; Nash and Sutcliffe, 1970) was used to measure the predictive skill of the model relative to the variance in the observed data (Equation 25):

(25) Ns=1−∑i=1n(yi−xi)2∑i=1n(yi−y¯)2

An Ns value of 1 indicates a perfect match, while values less than 0 indicate poor model performance. Generally, positive Ns values are considered acceptable (Gupta and Kling, 2011). Finally, the Pearson correlation coefficient (r, Equation 26) assessed the strength of linear association between modeled and observed values:

(26) r=∑i=1n(xi−x¯)(yi−y¯)(n−1)σxσy

Here, σ_x and σ_y are the standard deviations of the modeled and observed variables, respectively (Press et al., 2007). Statistical comparisons were based on paired model–observation values at the times of the accepted ADCP-HR bursts. The SE1 comparisons used n = 1,128 paired samples, whereas the SE3 comparisons used n = 1,325 paired samples.

Because consecutive observations were serially correlated, the 95% confidence intervals for the mean, MBE, RMSE, NSE, and r were calculated using a moving-block bootstrap with 1,000 realizations. Each block contained 124 consecutive observations, corresponding to approximately 12.4 hour at the 6-min sampling interval. This block length preserved the temporal covariance associated with the dominant semidiurnal tide. The lower and upper confidence limits were defined by the 2.5th and 97.5th percentiles of the bootstrap distributions, respectively.

Correlation significance was evaluated using a two-tailed Student’s t-test based on an effective sample size that accounts for serial correlation in the paired observed and simulated time series. The effective sample size was calculated using Equation 27:

(27) neff=n1+2∑j=1n−1(1−jn)ρx(j)ρy(j),

where ρx(j) and ρy(j) are the lag-j autocorrelation coefficients of the simulated and observed time series, respectively. The test statistic was then calculated using Equation 28:

(28) t=r(neff−21−r2)1/2,

with n_eff – 2 degrees of freedom. Values of p< 0.05 were considered statistically significant. These statistical metrics have been widely used in model-evaluation studies in hydrology and atmospheric sciences (e.g., Krause et al., 2005; Elsanabary and Gan, 2014; Jia et al., 2013; Pinker et al., 2014).

3 Results and discussion

Over the course of the 94-day deployment, four distinct Shamal wind events were identified, each lasting approximately 2 to 3 days. These Shamal events (SE) occurred on Feb 3–5 (SE1), Feb 12–14 (SE2), Mar 6–8 (SE3), and Mar 24–27 (SE4). Owing to sensor degradation from biofouling during the latter part of the study, SE4 was excluded from the analysis. Of the remaining events, SE1 and SE3 exhibited comparable duration (65 and 61 hours, respectively) and wind characteristics (average wind speeds of 8.64 and 9.11 ms^–1, and peak wind speeds of 13.97 and 15.12 ms^–1, respectively), with peak surface heat fluxes (Jq0) of 420 and 511 Wm^–2, respectively. As the longest lasting events with the most intense winds, and therefore the strongest expected influence on the hydrographic structure, SE1 and SE3 were selected for detailed analysis. All time references in the analysis correspond to UTC, as used in the figure panels.

3.1 First Shamal event (SE1, 3–5 February)

3.1.1 Atmospheric forcing and surface heat fluxes

The first Shamal event (SE1) occurred from 02:10 on Feb 3 to 19:20 on Feb 5 (Figure 4). Wind conditions prior to and following the event remained light and variable, with magnitudes below 7.5 ms^–1. In contrast, two bursts of intensified northwesterly winds were recorded during SE1. The first burst, spanning the evening (17:30) of Feb 3 through the morning (09:40) of Feb 4, reached maximum wind speeds of 14 ms^–1 and a peak surface stress of 0.58 Nm^–2. The second, more moderate peak occurred overnight between Feb 4 and Feb 5, with wind speeds reaching 9.95 ms^–1 and surface stress estimated at 0.21 Nm^–2 (Figures 4A, B).

Both air and SST displayed a pronounced diurnal cycle throughout the event. The most notable cooling occurred during the first wind burst, when air temperature decreased by approximately 4 °C and SST dropped by 1 °C (Figure 4C). Concurrently, relative humidity reduced by 35%, reflecting the influence of continental dry air advected by the northwesterly flow (Figure 4D). Atmospheric pressure rose by 11.4 mb as the frontal system advanced over the region (Figure 4E).

The net shortwave heat gain was reduced to a daytime magnitude of 642 Wm^–2 on Feb 3, compared with 707 Wm^–2 under clearer skies on Feb 6 (Figure 4F). This reduction was attributed to suspended dust, a common feature of Shamals originating over the Mesopotamian plain (Al-Senafi and Anis, 2015). Net longwave radiation losses peaked at 103 Wm^–2 during the clearest, driest period of the event, whereas the lowest loss of 26 Wm^–2 occurred under overcast conditions prior to the event onset (Figure 4G).

The sensible heat flux also exhibited increased losses during the cooler air periods, triggered by the Shamal winds, reaching 40 Wm^–2 (Figure 4H). The latent heat flux dominated the surface heat loss, driven by increased evaporation under dry and windy conditions, with peak losses of 334 Wm^–2 (Figure 4I). Collectively, these fluxes produced a maximum net surface heat loss of 420 Wm^–2 during SE1, followed by a recovery to a daytime net heat gain of 603 Wm^–2 in the post-event period (Figure 4J).

These flux characteristics place the present observations within the established range of cold, dry wind forcing in semi-enclosed and marginal seas. The dominance of latent heat in the surface loss during SE1 agrees with direct observations at the same site during the March 2013 Shamal, in which latent heat loss increased sixfold and dominated the net flux (Li et al., 2020b), and with the reanalysis assessment of Al Senafi et al. (2019), who attributed approximately 81% of the winter-Shamal net heat loss over the northern Gulf to latent heat. The maximum net surface heat loss of 420 Wm^–2 observed during SE1 is comparable to the 442 Wm^–2 reported by Li et al. (2020b) for a moderate March event and to the extreme turbulent heat losses exceeding ∼400 Wm^–2 documented for cold, dry northerly outbreaks over the northern Red Sea (Papadopoulos et al., 2013), another semi-enclosed, hypersaline basin, though it remains below the basin-scale Shamal extremes simulated by Thoppil and Hogan (2010) (net loss in excess of 1000 Wm^–2, with SST cooling up to 10 °C).

3.1.2 Hydrographic response of the water column

Light intensity within the water column was primarily driven by incoming solar radiation and decayed exponentially with depth. Peak values of 101 Wm^–2 at 2.2 m and 11 Wm^–2 at 22.0 m were recorded on the afternoon of the post-Shamal day (Figure 5B).

The thermal structure displayed three distinct patterns over the observation period. First, diurnal heating consistent with surface net heat flux (Figure 5A) produced warmer temperatures during daylight hours and cooler temperatures at night (Figure 5C). Second, the water column cooled gradually during the Shamal event and warmed afterward: daily average temperatures across the water column were 16.37 °C pre-Shamal (Feb 2), 16.35 °C during the Shamal (mid-Feb 3 to mid-Feb 5), and 16.67 °C post-Shamal (Feb 6). Third, horizontal advection appeared to influence the temperature structure, as seen in the alongshore and cross-shore currents (Figures 5E, F). A scatter plot of temperature versus horizontal currents at 22 m depth (Figure 6A1) – chosen to minimize surface forcing – showed an average temperature difference of 0.24 °C between currents from the western (northwesterly 16.23 °C, southeasterly 16.08 °C) and eastern quadrants (northeasterly 16.49 °C, southeasterly 16.29 °C), with the largest contrast (0.41 °C) between southeasterly and northwesterly flows. These patterns are consistent with advection of warmer water from the southeast during the flooding tide (northwesterly flow) and of cooler water from the northeast during the ebbing tide (southeasterly flow). Because these inferences are based on a single mooring, which cannot directly resolve horizontal property gradients, they indicate a likely advective contribution rather than a directly measured one.

Dissolved oxygen (DO) also exhibited two distinct patterns. The first was its correlation with the horizontal currents (Figure 6A2): eastern flows (northeast and southeast) were associated with higher DO (6.78 mg L¹), whereas western flows (northwest and southwest) coincided with lower DO (6.74 mg L^–1). This is consistent with the advection of comparatively DO-rich water from the southeast during the flooding tide. Second, DO increased gradually through the Shamal period, from 6.95 mgL^–1 pre-event to 7.02 mg L^–1 post-event (Figure 5D). This gradual increase is unlikely to reflect lateral advection alone. As the water column cooled during the event, the solubility of oxygen increased, raising the saturation concentration and favoring higher DO. The strong Shamal winds would also have enhanced air–sea gas exchange (Wanninkhof, 2014), promoting ventilation of the surface layer toward saturation. The observed DO rise is therefore consistent with the combined effects of increased oxygen solubility under cooling, enhanced wind-driven air–sea exchange, and vertical mixing, in addition to lateral advection.

The current structure followed a semidiurnal tidal pattern (Figures 5E, F), with flooding directed to the northwest and ebbing to the southeast. Current velocities ranged from –0.19 to 0.28 ms^–1 cross-shore and from –0.51 to 0.53 ms^–1 alongshore. Acoustic backscatter from the ADCP, used here as a qualitative proxy for suspended particulate matter (SPM) in the absence of an independent calibration, was highest at depths of 20–25 m, particularly during the first 4 hours of each Shamal peak, with echo intensities reaching up to 69 dB (Figure 5G).

Buoyancy frequency (N²) was highest (> 3.2 × 10–⁵ s^–2) near the mixed layer base, indicating strong stratification (Figure 5H). Stratification intensified before the Shamal as a warm layer developed between 5 and 17 m, increasing N² to 1.7 × 10–³ s^–2. During the first Shamal peak, wind-driven vertical mixing (evidenced by a deepening Monin–Obukhov length scale; Figure 4F) weakened the stratification. During the second peak, convective forcing (near-zero Monin–Obukhov scale) weakened it further.

Shear frequency (S²) was elevated (> 1.6 × 10–³ s^–2) near the mixed layer and seabed (Figures 5I, 7C), especially at 23–26 m depth, where shear followed the alongshore velocity trends (Figures 5E, F). At these depths, S² reached 0.1 s^–2, producing instability as the Ri fell below 0.25 (Figure 5J). Shear was relatively low during the Shamal peak periods (< 3.2 × 10–³ s^–2), likely due to vertically homogeneous velocities.

The hydrographic response during SE1 reflects the tidally dominated dynamics of this shallow, inverse-estuarine basin (Kämpf and Sadrinasab, 2006). Harmonic analysis of a three-month record at the same mooring showed that tides account for approximately 97% of the total current energy and 86% of the surface-elevation energy (Li et al., 2020b,a), consistent with the strong semidiurnal tidal circulation of the Gulf (Pous et al., 2012, 2015) and with the semidiurnal current structure and tidally modulated shear observed here. The weakening of stratification during the Shamal peaks, indicated by the reduction in buoyancy frequency, likewise mirrors the near-complete thermal homogenization of the water column reported by Li et al. (2020b) for the March 2013 Shamal at this site.

3.1.3 Bottom boundary layer turbulence: observations versus model

Turbulence characteristics within the BBL (approximately 23–27 m), including TKE, its dissipation rate, Reynolds stress, and vertical eddy diffusivity (Figures 7G–J), were primarily influenced by shear-induced instability associated with the semidiurnal tidal currents. These instabilities became pronounced when the Ri dropped below 0.25 (Figure 7F), often coinciding with elevated shear frequency (Figure 7C). The highest estimated TKE value recorded during the study period, 2.5 × 10–³ m² s^–2, occurred on Feb 2 at 08:37 UTC during a strong pre-Shamal current event (0.25 ms^–1; Figure 7G). During the two Shamal wind peaks, the average TKE decreased to 5.65 × 10–⁴ m² s^–2, below the five-day mean of 6.92 × 10–⁴ m² s^–2, and correlated more weakly with the tidal forcing.

Both the k–ϵ and k–kl models represented the reduction in TKE during the peak wind phases, with mean simulated values of 5.35 × 10–⁴ m² s^–2. Across the five-day period, the simulations were significantly correlated with the observed TKE (r = 0.68 and 0.66, respectively; p< 0.05). However, the corresponding NSE values were only 0.37 and 0.08, indicating limited ability to reproduce the point-by-point magnitude of the observed TKE variability. The MBE was ≤ 3.8×10–⁴ m² s^–2, and the RMSE was ≤ 6.06×10–⁴ m² s^–2 (Table 1). Thus, the closures underestimated the peak dissipation rate and Reynolds stress by factors of approximately five and three, respectively, demonstrating that they did not reproduce the magnitude of this short-duration turbulence event.

Table 1. Statistical comparison of observed and simulated turbulence parameters during SE1, including minimum, maximum, mean values, and corresponding 95% confidence intervals (in parentheses).

The TKE dissipation rate and the Reynolds stress (Figures 7H, I) followed similar trends, with period means of 1.35 × 10–⁵ m² s^–3 and 0.27 Pa, respectively. A marked increase in mixing occurred on Feb 4 at 02:07 UTC, driven by a cold, dense water intrusion associated with a strong Shamal wind pulse. This event, indicated by a vertical heat flux of –14.73 Wm^–2 in the BBL (Figure 7B), reduced the shear (Figure 7C) but elevated the dissipation and stress to 6.07 × 10–⁵ m² s^–3 and 0.82 Pa, respectively. Both models substantially underestimated this event, yielding dissipation rates of about 1.2 × 10–⁵ m² s^–3 and Reynolds stresses near 0.28–0.29 Pa.

The second Shamal peak also triggered convective mixing, reducing stratification and shear in the BBL, although the mixing was weaker than during the first episode, with a dissipation rate of 7.43 × 10–⁶ m² s^–3, below the five-day average. Despite underestimating the short-duration peaks, both closures represented the broader temporal variability associated with tidal shear and Shamal-induced cooling. Simulated and observed values remained significantly correlated (p< 0.05), with correlation coefficients up to 0.63, positive Nash–Sutcliffe efficiencies (NSE ≤ 0.33), and low error statistics (Table 1).

Vertical eddy diffusivity (Figure 7J) showed similar temporal behavior. The five-day mean observed diffusivity was 0.06 m² s^–1, closely matching the mean from the k–kl model. Observed and simulated diffusivities ranged from 7 × 10–⁴ to 0.48 m² s^–1, peaking during the first Shamal wind burst and reaching a minimum on Feb 3 at 12:55 UTC. For K_z, the model–observation comparisons yielded r = 0.79 and NSE = 0.60, with MBE< 5.8 × 10–³ m² s^–1 and RMSE< 5.11 × 10–² m² s^–1. The k–ϵ and k–kl simulations produced very similar estimates (r = 1.00 and NSE > 0.98). Collectively, these results show that both closures represented the general tidal modulation of the turbulence quantities but had limited ability to reproduce the magnitude of some intermittent turbulence events.

The magnitude and behaviour of the BBL turbulence during SE1 are consistent with observations from other energetic tidal environments. The dissipation rates (mean 1.35 × 10–⁵ m² s^–3, peaking at 6.07 × 10–⁵ m² s^–3) and near-bed Reynolds stresses (up to 0.82 Pa) are comparable to those reported for the tidally energetic bottom boundary layer of Red Wharf Bay, where near-bed dissipation peaks at ∼ 10–² W m^–3 (∼ 10–⁵ W kg^–1) and the bed stress reaches ∼ 1.1 Pa (Rippeth et al., 2003; Simpson et al., 1996). The strong tidal modulation of the turbulence, and its modulation during the Shamal, echo the response documented at this site by Li et al. (2020b) for the March 2013 Shamal. Moreover, the tendency of the closure schemes to reproduce the tidal modulation while underestimating the short-lived turbulence peaks is consistent with the benchmark GOTM evaluation of Burchard et al. (2002), who compared the same two-equation closures and Canuto stability functions against bottom-boundary-layer microstructure in the northern North Sea, using an equivalent relaxation of temperature toward observations, and likewise reported fair overall agreement but an under-representation of intermittent high-turbulence events. This comparison therefore places the observed BBL turbulence within the established range for energetic tidal shelf seas and indicates that the closure performance found here is typical of one-dimensional turbulence models rather than specific to this site.

3.2 Second Shamal event (SE3, 6–8 March)

3.2.1 Atmospheric forcing and surface heat fluxes

The third Shamal event (SE3) was observed from 01:30 UTC on Mar 6 to 15:10 UTC on Mar 8 (Figure 8). Preceding the event, southeasterly Kaus winds reached speeds of up to 10.14 ms^–1 (Figures 8A, B), similar in intensity to those recorded before SE1. As in the previous event, SE3 featured two distinct wind peaks. The first occurred between 04:40 and 09:10 UTC on Mar 6, reaching a maximum hourly wind speed of 11.38 ms^–1 and a surface stress of 0.21 N m^–2. The second, more intense peak, originating from the northwest, spanned Mar 6 21:40 to Mar 7 23:40, with winds peaking at 15.12 ms^–1 and surface stress reaching 0.67 Nm^–2 (Figures 8A, B). With an overall average wind speed of 9.11 ms^–1, SE3 was the strongest of the three events in terms of wind speed.

Air temperature followed a diurnal cycle, ranging from 14.3 °C to 22.4 °C (Figure 8C), with the lowest values during the second wind peak. As the Shamal intensified, the daily average air temperature dropped from 20.9 °C to 18.9 °C. SST varied between 19.6 °C at the onset of the Kaus wind and 17.6 °C at the midpoint of the second Shamal peak. Relative humidity decreased sharply from 94% to 33% under the dry Shamal winds (Figure 8D). Meanwhile, passage of the associated frontal system increased the mean atmospheric pressure by 10 mb, from 1007 to 1017 mb (Figure 8E).

Heat flux patterns also reflected the influence of synoptic forcing. The net shortwave heat gain followed a diurnal cycle, its daily-average magnitude reduced to 299 W m^–2 on Mar 6 by a dust storm accompanying the Shamal winds, and peaking at 817 W m^–2 in the post-Shamal afternoon (Figure 8F). Net longwave radiation losses were minimal (∼4 W m^–2) under the cloud cover and humidity associated with the pre-Shamal Kaus winds, whereas the maximum loss (111 W m^–2) occurred during the second Shamal peak (Figure 8G).

Sensible heat fluxes followed the air–sea temperature gradient. The ocean gained heat (a downward, negative sensible flux) when the air was warmer than the SST, reaching 43 W m^–2 during the Kaus phase when the air–sea temperature difference was 3.3 °C. Conversely, the ocean lost heat, up to 69 W m^–2, when the SST exceeded the air temperature by 4.3 °C during the second peak (Figure 8H).

Latent heat fluxes were driven by evaporation, with peak losses of 406 W m^–2 during the second wind burst (Figure 8I). The net surface heat flux showed strong diurnal variability, ranging from a net loss of 701 W m^–2 during the second Shamal peak to a net gain of 511 Wm^–2 in the post-Shamal afternoon (Figure 8J).

As in SE1, the surface heat loss during SE3 was dominated by latent heat, consistent with previous northern Gulf studies (Li et al., 2020b; Senafi et al., 2019). The maximum net surface heat loss of 511 W m^–2 during SE3 exceeded that of SE1 (420 W m^–2), reflecting the stronger winds of SE3, yet both events remained below the basin-scale Shamal extremes (> 1000 Wm^–2) simulated by Thoppil and Hogan (2010) and the Mistral-driven daily losses of 600–800 Wm^–2 that force open-ocean convection in the Gulf of Lions (Houpert et al., 2016). Importantly, whereas Mistral and Bora events drive convection through hundreds to thousands of metres to form deep and dense water (Schott et al., 1996; Houpert et al., 2016; Lee et al., 2005), at the shallow (∼30 m) study site, the same latent-heat-dominated cooling homogenized the entire water column to the seabed during both SE1 and SE3, coupling the surface forcing directly to the bottom boundary layer. These observations therefore confirm the cold, dry wind convective-mixing paradigm established for these deeper systems while documenting its distinct, full-depth expression at a shallow, tidally energetic site.

3.2.2 Hydrographic response of the water column

Hydrographic conditions during SE3 (Figure 9) were driven by both convective mixing and two synoptic wind events: the pre-Shamal Kaus and the Shamal itself. Episodes of convective forcing were identified on the nights of Mar 3, Mar 7, and Mar 8, characterized by near-zero Monin–Obukhov length scales. Light penetration into the water column decreased over time, with mean light intensity weakening from 52 W m^–2 on Mar 4 to 36 Wm^–2 on the afternoon of Mar 5 during the Kaus event (Figure 9B). On the first day of the Shamal, rising dust levels further reduced surface radiation, lowering the maximum observed light intensity to 43 W m^–2.

As in the earlier Shamal event, the temperature profile (Figure 9C) followed three dominant patterns. The first was a diurnal heating cycle that responded to solar input and surface heat flux. During the low-wind evening of Mar 3 (wind speed< 5 ms^–1), convective overturning extended the mixed layer to the seabed, cooling the entire water column. On Mar 4, surface warming driven by daytime heat fluxes shoaled the mixed layer from 26.8 to 17.2 m, with surface temperatures peaking at 19.23 °C. The second pattern corresponded to changes driven by the Kaus and Shamal events. Between Mar 4 and Mar 5, Kaus conditions deepened the mixed layer from 16.5 to 17.5 m and raised the mean water column temperature from 18.63 °C to 18.73 °C as surface heating penetrated to greater depths. This warming was followed by gradual cooling during the Shamal event, ending with a mean water column temperature of 18.58 °C. The third pattern was consistent with horizontal advection. To minimize advective effects, data at 22 m depth during calm pre-Kaus conditions were analyzed. As in SE1, easterly currents (northeasterly and southeasterly) were associated with warmer water (18.48 °C), whereas westerly currents (northwesterly and southwesterly) corresponded to cooler water (18.39 °C; Figure 6C1). This pattern is again consistent with advection of warmer water from the southeast during flood tides and of cooler water from the northeast during ebb tides; as for SE1, the single-mooring configuration permits only an inferred advective contribution rather than a direct measurement of horizontal gradients.

DO concentrations ranged from 6.58 to 7.04 mg L^–1 between Mar 3 and Mar 5 (Figure 9D); data after Mar 5 were excluded owing to sensor biofouling. Vertical mixing during the convective evenings of Mar 3 and Mar 4 enhanced near-bottom DO from 6.71 to 6.87 mg L^–1. This near-bottom oxygenation is the shallow-basin counterpart of the ventilation that winter convection supplies to the deep interior of other semi-enclosed seas, such as the Red Sea (Papadopoulos et al., 2015) and the northwestern Mediterranean (Ulses et al., 2021); here, the convective overturning ventilates the near-bottom layer directly. DO also showed a pattern consistent with lateral advection (Figure 6C2), with western currents associated with marginally higher DO (6.83 mg L^–1) than the eastern flows (6.81 mg L^–1). As in SE1, part of the DO variability may also reflect temperature-dependent oxygen solubility and enhanced air–sea exchange under strong winds (Wanninkhof, 2014), in addition to lateral advection and vertical mixing.

The current structure was dominated by semidiurnal tidal forcing (Figures 9E, F), as in SE1, with northwestward flow during flood phases and southeastward flow during ebb phases. Cross-shore currents ranged from –0.20 to 0.21 ms^–1, with the extremes on the afternoon of Mar 3, and alongshore currents from –0.53 to 0.56 ms^–1 over the same period.

Acoustic backscatter, again interpreted as a qualitative proxy for suspended particulate matter, peaked during two intervals, post-Kaus on the evening of Mar 5 and during the second Shamal wind peak, with average water-column backscatter reaching 69.65 dB.

Buoyancy frequency (Figure 9H) showed maximum stratification (> 3.2 × 10–⁵ s^–2) near the mixed layer base. On Mar 3 and Mar 4, surface warming before the Kaus event increased N² to 1.1 × 10–³ s^–2 near the surface and shoaled the mixed layer to 16.5 m. As the Kaus intensified, cooling of the upper layer deepened the mixed layer to 17.5 m and reduced N² to 3.98 × 10–⁵ s^–2. Buoyancy frequency continued to decline throughout the Shamal, falling below 7.94 × 10–⁶ s^–2 during the two wind peaks, before rising to a post-event average of 1.65 × 10–⁵ s^–2.

Shear frequency (Figure 9I) peaked near the seabed (Figure 10C), particularly within the bottom 3 m, with average values of 1.3 × 10–³ s^–2. Shear increased with rising alongshore velocity, reaching 5 × 10–³ s^–2, and promoted dynamic instability where Ri< 0.25 (Figure 9J). The average water column shear decreased from 5.91 × 10–⁴ s^–2 before the Shamal to 3.29 × 10–⁴ s^–2 during the second peak, reflecting a transition to more vertically uniform velocity profiles.

The SE3 observations provide the clearest illustration of the shallow-basin convective response and connect it to the wintertime overturning of the Gulf. During the low-wind, convective evening of Mar 3, overturning extended the mixed layer to the seabed and cooled the entire water column, and the buoyancy frequency fell below 7.94 × 10–⁶ s^–2 during the Shamal peaks, indicating near-complete destratification, consistent with the full-depth homogenization reported by Li et al. (2020b) at this site. This event-scale, full-depth convection may represent a local example of the wintertime buoyancy loss that, integrated over the shallow northern Gulf, forms the dense water documented by Swift and Bower (2003), whose densest water originates in winter at the northern end of the basin. It is also the shallow-water counterpart of the deep convection driven by Mistral and cold northerly outbreaks in the northwestern Mediterranean and northern Red Sea, where the same forcing deepens the mixed layer over hundreds of metres (Houpert et al., 2016; Papadopoulos et al., 2013); at the shallow (∼30 m) study site, however, the convection reached the seabed rather than forming a deep mixed layer within a much thicker column. As in SE1, the current field remained dominated by the semidiurnal tide (Li et al., 2020b), so that in both events the Shamal acted as an episodic convective perturbation superimposed on a persistently tide-dominated flow.

3.2.3 Bottom boundary layer turbulence: observations versus model

Turbulence characteristics within the BBL (approximately 23–27 m) during SE3, comprising TKE, its dissipation rate, Reynolds stress, and vertical eddy diffusivity, exhibited a systematic response to shear-induced instabilities associated with the semidiurnal tidal forcing. These responses were most evident when the Ri dropped below 0.25 (Figure 10F), with shear frequency modulated by the tidal cycle (Figure 10C), consistent with those observed during the preceding Shamal event. Observed TKE ranged from 1.14 × 10–⁴ to 25.01 × 10–⁴ m² s^–2 (Figure 10G). As in SE1, the highest TKE was recorded during the peak alongshore current speed (0.56 ms^–1) on the afternoon of Mar 3. TKE was lower during the Shamal phases (5.22 × 10–⁴ m² s^–2) than the period average (7.14 × 10–⁴ m² s^–2). Simulated TKE from both closures was significantly correlated with the observations (p< 0.05; r = 0.71–0.73). However, the NSE values were only 0.17–0.19, indicating limited agreement in point-by-point magnitude despite representation of the general temporal variability. The MBE was ≤ 2.36 × 10–⁴ m² s^–2, and the RMSE was ≤ 4.88 × 10–⁴ m² s^–2 (Table 2).

Table 2. Statistical comparison of observed and simulated turbulence parameters during SE3, including minimum, maximum, mean values, and corresponding 95% confidence intervals (in parentheses).

The TKE dissipation rate and Reynolds stress (Figures 10H, I) closely followed the TKE and were strongly modulated by the tidal cycle, with six-day means of 1.30 × 10–⁵ m² s^–3 and 0.26 Pa, respectively. Both peaked near the end of the second Shamal wind phase, reaching 6.45 × 10–⁵ m² s^–3 and 0.88 Pa. During the two wind peaks, strong Shamal winds generated cold plumes, evidenced by vertical heat fluxes exceeding 20 Wm^–2 in magnitude (Figure 10B). These plumes homogenized the upper water column, reducing stratification and shear frequency (Figure 10C) and thereby suppressing turbulent mixing; consequently, the dissipation rate and Reynolds stress were lower during these periods. Both closures represented the suppression of mixing during the cold-plume periods and the intensification of turbulence under stronger tidal shear. Agreement was stronger for the dissipation rate and Reynolds stress than for TKE, with significant correlations (p< 0.05; r ≥ 0.84), NSE ≥ 0.67, and relatively low MBE and RMSE (Table 2). Vertical eddy diffusivity (Figure 10J) followed the general pattern of TKE and ϵ. The six-day average observed diffusivity was 0.03 m² s^–1, closely matching the averages from both models. The maximum observed K_z during SE3 was 0.861 m² s^–1. This unusually large value reflected the combination of a relatively high TKE dissipation rate and a small but positive temperature-derived N². Peak diffusivity coincided with the afternoon of Mar 3, following the maximum observed TKE and alongshore current speed. Model outputs agreed well with the observations (r = 0.78, p< 0.05; NSE = 0.60), with low MBE (< 3.2 × 10–³ m² s^–1) and RMSE (≈ 0.01 m² s^–1).

Consistent with SE1, the k–kl and k–ϵ simulations produced very similar turbulence estimates (r = 1.00 and NSE ≥ 0.98) and represented the general variability associated with semidiurnal tidal forcing and the Shamal events. This similarity reflects agreement between the two closure schemes and should not be interpreted as evidence of equivalent agreement with the observations.

The SE3 turbulence both confirms and refines the response seen in SE1. As in SE1 and in energetic tidal shelf seas (Rippeth et al., 2003; Simpson et al., 1996), the six-day mean dissipation rate and Reynolds stress (1.30 × 10–⁵ m² s^–3 and 0.26 Pa) were modulated by the semidiurnal tide. A distinctive feature of both events, most clearly expressed in SE3, is that the Shamal peaks did not enhance the near-bottom turbulence but suppressed it: convective homogenization of the velocity field reduced the vertical shear, and hence the shear production, lowering the BBL dissipation and Reynolds stress during the wind maxima. This apparent contrast with the fivefold increase in event-mean dissipation reported at this site by Li et al. (2020b) reflects the different processes involved, in that the Shamal enhances full-depth convective mixing and the event-averaged dissipation while simultaneously eroding the tidal shear that sustains the near-bottom turbulence during the wind peaks. Such competing behaviour is a direct consequence of the shallow basin, where the convection reaches the seabed and vertically homogenizes the currents rather than deepening a mixed layer within a thicker, more strongly sheared water column. Both closures reproduced this suppression and the intervening tidal intensification, consistent with the benchmark GOTM evaluation of Burchard et al. (2002) for a comparable tidally forced bottom boundary layer, although, as in SE1, they underestimated the largest short-duration peaks.

Because temperature and velocity were relaxed toward observations, the simulations were designed to evaluate the response of the turbulence closures under the measured hydrographic and current conditions rather than to test the unconstrained predictive skill of GOTM. The turbulence quantities were not prescribed or assimilated and therefore provide an evaluation of how the two closure schemes respond to the observed shear and temperature-derived stratification.

With salinity assumed to be vertically uniform, N² represents temperature-driven stratification, while Ri and K_z were calculated using this temperature-based estimate of N², as discussed in Section 2. Sensitivity calculations indicate that unresolved salinity gradients may influence the absolute magnitudes of these quantities, particularly during the weakly stratified SE3 period. Consequently, quantitative interpretations based on the Ri = 0.25 threshold and K_z should be treated cautiously. Nevertheless, the qualitative distinction between tidal shear turbulence and Shamal-induced convective mixing remains independently supported by the observed surface heat loss, temperature and velocity homogenization, tidal modulation of vertical shear, and ADCP-derived TKE, dissipation rate, and Reynolds stress.

3.3 Tidal shear and surface-forcing sensitivity in the BBL

The hydrographic and turbulence analyses across both Shamal events indicate that tidal currents play a central role in modulating the thermal structure, dissolved oxygen distribution, and turbulence characteristics of the BBL, largely through current-induced shear. To isolate the contribution of tidal forcing, harmonic analysis following Pawlowicz et al. (2002) was applied to the observed current profiles (Figures 11C, F) to separate the tidal components (Figures 11D, G) from the non-tidal (residual) flow components (Figures 11E, H).

To examine the sensitivity of the simulated BBL turbulence to the forcing configuration, simulations using both closures were forced by the extracted tidal currents, while the explicit surface wind stress and heat flux were excluded. Temperature remained relaxed toward the observed profiles, which retain the effects of Shamal-induced cooling and convection. This experiment therefore represents a surface-forcing sensitivity test rather than a dynamically independent tide-only control and does not provide a complete separation of tidal and Shamal forcing.

The analysis focuses on the strongest wind phase of SE3, the ∼24-hour period spanning the second Shamal peak (Mar 7). During this interval, the surface-forced simulations gave a lower average shear frequency (5.12×10–⁵ s^–2; Figure 10C) than those without surface forcing (5.83×10–⁵ s^–2; Figure 12B). This result, together with the full record, supports the interpretation that Shamal winds reduce BBL shear by vertically homogenizing the water column.

The increased shear in the no-forcing simulations produced more unstable conditions, reflected in lower Richardson numbers (mean Ri = 0.07; Figures 12C, E) than in the surface-forced simulations (mean Ri = 0.12; Figure 10E). Consequently, the TKE was roughly twice as high and the dissipation rate elevated by 3.20 × 10–⁶ m² s^–3 (Tables 2, 3). The Reynolds stress remained largely unchanged between the two sets, whereas the vertical eddy diffusivity differed by approximately 0.01 m² s^–1.

Table 3. Statistical comparison of observed and simulated turbulence parameters during the surface-forcing sensitivity experiment for SE3.

The simulations excluded explicit surface wind stress and heat flux but retained forcing by the extracted tidal currents and relaxation toward the observed temperature profiles, and were compared against the same SE3 observations as in Table 2. MBE was calculated as model minus observation. Mean values and corresponding 95% confidence intervals are shown in parentheses.

Although the no-forcing simulations did not match the observations during the Shamal-active periods, they agreed much better under weaker winds (e.g., Mar 3; Figures 12F, I), where model–observation comparisons yielded significant correlations (p< 0.05; r > 0.6) and positive Nash–Sutcliffe efficiencies (> 0.1). These results confirm that tidal currents are a dominant driver of BBL turbulence under low-wind conditions, whereas during strong wind events the influence of tidal shear diminishes as wind-driven mixing reduces the vertical velocity gradients.

These results reinforce the paradigm, well established for tidally energetic shelf seas, that the barotropic tide is the primary source of bottom boundary layer turbulence (Simpson et al., 1996; Rippeth et al., 2003; Rippeth, 2005), and they are consistent with the tidal dominance of the current field at this site, where tides account for approximately 97% of the current energy (Li et al., 2020b). The present analysis adds an element specific to shallow, convectively forced basins: rather than augmenting the near-bottom turbulence, the strong Shamal winds reduce the tidal shear production by vertically homogenizing the velocity field, so that surface convection and tidal shear act in opposition within the BBL. This differs from the classical shelf-sea picture, in which wind and convective forcing generally augment the turbulent kinetic energy budget (Simpson et al., 1996; Rippeth, 2005); at the shallow study site, the full-depth convection instead erodes the vertical velocity gradients that sustain shear-driven BBL turbulence. The observations therefore identify tidal shear as the dominant control on BBL mixing under quiescent conditions, with episodic Shamal convection acting as a competing, shear-suppressing influence during strong wind events.

4 Conclusions

This study investigated the hydrographic structure and turbulence dynamics of the BBL in the northwestern Arabian/Persian Gulf during two winter Shamal events. Each event lasted more than 60 hours and produced abrupt changes in wind, temperature, humidity, and pressure. Shamal winds enhanced surface heat loss, primarily through latent heat, driving convective mixing that reached the BBL. The main findings are summarized below.

The Shamal winds peaked at 15.12 ms^–1 and were accompanied by sharp reductions in relative humidity (to 33%), air temperature (by ∼4 °C), and SST (by ∼2 °C), with latent heat losses reaching ∼406 Wm^–2. Under this forcing, the temperature field reflected three processes: diurnal heating, progressive cooling during the passage of the Shamal, and an advective signal associated with the semidiurnal tidal currents. Both the alongshore and cross-shore currents followed the semidiurnal tide, reaching maximum speeds of 0.67 and 0.37 ms^–1, respectively. Within the bottom boundary layer, both the k–kl and k–ϵ closures reproduced the general temporal variability and semidiurnal modulation of the turbulence quantities but underestimated some short-duration peaks; during SE1, the peak dissipation rate and Reynolds stress were underestimated by factors of approximately five and three, respectively, and the low TKE NSE values indicated limited agreement in point-by-point magnitude. Taken together, the observations show that BBL turbulence was governed mainly by tidal shear under weak winds, whereas convective overturning during the Shamal peaks homogenized the water column, reduced the tidal shear, and suppressed shear-driven mixing.

The interpretation of these results is subject to several limitations. The constant-salinity assumption means the buoyancy frequency, Richardson number, and vertical eddy diffusivity reflect temperature-driven stratification only. The model comparison is observation-constrained rather than fully independent. The one-dimensional framework cannot resolve horizontal advection, and the vertical eddy diffusivity is sensitive to weak stratification through the Osborn relation. These constraints, discussed in Sections 2 and 3, motivate future salinity-resolving observations and three-dimensional modeling.

The work reported here addressed winter Shamals only; summer events warrant investigation to assess seasonal variability. Extending observational coverage to other sectors of the Gulf would clarify the spatial variability of wind-driven mixing. Given the strong influence of Shamal winds on both the atmosphere and the ocean, three-dimensional hydrodynamic modeling is a natural next step for representing the regional response to these events more completely.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. ChatGPT (OpenAI) was used to assist with language editing (grammar and spelling) of the manuscript text. Gemini (Google) was used to help redraw the schematic of the meteorological station setup shown in Figure 3.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2026.1949202/full#supplementary-material

Summary of meteorological instruments deployed at the study site and their corresponding Supplementary material at the source

Hydrographic instrumentation deployed on the vertical mooring and bottom pod, including sensor type, number of units, and sampling intervals. The instruments were configured to resolve temperature, salinity, dissolved oxygen, light, and velocity structure throughout the water column with high temporal and vertical resolution. Supplementary material at the source

Configuration of the one-dimensional GOTM simulations. All numerical, forcing, initialcondition, and boundary-condition settings were identical for SE1 and SE3 and for the k–ϵ and k–kl simulations, except for the simulation period and turbulence closure. Supplementary material at the source

Footnotes

Edited by: Aboobacker Valliyil Mohammed, Qatar University, Qatar

Reviewed by: Abdul Azeez Saleem, King Fahd University of Petroleum & Minerals, Saudi Arabia

Wanwei Zhang, Hohai University, China