Dynamics of turbulent western-boundary currents at low latitude in a shallow-water model

. The dynamics of low latitude turbulent western-boundary currents (WBCs) crossing the Equator are considered using numerical results from integrations of a reduced-gravity shallow-water model. For viscosity values of 1000 m 2 s − 1 and greater, the boundary layer dynamics compares well to the analytical Munk-layer solution. When the viscosity is reduced, the boundary layer becomes turbulent and coherent structures in the form of anticyclonic eddies, bursts (violent detachments of the viscous sub-layer, VSL) and dipoles appear. Three distinct boundary layers emerge, the VSL, the advective boundary layer and the extended boundary layer. The ﬁrst is characterized by a dominant vorticity balance between the viscous transport and the advective transport of vorticity; the second by a balance between the advection of planetary vorticity and the advective transport of relative vorticity. The extended boundary layer is the area to which turbulent motion from the boundary extends. The scaling of the three boundary layer thicknesses with viscosity is evaluated. Characteristic scales of the dynamics and dissipation are determined.


Introduction
Strong western-boundary currents (WBCs) are a dominant feature of the world's oceans.They are also present at low latitudes in the Atlantic and the Indian oceans, where they are called the North Brazil Current (NBC) and the Somali Current (SC), respectively.We refer the reader to Richardson et al. (1994), Garzoli et al. (2003) and Fratantoni and Richardson (2006) for a detailed discussion of the NBC and the subtropical gyre in the Atlantic Ocean.A detailed description of the circulation of the northern Indian Ocean, of which the SC is the most energetic part, is given by Schott and McCreary (2001), Schott et al. (2009), Beal and Donohue (2013) and Beal et al. (2013).These currents are variable in time.Part of this time dependence is due to the timedependent wind forcing and the other part is due to the internal dynamics of the ocean.In the present study we completely neglect the former by using time-independent forcing, and the latter is the subject of the present publication.Even when subject to time-independent forcing, low latitude western-boundary currents retroflect (i.e., separate away from the boundary and turn anticyclonically for more than 90 • ) and form anticyclonic eddies.The NBC retroflect near 6-8 • N and sheds eddies exceeding 450 km in overall diameter (see, e.g., Richardson et al., 1994;Garzoli et al., 2003;Fratantoni and Richardson, 2006).The SC and the East African Coastal Current retroflect to form eddies called the Great Whirl (GW) between 5 and 10 • N and the Southern Gyre (SG) near the Equator with overall diameter between 350 and 540 km (see, e.g., Schott and McCreary, 2001;Beal and Donohue, 2013;Beal et al., 2013).There are substantial differences between the near-surface circulation in the low latitude Atlantic and Indian Ocean, mainly due to the distinct wind forcing and coastline geometry.
There is a large number of numerical work on the dynamics of the SC and the NBC with a realistic coast line and topography (see, e.g., Fratantoni et al., 1995;Barnier et al., Published by Copernicus Publications on behalf of the European Geosciences Union.2001; Wirth et al., 2002;Garraffo et al., 2003).Although realistic models permit representation of the observed features of the world's oceans, it is difficult to learn about isolated processes because all the phenomena take part simultaneously in the dynamics and interact non-linearly.The only way to guarantee our understanding of the ocean dynamics is to decompose it into processes.
Idealized rectangular basin studies which address the formation and dynamics of the large anticyclones have been previously performed (see, e.g., Cox, 1979;Lin and Hurlburt, 1981;Philander and Pacanowski, 1981;Philander and Delecluse, 1983;McCreary and Kundu, 1988).A detailed determination of the vorticity balances, fluxes and stability of idealized low latitude turbulent WBCs has been accomplished by Edwards and Pedlosky (1998a) and Edwards and Pedlosky (1998b) on the WBC, and by Fox-Kemper (2005) on the dynamics of single and multiple gyres in a barotropic constant depth β-plane model.All the above works proved with no doubt that the large anticyclones and their nonstationary dynamics have a crucial impact on the mean circulation in the boundary regions.All these studies focused on the larger-scale features of the WBCs such as the large anticyclones.
It has been shown in engineering fluid dynamics that boundary-layer dynamics involves different types of coherent structures (see, e.g., Aubry et al., 1988;Robinson, 1991).In satellite observations of the SC, small flanking cyclones correlated with the large anticyclone tend to circulate clockwise around it (Beal and Donohue, 2013).The coarse resolution in space and time of satellite data do not allow for a detailed study of these small-scale structures.Such flanking vortices are also present in the laboratory experiments of geophysical fluid dynamics (see, e.g., Van Heijst and Flor, 1989) and are also clearly visible in fine-resolution realistic simulations of the ocean dynamics.
The purpose of the present work is the identification and the study of the smaller-scale coherent structures, their interaction and their influence on the large-scale circulation.Indeed, to the best of our knowledge, there is so far no description or theory of near-coastal turbulence in the western-boundary current that goes beyond the large anticyclonic eddies.For oceanic WBCs in general, the quantitative description is mainly based on Munk-layer theory (Munk, 1950) or inertial-layer theory (Stommel, 1995;Fofonoff, 1954;Charney, 1955) and the analysis of their stability (see, e.g., Edwards and Pedlosky, 1998b;Ierly and Young, 1991).This is in stark contrast to engineering fluid dynamics, where the turbulent boundary-layer theory is the leading domain since its birth in the beginning of the 20th century (Prandtl, 1904).In this article we study the dynamics of low latitude turbulent WBCs in the viewpoint of boundary-layer theory with emphasis on coherent structures.
In a highly idealized equatorial single gyre configuration, we focus on the dynamics of low latitude turbulent WBCs to determine its structure, its dependence on the Reynolds number, by varying the viscosity between experiments.The experimental setup comprises essential prerequisites such as a fine resolution throughout the entire domain in both horizontal directions and long-time integrations to obtain statistically converged results.
The physical situation considered, the mathematical model to study its dynamics and its numerical implementation are detailed in the next section.Results on the taxonomy of the coherent structures, the turbulent fluxes, their parameterization and the vorticity balance are given in Sect. 4 and discussed in Sect. 5.

