Dynamics of Upwelling and Downwelling in a Channel Basin of the Baltic Sea

2.1 Basic model configuration

We use the Massachusetts Institute of Technology general circulation model (MITgcm; http://mitgcm.org), an open-source FORTRAN solver for Primitive Equations (the Boussinesq form of the Navier–Stokes equations). MITgcm employs a finite-volume discretization rendered on a horizontal Arakawa C-grid, with vertical z-levels employing partial cells (Adcroft et al., 2004; Marshall et al., 1997). The MITgcm has been used previously for idealized regional circulation studies including upwelling (Dawe and Allen, 2010; Howatt and Allen, 2013; Thomsen et al., 2021).

In our study, we consider a high-resolution, idealized, hydrostatic setup of the WGB. A Cartesian coordinate system (${\displaystyle x,y,z}$) has horizontal axes aligned with the cross-shore direction (x; width ${\displaystyle L_x=130 \mathrm{k}\mathrm{m}}$) and the alongshore direction (y; width ${\displaystyle L_y=100 \mathrm{k}\mathrm{m}}$). The three-dimensional velocity vector in this coordinate system is ${\displaystyle \overrightarrow{u}=(u,v,w)}$. Periodic boundary conditions are applied in the alongshore direction, and closed boundaries in the cross-shore direction at the Swedish east coast and Gotland (Figure 1e). On these scales, the latitudinal variation of the Coriolis parameter is negligible and we apply a f-plane for 58° N (${\displaystyle f=1.23\cdot {10}^{-4} {\mathrm{s}}^{-1}}$), corresponding to one inertial period of ${\displaystyle T_i=2\pi /f\approx 14 \mathrm{h}}$.

Our idealized model configuration is focused on the dynamics of sub-inertial motions (inertial waves are merely transient phenomena given the applied forcing and wind waves are not resolved here). We apply a fully implicit time stepping (https://mitgcm.readthedocs.io/en/latest/algorithm/crank-nicol.html), thus neglecting the time variations of free surface when computing the horizontal transports in the continuity equation and in tracer and momentum advection terms. This method is unconditionally stable and avoids computational effort related to the implementation of the non-linear free-surface formulation, but it damps the fast gravity waves. The horizontal resolution is ${\displaystyle \mathrm{\Delta }x=1 \mathrm{k}\mathrm{m}}$ and the vertical resolution is ${\displaystyle \mathrm{\Delta }z=1 \mathrm{m}}$. In the horizontal, the viscosity and diffusivity are parameterized with a Laplacian operator with a viscosity of 10^–5 m² s^–1 and temperature and salinity diffusivities are kept at lowest possible, i.e., default molecular diffusivities of 10^–7 m² s^–1 and 10^–9 m² s^–1, respectively, to avoid excessive smoothing of fronts (so in practice the horizontal diffusivities are negligibly small). In the vertical, a non-local K-Profile Parameterization (KPP) scheme (Large et al., 1994) is used to parameterize vertical mixing: the background value of vertical viscosity and tracer diffusivity (both equal to 10^–6 m² s^–1) are modified in the surface boundary layer by employing the Monin-Obukhov similarity theory, accounting for non-locality of the turbulent fluxes, which depend on the boundary layer depth and surface fluxes, in addition to local properties and gradients. In the interior, the vertical viscosity and diffusivity depend on the velocity shear and stratification to represent shear instability and internal wave breaking. The vertical dissipation is thus a function of the system state and will vary between the configurations. The use of KPP circumvents the necessity of setting high values of vertical viscosity and diffusivity coefficients everywhere in the domain to ensure the numerical stability in the more vertically diffusive boundary layers. It also allows us to quantify vertical turbulent transport associated with the resolved dynamics (see Figure 10a–f and section 3e.).

The bottom stress is given by the linear friction law ${\displaystyle {\tau }_b^{(x,y)}=r{\rho }_0(u_b,v_b)}$, where u_b and v_b are horizontal bottom velocity components and the bottom drag coefficient is r = 10^–4 m s^–1 as used in Brink (2016), and with ${\rho }_0=998\hspace{0.17em}\text{kg}\hspace{0.17em}{\text{m}}^{-3}$ as the model reference water density. Heat fluxes are set to zero.

2.2 Wind forcing

The model is forced with an alongshore wind stress ${\displaystyle {\tau }_s^y}$ variable in time following:

This corresponds to a single wind pulse of the duration of 12 inertial periods where the period of maximum wind amplitude, ${\displaystyle \mathrm{\Delta }t}$, is equal to 8.45 inertial periods (five days), which is approximately the center of distribution of upwelling event duration (Figure 1d) and long enough for the instabilities to develop. The wind stress is initialized with ${\displaystyle t_r=1.7T_i}$ (1 day), before reaching the constant amplitude ${\displaystyle {\tau }_0=\pm 0.1}$ N m² (negative for upwelling on the Swedish east coast side and for downwelling on the Gotland side), corresponding to gentle breeze conditions on the Beaufort scale (7–10 m s^–1), typical to WGB (Figure 1c).

The wind is spatially uniform, which is a reasonable assumption based on wind observations in WGB (not shown), see also Lehmann and Myrberg (2008). The shape of the wind stress function ${\displaystyle {\tau }_s^y}$ (for positive values) is superimposed on Figure 1d. The wind impulse is defined as:

and quantifies the input of the wind momentum due to a given wind stress and duration of the wind event. The total time for model experiments and analysis is ${\displaystyle T=1}$ month = 28 days, i.e., longer than the transient wind time scale in the region (Figure 1c–d).

2.3 Bathymetry

The bathymetry setting is a channel with a sloping bottom on the west (Swedish coast) side and a vertical wall on the eastern (Gotland) side (Figure 1e):

where ${\displaystyle H_0=-160 \mathrm{m}}$ is the basin depth and the slope parameter α is equal to the average slope south of the Landsort Deep, ${\displaystyle \alpha =0.002}$. The bending point of the slope is at ${\displaystyle x_b=80 \mathrm{k}\mathrm{m}}$ where the bottom becomes flat until the Gotland coast at ${\displaystyle L_x=130 \mathrm{k}\mathrm{m}}$.

Two model runs using this bathymetric setup with alongshore winds of opposite directions result in four flow regimes: UPSLOPE (upwelling on the slope side) with DOWNWALL (downwelling on the wall side), and DOWNSLOPE (downwelling on the slope side) with UPWALL (upwelling on the wall side). The regimes are demarcated by the slope bending point at ${\displaystyle x=80 \mathrm{k}\mathrm{m}}$.

We complement our study with the simulation LANDSORT, including a sinusoidal alongshore bathymetric variation (alongshore wavenumber ${\displaystyle l=1}$) mimicking the basin-scale feature in the WGB (a trough at the Landsort Deep). We do not consider the deepest part of the Landsort Deep below –160 m, but since the density gradients are very weak below –100 m (Figure 1f and Section 2.4), we assume that it does not impact the dynamics of the coastal upwelling at the basin scale. We also run a simulation with no slope vertical walls on both sides) to compare with the analytical two-layer model. We evaluated the sensitivity to horizontal resolution by running a higher resolution simulation (${\displaystyle \mathrm{\Delta }x=500 \mathrm{m}}$) and also to the bottom drag formulation (by running a simulation with quadratic bottom drag). Both sensitivity runs have been initialized with winter stratification.

2.4 Stratification

The simulations are initialized with continuous stratification profiles representative of summer (August) and winter (March) conditions taken from the Copernicus Marine Service (CMEMS) reanalysis and forecast [https://data.marine.copernicus.eu/product/BALTICSEA_ANALYSISFORECAST_PHY_003_006/. Product DOI: https://doi.org/10.48670/moi-00010. The reanalysis profiles were averaged over a transect parallel to the coast and idealized by the following analytical expressions:

where S₀ and T₀ are, respectively, the salinity and temperature at the surface, while S₁ and T₁ represent, respectively, the bottom values and z₀ and z₁ are fitting parameters. The idealized profiles smooth over noise and neglect the deep temperature minimum in the summer profiles, which is of minor importance for the dynamical processes studied here, but it might be important for more realistic studies of WGB circulation in the future. The resulting profiles used to initialize the model are shown in Figure 1f. During the month of August, the summer thermocline is well developed at 20 m with surface waters up to 10°C warmer than deeper layers. Example upwelling and downwelling expressions on the summer SST field are shown in Figure 1a–b. During winter (March), the upper thermocline is reversed (colder water at the surface) and weaker (3°C difference), with a more distinct mixed layer.

The salinity and halocline do not show a pronounced seasonal cycle, as it is set by the estuarine circulation in the Baltic Sea dominated by interannual variability (Kullenberg, 1981; Lehmann et al., 2022). Figure 1f includes the corresponding profiles of the squared buoyancy (Brunt–Väisäla) frequency computed from the potential density referenced to the surface:

This displays one, and then two, distinct pycnoclines in summer and winter, respectively: the summer thermocline forms a distinct summer pycnocline above the main one, strenghtening the stratification; while the winter thermocline is located only slightly above the main pycnocline (halocline) and has a weakening effect on the stratification.

2.5 Scales of motion

For ${\displaystyle H_0=-160 \mathrm{m}}$ (flat bottom), the barotropic Rossby radius of deformation is ${\displaystyle R_{bt}=c_{bt}/f=322 \mathrm{k}\mathrm{m}}$ (the gravity wave speed ${\displaystyle c_{bt}=\sqrt{g|H_0|}=40}$ ${\displaystyle m{s}^{-1}}$). We can estimate the baroclinic Rossby radius of deformation, R_bc, from the two-layer model which is representative for the Baltic winter stratification:

where ${\displaystyle ϵ=({\rho }_1-{\rho }_2)/{\rho }_1}$, and ${\displaystyle h_1,{\rho }_1}$ are the thickness and density of the layer above a pycnocline, and ${\displaystyle h_2,{\rho }_2}$ are the thickness and density of the layer below, so that ${\displaystyle H_0+h_1+h_2=0}$. The value of R_bc given by equation 6 is comparable with the estimate of the baroclinic Rossby radius of deformation using the first eigenvalue of the vertical normal mode decomposition Fennel et al. (1991) and Zhurbas et al. (2006), which yields about 8 km for summer and 7 km for winter stratification. Throughout the manuscript, we only use R_bc as an estimate of the baroclinic Rossby radius. The baroclinic gravity wave speed estimate yields ${\displaystyle c_{bc}=R_{bc}\cdot f \approx 1 m{s}^{-1}}$.

A non-dimensional number quantifying the effects of sloping topography and stratification is the topographic Burger number (Lentz and Chapman, 2004):

This number gives an estimate of the proportion of the volume transport confined to the bottom boundary layer (Jacox and Edwards, 2011): lower values indicate that the surface wind stress is compensated by the bottom stress, indicating strong cross-shore volume transport in the boundary layer. As the Burger number increases, the contribution from cross-shelf momentum divergences increases, together with the volume transport in the interior. In our configuration, given the varying stratification, the Burger number can be estimated to be in the interval ${\displaystyle Bu=0.23-0.25}$ for all simulations.

To interpret the results related to sub-inertial dynamics, the time is non-dimensionalized by the local inertial period:

The cross-shelf stream function ${\displaystyle \mathrm{\Phi }}$ is defined in such a way that positive values describe an anti-clockwise overturning cell in an x–z plane and is non-dimensionalized as:

together with the cross-shore velocity field u_w:

where

is an estimate of the Ekman layer thickness calculated using the value for the vertical viscosity computed by KPP scheme in the surface or bottom layer. The non-dimensionalization of Equations 9–11 sets the reference for the interpretation of the results.

2.6 Energy and instability diagnostics

We use a set of diagnostics introduced by Orlanski and Cox (1973) and applied in Brink (2016) and Durski and Allen (2005). Starting by defining the alongshore mean ${\displaystyle \{\varphi \}=\{\varphi \}(x,z,t)}$ for a generic three-dimensional quantity ${\displaystyle \varphi =\varphi (x,y,z,t)}$:

together with its fluctuation from the alongshore mean:

with the spatially uniform wind forcing and bathymetry defined in equation 3, the deviations from the alongshore mean diagnose the instability processes in the domain.

The kinetic energy associated with the alongshore mean flow (MKE) is:

and the eddy kinetic energy (or EKE) associated with its perturbation:

where (${\displaystyle W_1,W_2}$) are set respectively to (0,80) km and (80,130) km for the flow regimes on the sloping (subscript 1) and wall (subscript 2) sides separated by the slope bending point, and A is the averaging area.

Applying the averaging procedure to the equations of motion, the time variation of the total (summed slope- and wall side contributions defined in equation 15) eddy kinetic energy can be expressed as (Orlanski and Cox, 1973):

Note also the cross-shelf energy flux in equation 16 must sum to zero. In the following, we focus on the energy conversion terms diagnosing instability. The conversion from potential to eddy kinetic energy (cpe) can be split to the sum of a term associated with geostrophic adjustment ${\displaystyle \{\rho w\}}$ and another term associated with baroclinic instability:

The conversion from mean to eddy kinetic energy (cmke) can be split into the component associated with the horizontal shear and another component with the vertical shear. The former is related to barotropic instability, while the latter is related to Kelvin-Helmholtz instability:

The onset of instability is estimated as the time when either of the energy diagnostics becomes different from zero. We estimate the instability growth rate as:

where ${\displaystyle t_{1/2}}$ and ${\displaystyle t_{1/4}}$ are the time at which the EKE reaches half and one-quarter of its maximum (Brink, 2016). To gauge the mechanism allowing for the development of instabilities we complement our analysis with Ertel potential vorticity (PV) computed following Morel et al. (2019), their Equations 11–14:

where the relative vorticity:

and ${\displaystyle \overrightarrow{f}=(0,0,f)}$; together with thermal wind and the associated horizontal density gradients.

3.1 Analytical linear solution for coastal jets in a channel basin

To determine the major features of upwelling and downwelling in an idealized WGB channel, and to assess our numerical model simulations, we derive an analytical solution for two-layer inviscid linear Ekman flow in a channel following the approach of Csanady (1982) however, we apply a wall boundary on both sides of the domain instead of the semi-infinite (open ocean) boundary used in previous works. For an alongshore wind forcing switched at ${\displaystyle t=0}$ and staying constant thereafter, the solution valid for small displacements of the interfaces (the free surface and the pycnocline) gives the following expression for a non-oscillatory component of the free surface displacement ${\displaystyle \eta }$ that grows linearly in time and is negative (positive) at the coast for upwelling (downwelling):

The constants ${\displaystyle (A_1,A_2)}$ and ${\displaystyle (B_1,B_2)}$ are given by the:

The free surface displacement is decomposed into barotropic (${\displaystyle {\eta }_{bt}}$) and baroclinic (${\displaystyle {\eta }_{bc}}$) modes, with the latter giving a minor contribution ${\displaystyle ϵ=}$${\displaystyle ({\rho }_1-{\rho }_2)/{\rho }_1}$ ${\displaystyle \approx 0.002}$ to the free surface displacement. The integrated velocity in the upper layer ${\displaystyle (U,V)}$ is also decomposed into barotropic and baroclinic components (Csanady, 1982):

The amplitude of the interface displacement ${\displaystyle {\eta }^{\prime }}$ is also evolving linearly in time and the transport in the lower layer ${\displaystyle (U^{\prime },V^{\prime })}$ yields (positive for upwelling):

Note that cross-shore velocities are independent of time while alongshore velocities grow linearly in time and form coastal jets trapped to the coast within a distance determined by channel width L_x and stratification (${\displaystyle R_{bt},R_{bc}}$, Equations 22–24), anti-symmetrically about the zero-crossing point of the hyperbolical function at approximately ${\displaystyle L_x/2}$. The closed boundary conditions in the cross-shore directions within the channel require that (${\displaystyle U_1+U_2=0}$ at ${\displaystyle x=0}$ and ${\displaystyle x=L_x}$).

The above solution describes wind-driven flows in a two-layer channel with vertical walls (no slope, no friction effects) and is valid only for the time period before the interface displacements become large and nonlinear effects kick in. Analyzing differences between the analytical and numerical solutions allows us to assess the role of friction, slope, the seasonal variations in stratification and the nonlinear effects. Figure 2a–b show the free surface profile (related to the barotropic component of the coastal jet) in the winter simulation with walls on both sides as well as simulations with slope on one side (scaled by the difference in domain volumes) for summer (three-layer-like continuous) and winter (two-layer-like continuous) stratification, and for both, constant wind stress and a wind impulse stress (see Equation 1) applied in our numerical simulations (only the simulations with negative wind direction are shown).

Figure 2

After a spinup of t_i = 2, the free surface displacement in the winter wall-wall simulation follows the analytical solution nearly exactly during the active phase. The free surface displacement over the slope is larger by 50% compared to the wall side due to higher wind momentum transfer per unit depth and the zero crossing is moved by about 10 km towards the coast. As a result, the free surface gradient (${\displaystyle \mathrm{\partial }\eta /\mathrm{\partial }x}$) is stronger on the slope side, leading to a stronger barotropic flow there (over the course of time, this excess of momentum transfer in simulations is redistributed in the domain to the wall side leading to a stronger flow than the analytical solution there as well).

Notably, the distance over which the coastal flows on both sides are confined to the coast is set primarily by the channel geometry (width and slope). Seasonal variations of stratification (resulting in only a small range of R_bc) have a negligible effect. On the wall side, the free surface displacement continues to grow linearly and anti-symmetrically for upwelling and downwelling reaching ~0.45 m at the end of the active phase in agreement with the analytical solution. On the slope side, it grows at a higher rate albeit in a sub-linear manner (slower in the beginning as it takes time to develop boundary layer circulation). A slow relaxation during the passive phase follows reducing the displacement by half after one month.

Because the simulated stratification is continuous and the definition of interfaces is thus arbitrary, we can only gauge a qualitative comparison with the analytical solution. The model thermocline (halocline) are defined as isotherms (isohalines) corresponding to the maximum vertical gradient of the initial temperature (salinity) profiles in Figure 1f. Their maximum displacements also grow linearly in time during the active phase in agreement with the two-layer theory (Equation 27) but only for a period before the outcrop of the summer thermocline for UPWALL-S and UPSLOPE-S (about t_i = 4) preceding onset of the instability (see Section 3.3), and similarly before shoaling of the winter thermocline for UPWALL-W after t_i = 10 (for UPSLOPE-W the winter thermocline stays below 10 m, not shown). To further quantify this, we plot the simulated interface displacements scaled with the analytical solution (Equation 27) in function of the wind impulse in Figure 2d–e. The agreement between the simulated interface displacements and the analytical solution (derived for a channel with two walls) is within ±50% and is, as expected, better on the wall side (half of the values fall within ±25%). On the slope side, the analytical solution tends to be underestimated, in particular for the summer halocline. Note that there is uncertainty related to the definition of the interface and continuous stratification in the simulations). For stronger wind impulses, the interface outcrops and instabilities make the interface displacements deviate from analytical solution for the UPWALL (by a factor of >1.25) where instabilities set on early during the active phase (see Section 3.3). The scaling with wind impulse suggests that although neither the winter thermocline nor halocline outcrop in simulations presented here, they would do so if the wind impulse (strength and/or duration) would be increased by about 15% and 30% respectively (these however correspond to quite rare wind conditions in the WGB, see Figure 1c).

3.2 Coastal jets in a channel basin

The evolution of the flow for positive wind stress (UPSLOPE/DOWNWALL) and negative wind stress (DOWNSLOPE/UPWALL), summer and winter, is further analyzed in Figure 3 for selected time instances t_i = 8, 14, 22. The first time instance is representative of the flow when the coastal jet is fully developed in the active phase, the second and the third represent, respectively, early and late stages of the passive phase when instabilities are present. Shown are the density anomaly, cross-shore, alongshore velocity, and streamfunction.

Figure 3

In the first 20 km from the coast, the outcropped summer thermocline forms a surface mixed layer of about 20 m. The density difference between the upper and lower layer due to summer thermocline and permanent halocline reaches 7 kg m^–3. The outcropping summer thermocline (Section 3.1) reveals upwelling at the surface (Figure 1a–b) enabling statistical quantification and impact studies using SST imagery (Gurova et al., 2013; Lehmann et al., 2012; Norman et al., 2013; Sproson and Sahlèe, 2014).

In winter, the reversed thermocline (warmer water at the surface) lies deeper (55 m) i.e., only a few meters above the halocline marking a deep mixed layer and partly compensating the salinity effect on density in reducing its density difference to 3 kg m^–3. The shoaling of the winter thermocline has a SST signature in our simulations (not shown). Notably, a ‘warm’ winter upwelling has only recently been quantified in the Baltic Sea using in-situ measurements and few available SST pictures in a channel basin of the Gulf of Finland (Suursaar, 2021), though the local conditions there feature higher run-off and presence of sea ice and are thus not exactly comparable to these in the WGB; The largely overlooked ‘warm’ winter upwelling in the Baltic Sea might have important consequences for triggering the early spring algae blooms (e.g. (Beltran-Perez and Waniek, 2022)).

During the active phase (${\displaystyle t_i=8}$), the depth of the model surface Ekman boundary layer marked by the positive cross-shore transport u_w (second row of Figure 3) is about 15 m in summer (a few meters above the summer thermocline) and 20 m in winter (much shallower than the winter mixed layer). These values are consistent with the theoretical prediction for D_e (Equation 11) using mean values of simulated vertical viscosity values averaged over the upper 70 meters (K_V = 0.0015 m² s^–1 in summer and 0.002 m² s^–1 in winter; see Section 2e).

One estimate for the density front location proposed by Austin and Lentz (2002) reads respectively for upwelling and downwelling:

which yields around 20 km after ${\displaystyle t_i=8}$ both summer and winter and both sides of the domain, consistently with our simulations (Figure 3, first and third row) and reanalysis of observations (Figure 1a–b). The front location in the WGB is thus primarily determined by the wind impulse. The location of the frontal jet determines, in turn, the extent of the near-shore recirculation on the slope so-called inner shelf (Austin and Lentz, 2002) for the summer and winter upwelling (Figure 3), v and streamfunction). The inner shelf will in turn affect the circulation and oxygen levels in the nearshore coves and bays (Fredriksson et al., 2024). For DOWNSLOPE, the inner circulation develops during the passive phase related to the gradual offshore migration of the frontal jet (approximately 10 km over the month-long simulation). Kämpf (2019) showed that for stronger and persistent wind events, extreme bed stresses (~0.6 Pa) arising due to downwelling jet on a slope can lead to a shutdown of the cross-shelf circulation, which we do not see in our simulations forced with impulse-like and weaker wind forcing. The fronts detach in all cases except for DOWNWALL, with implications for jet instability (see below).

