The hydrographic data used in this study were obtained from a survey of the Cape Ghir region in northwestern Africa, which formed part of the project “Mixing Processes in the Canary Basin (PROMECA)” ( Figure 1 ). The survey was conducted aboard the R/V García del Cid in early fall (18–29 October 2010) when the Trade Winds typically weaken. Conductivity-temperature-depth (CTD), expendable bathythermograph (XBT), and microstructure data were collected at stations located approximately 10 km apart from one another along three transects (A, B, and C in Figure 1).

Wind speed, instantaneous wind speed, wind direction, air temperature, water temperature, relative humidity, air pressure, and incoming solar radiation were recorded at ~2-min intervals by a meteorological station installed onboard the ship. To remain consistent with hydrographic station measurements, the meteorological data were averaged over 2-h intervals (i.e., the time resolution of the CTD profiles).

During the cruise, current velocity data were also collected nearly continuously with a vessel-mounted 75 KHz Ocean Surveyor Acoustic Doppler Current Profiler (SADCP; Teledyne Technologies, Thousand Oaks, CA, USA). The data were processed with Common Ocean Data Access System (CODAS) software (Firing et al., 1995) to obtain vertical bin sizes of 10 m averaged over 2-h periods. An SBE911 plus CTD (Sea-Bird Scientific, Bellevue, WA, USA) was used to produce temperature, salinity, and density profiles (1 dbar vertical resolution). Between each CTD station, vertical temperature profiles were also obtained with XBT Sippican T5 probes (Lockheed Martin, Bethesda, MD, USA) that transmitted to 2000 m depth. The temperature profiles obtained with the XBT probes were smoothed using a classic, low-pass Butterworth filter. A comparison of the CTD and XBT temperature profiles revealed a discrepancy of approximately 10 m. Although uncommon, this offset has been observed in similar studies of the Cape Ghir filament (Pelegrí et al., 2005b). The results of this study are based on CTD rather than XBT vertical profiles.

A TurboMAP-L microstructure profiler Wolk et al. (2002) was used to obtain profiles of ε o . The TurboMAP-L is a vertical free-fall profiler that carries microstructure sensors, including two shear probes and an FP07 thermistor, CTD sensors, and internally mounted accelerometers. The profiler freely falls at a speed of ~0.7 m s^-1 while sampling at a rate of 512 Hz. All data were binned at 2-m intervals down to ~470 m depth and processed using TMTools v. 3.04 A.