The physical problem considered
The basin is a rectangular box that straddles the Equator with dimensions L x ×L y (zonal width and meridional width, respectively; values are listed in Table 1).It spans from 1000 km south of the Equator to 3000 km north.The domain extends further northward than southward, as our research is directed towards studying the SC and NBC, ranging within the most energetic structures in the world's oceans and both occurring north of the Equator.The model is comparable to those used in idealized configuration to study mid-latitude gyres (see, e.g., Jiang et al., 1995;Sushama et al., 2007;Speich et al., 1995) and low latitude WBCs (see, e.g., Edwards and Pedlosky, 1998a;Edwards and Pedlosky, 1998b;Fox-Kemper, 2005).A Cartesian grid is used and the Coriolis parameter varies linearly with latitude; this geometry is usually referred to as the equatorial β-plane.We further suppose that the dynamics considered is this of an homogeneous fluid layer of an average thickness H which superposes a constant density motion-less fluid layer of infinite depth.The density difference between the layers is expressed by the reduced gravity g .The values of these parameters are inspired by the water-mass properties in the Indian Ocean.The layer is forced by a wind shear at its surface.

The mathematical model
The governing reduced-gravity shallow-water equations are where u and v represent the zonal and meridional velocities, respectively, and η is the variation of the layer thickness.The Coriolis parameter is given by is the reduced gravity, and g the acceleration of gravity.The parameters for the experiments performed here are listed in Table 1.The system is subject to a meridional wind-stress forcing τ y and no-slip lateral boundary conditions.A Laplacian lateral diffusion with a viscosity ν is used.This is necessary to insure the no-slip lateral boundary condition and its role is also to prevent the accumulation of energy/enstrophy at the smallest scales that are resolved numerically.Please see Frisch et al. (2008) for a detailed discussion of this bottleneck phenomena.
The associated equation for vorticity is or its conservative form: where F is the curl of the forcing.The wind stress in Eq. ( 2) is The form of the wind stress is chosen to have a strong alongshore wind stress at the western boundary (please see Mc-Creary and Kundu (1988) for more detailed explanation) and an inversion of the wind stress at about 3500 km from the western boundary.It corresponds roughly to the Indian Ocean during summer monsoon wind forcing.To avoid the initial shock, the wind stress increases exponentially from zero with spin-up time of t c = 180 days.This wind forcing leads not only to a single gyre extending over the entire domain but also to an almost vanishing zonal velocities.The zonal velocities lead to an inertial boundary current (Charney (1955); see also Pedlosky (1979)) and have a stabilizing (when westward) or destabilizing (when eastward) effect on the western-boundary current.This behavior is subject of a future publication.