The Ekman flows on both sides of the channel are connected by closed streamlines (Figure 3, last row), with a stronger surface-intensified flow and an inner shelf expression on the slope side. On the wall side, the Ekman transport in the surface boundary layer is compensated by upwelling-downwelling in the interior (below the surface layer). On the sloping bottom, in a low Burger number regime characterizing our simulations (${\displaystyle Bu\sim 0.25}$), most of the compensating cross-shore transport takes place in the bottom boundary layer (Lentz and Chapman, 2004). This result is valid within the entire range of slope values across the WGB estimated from the bathymetry (${\displaystyle Bu\sim 0.20-0.50}$) and could be used in models for the nutrient fluxes (Jacox and Edwards, 2011).

The bottom boundary layer depth is about 6 m for UPSLOPE and twice as thick for DOWNSLOPE. These values are consistent with the theoretical prediction for Ekman boundary layer thickness (Equation 11) using local interior values of the vertical viscosity output from KPP, which are also higher for DOWNSLOPE due to down-welled fresher (and warmer in the summer) water weakening the stratification.

The deeper bottom boundary layer in DOWNSLOPE compared to UPSLOPE is due to by the downward bending of isopycnals implying weaker hydrostatic stability (Brink, 2016). The Ekman flow on the slope generates horizontal density gradients resulting in thermal wind shear that is expected to eventually bring the bottom alongshore velocity to zero neutralizing the bottom stress and Ekman cross-shore transport, so-called buoyancy arrest (Lentz and Chapman, 2004). The timescale for the buoyancy arrest depends on the Burger number and the friction and for the downslope flow can be estimated as (Brink, 2016):