The location, evolution, and coverage of the Cape Ghir upwelling filament during the survey were determined via sea surface temperature (SST) and chlorophyll-a (Chl) satellite images ( Figure 2 ) obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors of the AQUA and TERRA satellites. Satellite images were also downloaded from Ocean Color Web (http://oceancolor.gsfc.nasa.gov). Geostrophic surface currents were derived from the sea level anomaly provided by the AVISO altimeter products at a spatial resolution of 1/4° × 1/4°, which were downloaded with OpenDAP from the AVISO Website (http://www.aviso.altimetry.fr/en/home.html).

To obtain ε o , we first removed spikes from the turbulence profiles to obtain the shear fluctuation power spectrum, ψ ( k ) , from the fall speed of the profiler, where k is the wavenumber determined from its fall speed. Assuming isotropic turbulence Hinze, 1979, observed ε o was estimated by integrating ψ ( k ) within an appropriate wavenumber range following the methods of Oakey (1982):

where v is the molecular viscosity coefficient with a value of 1 × 10^-6 m²s^-1 and k 1 is the lowest cutoff wavenumber, which was set to 1 cpm given the physical scale of the microstructure shear probe. In contrast, the upper limit of integration, k 2 , is the highest vertical wavenumber free of noise and is usually taken to be the Kolmogorov wavenumber k c = ( 2 π ) − 1 ( ε ν − 3 ) 1 / 4 . After segmenting each microstructure profile within vertical bin sizes of 2 m, the observed shear spectrum was fitted to the Nasmyth empirical universal spectra for turbulence (Oakey, 1982; Wolk et al., 2002). This spectrum is considered to be representative of the spectral form of oceanic turbulence and is commonly used to verify shear spectrum measurements. An example of micro-shear power spectral density ψ ( k ) and Nasmyth spectra can be seen in Figure 3 for a depth range of 43 to 45 m at station 24. Given that records of micro-scale velocity shear can be contaminated by noise (e.g., mechanical vibration of the instrument and the influence of the rope at the surface) until the probe reaches a quasi-constant free-falling velocity, the data within the first 16 m of the water column were removed. In addition, the vertical profiles of ε o can exhibit large variability over short time periods. As such, the casts were repeated at least twice at each station.

To assess the relationship between atmospheric forcing and upper ocean turbulence, meteorological data were used to compute the net surface heat flux ( J q o ) as the sum of four individual components J q o = J q S w + J q L w + J q S e + J q L a , where J q S w is the net shortwave radiation flux, which is the main contributor to J q o during daytime and was directly obtained from the on-board meteorological station. Net longwave radiation ( J q L w ) was determined with the equation of Berliand and Berliand (1952):

where E L w = 0.985 is the longwave emissivity from from Dickey et al. (1994), σ = 5.67 × 10 − 8 m 2 K − 4 is the Stefan-Boltzman constant, T a i r is the air temperature measured at a height of 10 m above the sea surface, T s w is the surface water temperature, e a is vapor pressure, and F c is the cloud correction factor with values that range from 0.4 to 1 during daytime, as the sky was mostly clear throughout the study. The sensible heat flux ( J q S e ) and latent heat flux ( J q L a ) were determined by using the Tropical Ocean Global Atmosphere Coupled Ocean-Atmosphere Response Experiment TOGA-COARE code available in the Matlab Air–Sea toolbox (version 3.0; http://sea-mat.whoi.edu) developed by the air-sea fluxes science group of the TOGA COARE project, which is a version of the bulk flux described in Fairall et al. (1996).

Once the net surface heat fluxes were obtained, we computed the net surface buoyancy flux ( J b o ) as:

where g is acceleration due to gravity, C p is the specific heat capacity of seawater ( C p ¯ = 3.98 × 10 3 ± 0.42 J kg − 1 ° C − 1 ), α is the thermal expansion coefficient of seawater ( α ¯ = 2.58 × 10^-4 ± 6.92 × 10^-6°C^-1), β is the haline contraction coefficient of seawater ( β ¯ = 7.44 × 10 − 4 ± 1.08 × 10 − 6 ° C − 1 ), and ρ o is a density reference value (1026 kg m^-3). The first term corresponds to the thermal surface buoyancy flux ( J b T ), and the second term is related to the contribution of the haline surface buoyancy flux ( J b S ), as suggested by Dorrestein, 1979, where So is surface salinity and (E - P) is the difference between the evaporation and precipitation rates. Surface wind stress ( τ o ) was computed as τ o   = ρ a C d   u r 2 , where ρ a is the reference density of the air 10 m above the sea surface, and u r is the wind speed. The drag coefficient ( C d ) was derived from the equation of Large and Pond, 1981 with C d = 1.14 × 10 − 1 for u r < 10   m s − 1 and C d = ( 0.49 + 0.065 u r ) × 10 − 3 for u r   > 10   m s − 1 . Note that only the meteorological-related quantities averaged over a window of 2 h were used for calculations in this study (denoted by superscript s in Figure 4 ). The window size of 2 h was chosen to include the hour before and the hour after the release of the CTD and reflects the time usually required to complete a CTD cast.

Entrainment rates can be parameterized by solving the turbulent closure scheme of the 1D TKE budget (Niiler and Kraus, 1977). Therefore, depending on the sources and sinks that balance the TKE budget of the mixed layer, an assessment of the entrainment parameterization can also help to elucidate which sources of energy drive mixing in the upper boundary layer.

The 1D TKE budget primarily depends on two sinks (i.e., the dissipation term and buoyancy fluxes during daytime) and three sources, namely (i) the production of TKE from mechanical stirring induced by wind stress whose velocity scale is the friction velocity u * = ( τ o / ρ o ) 1 / 2 ; (ii) the generation of TKE by buoyancy forces during free-convection in which the velocity scale is the free-convection velocity w * = ( J b o h ρ ) 1 / 3 (Deardorff, 1970); and (iii) TKE produced by shear, which is parameterized via the square of the vertical shear S h 2 = [ ( ∂ u / ∂ z ) 2 + ( ∂ v / ∂ z ) 2 ] at the base of the mixed layer. In this study, this is represented by the Dynamic Instability Term, D I T = 0.5 w e ( δ u 2 + δ v 2 ) , where δ u and δ v represent velocity jumps at the base of the mixed layer calculated as the difference in velocity just below and above the mixed layer depth. All variables and parameters used in this study are summarized in Table 1 . It should be noted that other entrainment sources can be included, such as internal waves breaking at the interface (Strang and Fernando, 2001) or Langmuir vortices within the mixed layer (Flór et al., 2010), although further research is required to elucidate the roles of additional sources in entrainment phenomena.

Several equations have been proposed to compute the entrainment rate as the effective vertical velocity that transports buoyancy across the base of the mixed layer. This study focuses on the parameterization of w e following the widely used mixed layer bulk model of Gaspar (1988) ( w e G , hereinafter), in which the evolution of the mixed layer can be computed as follows:

where m _{1 =} 2.6 and m _{2 =} 1.9 are empirical constants, δ b is the buoyancy jump at the base of the mixed layer calculated as the difference in buoyancy just below and above the mixed layer depth, and ε p is the parameterized TKE dissipation rate. This equation is similar to that of Niiler and Kraus (1977). The main difference between these expressions lies in how ε p is parameterized. Equation (4) does not consider the effect of vertical shear at the base of the mixed layer as a relevant source of TKE. To take into consideration shear-driven turbulence (hereinafter w e D I T ), we also employed the equation used in the modeling study of Samson et al. (2009):

These authors included the DIT term to evaluate the response of the oceanic mixed layer during a hurricane. In this study, the effect of the DIT term was analyzed under normal wind conditions. In this case, vertical shear will be generated by other sources such as the baroclinicity of the filament front and/or breaking internal waves.

Upwelling areas are highly dynamic systems (e.g., Barton et al., 1998; Hagen, 2001; Pelegrí et al., 2005b), and studying turbulent processes within these systems is understandably challenging. Unlike in laboratory or modeling experiments, the interactions between processes that operate on different spatiotemporal scales within upwelling systems makes isolating and analyzing a single process quite difficult. In order to describe the physical processes involved in the dynamics of an area, such as within frontal systems, mesoscale eddies, and submesoscale structures, a hydrographical description of the area is required.

The SST and Chl satellite images in Figure 2 show the extension and location of the Cape Ghir upwelling filament over a two-day period (20–21 October 2010). The meander-like structure is identifiable as a cool feature with a high Chl concentration that crosses the sampling stations in transects A and B ( Figure 2 ). The filament waters exhibited a difference of ~2°C from the surrounding waters and were highly conditioned by the mesoscale dynamics of the area, as observed in the superimposed geostrophic currents from the AVISO product in Figure 2 . Although representative images with low cloud coverage were chosen as examples, additional images indicated that the filament was present during the sampling period. The filament signal intensified after the last CTD station (i.e., station 50) when stronger wind speeds were recorded. The general description of this filament agrees with those of previous observational studies of African upwelling filaments (e.g., Barton et al., 2001; Hagen, 2001), especially that of Pelegrí et al. (2005b) for the same filament and season in 1995 and 1997.

The vertical sections of salinity (S) and the potential density anomaly ( σ θ ) in Figure 5 show doming of isohalines and isopycnals in the northeastern stations of transect A ( Figures 5A, B ). In addition, the combined vertical CTD-XBT sections also show an elevation of isotherms at the same locations ( Figure 6A ). This dome-like structure seems to be well correlated with stations that cross the filament front according to satellite imagery ( Figure 2 ).

Although it was not the aim of this study, it was interesting that the temperature-salinity ( T θ − S ) relationship (not shown) indicated the presence of a subsurface salinity maximum below the mixed layer (~50 m), with mean values of S = 36.55, T_θ = 19°C, and σ_θ = [26-26.5] kg m^-3 in stations 38 and 40 of transect B ( Figures 5C, D , 6B ). This subsurface salinity maximum agrees with the one found by Pelegri et al. (2005b). In their study, these authors suggested that the presence of this subsurface maximum may indicate an interconnected horizontal recirculation cell in the Cape Ghir region. A T θ − S diagram with the same data set used in this study can be found in Arcos-Pulido et al. (2014), which shows the subsurface salinity maximum (their Figure 1 ).

An abrupt upwelling of isohalines and isopycnals ( Figures 5C, D ) and isotherms ( Figure 6B ) from station 44 to 48 was apparent in transect B, which was located mostly onshore. The combined CTD-XBT vertical section allowed for the creation of an additional transect, C ( Figure 6C ), in which isotherm outcropping was smoother than that of transect B. This may have been due to the proximity of transect B to the filament.

In general, it appeared that the rising of isotherms and isopycnals was concentrated in the first 200 m of the water column, particularly in the stations located within the upwelling filament. This finding agrees with what was reported by Pelegrí et al. (2005b), who suggested that this filament originates at shallow depths within the water column. On the other hand, the area of interest was strongly influenced by strong mesoscale dynamics. In order to effectively show these structures, geostrophic velocities derived from AVISO altimeter products were superimposed on the satellite images ( Figure 2 ). The current velocity revealed a cyclonic mesoscale structure, which seemed to induce northward along-front flow in transect A. A poleward geostrophic current crossed transect B in opposite direction to the filament flow as the result of an anticyclonic structure. The geostrophic currents revealed a cyclonic mesoscale structure, which seemed to induce southward along-front flow in transect A and northward along-front flow due to the interaction between a mesoscale cyclone and an anticyclone that crossed transect B.

Based on SST and Chl satellite imagery and its agreement with isopycnal shoaling, we concluded that the hydrographic stations could be separated into two groups: group-F, which was influenced by the upwelling filament and included stations 26–32 and 38–46, and group-nF, which included the remaining stations (24, 34–36, and 48–50). In a related study, Arcos-Pulido et al. (2014) indicated that group-F experienced enhanced nutrient concentrations below the mixed layer.

Surface fluxes exhibited a regular diurnal cycle in which both J q o and J b o were positive when upward fluxes occurred during the night ( Figures 4A, B ). The J q o values varied between -3 W m^-2 and -731 W m^-2 during daytime and between 2 W m^-2 to 309 W m^-2 at night ( Figure 4A ), which resulted in convective conditions with mean values of 1.18 × 10^-8 W kg^-1 during the night ( Figure 4B ). In this regard, fluxes of J b S were always lower than those of J b T and differed by one order of magnitude. The main contributors to the net total heat flux during the night were J q L a and J q L w ( Figure 4A ).

The air temperature was generally higher than the sea water temperature except in CTD stations 30–32 and 42–44, which had comparable T a i r and SST values ( Figure 4C ). This decrease in seawater temperature may be related to these stations being located within cool upwelling filament waters. It is noteworthy that the mean difference in temperature between waters within and outside of the filament was 1.09°C.

Weak northeasterly winds prevailed throughout the cruise ( Figure 4D ), with a mean speed of u _r = 3.78 m s^-1 and τ o values that varied from 0.0026 N m^-2 to 0.065 N m^-2 on the last day of sampling. This wind pattern was characteristic of early fall when trade winds are relatively less intense in this region. A sudden increase in wind speed was observed on 24 October ( Figure 4D ), which interfered with observations. Thus, the measurements concluded on 23 October.

Given that both wind-stress and net surface heat fluxes can act as forcing mechanisms to enhance mixing at the base of the mixed layer and increase w e , a meteorological-related quantity, such as the Monin–Obukhov length scale ( L m o ) ( Figure 4E ), is useful for determining which forcing process dominates the dynamics of the upper ocean layer. To indicate the depth at which both mechanical and convection forces are comparable, this scale can take the following form:

where κ = 0.4 is the von Kármán constant (e.g., Huffman and Radshaw, 1972). Values of L_mo > 0 ( Figure 4E ) indicate that the turbulence generated by wind stirring is suppressed by the stable stratification present during daytime, while L_mo < 0 occurs under unstable conditions at night. Moreover, in our study, |L_mo | was generally shallower than h ρ ( Figure 4E ), which suggests that wind-induced mixing does not notably control the depth of the mixed layer.

Bulk mixed layer models and consequently entrainment rates heavily depend on the depth of the mixed layer. This dependence has been attributed to velocity and buoyancy values across the mixed layer layer (Ravindran et al., 1999), which vary greatly with h ρ . Therefore, we carefully determined h ρ by comparing different methods ( Figure 7 ).

The algorithm developed by Holte and Talley (2009) can be based on the shape of either the potential density or potential temperature profile of each hydrographic station. With the exception of station 50, we found that the σ θ method from Holte and Talley (2009) resulted in an adequate fit to the observed h ρ . For station 50, the T θ method of Holte and Talley (2009) was more appropriate, which was largely due to the existence of salinity barriers around the area caused by the subsurface salinity maximum (Pelegrí et al., 2005b). This algorithm agrees rather well with the depth of the maximum gradient of the buoyancy frequency squared (N ², rad²s^-2) ( Figure 7 ) and with the threshold method of Δσ_θ = 0.03 kg m^-3 (de Boyer Montégut et al., 2004).

Further, we found that the algorithm given by Kara et al. (2000) tended to overestimate h ρ , whereas the gradient methods of Dong et al. (2008) tended to underestimate h ρ . These results agree with the findings of Holte and Talley (2009) for large databases of hydrographic profiles.

Several authors have emphasized the need to consider a boundary region of elevated and active mixing in the pycnocline at the base of the mixed layer (e.g., Dewey and Moum, 1990; Brainerd and Gregg, 1995; Nagai et al., 2005; Inoue et al., 2010; Sutherland et al., 2014), which is usually referred to as the mixing layer depth ( h ε ). The difference between h ρ and h ε is that the latter is the depth at which turbulent processes are active and maintain the homogeneity of the mixed layer and entrain buoyancy across the pycnocline. In contrast, h ρ represents the depth at which these surface fluxes have been mixed in the recent past (i.e., in a daily cycle or longer), thus h ρ represents the history of several mixing events. Depending on the time scale and spatial resolution, the use of h ε could be more appropriate when evaluating entrainment than h ρ (Brainerd and Gregg, 1995). For example, fine-scale processes forced by transient mixing events in which flux conditions can change at any given moment of the day require measurements of the active mixing layer. As such, the definition of h ε has important repercussions for Sverdrup’s critical depth theory of phytoplankton blooms (Sverdrup, 1953). Franks (2014) revisited Sverdrup’s critical depth theory and argued that h ρ , defined by density gradients, may not always reflect the intensity of active turbulence or its vertical extension. In this sense, the critical turbulence hypothesis (Huisman and Weissing, 1999) suggests that active mixing described by the definition of h ε must be considered, as it can act to reduce stratification in the upper layer, which may lead to favorable conditions that allow phytoplankton blooms to develop (e.g., Ferrari et al., 2015).

In general, h ε can be directly obtained with the ε o of microstructure data (Brainerd and Gregg, 1995; Inoue et al., 2010; Franks, 2014; Sutherland et al., 2014), microstructure temperature gradients (Nagai et al., 2005), or measurements of the largest turbulent overturning length scales such as Thorpe scales (Brainerd and Gregg, 1995). In this study, we determined the base of the mixing layer when the dissipation rate decreased by two orders of magnitude from the surface, which was similarly undertaken by Inoue et al. (2010). However, not all ε o profiles exhibited clear transitions. Therefore, a turbulent length scale, such as the Ozmidov scale L o (Ozmidov, 1965), could help clarify the extension of the active mixing or h ε , where L o takes the following form:

In such a case, L o is interpreted as the size of the largest turbulent eddy in a region of stable stratification. An example of h ε and L o can be viewed in Figure 8 , in which the maximum vertical gradient of L o agrees with a decrease in ε o of one order of magnitude within the upper boundary layer.

Dewey and Moum (1990), in their examination of microstructure data in upwelling filaments off northern California, noted that entrainment can only take place when h ε exceeds h ρ ; otherwise, when mixing is not strong enough to overcome the effects of stratification in the pycnocline, entrainment cannot occur. Inoue et al. (2010) also distinguished between h ε and h ρ to measure entrainment heat fluxes based on dissipation profiles across the Gulf Stream provided that h ε was deeper than h ρ . The same approach was considered in this study by defining the entrainment zone or interface as Δ h = h ε − h ρ .

When the active mixing layer extends deeper than the mixed layer, that is h ε > h ρ (Dewey and Moum, 1990; Brainerd and Gregg, 1995; Inoue et al., 2010), non-turbulent fluid from the pycnocline is subject to sufficient mixing to break stratification at the pycnocline and enters the turbulent fluid.

Some authors have speculated that enhanced turbulence below h ρ could be related to shear instabilities produced by internal waves as a result of impinging convective plumes from beneath the mixed layer (Moum et al., 1989). MacKinnon and Gregg (2005) found higher correlations between ε o and internal wave energy than between ε o and τ o below the mixed layer. Here, an example of this can be seen in Figure 8 , in which h ε is deeper than h ρ for a station in group-F ( Figures 8A, B ). In contrast, h ε and h ρ exhibited the same depth in a station of group-nF ( Figures 8C, D ). Both stations were measured during daytime conditions (i.e., during the restratification of the mixed layer when convection forcing is absent).

A complete view of active turbulent mixing is given in Figure 9 , which shows vertical profiles of ε o and superimposed profiles of L o scales computed for each microstructure station with their respective h ρ and h ε . Interestingly, h ε was deeper than h ρ (i.e., Δ h > 0 ) in filament stations (thick profiles in Figures 9 , 10A, B ). This suggests that only filament stations are susceptible to diapycnal entrainment by surface-forced processes (Dewey and Moum, 1990; Fer and Sundfjord, 2007; Inoue et al., 2010). This may be the result of elevated ε o at group-F stations and the shallow mixed layers produced by isopycnal outcropping due to the isopycnal morphology of the frontal filament system.

It is important to note the enhancement of ε o below the mixed layer of filament waters, which can be seen in Figure 9 . A weak wind regime during sampling days, as evidenced by the L m o scales ( Figure 4E ), indicates that the turbulence generated within the upwelling filament by wind stress was likely small when compared to turbulent diapycnal buoyancy fluxes or vertical shear at the base of the mixed layer. However, this increase in turbulence has also often been linked with elevated vertical shear levels ( Figure 11 ).

Vertical current shear was relatively high in areas associated with isopycnal outcropping, such as in group-F stations ( Figure 11 ), thus we examined the nature of these shear values. To ensure consistency when comparing ADCP- S h 2 with geostrophic shear ( S h g 2 ), the ADCP velocities were rotated to the same reference system and interpolated between pairs of stations. Geostrophic shear was computed through the thermal wind relation by setting a reference level of no motion at 650 m (see Pelegrí et al., 2005b). In Figures 11A, C , it can be seen that the S h 2 in group-F stations is mainly driven by the resolved component of the geostrophic velocity field in transect A. In transect B, the isopycnals outcropped towards the coast due to the presence of the filament and proximity to the coastal transition zone ( Figure 11D ). Here, the resolved geostrophic component became less relevant and ageostrophic effects arose to force total vertical shear in the stations of group-F ( Figures 11B, D ). These results were also supported by the geostrophic currents of altimetry data (see Figure 2 ). Further, the presence of a southward geostrophic current agreed with the predominance of geostrophic shear in transect A. In contrast, a departure from geostrophy appeared to take place in the coastal transition zone of transect B, in which the filament flowed southward in the opposite direction of geostrophic flow ( Figure 11D ), which was also observed by Dewey et al. (1993). As the filament moves offshore, it is balanced geostrophically, as was observed in transect A ( Figure 11C ). The enhancement of ε_o within the filament waters of group-F, the increase of vertical shear at the base of the mixed layer, and the resulting Δh > 0 are summarized in Figure 12 .

In a modeling study, Small et al. (2012) argued that shear instability below the mixed layer may be generated in large part by propagating inertial waves, which ultimately increase the entrainment of cold waters from below the thermocline. Although it was not the aim of the present study to evaluate inertial wave activity, the lack of a time series of velocity currents did not allow for the presence of such near-inertial oscillations to be corroborated.

The entrainment of cold, nutrient-rich waters from below the pycnocline to the surface layers and the subsequent deepening of the mixed layer are sensitive to turbulent mixing and the depth of the upper boundary layer. However, most bulk mixed layer models (Kraus and Turner, 1967; Niiler and Kraus, 1977; Gaspar, 1988) do not account for either (a) vertical shear as a source of TKE or (b) variable h ε . Thus, such models are inappropriate for frontal areas with baroclinicity in which shallow h ρ values are driven by the rising of the pycnocline and in which breaking internal waves may be important for the generation of strong sheared flows.

Accordingly, in the case of (a), we took into account the DIT term, which is described in Eq. (5). It is assumed that turbulent entrainment can occur only when active mixing in a stratified fluid is strong enough to overcome buoyancy effects and erode the pycnocline, that is when Δ h > 0 (Dewey and Moum, 1990; Inoue et al., 2010). Given that we had direct measurements of h ε , we restricted our analysis to stations in which ε o was available and Δ h > 0 . In the case of (b), we considered the diapycnal turbulent buoyancy flux integrated over the entrainment zone, ( J b Δ h ) I . This term represents a TKE sink in the entrainment zone when Δ h > 0 , as the energy is consumed when stratification breaks down. The procedure to obtain ( J b Δ h ) I is described in Appendix A. Thus, the integration of the dissipation term over the entire active mixing layer ∫ − h ε 0 ε o ( d z ) = ε o ¯ h ε should also be considered, similarly to what can be observed in Dewey and Moum (1990). By including these considerations in Eq. (5), the TKE-based entrainment parameterization can be described by the following equation:

Unlike other entrainment parameterizations, Eq. (8) considers TKE sources and sinks due to mechanical stirring by means of u * , vertical shear at the base of h ρ through the DIT term, buoyancy fluxes at the surface and turbulent buoyancy fluxes at the entrainment zone (i.e., J b o and ( J b Δ h ) I ), and TKE dissipation across the layer in which active mixing operates. The results from Eq. (8) yield mean entrainment rates of w e Δ h ¯ = 6.99 × 10 − 5 ± 5 × 10 − 5  m s − 1 (~6 m d^-1), which differ by a factor of 6 when compared to those of w e G ¯ = 1.51 × 10 − 5 ± 5 × 10 − 5 m s − 1 (1 m d^-1) ( Figure 10 ). The addition of the DIT term resulted in a slight increase in the entrainment rates, with a mean value of w e D I T ¯ = 1.92 × 10 − 5 ± 1 × 10 − 5 ± m s − 1 (~2 m d^-1).

It is important to note that although this analysis was based on a 1D approach, several scales can coexist within ocean fronts. Two-dimensional (2D) submesoscale processes in an upwelling filament may play major roles in mixing by serving as pathways that convey the kinetic energy extracted from the mean flow to the final dissipation range throughout the secondary instabilities of the submesoscale flow (Capet et al., 2008; D’Asaro et al., 2011; Thomas et al., 2013). To determine whether this type of process affected the front in this study, the vertical distributions of potential vorticity and submesoscale energy sources were evaluated (Appendix B). As a first approximation, our results indicate that at the moment of sampling, submesoscale 2D instabilities, such as symmetric instabilities, appeared to weakly affect the frontal structure. The cross-front wind in transect A (see Figure B.1 in Appendix B) was likely insufficient to induce Ekman transport, which would cause the front to be susceptible to symmetric instabilities (e.g., D’Asaro et al., 2011). Thus, a 1D perspective based on the TKE balance is appropriate in this study. Furthermore, D’Asaro et al. (2011) indicated that although turbulent mixing near fronts can arise from these instabilities, this does not seem to be responsible for the observed elevated mixing in the Cape Ghir upwelling filament front.

The TKE budget and consequently entrainment are highly sensitive to the parameterization of the dissipation term, which has proved to be a dominant sink in the energy balance (e.g., Deardorff, 1983; Gaspar, 1988; Wada et al., 2009). Our results highlight that the role of the dissipation term may be particularly important in frontal areas such as upwelling filaments. However, a dependence on observational ε o values may make it difficult to include parameterizations of the dissipation term in bulk mixed-layer models.

We assumed a well-mixed upper ocean layer and averaged ε o from the near-surface (far enough away from the effects of breaking waves at the surface ~16 m) to h ε and compared this to the parameterized ε G included in w e G following the methods of Gaspar (1988). Our results, which are shown in Table 2 , indicate that the normalized parameterized dissipation term ( ε G ¯ | h ρ ) can be overestimated by two orders of magnitude with respect to the observed dissipation term ( ε o ¯ | h ε ).

In the context of a 1D TKE balance, most of the energy transferred to the interface is dissipated by viscosity. The remaining energy is stored in the form of potential energy that must work against buoyancy forces or as energy that will be subsequently converted to kinetic energy by entrainment. As a result of this 1D balance, the dissipation term may be overestimated due to entrainment velocities that are lower than expected, which suggests that the parameterized dissipation term extracts more energy from the system and leaves little energy for entrainment. Thus, decreasing the dissipation term in entrainment parameterizations to match the observed dissipation term will lead to an enhancement of w e . This can be seen in Figure 10 , in which w e Δ h exhibits velocities that are higher than those given by w e G and w e D I T . This suggests that the dissipation term within the entrainment parameterization of Gaspar (1988) might not be well resolved in frontal systems.

We also compared ε o to the internal wave scaling ( ε G 89 ) proposed by Gregg (1989), which takes the internal wave spectrum into account:

where No = 5.2 × 10^-3 (rad s^-1), and S_GM = 1.91 × 10^-5 (s^-1) is the Garret-Munk shear spectrum (Gregg, 1989).

Although the fine-scale parameterizations of ε G 89 differ from the observed values of ε o , when overestimating or underestimating ε G 89 along the water column ( Figure 13 ), the averaged dissipation term ε G 89 ¯ | h ε 10 G 89 resembles the observed dissipation term ε o ¯ | h ε , particularly when this is averaged within 10 m depth intervals ( ε o 10 ¯ | h ε 10 ) ( Table 2 ). It is worth noting that the mixing layer depths obtained from ε G 89 are quantitatively similar to those given by the observed ε o 10 ( Figure 10A ). A comparison between mixing layer depths and the resulting entrainment zones is shown in Figures 10A, B . Note that for coarser resolutions, h ε 10 is always deeper than the mixed layer depth, resulting in Δ h > 0 for all stations. Moreover, the results also show agreement between ε o and ε G 89 in that both show enhanced turbulent mixing in stations of group-F ( Figure 13 and Table 2 ).

These results indicate that the observed dissipation term in the proposed entrainment parameterization of Eq. (8) could be replaced by the scaling of Gregg (1989) (hereinafter w e G 89 ). The mean value of w e G 89 was 6.91 × 10^-5 ± 5 × 10^-5 m s^-1 (6 m d^-1), which is similar to w e Δ h , as indicated in Figure 10B . This result also suggests that the role of internal waves acting at the base of the mixed layer should not be ignored given the agreement between w e G 89 and w e Δ h .

It is to be expected that the enhancement of diapycnal entrainment rates observed in group-F stations would favor the elevation of deep nutrient-rich water to the photic layer. In the same study region, Arcos-Pulido et al. (2014) observed that onshore stations, which belong to group-F, exhibited diapycnal nutrient fluxes, F z ( N u t ) , which were one order of magnitude greater than those observed in offshore stations calculated at a reference layer located 20 to 30 m below the mixed layer.

To verify that the enhancement of entrainment rates in the upwelling filament contributed to increasing nutrient availability below the mixed layer, we computed F z ( N u t ) immediately below h ρ for each microstructure station with the following equation:

where K z is the diapycnal diffusivity coefficient given by the classic expression ( K z = Γ ε o / N 2 ) of Osborn (1980). Buoyancy diffusivity and heat were assumed to be of the same order. The mixing efficiency (Г) is the ratio between the vertical buoyancy flux and ε o set to 0.2 (Oakey, 1982). The nutrient concentration, which consists of the sum of nitrites and nitrates (Nut), was linearly interpolated to a vertical 2-m grid given the microstructure vertical resolution.

Nutrient gradients ( ∂ N u t / ∂ z ) were calculated over depths ranging from 5 m to 20–30 m below h ρ depending on the availability of nutrient (Schafstall et al., 2010). This differs from what was reported by Arcos-Pulido et al. (2014), who calculated fluxes for a reference layer far below h ρ . Then, K z was averaged for the same depth range over which ∂ N / ∂ z was calculated. Positive values of F z ( N u t ) imply an upward transport of nutrients from the pycnocline to the surface layers, where they may then be assimilated by phytoplankton.

The diapycnal nutrient fluxes ( Table 3 ) indicated that stations of group-F exhibited large F z ( N u t ) immediately below the mixed layer associated with maximal diapycnal mixing, which was three times larger than those of group-nF stations. However, station 34 was an exception to this pattern. Station 34 was located adjacent to the filament edge (see Figure 6A ) and was surveyed under weak wind, non-favorable convection, and low δ V conditions. Station 34 also exhibited more intense Lo at the mixed layer depth ( Figure 9 ) and high nutrient gradients ( Table 3 ) when compared to the other stations in its group. In this case, Δ h was nearly zero at station 34, indicating that both h ρ and h ε were similar. Stations of group-F that were near station 34 (i.e., 30, 32, and 36) also exhibited high vertical nutrient gradients (not shown). This suggests that it cannot be ruled out that station 34 was influenced by either the upwelling filament or other processes, such as breaking internal waves or Langmuir cells, capable of promoting intense vertical mixing at the base of the mixed layer (Flór et al., 2010; Schafstall et al., 2010). However, the lack of microstructure observations and time series of physical variables did not allow for an evaluation of the underlying causes. Stations 48 and 50 also exhibited positive F z ( N u t ) . Despite forming part of group-nF, these stations exhibited Δh > 0, indicating that favorable entrainment conditions were present, yet their fluxes were lower than those of group-F stations

In general, stations of group-F exhibited diapycnal nutrient fluxes that were two orders of magnitude larger than those of group-nF stations. The values of F z ( N u t ) for stations of group-F were at the upper end of those reported for the North Atlantic (Bahamón et al., 2003; Dietze et al., 2004; González-Dávila et al., 2006; Mouriño-Carballido et al., 2011) but of the same order of magnitude as those of other EBUS sites (Hales et al., 2005; Hales et al., 2009; Li et al., 2012). Most studies conducted in the North Atlantic have not included calculations of nutrient fluxes immediately below the mixed layer using direct microstructure observations nor in highly turbulent environments such as frontal systems. The values reported for frontal systems in other EBUS range from 5.2 (mmol m^-2 d^-1) in the California Current EBUS to 19.11 (mmol m^-2 d^-1) in the Chilean EBUS (Hales et al., 2005; Hales et al., 2009; Li et al., 2012; Corredor-Acosta et al., 2020). The latter values agree with our maxima of 8.34 (mmol m^-2 d^-1) for station 42 (group-F) and average of 1.35 (mmol m^-2 d^-1) for all stations influenced by the upwelling filament front.