The numerical implementation
The numerical method used to solve the Eqs.( 1)-( 3) is a centered, second-order finite difference scheme in space, using an Arakawa A-grid, and a second-order Runge-Kutta scheme is used for time stepping.A fine numerical resolution of square geometry ( x = y = 2.5km) is employed throughout the entire domain.The scheme was successfully tested changing resolution in space and time in Wirth (2013).This uncommon choice, of not using grid refinement at the boundary, is justified by the results presented in Sect.4, where it is clearly seen that for high Reynolds number flow, parts of the viscous sub-layer (VSL) are torn off the wall and transported away from it by the surrounding turbulent flow.This leads to small-scale structures, also far from the boundary.Such a kind of process can only be represented when there is fine resolution in both horizontal directions throughout the extended boundary layer (to be defined in Sect.4.4).Please note that the resolution is well below the Munk scale δ M = (ν/β) 1/3 , which is around 18 km in the lowest viscosity experiment.We favor fine-resolution rather than highorder schemes as it also insures an isotropic representation of the smaller-scale structures.The time step is 90 s, which is almost 5 times shorter than the Courant-Friedrichs-Lewy (CFL) time step imposed by the speed of the flow and the gravity waves.In the nonlinear boundary layer the high vorticity in the boundary layer (VSL, detailed in Sect.4.5) is intermittently torn off the boundary.This process is the equivalent of bursts in 3-D boundary layers (see, e.g., Robinson, 1991).To insure a correct representation of this intermittent, rapid and violent process and its nonlinear evolution, a short time step was used.The physical parameters are such that in the present dynamics a vanishing of the fluid layer (outcrop) does not occur.

Experiments
The system is forced by the prescribed wind field and energy is dissipated by viscosity ν.The spatiotemporal complexity of the system's behavior increases with decreasing viscosity ν.The other model parameters are kept constant and are given in Table 1.Experiments with different values of the viscosity were performed.The name of experiments is referred to by the abbreviation "EXP" followed by its viscosity value: EXP1000 is an experiment with a viscosity ν = 1000 m 2 s −1 .The highest viscosity experiments with ν = 1000 m 2 s −1 converged with time towards a laminar dynamics, the corresponding Reynolds number based on the maximal value v 0 of the time-averaged meridional velocity in the boundary current and the Munk-layer thickness δ M = (ν/β) 1/3 (where ν and β are model parameters and thus obtained before running EXP1000) at y = +1500 km is Re = v 0 δ M /ν = 42 .The transport in the boundary layer is, to leading order, imposed by the wind forcing over the entire basin.This leads to the fact that velocity times the boundary layer thickness is constant and the Reynolds number scales as the inverse of the viscosity.The numerical resolution, the numerical scheme and the model parameters allowed one to perform calculations with viscosities down to ν = 300 m 2 s −1 .
In the high viscosity experiment the boundary layer dynamics converges towards a stationary state in about 3000 days of the dynamics.Lower viscosity experiments converge to a statistically stationary state.To increase the significance of the statistics, experiments were performed for 5000 days of the dynamics and averages used herein were calculated over the last 2000 days.

Large-scale circulation
Figure 1 shows the layer thickness variation (η) contours and horizontal velocity arrows from experiment EXP1000, at t = 2000 days, when the dynamics has converged to a stationary state.This figure shows the classic Sverdrup interior solution with a Munk boundary layer.
For all the experiments, strong western-boundary currents crossing the Equator in the northward direction, with a recirculation in the rest of the domain, were observed, forming a single gyre circulation.
In experiments with lower viscosity, time dependence arises in the form of coherent anticyclones moving northward along the western boundary.For the lowest viscosity experiments the dynamics is fully turbulent in the vicinity of the western boundary, with chaotic motion over a range of spatial scales (see Sect. 4.4).The time-averaged large-scale circulation of the low viscosity experiments is qualitatively similar to the stationary flow of higher viscosity (ν ≥ 1000 m 2 s −1 ) experiments.

Laminar western-boundary layer
For the higher value of the viscosity, the stationary solutions of the boundary layer are, to leading order, given by a balance of the meridional transport of planetary vorticity (4th term in Eq. 6) and the viscous dissipation (last term on the left-hand side of Eq. 6).This dynamics is described by the Munk-layer theory (Munk, 1950;Pedlosky, 1990) and the solution is where δ M = (ν/β) 1/3 is the characteristic boundary layer thickness of the Munk layer and v 0 M is a velocity scale.The vortex stretching is given by the fifth term in Eq. ( 6).We found the vortex stretching to be important only very close to the boundary.It decreases rapidly with the distance from the boundary (∼ 20 km, not shown), before the meridional velocity reaches its maximum.As zonal velocities and meridional vorticity gradients are small so is the advection of relative vorticity (second and third term in Eq. 6).The analytic solution of the Munk theory, for ν = 1000 m 2 s −1 vanishes at 2π √ 3 δ M = 133 km.The laminar experiment EXP1000 has a vanishing meridional velocity (width of the boundary current) of around 150 km (not shown).The meridional velocity is thus close to the Munk-layer solution.This is in agreement with the results of Edwards and Pedlosky (1998a).It is noteworthy to mention that the zonal velocity vanishes almost completely, except near the southern and northern boundaries of the domain and thus inertial effects (the zonal transport of vorticity, see Fofonoff, 1954;Charney, 1955;Stommel, 1995;Pedlosky, 1990) can be neglected in the present experiments.