where ${\displaystyle d=C_{d}N/f}$, where C_d is the quadratic bottom drag coefficient, which we substituted with our linear drag coefficient r divided by an indicative value for cross-shore velocity ${\displaystyle \tau /(D_E\rho f)\sim 0.06 \mathrm{m}/\mathrm{s}}$) (Arbic and Scott, 2008). For upwelling the buoyancy arrest timescale is one order of magnitude smaller (Brink and Lentz, 2010):

where in both cases we used the mean value of N and noted that the result is sensitive to the value of N (using the bottom values instead of the mean makes the time scales longer of a factor of ~1.8 for downwelling and 1.5 for upwelling). Because of the relatively long scales of the buoyancy arrest that are comparable to, or longer, than the active phase and the onset of instability (see below), we do not detect clear buoyancy arrest in our simulations: the bottom cross-shore Ekman transport persists throughout the active phase (Figure 3) and is subject to relaxation and instabilities in the passive phase.

3.3 Development of instabilities

The summary of instability analysis is given in Figure 4. Time series of the conversion rates for the flow regimes at the slope and wall sides (Section 2.6) are shown in Figure 5a–d. The onset time of baroclinic instability is identified by a non-zero value of cpe (and EKE, not shown).

Figure 4

Figure 5

The conversion of PE to EKE (cpe), associated with the baroclinic instability, dominates the energetics for UP-, and DOWNSLOPE as well as UPWALL, while DOWNWALL does not exhibit instability. The onset of instability is approximately the same for UPSLOPE summer and winter ($t_i^{\mathit{onset}}\approx 3$, corresponding to a wind impulse of ~12000 Pa s, see Figure 5e) the onset is approximately at the same time for UPWALL summer ($t_i^{\mathit{onset}}\approx 4$) and for UPWALL winter the onset is postponed to $t_i^{\mathit{onset}}\approx 5$, corresponding to a wind impulse of ~22 500 Pa s. In all upwelling simulations the onset occurs already during the active phase.