Anticyclones
In experiments with higher Reynolds number (lower viscosity), the flow becomes time dependent and coherent structures appear.The most conspicuous coherent structures are the anticyclonic eddies along the western boundary (Ierly and Young, 1991;Edwards and Pedlosky, 1998a;Schott and Mc-Creary, 2001;Wirth et al., 2001;Richardson and Schmitz, 1993).In our experiments they start to appear at viscosity values of ν = 1000 m 2 s −1 as transient features during the spin-up in the form of poleward traveling waves in the boundary layer.They travel northward along the boundary at a speed of V eddy ≈ 2.3 • 10 −1 m s −1 .In experiments with lower viscosity values their size increases.At a viscosity of ν ≈ 500 m 2 s −1 , they are coherent regular vortices.Their diameter is then around the equatorial Rossby radius of deformation L β = g H /β = 350 km, a size that compares well to the size of the eddies in the SC (Schott and Mc-Creary, 2001;Wirth et al., 2001) and to the eddies of the NBC (Richardson and Schmitz, 1993).When inspecting the potential vorticity (PV) they appear as negative PV anomalies that move poleward with an average speed of V eddy ≈ 1•10 −1 m s −1 , while the fluid velocity in their interior reaches a speed of v eddy = 2 m s −1 .This demonstrates that the eddies are advected water masses and not a wave-like phenomena.In the literature, eddy or ring are often used interchangeably to denote the same object.A closer inspection (not shown) of the velocity field shows that they are eddies in almost perfect solid-body rotation.They are not vortex rings with an almost motionless core (eye).With decreasing viscosity their shape and poleward displacement exhibit a random-like behavior (Wirth et al., 2001) as can be verified analyzing Hovmöller diagrams (not shown) indicating a chaotic dynamics.For the lower viscosity values the eddy dynamics becomes more chaotic, some of the eddies migrate into the interior of the basin, merge with other eddies or are disintegrated by them in a 2-D-turbulent eddy dynamics.At the lowest viscosity value of ν = 300 m 2 s −1 , the average northward displacement velocity is around V eddy ≈ 6 • 10 −2 m s −1 , while the fluid velocity in their interior reaches a speed of v eddy ≈ 2.4 m s −1 .

Bursts
In all calculations a VSL is present at the boundary, where the vorticity has large positive values and the viscous dissipation is a dominant process.It is a thin layer of a few tens of kilometers thickness, for the lower viscosity values.It is discussed in detail in Sect.4.5.For the lowest values of the viscosity, intermittent detachments of the VSL just northward of the eddy center are observed at the boundary (see Fig. 2).The detachments are the most violent phenomena in the simulations (with viscosities ν = 500 m 2 s −1 and lower) with the strongest velocity and vorticity gradients.When the sheet of positive vorticity (the VSL) along the western boundary in the Munk layer breaks due to the action of an anticyclone, the southern part of the VSL detaches and is torn off the boundary by the anticyclone and accelerates away from the boundary (see Fig. 2).North of the detachment, the vorticity anomaly and the meridional velocity are negative.The northern part of the VSL continues to flow northward along the boundary.These events are the analog to bursts or ejections in the classical boundary layer (Robinson, 1991) and are thus given the same name here.They are strongly spatially localized and temporally intermittent ejections of fluid and vorticity away from the wall, initiated by the large anticyclonic eddies.The separation of the boundary layer plays a key role in boundary layer dynamics; see, e.g., Prandtl (1904) and also Schlichting and Gertsen (2000).
The ejection of the VSL and its offshore transport asks for fine resolution in both horizontal directions not only in the vicinity of the VSL but also in areas to which the boundary layer fragment is transported.
In our analysis we identify bursts as events of negative meridional velocity in the VSL.Please note that the dynamics in the VSL is not laminar, as bursts are intermittent and violent detachments of the VSL.This feature is also found in turbulent wall-bounded flows in engineering applications (Robinson, 1991).To quantify the occurrence of bursts, the fraction along the western boundary in the interval y ∈ [+125, +2250km] at which a flow reversal occurs is calculated and then average over time, to obtain the value T presented in Table 2.For viscosities ν = 1000 m 2 s −1 or greater there are no bursts.Bursts are observed for ν = 500 m 2 s −1 and lower.The fraction of time with flow reversal strictly increases with decreasing viscosity in all the experiments performed and reaches values of around 19 % for the lowest values of the viscosity, showing that they are a recurrent dominant feature of low viscosity boundary currents when inertial effects are absent.

Dipoles
In many instances the positive vorticity anomalies, ejected from the boundary during bursts (fractions of the VSL), pair with negative vorticity anomalies from within the anticyclones and form asymmetric dipoles (see Fig. 2) which then travel ballistically (at almost constant velocity) over distances of several eddy diameters.The size of the dipoles measured by the distance of the vorticity minima and maxima spans between the thickness of the viscous boundary layer δ ν (see below) and the size of the coherent anticyclones.Two other scenarios (not shown) are observed for the smaller cyclonic vortex.It can separate from the large anticyclonic eddy and drift offshore.It can also remain attached to the large anticyclonic eddy, circle around it and return to the western boundary, where it collides with the boundary current, before being sucked up into the large anticyclonic eddy.

Scales of motion
As an example let us consider the velocity field u(y) = A sin(y/L+ωt) which has the vorticity ζ (y) = A cos(y/L+ ωt)/L.
The formula gives L, the scale of the velocity field.The following analysis is based on this formula.For an understanding of the dynamics it is essential to determine the spatial scales of the turbulent motion.We consider two key quantities.The first is twice the time-averaged kinetic energy (per unit mass) divided by the time-averaged enstrophy (square of vorticity): This quantity is shown in Figs. 3 and 4. In 3-D turbulence it is the Taylor-scale divided by √ 5 (see Frisch (1995)).This length scale characterizes the size of the velocity gradients.
The second length scale is the time-averaged enstrophy divided by the time-averaged palinstrophy (square of the vorticity gradient): This quantity is shown in Fig. 4. It is characteristic of the viscous dissipation length-scale in the enstrophy cascade  Eq. 11) and the viscous dissipation length-scale λ 2 (dashed lines, Eq. 12) at latitude y = +1500 km for the most turbulent experiments EXP300 and EXP400.(Bofetta and Ecke, 2012), the smallest scales in the vortical dynamics.The separation between the two scales gives an idea of the scale range over which turbulence is active.These scales are instructive in a turbulent environment but in the boundary layer dominated by viscosity their significance is limited.At the boundary λ 1 = 0 as the energy vanishes, which does not mean that we have infinitely small scales there.At high viscosity the smallest scale is given by the Munk scale δ M even through the analytic solutions for the laminar Munk layer are (with x = √ 3x/(2δ M )) 2 and ( 13) which oscillate between zero and infinity.This shows that the above scales are not useful for analyzing a time-independent flow.Note that traces of these oscillations remain in the lower viscosity experiments, as can be seen in Figs. 3 and 4.
Figure 3 shows the spatial distribution of the Taylor scale in the highest Reynolds number experiment EXP300.A striking feature is the wide extension of the small-scale values into the interior of the domain in both cases, the feeble variation within this domain and the sudden jump to high values clearly defined boundary as can be seen in Figs. 3 and 4. A clear plateau at around a scale of 60 km which extends over 1000 km into the interior of the domain is observed.We call the area of the plateau the extended boundary layer (EBL).The zonal extend of the extended boundary layer increases during the first part of each experiment but then stabilizes.The scale of 60 km is easily explained by the eddy size of 400 km ≈ 2π 60 km. Figure 4 shows that the width of the extended boundary layer increases with decreasing viscosity.The dissipation length scale λ 2 is smallest near the bound-ary and increases slowly there after, approaching the Taylor scale.When λ 2 reaches the eddy scale λ 1 , the velocity gradients are dissipated and turbulence disappears.A sharp boundary between turbulent areas and a laminar environment is observed in many instances when turbulence arises from a single process such as turbulent jets, planetary boundary layers, gravity currents and stratified turbulence.The behavior of both scales, λ 1 being constant and λ 2 increasing by barely a factor of 2 through the extended boundary layer, shows that grid refinement near the boundary might be useful in laminar, low Reynolds number simulations, but is not adapted for the fully turbulent case where small-scale structures dominate throughout the extended boundary layer.The zonal extension of the extended boundary layer increases with decreasing viscosity as shown in Figs. 4 and 6 and quantified in Sect.4.5.A striking feature is that, although the zonal extension of the extended boundary layer depends on viscosity, the scales within it appear almost independent of it, once the viscosity is low enough to allow for turbulent motion.Turbulent motion in the extended boundary layer is likely to include the range of scales from λ 1 down to λ 2 .
It is important to notice that in our calculations λ 2 is always more than 5 times the grid size showing that the dynamics is numerically well resolved in our calculations.