The onset for DOWNSLOPE occurs later during the passive phase ($t_i^{\mathit{onset}}\approx 14$ in summer, $t_i^{\mathit{onset}}\approx 19$ in winter) because of the weaker stratification and wind shear during the active phase. The longer onset makes DOWNSLOPE instability less likely to occur in the WGB than UPSLOPE as the northeasterly winds are less frequent (~10% time) than southwesterly winds (~30% time) while the winds tend to change direction often on short time scales (Figure 1c–d). The onset of instability for upwelling is directly related to the wind impulse during the active phase: the flow gets unstable for ${\displaystyle I_{\tau }>12000}$ Pa s for all flow regimes for the given geometry and stratification (Figure 5e).

The onset of baroclinic instability is associated with the off-shore migration of the coastal jets for all unstable regimes (UPSLOPE, UPWALL and DOWNSLOPE; Figure 3, middle panels) and development of frontal recirculations. The contribution of the baroclinic component to the total flow at the later time during the active phase is shown in Figure 6.

Figure 6

Clearly, while the circulation on the slope is more barotropic (see also Section 3.1), the circulation at the wall is more baroclinic. The baroclinic flows on both sides are coupled to both, thermocline and halocline, in summer (in winter both interfaces nearly coincide). Thus, the presence of the displacement of the summer surface thermocline intensifies the instability (by increasing the ${\displaystyle {\mathrm{\partial }}_x\rho }$), but it is not necessary for the development of instabilities due to the halocline. The maximum of EKE decreases by almost a factor of 5 from WALL to SLOPE and by a factor 2 from summer to winter because of the weaker winter stratification and thermal wind (not shown).