Multiple boundary layers
The vorticity balance in the laminar, time-independent, boundary layer is described in Sect.4.2.In the unstable boundary layer the vorticity balance changes.When time averaging (denoted by .and the fluctuations by a prime: a = a + a ) is applied to Eq. ( 7), it becomes In a statistically stationary state, the time average of an integration of the advection of vorticity over a closed basin vanishes and the integral balance is between the forcing (righthand side of Eq. 14) and the viscous vorticity flux through the boundary (last term on the left-hand side of Eq. 14).Within the basin the advection of vorticity can connect the (basinwide) source to the sink.The different terms on the left-hand side of Eq. ( 14) correspond to the relative vorticity advection (RVA; terms 1 and 2), turbulent relative vorticity advection (TRVA; terms 3 and 4), planetary vorticity advection (PVA; term 5), stretching (STR; term 6) and friction (FRIC; term 7).The stretching term is negligible and does not contribute significantly to the vorticity balance (see Fig. 5).For higher viscosity (ν ≥ 1000 m 2 s −1 ) the local vorticity balance in the boundary layer is, to leading order, between the planetary vorticity advection (term 5) and the vorticity dissipation (term 7), leading to a Munk layer as discussed in Sect.4.2.When the viscosity is reduced, the relative vorticity advection term and its turbulent part play an increasing  14) are plotted for the lowest viscosity experiment EXP300 at y = +750 km.The different terms of Eq.( 14) plotted correspond to the relative vorticity advection (RVA; terms 1 and 2), turbulent relative vorticity advection (TRVA; terms 3 and 4), planetary vorticity advection (PVA; term 5), friction (FRIC; term 7) and S comprises forcing, stretching and residual time dependence.role in the vorticity balance.The advection of relative vorticity spatially connects the transport of planetary vorticity and the viscous dissipation and both can exhibit a different zonal length scale.This is clearly visible in Fig. 5: the friction dominates in a narrow region near the boundary, whereas the planetary vorticity advection extends further from the boundary.We call the area of the viscous dissipation the VSL while we choose the expression "advective boundary layer" (ABL) for the wider area of large average meridional velocity.The thickness of the former is denoted by δ ν while the thickness of the latter is given by the symbol δ A .In the Munk-layer theory they both coincide δ ν = δ A = δ M .According to the shape of the different terms in Eq. ( 14) (shown in the Fig. 5), we estimate the thickness of the VSL by the distance from the boundary at which the absolute value of the Laplacian of the average vorticity has reduced to a third of its maximal value.The same criterion was applied to the average meridional velocity to obtain δ A .Results for the corresponding boundary layer scales at different latitudes as a function of viscosity are assembled in Fig. 6.For the VSL results show that its thickness drops well below the Munk-scale for the lower viscosities, while the thickness of the advective boundary layer is always above.
The scaling of the advective boundary layer thickness δ A shows a slight increase with decreasing viscosity (see Fig. 6) and a possible saturation around 200 km.