The baroclinic instability on the slope is further assessed using the ratio of the bottom slope to the isopycnal slope (${\displaystyle \mathrm{\Delta }=\frac{\mathrm{\partial }H/\mathrm{\partial }x}{\mathrm{\partial }\rho /\mathrm{\partial }x}}$) as in the Eady problem with a sloping bottom (Isachsen, 2011; Stipa, 2004; Zhurbas et al., 2006). For ${\displaystyle 0<\mathrm{\Delta }<1}$ the flow is unstable with a fastest growing mode of ${\displaystyle R_{bc}{k}_m\approx 1.14}$ (1.61 for the Eady model). For ${\displaystyle \mathrm{\Delta }>1}$ no unstable solutions can be found, for ${\displaystyle \mathrm{\Delta }<0}$ (opposing isopycnal and bottom slopes) the flow is unstable and the growth rate decreases with ${\displaystyle |\mathrm{\Delta }|}$, while there is a strong stabilization of long waves and the short-wave cut-off moves to even shorter waves (see Stipa (2004), their fig. 2 and Isachsen (2011), their fig. 1]. In our simulations, ${\displaystyle \mathrm{\partial }H/\mathrm{\partial }x<0}$. For UPSLOPE, ${\displaystyle \mathrm{\partial }\rho /\mathrm{\partial }x<0}$, so ${\displaystyle \mathrm{\Delta }>0}$ and the flow can be stable or unstable depending primarily on the strength of the upwelling (wind forcing) and stratification; since ${\displaystyle (-\mathrm{\partial }T/\mathrm{\partial }x<0,\mathrm{\partial }S/\mathrm{\partial }x<0)}$ in summer and ${\displaystyle (-\mathrm{\partial }T/\mathrm{\partial }x>0,\mathrm{\partial }S/\mathrm{\partial }x<0)}$ in winter, ${\displaystyle {\mathrm{\Delta }}_{UP}^s<{\mathrm{\Delta }}_{UP}^w}$ implying that upwelling on the slope is more prone to instability in summer than in winter for a given wind magnitude due to the outcropping summer thermocline and the unstable wavelength compares to R_bc, and the estimated growth rates are about 4 day (Figures 4, 5). For DOWNSLOPE, ${\displaystyle \mathrm{\partial }\rho /\mathrm{\partial }x>0}$ so ${\displaystyle \mathrm{\Delta }<0}$ and the system is always unstable, but because of the seasonality in temperature stratification, ${\displaystyle |{\mathrm{\Delta }}_{DOWN}^s|<|{\mathrm{\Delta }}_{DOWN}^w|}$, the growth rate (onset) of instability is larger (faster) for the summer, consistently with our simulations (Figures 4, 5), and the short waves are more unstable.

We zoom into the flow under baroclinic instability on the slope in summer in Figure 7 (UPSLOPE) and Figure 8 (DOWNSLOPE) showing the cross-shore density gradients (thermal wind, panels a) and the Ertel PV field (Equation 20, panel b) a couple of inertial periods after the onset of instability.

Figure 7

Figure 8

The horizontal density gradient (thermal wind) is negative (positive) for UPSLOPE vs. DOWNSLOPE and since the Ekman circulation is reversed, the PV changes sign (a necessary condition for the baroclinic instability) from positive to negative within a broader region along the slope bounded by the jet (front) maximum and the halocline (a region of baroclinic flow, see Figure 6). The maximum cpe is localized at the jet location due to the strong thermal wind there (sufficient condition) in both cases, however, while for the UPSLOPE this maximum is above the boundary layer, while for DOWNSLOPE it is within. In the result, the instabilities marked by the cross-shore velocity changes (panels c) develop above the bottom boundary layer for UPSLOPE already during the active phase, while the bottom boundary layer upwelling flow continuously refurnishes the cross-shore gradient and the thermal wind; this boundary flow diminishes in the passive phase making the instability waves extend to the bottom (not shown). For DOWNSLOPE the instability develops during the passive phase when three-layer structure associated with a detachment of the jet spans the bottom boundary layer. This is consistent with Stipa (2004) who found that in case of ${\displaystyle \mathrm{\Delta }<0}$ (DOWNSLOPE) the long waves are confined to the lower boundary. For increasing wavenumbers (shorter waves), the upper boundary comes into play, triggering interaction between the perturbations on both boundaries and their growth in time, because the shear maintains their phase difference while the mutually induced velocity field tends to diminish it. If bottom Ekman pumping is applied in this case, the short-wave bottom-trapped Eady solution is barotropised and becomes strongly decaying by increased Ekman pumping, explaining why the onset of instability for DOWNSLOPE in our simulations is postponed to the passive phase. The DOWNSLOPE pattern is similar in winter albeit with weaker thermal wind and PV gradients since the winter thermocline (colder water at the surface) weakens the salinity stratification and is also located close to the halocline, resulting in a deeper mixed layer (not shown).

For the wall side, the baroclinic instability is referenced directly to the PV gradient (Vallis, 2019). The vertical component of relative vorticity in Ertel PV (Equation 20) dominates and thus the cross-shore gradient yields ${\displaystyle {\mathrm{\partial }}_x({\mathrm{\partial }}_z({\zeta }^z+f)\rho )}$. Since ${\displaystyle {\mathrm{\partial }}_{xz}\rho }$ is not changing sign close to the wall during the period when the instabilities’ set on (Figure 3, first row), the necessary condition for the instability can only be satisfied for UPWALL if and when ${\displaystyle {\mathrm{\partial }}_{x}v}$ changes sign when the jet detaches from the wall (Figure 3, third row). The signs are reversed in DOWNWALL, but the same argument holds; DOWNWALL does not exhibit detachment (Figure 3, third row) which implies that it is always stable, consistently with energy diagnostics (Figure 5e–h). See Figure 4 for the summary.

3.4 Coastal eddy field

Baroclinic instability results in meandering of the coastal jet and formation of an eddy field. Because instabilities associated with the relaxation of the downwelling front occur without any signature on the surface density, following Durski and Allen (2005) we choose to use the depth-averaged cross-shore velocity (transport per unit length in along-shore direction) ${\displaystyle U_w=\int u_{w}dz}$ (Figure 9a–f). In a uniform alongshore channel, the depth-averaged cross-shore velocity does not depend on time, any deviation in time would mark the presence of instabilities (see sum of equations 25 and 28). Note that the cross-shore velocity shows a vertical structure (Figures 7–8) due to growing perturbations (Stipa, 2004) as well as detachment of coastal jets in the relaxation phase (Figure 3). This vertical structure is averaged out in Figure 9a–f, which are used here to quantify the integrated cross-shore transport due to turbulent features generated by the instability.

Figure 9

The dominant wavelength is estimated from the first zero crossing of the autocorrelation function of U_w with distance off-shore. For the UPWALL, the unstable wavelength ${\displaystyle {\lambda }_{UPWALL}^{(S,W)}}$ growths from about R_bc to 50 km. The presence of the slope reduces the unstable wavelength (${\displaystyle {\lambda }_{UPSLOPE}^{(S,W)}\approx {\lambda }_{DOWNSLOPE}^{(S,W)}\approx 20 \mathrm{k}\mathrm{m}}$).

Those differences are qualitatively consistent with the results from the two-layer model of baroclinic instability. The absence of the sloping bottom maps onto zero topographic ${\displaystyle \beta =0}$, implying that there is only a high wavenumber cutoff ${\displaystyle k_c{R}_{bc}=2.8}$ [(Philips model, (Vallis, 2019))] and the size of unstable waves is restricted solely by the geometry of the channel (about 50 km in our simulations). On the sloping bottom, the most unstable wavelength ${\displaystyle k_m{R}_{bc}\approx 1.14}$ for UPSLOPE and for the DOWNSLOPE the most unstable wavelength is in the short wave range and depends on exact value of ${\displaystyle |\mathrm{\Delta }|}$ (Isachsen, 2011; Stipa, 2004).

Hovmöller diagrams (Figure 9g–i) show propagation of disturbances alongshore in time while they grow in size. In the absence of the mean flow, unstable waves on the slope would propagate with the shallow side to the right (negative velocity alongshore) due to the topographic beta effect which superimposes on the advection by coastal jets, which (vertically averaged) velocity drops by about half over the passive phase. For UPSLOPE (Figure 9g), the propagation is initially dominated by the advection by the positive coastal jet resulting in a positive alongshore propagation speed 100 km/3Ti ≈ 0.6 m/s, which implies a fast transport due to upwelling jet out of the WGB. During the passive phase the jet fades away and the slope effect yields a slow net positive propagation (50 km/5Ti ≈ 0.2 m/s). For the DOWNSLOPE, the propagation on the slope has same direction as the coastal jet (negative) but the instability sets in well into the passive phase when the jet is already weakened resulting in net negative propagation of about –0.4 m/s at ti = 15 and dropping to –0.3 m/s at the end of the simulation, i.e., continued significant transport southward of the domain long after the wind has ceased. Note also that there is a signature of small scale, weaker instabilities on the inner shelf for both DOWNSLOPE and UPSLOPE. For UPWALL, there is no topographic beta effect and the instabilities are advected by the negative jet alongshore with a of 0.8 m/s at the beginning of the passive phase and then slow down toward the end of the passive phase to 0.3 m/s as they also interact with disturbances from the slope side propagating offshore.

3.5 Implications for mixing