Estimation of the eddy viscosity via the Munk formula
When turbulence is present, the shape of the time-averaged meridional velocity still somehow resembles the Munk-layer solution with the meridional velocity vanishing at a dis- 3)δ M (zero of Eq. 9).The meridional gradient in layer thickness (s) imposed by the large-scale circulation adds a topographic β topo = −f s/H to the planetary value.Its value depends only weakly on the viscosity.When the effective β-term, composed of the planetary and topographic part, is constant, the Munk-layer scale is proportional to the cubic root of the (eddy) viscosity and so is x 0 .The idea is now to calculate an eddy viscosity ν eddy based on x 0 .To this end we measure the value x 0 in an experiment with high viscosity ν stat = 1000 m 2 s −1 that has a time-independent dynamics and compare it to the value obtained from the average of a turbulent experiment at the same latitude.The eddy viscosity can then be obtained by using the proportionality: A clear scaling for ν eddy = ν eddy −ν as a function of the zonal maximum of the root mean squared (r.m.s.) velocity fluctuations u r.m.s. is shown in Fig. 7 at latitudes y = +1500 km.The scatter plot is well fitted by an affine regression line of equation which means that whatever the forcing and the viscosity, there is a correlation between the eddy viscosity and the fluctuating velocity.The correlation of the best fit linear regression is R = 0.97.The finding that for small values of u r.m.s.there is no turbulent contribution to the eddy viscosity is explained by the fact that the small perturbations have a wavelike structure which does not lead to turbulent fluxes.The simplest way to estimate an eddy viscosity proposed by Prandtl (1925)'s Mischungsweg (mixing length) λ and the fluctuating velocity u r.m.s. is ν eddy = ν eddy − ν = αλ 1 u r.m.s.
(17) The results of the nonlinear experiments confirm this proportionality.For our data and λ 1 = L eddy /(2π ) = 60 km calculated previously, we obtain α ≈ 0.1.If we suppose that the eddy viscosity is due to the anticyclones this value of α is within the range proposed by Smagorinsky (1993).The values of λ 1 and u r.m.s.can not be obtained from external parameters but are a result from the numerical experiment.In concrete cases, they can often be obtained from observation or fine-resolution numerical simulations.
Using α = 0.1 and the typical values for the SC of L eddy = 400 km and u r.m.s.= 1 m s −1 leads to ν eddy ≈ 6000 m 2 s −1 and a δ Munk ≈ 70 km.A consequence of this is that even a non-eddy permitting ocean model should have a grid size not exceeding 50 km to capture the boundary layer dynamics and the associated meridional heat transport at least in an average sense and no value of the eddy viscosity larger than 6000 m 2 s −1 should be used.
This pragmatic approach leads to a viscosity and a boundary layer thickness that compares well to average values in the turbulent boundary current.This approach is of course questionable as the eddy size is larger than the mean current, that is, the scale separation is smaller than unity and the eddy viscosity approach asks for large-scale separations.This problematic was already noticed by Charney (1955) who states: "In order to account for the observed width of the current, Munk was forced to postulate an eddy viscosity so large that the eddy sizes were themselves comparable to the width." We have estimated the eddy viscosity based on the average meridional velocity and have shown that it can be connected via Prandtl (1925)'s formula to the velocity fluctuations.

Discussion and conclusions
The western boundary is a turbulent region with interacting eddies, bursts and dipoles and frequent velocity inversions.A laminar boundary layer structure can be recovered in an average sense.The turbulent dynamics leads to a split up of the boundary layer into three layers: a VSL, an advective boundary layer and an extended boundary layer.The thickness of the VSL increases with viscosity, the thickness of the extended region (EBL) decreases and the advective region (ABL) stays essentially unchanged, once the viscosity drops below values that allow for turbulent motion.
For the lower viscosity experiments, we identified a sequence in the evolution of the dynamics of the coherent structures: anticyclones are generated by instability, during their northward migration they intermittently detach parts of the VSL containing strong positive vorticity called bursts.These bursts pair with negative vorticity from within the anticyclones and form dipoles which then travel ballistically (at almost constant velocity) over distances of several eddy diameters.Two other trajectories are possible as said in the Sect.4.3.3.Large anticyclonic eddies creating bursts and a strong dipole are also clearly visible in numerical simulation of Spall (2014, his Fig. 9).In observations (called "flanking cyclones" by Beal and Donohue (2013)) and a fineresolution ocean general circulation model (Akuetevi et al., 2015) bursts are seen to lead to the substantial upwelling of cold and nutrient rich water masses from the deep.The dipoles transport these water masses offshore, leading to an increased biological production several hundreds of kilometers from the coast (Kawamiya and Oschlies, 2003;Wirth et al., 2001).The above is an example of how meso-and sub-meso-scale activity can increase biological activity.
We showed that the turbulent eddy dynamics is the natural state of the high Reynolds number low latitude westernboundary current, when the stabilizing effect of inertial effects is absent.