The role of coastal upwelling for mixing in the Baltic Sea have been studied previously in an idealized configuration (Zhurbas et al., 2006) and its sub-basins [e.g., Gulf of Finland, (U. Lips et al., 2017; Väli et al., 2017; Zhurbas et al., 2008), but no studies have considered WGB so far. We quantify the implications for momentum and tracer mixing using the vertical eddy mixing coefficient K_V (state-dependent output from the model KPP parameterization scheme) and the surface cross-shore eddy salinity flux ${\displaystyle \{(u-\{u{\}}_x)(s-\{s{\}}_x){\}}_x}$, where ${\displaystyle \{\dots {\}}_x}$ indicates the cross-shore mean (Figure 10). Large values of K_V ($>{10}^{-3}{\text{m}}^2\hspace{0.17em}{\text{s}}^{-1}$) prevail in the surface and bottom Ekman boundary layers during the active phase, with a thicker bottom layer and higher values for the DOWNSLOPE (due to a lower stratification there) compared to UPSLOPE (Figure 10a and d). While for the UPSLOPE the high vertical mixing layers first thicken at the unstable frontal jet and then successively fade away during the passive phase (Figure 10b–c), they persist for DOWNSLOPE throughout the simulation time (Figure 10e–f). A striking result is that for DOWNSLOPE, the high K_V values are found on the inner shelf, suggesting that the downwelling jet on the sloping bottom may play a prominent role in determining the vertical mixing of the broad inner shelf parts of the WGB. The DOWNWALL has no clear signal in K_V (Figure 10a–c) while baroclinically unstable jet and eddies for the UPWALL come across as centres of elevated vertical mixing K_V ($>{10}^{-2}{\text{m}}^2\hspace{0.17em}{\text{s}}^{-1}$) consistent with the ‘chimney effect’ of baroclinic eddies for vertical mixing found in earlier studies (Koszalka et al., 2010).

Figure 10

On the contrary, the surface horizontal mixing is generally low during the active phase (not shown) and gets successively elevated due to instabilities and baroclinic eddies (Figure 10g). This is consistent with observations and modeling results about the important role of unstable jets and baroclinic eddies for cross-shore transport in other regions (Isachsen et al., 2012). The cross-shore eddy fluxes are one order of magnitude larger for UPSLOPE and UPWALL (due to baroclinic instability with larger unstable wavelengths) compared to DOWNSLOPE (baroclinic instability with smaller wavelengths, Figure 9) and DOWNWALL (no instability) and concentrated within the first 20 km from the coast, i.e., marked by the location of the coastal jet (Section 3.2). Thus, while the northeasterly wind regime with downwelling on the Swedish coast and upwelling on Gotland side is a much more rare event in the WGB (10% time compared to about 30% for upwelling; Figure 1c), it is a more effective ‘bottom vertical mixer’ (with implications for the bottom boundary layer biogeochemistry and coastal hypoxia (Fredriksson et al., 2024))]. On the other hand, the southwesterly wind regime (upwelling on the Swedish coast and downwelling on the Gotland side) are more effective ‘surface horizontal mixers’ [with implications for, e.g., algae blooms and air-sea exchange; (I. Lips and Lips, 2008); (Dabuleviciene et al., 2020); (Norman et al., 2013)).

3.6 Topographic effects and relevance of idealized simulations

So far we focused on simulations with a constant slope on the Swedish coast. In reality, the slope is not uniform but steepening from south towards the Landsort Deep. However, for the range of observed slopes Bu is always small (between 0.2 and 0.4), barely modulating the contribution of the boundary layer to the offshore transport, so our results about the main features of the coastal circulation hold over.

The bathymetry of the WGB is complex, resembling a sinusoidal shape but with a deep depression of the Landsort Deep in the northern part of the domain and a few larger bays and gulfs (Figure 1c). As a next step in refining the representation, we considered a sinusoidally varying slope topography (Section 2.3). It resulted in a meandering of the coastal jets due to squeezing and stretching of the water column, but the qualitative picture of baroclinic instability did not change (Figure A2). Several other studies addressed the interaction of locally varying bathymetry with upwelling in similar idealized configurations. For example, Saldías and Allen (2020) studied the upwelling on the shelf with a deep canyon and open sea boundary. The presence of a submarine canyon imprinted on the circulation and coastal upwelling structure close to the canyon; part of which can be associated with squeezing and stretching resulting in jet meandering. The baroclinic instability was still the main instability mechanism while the presence of the canyon increased the wavelengths downstream. They attributed part of the observed variability to the presence of coastal trapped waves (CTWs) (i.e., Kelvin-like waves propagating along a sloping coast) (Walin, 1972a, b). Also Fennel and Seifert (1995) addressed the interaction of such waves with upwelling in the southern Baltic Sea with idealized model configuration but did not quantify them in detail. The variability in hydrographic observations on western Gotland’s coast have been traced back to propagation of CTWs excited at the southern tip of Gotland (Pizarro and Shaffer, 1998). Fennel et al. (2010) studied the interaction of Kelvin waves generated at the northern and southern tip with upwelling around Gotland. The upwelling signal is amplified or reduced by Kelvin waves at different locations around the island, depending on local stratification and water depth which determines the wave dispersion and damping. The (frequent) downwelling signal from the western part of the island is propagated by the waves to and along the eastern side reducing the upwelling signal there, but the opposite is not occurring on the western side due to stratification-depth conditions there. Our idealized simulations are not designed to study the dynamics of CTWs in the WGB. A follow-up study should address the impact of locally complex bathymetry and the presence of Landsort Deep in the WGB, including the interaction with CTWs to identify locations of persistent and wave-modified upwelling/downwelling patterns and their impacts on biogeochemistry in the WGB. Our idealized configuration does not fully represent the role of inertial and inertial-internal waves in the dynamics because of the lack of continuous variations in wind forcing. These waves are transient phenomena in our simulations and the inertial frequency band (corresponding to periods of 12–16 hours) makes up to 13% of the total variance for the cross-shore velocity component at the depth of the summer thermocline, and below 5% for other components). Inertial and internal waves were however found significant in the observations of upwelling (Walter et al., 2016), idealized upwelling simulations with forcing oscillating in time (Brink and Lentz, 2010) and idealized simulations of ocean eddy field forced with spatially-varying winds (Koszalka et al., 2010). Recent observational and modelling studies suggest that horizontal transport due to the wind waves is important in coastal seas (Stanev et al., 2019; Staneva et al., 2021). Future more realistic and focused numerical studies are needed to assess the wave dynamics and wave-induced transport in the WGB.

Another important aspect not considered here is the modulation of Ekman dynamics by surface buoyancy fluxes on daily to sub-seasonal time scales. A modeling study addressing this issue was done by Nerheim and Stigebrandt (2006) to explain variability in observed buoy data. Including buoyancy flux would be a necessary step towards more realistic studies of the dynamics in the WGB.