Conclusions concerning numerical simulation of turbulent boundary layers
It is the thickness of the VSL that imposes the spatial resolution of a numerical model.The thickness of the turbulent VSL decreases faster with decreasing viscosity than the prominent one-third scaling from Munk-layer theory, in all our experiments performed and at all latitudes considered as demonstrated in Fig. 6.The laminar Munk-layer theory is however commonly used to determine the (hyper-)viscosity for a given spatial resolution in today's simulations of ocean dynamics.The results presented here prove that for the turbulent boundary layer the choice of spatial resolution based on the Munk-layer theory is far from being sufficient.From Fig. 6 it is clear that the difference between the thickness of the extended boundary layer and the VSL widens with increasing the Reynolds number.The difference is a measure-C.Q. C. Akuetevi and A. Wirth: Dynamics of turbulent western-boundary currents at low latitude ment of the complexity of the numerical calculations of low latitude turbulent WBCs as the finest scale δ ν of the VSL has to be resolved throughout the extended boundary layer δ ext in both horizontal directions.This shows that grid refinement near the boundary that is using a finer grid closer to the boundary than further away, has no place in simulations of the turbulent boundary layer as (i) the structures are almost isotropic and (ii) the small scales extend far from the boundary.

Conclusions concerning the parameterization of the turbulent boundary layers
One of the major challenges in the numerical simulation of the ocean dynamics is to parameterize the effect of the smallscale dynamics not explicitly resolved on the explicitly resolved large-scale flow.
Our determination of the eddy viscosity in Sect.4.6 via the Munk formula is a parameterization as we related eddy viscosity to the maximum fluctuating velocity.These show that for the lowest viscosities, δ A saturates at a value corresponding to ν ≈ 6000 m 2 s −1 .A consequence of the results presented above (Sect.4.6) is that choosing viscosity values lower than ν ≈ 6000 m 2 s −1 but above the threshold for fully turbulent boundary layers ν ≈ 300 m 2 s −1 leads to an unrealistically thin average boundary layer, worsening the representation of the advective boundary layer dynamics.In numerical simulations of the boundary layer dynamics, one should either simulate the turbulent dynamics or parameterize it.In other words, our findings discussed above suggest that one can either use fine resolution and a low viscosity (ν <≈ 300 m 2 s −1 ) to simulate the turbulent boundary, or one can use coarse resolution and a high viscosity (ν ≈ 6000 m 2 s −1 ) and recover the time-averaged boundary layer dynamics.Using viscosities in the interval 300 m 2 s −1 < ν < 6000 m 2 s −1 or values greater than 6000 m 2 s −1 leads to a wrong time-averaged boundary-layer dynamics.
In our simulations we varied the viscosity parameter by roughly a factor of 3. The corresponding necessary spatial resolutions vary from those of today's coarse-resolution climate models down to those of fine-resolution regional models.Our calculations suggest that even lower viscosity values lead to smaller boundary layer scales and higher velocities.At smaller scales the hydrostatic approximation, on which the shallow-water equations are based is no longer valid as the dynamics becomes truly three-dimensional.Higher velocities lead to Froude numbers exceeding unity, that is, the fluid velocity is higher than the speed of the gravity waves.In this case, hydraulic jumps occur and the flow becomes fully three-dimensional, such phenomena can not be explicitly resolved by the two-dimensional shallow-water equations.In Fox-Kemper and Pedlosky (2004) and Fox-Kemper (2004) these problems are bypassed by using a constant depth model, where Froude number vanishes and by increasing the viscosity in the vicinity of the boundary.
We did not consider the more involved behavior of hyperdissipation operators (hyper-viscosity, powers of the Laplacian) which ask for boundary conditions for derivatives of the velocity field and which lead to thermalization at small scales of the turbulent dynamics as explained by Frisch et al. (2008).

Figure 1 .
Figure 1.Instantaneous velocity arrows (plotted every 20 × 50 grid points in the x and y directions) superimposed on layer thickness variation (η) at time t = 2000 days, when the dynamics has converged to a stationary state, for the high viscosity experiment EXP1000.

Figure 2 .
Figure2.Sequence of potential vorticity (m −1 s −1 ) snapshots showing bursts and its subsequent development into a dipole for the lowest viscosity experiment (EXP300).From bottom to top, the snapshots were taken at t = 4023, 4026 and 4034 days.

Figure 3 .
Figure 3. Taylor scale λ 1 (Eq.11) for the lowest viscosity experiment EXP300.Note that the color bar stops at 100 km to emphasize the behavior in the extended boundary layer (region of λ 1 ∼ 60 km).

Figure 6 .
Figure 6.Thickness of the viscous sub-layer (VSL), the advective boundary layer (ABL) and the extended boundary layer (EBL) for all the experiments at different latitudes y.

Figure 7 .
Figure 7. Scatter plot diagram of eddy viscosity ν eddy = ν eddy − ν computed from the data using the Munk formula approach of Eq. (15), as function of the maximum fluctuating velocity for all the experiments at latitudes y = 250, 750, 1500, 1750 and 2000 km and the red line is the best fit affine regression line.