The idealized simulations can help to isolate the underlying instability dynamics and associated time and space scales. It is not easy to isolate the upwelling and downwelling in observations because of the transient nature of the wind events and interaction with complex topography, coastally trapped waves, river run-off and air-sea interactions. We downloaded and inspected GHRSST (https://www.ghrsst.org/), and MODIS data (https://oceandata.sci.gsfc.nasa.gov), but the upwelling and downwelling development on two sides of the WGB is hard to study systematically using satellite products because of the frequent cloud coverage and few available swaths over the area. Moreover, downwelling is not easily detectable in SST images, and upwelling is hard to detect during winter. Even during the summer, the detectability of the upwelling may be limited, depending on the intensity and duration of the wind event that determines the outcrop of the summer thermocline (Section 3a–b). Using reanalysis products circumvents the issue with the cloud cover and provides velocity fields, but is instead prone to errors due to limited model resolution (nominal resolution of 1 nm = 1,852 km) and data assimilation.

Nevertheless, in lieu of satellite data of sufficient coverage and resolution, the reanalysis products can help to analyze velocity- and derived quantities and regionalize our idealized results as well as to discuss their limitations in representing the more realistic representation offered by the reanalysis. To exemplify the representativeness of our idealized simulations, we choose a north-easterly wind event in August 2021 evidenced in wind observations (Figure 11a) and show results from reanalysis at about t_i = 16 duration of the wind event (SST and surface relative vorticity; Figure 11b and e). This meteorological situation is less frequent in WGB (north-easterly winds occur with about 10% frequency, Figure 1c) and corresponds to our idealized DOWNSLOPE-UPWALL configuration where the jets on both basin sides are unstable (with instability onset times of about 14 and 4 inertial periods for DOWNSLOPE and UPWALL, respectively, Figure 4) and an eddy field spanning the whole basin develops from the unstable UPWALL coastal jet (Figure 9).

Figure 11

Figure 11c and f show, respectively, the SST anomaly and the relative vorticity from our idealized simulations with summer stratification for uniform bathymetry, and panels d and g show the same for the sinusoidally varying bathymetry (LANDSORT). The SST anomaly is computed from the mean in the domain before the onset of upwelling (26th August 2021). Consistent features between the reanalysis and idealized simulations are: 1) a downwelling jet on the eastern (DOWNSLOPE) side indiscernible in SST field but marked by negative values in the relative vorticity because of squeezing over the sloping bottom [e.g., (Zhurbas et al., 2006) and showing signs for onset of instability after about two weeks with a few kilometer-scale meanders, 2) unstable eddy visible in SST and relative vorticity field with about 25 km length scale formed due to instability of upwelling jet on the western (UPWALL) side (for a wind impulse of about 11,000 Pa s; in accordance to our results in Figure 5e).

Note that the observed wind is varying during the event duration so the observed wind impulse is lower than in our idealized simulations (Figure 1d), but it is sufficiently large to generate instabilities in accordance to our simulations (EKE/MKE >0, see Figure 5e) while they are, as expected, less intense (EKE is smaller in reanalysis as only one eddy feature is observed). There are also differences, however. In the reanalysis, the downwelling jet seems to be advecting a cold water anomaly into WGB from the northern direction disrupting the downwelling signal inside the basin, there is signature of runoff and small scale flow features at the eastern coast, and the coastal jets on both sides of WGB interact with the flow associated to the Oland Island located in the south. These interactions cannot be captured in our idealized simulations with periodic north-south boundaries and idealized coasts; however, this comparison shows that the salient dynamics of coastal jets and related instabilities is well represented in idealized simulations on length scales comparable to, and larger than, the baroclinic Rossby radius and time scales longer than inertial.

To address the relevance of an idealized simulations of the WGB channel of a finite width, we investigate the time evolution (Hovmöller diagrams) of the velocity components (Appendix Figure A3, top three rows). These show evidence for interaction of the coastal flows across the domain after the onset of instabilities, involving both the subinertial dynamics and inertial waves (evident mostly in the cross-shore velocity component). We quantify further the interaction of coastal jets due to nonlinear dynamics by alongshore- and vertically averaged alongshore velocity tendency due to cross-shore advection during one inertial period (${\displaystyle u\mathrm{\Delta }v/\mathrm{\Delta }x\mathrm{\Delta }t}$, ${\displaystyle \mathrm{\Delta }t=t_i}$, bottom row in Figure A3). This interaction leads to a change in alongshore velocity of the order of a few cm/s during one inertial period and for DOWNSLOPE-UPWALL, it extends across the domain after the instabilities develop. Note that the interaction between the flows on two sides of the domain implies that it is not straightforward to replace our channel configuration by separate model configurations of coastal jets on the two sides of the channel as in e.g., Brink (2016). Appropriate open boundary conditions would be needed not only to cope with the artificial reflection of outgoing waves offshore but also tackle the non-zero geostrophic flow in the middle of the domain (see also numerical experiments and discussion in Nycander and Doos (2003)). The MITgcm has a choice for prescribing the flow at the open boundary following Stevens (1990) that could possibly be adopted to this case in a future work.

3.7 Sensitivity studies

We ran a simulation with a higher horizontal resolution (500 m), however, this showed no notable changes in the coastal jet and instability dynamics, implying that the nominal 1 km resolution was sufficient to represent the baroclinic flows given the regionally typical winds and stratification (${\displaystyle R_{bc}=7-8 \mathrm{k}\mathrm{m}}$). We ran simulations with the quadratic bottom drag coefficient (${\displaystyle C_d={10}^{-4}}$) to assess the sensitivity to the drag parameterization. The formation of the boundary layer and pycnic displacements showed no significant change, nor did the instability dynamics. The bottom velocity however, became stronger at the end of the active phase and the decay of the mean flow was slower in the passive phase for the quadratic drag case compared to the linear one. The estimated buoyancy arrest time scales (Equations 31–32) are even longer (${\displaystyle t_i\sim 13}$ for UPSLOPE, ${\displaystyle t_i\sim 500}$ for DOWNSLOPE) for the quadratic bottom drag. It is to be noted that the buoyancy arrest time scales depend on the product of the drag coefficient and stratification ${\displaystyle d=C_{d}N}$ (in a square-root manner for upwelling, and linearly for downwelling; Equations 31–32), which implies a sensitivity to the value of C_d, although one can expect that the C_d and N are not independent (stronger friction would lead to a weaker bottom flow and more quiescent conditions and thus stronger stratification). To test this, we also run a simulation with much higher $C_d={10}^{-3}$, and we did not observe a substantial change in the dynamics, consistently with Brink (2016).

A sensitivity study to the vertical mixing scheme for coastal upwelling has been done by Durski and Allen (2005) who did not observe remarkable differences between simulations run with the Mellor–Yamada and the KPP schemes. We expect that their result would carry to our configuration where most of the cross-shelf transport occurs in the bottom boundary layer (${\displaystyle Bu\approx 0.25}$), with only a minor contribution from the interior and surface mixed layer (where one would expect largest differences due to vertical mixing schemes).