Mechanisms of bottom boundary fluxes in a numerical model of the Shetland shelf

Across-slope bottom boundary layer (BBL) fluxes on the shelf-edge connect this region to deeper waters. Two proposed ways in which across-slope BBL fluxes can occur, in regions that have a slope current aligned to the bathymetry, are the frictional veering of bottom currents termed the ‘Ekman drain’ and through local wind-forced downwelling (wind-driven surface Ekman flow with an associated bottom flow). We investigate the variability, magnitude and spatial scale of BBL fluxes on the Shetland shelf, which has a prominent slope current, using a high-resolution (∼2 km) configuration of the MITgcm model. Fluxes are analysed in the BBL at the shelf break near the 200 m isobath and are found to have a seasonal variability with high/low volume transport in winter/summer respectively. By using a multivariate regression approach, we find that the locally wind-driven Ekman transport plays no explicit role in explaining daily bottom fluxes. We can better explain the variability of the across-slope BBL flux as a linear function of the speed and across-slope component of the interior flow, corresponding to an Ekman plus mean-flow flux. We estimate that the mean-flow is a greater contributor than the Ekman flux to the BBL flux. The spatial heterogeneity of the BBL fluxes can be attributed to the mean-flow, which has a much shorter decorrelation length compared to the Ekman flux. We conclude that both the speed and direction of the interior current determines the daily BBL flux. The wind does not explicitly contribute through local downwelling, but may influence the interior current and therefore implicitly the BBL fluxes on longer timescales.


Introduction
Holt et al (2009) modelled the entire north-western European continental shelf to study bottom boundary layer (BBL) fluxes on the shelf-edge. They identified a regional downwelling circulation with two parts: wind-driven, on-shelf transport leading to downwelling; and local, off-shelf transport in a near-bed Ekman layer, termed the 'Ekman drain' (Souza et al, 2001). Observations by Simpson and McCandliss (2013) on the Hebridean shelf edge have contributed towards evidence of an 'Ekman drain'. Meanders of the regional slope current (Sherwin et al, 1999(Sherwin et al, , 2006 and eddies are considered elements of across-slope transport ) over the entire shelf-sea depth, but not specifically in the BBL.
Fluxes near the shelf break are important as they connect shallower shelf seas to deeper waters. Shelf seas play a significant role in the uptake of atmospheric CO 2 (Thomas et al, 2004;Ryther, 1969). Tsunogai et al (1999) proposed the continental shelf pump mechanism in an attempt to explain why this is the case. As part of it, the annual absorption of atmospheric CO 2 into shallow shelf seas requires it to be transported off-shelf to deeper waters. The development of a seasonal pycnocline inhibits this off-shelf transport in the ocean layers above the pycnocline, but not below. So in a near-bed bottom boundary layer, export of carbon off-shelf can occur year-round. The hydrodynamic processes responsible for BBL volume fluxes contribute to the sustained off-shelf export of carbon near the seabed.
Observations of carbon uptake on sections of the north-western European shelf (such as the North Sea) find it to be a carbon sink (Frankignoulle and Borges, 2001;Bozec et al, 2005). There is also evidence that the carbon on this shelf region is exported off-shelf, consistent with the continental shelf pump hypothesis (Thomas et al, 2004;Bozec et al, 2005). The Shetland shelf, part of the larger European continental shelf, is therefore an interesting and relevant region to study hydrodynamic shelf-edge BBL fluxes (which impact on the carbon fluxes).
The Shetland shelf is one side of the Faroe-Shetland Channel (FSC): a deep bathymetric channel extending north-eastwards between Scotland and the Faroe Islands ( Figure 1). A prominent feature is the Continental Slope Current. It has a high-speed core (Hansen and Østerhus, 2000;Berx et al, 2013) with current speeds increasing from 10 cm s −1 to 20 cm s −1 at the shelf-break (around 200 m; Dooley et al, 1976;Turrell et al, 1992) to 40 cm s −1 over the 1000 m contour (Dooley et al, 1976). It is vertically coherent being predominantly barotropic to depths of 500 m and is composed of North Atlantic Water (Sherwin et al, 2008). It is continuous along the shelf-edge, existing from the Hebridean shelf-edge (Souza et al, 2001) with its origins as far south as the Brittany peninsula (Pingree and Le Cann, 1989). Observations (e.g. Sherwin et al, 1999Sherwin et al, , 2006 and model studies (Oey, 1998) of the slope current have identified mesoscale meanders and eddies, identifying them as important in the across-shelf transport and mixing of slope current water with Modified North Atlantic Water located on the Faroe Plateau (Sherwin et al, 1999). Meanders on short timescales (∼ days) are theorised to be caused by baroclinic instability (Sherwin et al, 2006). On interannual timescales increased cross-shelf transport is associated with the negative mode of the North Atlantic Oscillation (Chafik, 2012).
There are several gaps in our knowledge. Firstly, we do not know the seasonal variability of BBL transport on the shelf-edge. Secondly, the hydrodynamic processes responsible for BBL fluxes have been identified, such as the Ekman drain and downwelling-however in the simple Ekman drain model (Souza et al, 2001), poor correlations and a large variability have been noted between across-slope BBL and calculated Ekman transports Simpson and McCandliss, 2013) despite an apparently consistent slope current flowing parallel to the slope. What causes the short-term variability in the Ekman-drain model? Furthermore, Holt et al (2009) propose that the contribution of the wind to BBL fluxes should be valid across many shelf seas. However in the south-eastern Australian shelf Schaeffer et al (2014) have rejected the wind as a contributor. Can a wind-driven Ekman term, potentially driving a classical 2-D coastal downwelling, help to understand BBL fluxes on the Shetland shelf-edge which is part of the north-western European continental shelf? Finally, we do not know the spatial scales of the hydrodynamic processes contributing to BBL fluxes-can they be understood locally, or only in the context of shelf-integrated transport?
Using a high-resolution regional ocean model of the Faroe-Shetland Channel our study investigates the variability, magnitude and spatial scale of across-slope BBL transport on the Shetland shelf with three aims. Firstly, what is the seasonal variability of BBL transport on the Shetland shelf? Secondly, building on the previous work of Holt et al (2009) can we confirm the Ekman-drain model in a higher resolution model and, for the first time, identify and quantify the cause of the observed variability (Simpson and McCandliss, 2013) of across-slope bottom transport? Finally can we, for the first time, demonstrate the contribution of local wind-forced downwelling on the across-slope BBL flux for the Shetland shelf?
The model was initialised and forced at the lateral boundaries by daily oceanic fields (temperature, salinity, north/eastward currents) provided by the HYbrid Coordinate Ocean Model (HYCOM) reanalysis (Chassignet et al, 2007). Initial and hourly atmospheric forcing was provided by Climate Forecast System Reanalysis (CFSR) (Saha et al, 2010). Our simulation ran from 1 st January 2002 to 31 st December 2006; this period was chosen as it had good overlap with our current observations. To miti-  (Sandwell and Smith, 1997)  We did not explicitly force with tides, so to mitigate major tidal effects we averaged our currents daily. Our horizontal eddy viscosity was 5 × 10 −4 m 2 s −1 . Larger values resulted in poorer validation. Vertical eddy viscosities were parameterised using the KPP vertical mixing scheme (Large et al, 1994) with a background viscosity of 1 × 10 −4 m 2 s −1 . Our simulations ran in hydrostatic mode because non-hydrostatic runs had no appreciable difference in validation. For the surface boundary condition we used an implicit linear free-surface. For the bottom boundary we used a free-slip condition with an explicit quadratic drag coefficient of 2.5 × 10 −3 .

Orientation
We focus on fluxes crossing the 200 m isobath and split the isobath into many 'stations' interpolated from the model grid points that effectively act as proxies for observation sites. We define our BBL as the 180 m to 200 m layer. At a given station the across-slope BBL flux, Q BBL = V BBL × h BBL , where V BBL is the across-slope BBL Fig. 2 Diagram of currents centred on a section of the 200 m isobath. Our study is 2-dimensional locally (no arrows have a vertical component) but we project into 3-dimensions for context: the 200 m isobath divides on-and off-shelf waters. The interior slope current, V int (thick light grey arrow), and BBL current, V BBL (thick dark grey arrow), are decomposed into their along-and across-slope components (thin dashed arrows). V BBL corresponds to the across-slope BBL flux. U int and V int are the along-and across-slope components of the interior flow. The slope angle is given from the bathymetry (θ LG ) or using the Taylor-Proudman Theorem (calculated from the time-mean of θ T PT ). A bottom stress τ (thick light grey dashed arrow) opposes the interior current and induces a perpendicular Ekman flux Q E (thick pale dashed arrow). current ( Figure 2) and h BBL is the BBL height (20 m). To calculate V BBL we need to know the slope angle, θ s , so we can orientate our eastward and northward currents to the local along-/across-slope directions. We used two methods to determine θ s at a given station. For the Local Gradient (LG) slope angle (θ LG s ), we interpolate latitude and longitude midway between stations (black dots in Figure 2) to calculate the angle, relative to east, between these points: so θ LG s = θ LG (Figure 2). To calculate the Taylor-Proudman Theorem (TPT) slope angle (θ T PT s ) we use the result that geostrophic flow on an f -plane follows isobaths. We assume the interior current V int is geostrophic. The angle it makes relative to east is θ T PT , so we set the slope angle θ T PT s to be the time-mean of θ T PT ; or θ T PT s = θ T PT where . . . denotes timeaveraging. This method is also an approximation to the LG method that we consider exact (given the model resolution), and is used by Simpson and McCandliss (2013) in their study of BBL fluxes. We aim to understand how the flux estimates depend on the orientation method.

Fluxes
We model Q BBL (the BBL flux over the 180 m to 200 m layer) as a function of the current above the 20 m BBL layer. This current is termed the interior current or V int , and is calculated as the depth-averaged current in an 80 m layer (100 m to 180 m) above the BBL. We chose an 80 m layer because: we want to sufficiently capture the mean flow in the ocean interior; be far away from the surface Ekman layer; and to minimise against any bottom boundary effects that may pervade into the next vertical cell.
In the absence of friction, and away from lateral boundaries, the interior current V int should extend to the seabed. Then the across-slope BBL flux, Q BBL , can be estimated from the mean flow, that is, the across-slope component of the interior current V int (section 3.4.3).
In the presence of friction an Ekman spiral develops where the current is deflected due to the effect of the bottom stress. The vertical integral of the Ekman spiral gives a transport that is perpendicular to the mean-flow. For brevity we will call this vertically integrated layer the veering layer. For this layer a bottom stress, τ, opposes the direction of V int , and so veers the BBL current to the left of the interior current in the Northern hemisphere. This is quantified by a veering angle, θ veer . The veering also corresponds to an Ekman flux that can be modelled as, where Q E denotes the Ekman flux (dimensions L 2 T −1 ) perpendicular to the bottom stress τ, with ρ the density and f the Coriolis parameter. As we would like a simpler representation of the across-slope BBL flux as a function of interior current components, we model τ as a quadratic function of the interior velocity, where γ 2 is a dimensionless quadratic drag coefficient. Then the across-slope Ekman flux, Q ⊥ E , is caused by τ , the component of the bottom stress parallel to the slope. Taking θ to be the angle between V int and the slope, τ is then where U int is the along-slope component of the interior current. Alternatively we can model τ as a linear function of the interior velocity, so where γ 1 is a linear drag coefficient with dimensions L T −1 . In the presence of tides the bottom stress can be parameterised linearly (Bowden, 1953;Hunter, 1975) though we have neglected them here. We investigate both linear and quadratic τ parameterisations for our model of BBL fluxes (Section 3.4.1) to see if there are any differences in approach. An additional flux considered by Holt et al (2009) andSchaeffer et al (2014) is a wind-driven two-dimensional downwelling circulation. Here a local wind-driven surface Ekman transport (that is perpendicular to the 200 m isobath), Q ⊥ W , may lead to an additional BBL flux below. We model this additional flux similarly to Eq. 3, with the interior current replaced by the 10 m surface wind components parallel to the 200 m isobath. We use the same CFSR surface wind fields as used previously for model atmospheric forcing (section 2.1) and consider separately the contribution of wind to the BBL flux (section 3.4.5).

Regression model
To test the simple Ekman drain model, we perform a univariate linear regression between BBL and Ekman fluxes (section 3.4.1). However, we wish to quantify for the first time the effect the across-slope interior flow has on the BBL flux (section 3.4.4). We propose a multivariate linear regression model where we decompose Q BBL into an Ekman variable (Q ⊥ E ) and a variable which characterises the mean flow, V int , which corresponds to the across-slope component of the interior flow. These variables have regression coefficients Γ and α. The unaccounted physics in the regression model can overall be characterised by an intercept, C. The first three terms provide a best fit of the data; the additional variability of any individual datum around this best fit is given by the ε term, or residual. In section 3.4.5 a wind-driven Ekman term, , is added to the multivariate regression of Eq. 5.
Performing this regression for each station yields local coefficients for the Ekman/interior flux variables (section 3.4.4). Alternatively, we can integrate Q BBL , Q ⊥ E and V int along the entire shelf to calculate integrated transport terms, then perform a regression that yields global coefficients for a single shelf-integrated transport model (section 3.5). The R 2 -value of the regression model indicates whether the independent variables can explain the variability in the dependent variable. As increasing numbers of parameters may artificially inflate the ordinary-R 2 in the regression model, we quote the adjusted-R 2 which accounts for different numbers of parameters.
In this study we quote sample estimates as the mean ± 1 standard deviation, or x ± s. However some distributions (e.g. R 2 for different locations) are not normal-so we quote the sample mean with interquartile range (IQR), orx[IQR 25 % , IQR 75 % ]. We perform hypothesis tests at the 1 % significance level and where appropriate provide 99 % confidence intervals (CI).

Decorrelation length scale
The decorrelation length scale is the distance over which two time-series are sufficiently far apart so that they are independent. This distance can help identify processes and instruct how far apart measurements have to be to optimise data collection to prevent spatial aliasing (Brink and Robinson, 2005). We use this length scale to diagnose the horizontal scales of variability for Ekman, BBL, V int and residual transports. We performed two methods of decorrelation for variables in our regression model (Eq. 5). Our first method finds the correlation coefficient, r, of the flux (BBL, Ekman, V int , residual) between all pairs of stations and bins this by the distance between the stations. For each bin we calculate the mean correlation (plus/minus one standard deviation). The mean (plus/minus one standard deviation) correlation per bin is a function of distance, and we infer the decorrelation length scale from the e-folding distance. The second method calculates the normalised root-mean square difference (RMSD) of the Ekman time-series between all pairs of stations. Normalisation is given by dividing the RMSD by the mean of the first time-series. For a given distance, the minimal-RMSD quantifies the maximal similarity of any two time-series along the shelf. The minimal-RMSD starts small and grows with distance. Where this growth stops (or continues to rise to a far-field value slowly) is the decorrelation length scale (Böhme and Send, 2005).

Mean currents
A slope current directed north-eastwards can be seen at 190 m ( Figure 3). Sherwin et al (2006) observe the fastest currents on average are near 61.25 • N,−2 • E from archive drifter data. In our model the fastest section of the slope current (average speed > 0.40 m s −1 ) is further downstream. The slope current has two fast sections (average speed > 0.30 m s −1 ) split around 60.5 • N,−3.5 • E. The circulation is also concentrated near 60 • N,−6 • E; this flow comes from the Wyville-Thompson Ridge (WTR) region and from the Faroe-Bank Channel (FBC). There is also a bifurcation (or meander) in the mean circulation at 60.5 • N,−5 • E. This meander has been observed by Sherwin et al (2006) in the mean current flow of drifters and from sea surface temperature (SST) snapshots, and by Sherwin et al (2008) over a week-long composite of SST fronts. East of the Faroe plateau there is a clockwise circulation of flow south-westwards, a well-known feature in the region (e.g. Hansen and Østerhus, 2000). A persistent eddy appears on the northern boundary at 62.75 • N,−1.5 • E. This could be caused by the HYCOM velocity boundary conditions on the northern boundary interacting with the strong slope current.

Sea Surface Temperature
Mean model SST is compared with observations ( Figure 4) from the Group for High Resolution Sea-Surface Temperature (GHRSST; Donlon et al, 2009). The model has a warm bias with a mean anomaly of (0.09 ± 0.67) • C. In a previous model study using the Regional Ocean Modelling System (ROMS), Broadbridge and Toumi (2015) show a cold bias for 2005 of −0.77 • C (−1.93 • C, −0.15 • C) where there they give the range. Some basic spatial structures are captured by the model: (i) a strong meridional temperature gradient between the northern Hebrides and southern Faroes; (ii) the north-eastward extension of the 9 • C to 10 • C isotherms; (iii) the warmest water (T > 11 • C) located west of the Hebrides and coldest water (T < 9 • C) north of the Faroes. The model is warmer compared to observations towards the northern boundary, with the largest warm anomaly north-east of the Faroe plateau (∆ T > 0.6 • C). The largest cold anomalies are relatively smaller (∆ T < −0.3 • C) and mainly towards the eastern boundary. We presume these warm/cold anomalous regions are due to our HYCOM boundary conditions which may be warmer/colder than they should be.

Current Profiles
We compare snapshot currents from our model and the HYCOM global ocean reanalysis (also used for our initialisation/boundary forcing) with current observations. Observational datasets were collected for British Petroleum (BP) and their partners by Fugro GEOS at two locations, Foinhaven (Foi) and Schiehallion (Sch), each at three depths. Current meter moorings and platform mounted Acoustic Doppler Current Profilers (ADCPs) were used to collect data over the period 17 September 1992 to 9 September 2007 (Foinhaven) and 21 September 1993 to 30 July 2007 (Schiehallion). Current datasets consist of 10-minute mean velocities. It was assumed that current data was fully screened for errors and spikes before archival in accordance with Fugro GEOS quality control. Quantile-Quantile (Q-Q) plots are used to scatter the percentiles of different current speed datasets against each other, and so compare whether the two datasets have the same underlying distribution ( Figure 5). At both sites HYCOM currents overestimate observational currents for the small and mid-range values, and underestimate the extremes. In comparison the MITgcm model currents slightly underestimate the observations for small and mid-range values but fit the extremes better than HYCOM. The MITgcm performs better than HYCOM at larger depths.
At the ADCP observational sites current direction is mostly north-eastwards (Figure 6 and 7). At Foinhaven, HYCOM currents are directed predominantly on a single bearing (60 • T) whereas the observations and our model have two main bearings (30 • T and 60 • T). Also, the size of the speed bins (colours in Figure 6) for our model match the observations better than for HYCOM. So the directional distribution of MITgcm currents are better than HYCOM-this may be due to the smaller  grid spacing (2 km vs. 10 km). At Schiehallion there is a clockwise bias of currents compared to observations, which is the same for HYCOM. We do not know why this is the case. Despite this bias our speed distribution for each direction remains superior to HYCOM.
The temporal correlations of both HYCOM and MITgcm with observations was poor (not shown). We note that the observations are located near the slope where there is potential for baroclinic instability (Sherwin et al, 2006). Instabilities create random variability in currents, which are extremely difficult to capture in models and may explain why the correlations were poor.

Hydrographic observations
We present model data at the Fair-Isle Munken (FIM) cross-section and compare it against a 14 year mean climatology from Berx et al (2013) (Figure 8; for location of the cross-section see Figure 1). We have not validated the Nolso Flugga (NOL) cross-section as we do not have corresponding observations, but we show the results here for additional visualisation and context of the regional dynamics. Both FIM and NOL cross-sections show downwelling of isotherms towards the Scottish continental shelf (Figure 8a,b). The largest temperatures (T > 10 • C) are also found here. At FIM, a clear thermocline exists near 500 m. However at NOL, the thermocline is less pronounced and the stratification is fairly uniform. Salinity profiles (Figures 8c,d) show much weaker stratification of isohalines, especially towards the Faroe plateau. In near-surface and mid-depth there is a meridional salinity gradient with fresher water towards the Faroe plateau and a high salinity core, indicative of the Continental Slope Current, located on the Scottish shelf. At greater depths salinity is more uniform across the channel.
We replicate the downwelling structure of isotherms well (Figure 8a) compared to observations. The model matches the depth of the 5 • C isotherm (400 m to 500 m). For salinity, model isohaline structure is not as well-matched (Figure 8c). In the model there is a freshwater bias (∆ S ∼ −0.1) in the upper 500 m and the absence of a clear halocline. We ascribe this bias to the HYCOM initial conditions, which are also too fresh in the upper layers compared to the observations (not shown). Our results on bottom fluxes (section 3.2 onwards) should not be significantly affected by this bias.
Mean flow along the channel is marked by the presence of the slope current centred on the 500 m isobath (Figure 9a,b). The largest slope current speeds (max. > 0.4 m s −1 ) are in the NOL section. Deep overflow waters do not reach as high speeds. Both sections also show a small south-westward return flow at depth and near the Faroe plateau. The across-channel mean circulation at FIM shows the slope current is directed towards the shelf (< 0 m s −1 ) (Figures 9c). At NOL its direction is away from the shelf (> 0.1 m s −1 ).
Comparing Figure 9a against mean along-slope currents from Berx et al (2013), we have replicated the well-known high-speed slope current structure (see e.g. Hansen and Østerhus, 2000). In the model the 0 m s −1 delimiting contour extends too far north. Deep current structures support the results of Broadbridge and Toumi (2015) who reported a complex flow field at the bottom of the channel. For the NOL section currents were rotated perpendicular/parallel to the cross-section (see Figure 1).

Shelf-edge exchange
From the model, the volume transport across the shelf-break (the shelf-break is defined as the 200 m isobath) is 3.7 Sv directed off-shelf. Averaging the volume transport quantity horizontally along the shelf (so dividing by the length of the shelf) gives a flux of 4.2 m 2 s −1 . Averaging this new quantity over the 200 m isobath (so dividing by 200 m) gives 21 × 10 −3 m s −1 ; this value is equivalent to the average off-shelf velocity for all model grid cells. Burrows and Thorpe (1999 ; Table 4) estimate the across-shelf mass flux to be −10 × 10 −3 m s −1 (summer) and 15 × 10 −3 m s −1 (winter) over 200 m; so combined it is off-shelf, about 5 × 10 −3 m s −1 . There is large uncertainty on the summer estimate as only one drifter crossed the 200 m isobath. Huthnance (1995) estimate the same quantity to be 6 m 2 s −1 to 7 m 2 s −1 over 500 m, or 12 ×10 −3 m s −1 to 14 ×10 −3 m s −1 vertically averaged. Turrell et al (1992) used a single current mooring near the 200 m shelf-edge, and found for a 3 month summer period the acrossshelf current was on-shelf at both 30 m and 187 m (−25 × 10 −3 m s −1 and −6 × 10 −3 m s −1 ). They defined the slope current direction inferring bathymetry from Admiralty charts. Comparing the previous vertically averaged estimates with our higher resolution estimate (21 × 10 −3 m s −1 ), we have off-shelf transport though it is not unreasonable given the large uncertainty and sparsity of observations.

BBL veering
Is the slope current topographically locked to the bathymetry? We can answer this by considering the difference between θ T PT s and θ LG s orientations (section 2.2; recall θ T PT s is defined from the time-averaged interior current direction whereas θ LG s is defined from bathymetric data). The shelf-averaged difference of θ T PT s − θ LG s is (5 ± 32) • . This is significantly greater than 0 (p < 0.01, 1-tailed t-test; 99 % CI 2 • to 9 • ). So on average, along the shelf, the interior current is not topographically locked but directed slightly off-shelf. This indicates the mean-flow is important to understanding the BBL fluxes.
Are BBL currents also directed off-shelf? θ LG veer is the time-mean veering of currents in the BBL with respect to the θ LG s orientation (V BBL in Figure 2). Averaged along the shelf, θ LG veer is (14 ± 33) • . This is statistically greater than 0 (p < 0.01, 1-tailed t-test; 99 % CI 10 • to 18 • ). So on average the BBL currents are directed off-shelf (as is the interior slope current).
However do the BBL currents veer with respect to the interior current, potentially due to Ekman dynamics? Under theoretical assumptions (e.g. constant flow field) the integrated transport in a bottom Ekman layer would be directed 90 • left of the mean-flow (northern hemisphere). Now θ T PT veer is the time-mean veering of currents in the BBL with respect to the θ T PT s orientation. Averaged along the shelf, θ T PT veer is (9 ± 8) • respectively. Though this result is not a 90 • veering (as would be the case for a theoretical Ekman layer) it is nonetheless significantly greater than 0 (p < 0.01, 1-tailed t-test; 99 % CI 8 • to 10 • ). So on average the BBL currents additionally veer off-shelf with respect to the interior current (which we presume is due to Ekman dynamics and aim to show in later sections). This indicates towards BBL fluxes being a combination of both mean-flow and Ekman dynamics, with the latter in spite of highly non-idealised conditions. We integrated the model BBL volume transport along the shelf (i.e. ∑VBBL × h BBL × L; where L is the mid-point distance between stations) into a time-series to investigate its monthly variation ( Figure 10). This shelf-integrated volume transport was also calculated for both orientations (θ LG s and θ T PT s ). Both orientations show similar sinusoidal variability though the θ T PT s has a smaller transport (by about 0.1 Sv) compared to θ LG s . There is a seasonality to the BBL volume transport, with a maximum in March and a minimum in August. The BBL transport is high (> 0.6 Sv) in the winter (Dec-Jan-Feb) and March. There is a rapid decrease in the latter spring months, to the lowest annual levels in summer (∼ 0.4 Sv). In September, the transport rapidly increases to the high winter levels. There are deviations from this general trend in June/October which have higher/lower transports when compared to their respective seasons.

Seasonal variation of the BBL fluxes
In the following sections we will decompose the BBL flux variability into the interior current speed (i.e. the Ekman term) and direction (i.e. on-/off-shelf given by V int ), and in doing so connect and quantify the seasonal variability of the BBL fluxes to the interior slope current dynamics. For now, an initial decomposition of the BBL transport into Ekman plus mean-flow transport ( Figure 10) shows a seasonality for both terms (irrespective of orientation), matching the seasonality of the BBL transport. The contribution of the mean-flow compared to the Ekman-transport is greater when using the θ LG s rather than the θ T PT s orientation. This is because the θ T PT s orientation is calculated from the time-averaged direction of the interior flow (θ T PT s = θ T PT from Figure 2). This results in some of the across-slope transport being absorbed into the definition of the TPT slope angle, that otherwise would have been ascribed to the mean-flow in the LG case.

Ekman flux
The basic Ekman-drain model linearly relates the across-slope BBL flux with the Ekman flux. Scatter plots between these fluxes show large variability with R 2 of about 0.3 and 0.1 for two locations (Figure 11). There is similar correlation when the bottom stress is parameterised as either a quadratic or a linear function of the interior current.
For the shelf-mean the correlation is also similar for linear or quadratic formulation (R 2 = 0.27[0.08, 0.43] and 0.25[0.09, 0.40] respectively). By definition the Ekman-drain model should have zero intercept. We find that the shelf-mean C (intercept) for the linear drag is (0.07 ± 0.48) m 2 s −1 and for the quadratic drag it is (0.30 ± 0.55) m 2 s −1 . Shelf-mean intercepts for the linear/quadratic drag are significantly different from 0, and also from each other (p < 0.01; 2-tailed t-test). We proceed with the quadratic drag as this parameterisation is used in the simulation (and has an explicit drag coefficient we can use for consistency; see Eq. 3). Overall, the low shelf-mean R 2 indicates that the basic Ekman-drain model does not explain much of the variability in the across-slope BBL flux.

Variation of the BBL flux with the interior current
To try and improve the basic Ekman-drain model we first investigate the variation around the best fit (Figure 11), and see if it can be explained by the across-slope component of the interior current, V int . Points above/below the line of best fit (i.e. residuals of the Ekman-drain regression fit) are frequently associated with the interior current directed off/on-slope. A method of examining this for all the locations is outlined: a regression line is fitted between across-slope BBL and calculated Ekman fluxes, but for the data subset where V int = 0. We count when a point above (or below) this line is also a point when V int is greater (or less) than zero. If V int is not related to the residual, we would expect this condition to be satisfied 50 % of the time. Averaged along the shelf, the percentage of time when this condition is satisfied is (84 ± 10) % (here the error is between shelf locations, not in time). This is significantly greater than 50 % (p < 0.01; 1-tailed t-test). This indicates that locations do have an additional transport, due to V int , from that predicted by the simple Ekman-drain model. So additional offshore/onshore fluxes are associated with an offshore (V int > 0) / onshore (V int < 0) interior flow.

The across-slope interior current, V int
Analogous to section 3.4.1, we also investigate an interior current-only model relating the across-slope BBL flux and V int where we find a shelf-mean R 2 = 0.68[0.58, 0.84] (not shown). These shelf-mean correlations are significantly higher (p < 0.01; 1tailed t-test) than in the case of Ekman-only flux (section 3.4.1).   Fig. 12 Regression plots of across-slope BBL flux against Ekman + V int fluxes (Eq. 5) for the same sites as in Figure 11. The line of best fit is provided for reference.

Multivariate regression of the BBL flux without wind
To further probe the relationships between the across-slope BBL, Ekman and V int fluxes we perform a multivariate linear regression (Eq. 5), using the quadratic drag formulation. The Ekman flux and V int are effectively independent from each other based on an analysis of Pearson's r: 75 % of locations are significantly independent (p < 0.01), with their observed t-statistic (t obs = r d f −2 1−r 2 ) smaller than the estimated t-statistic. For correlation the estimated t-statistic threshold is equivalent to an |r| > 0.06 because of the large sample size, and the shelf-mean r = 0.26[0.07, 0.54]. As such, for the locations with statistical 'significance' of correlation, the r-values are poor in general and do not equate to practical significance. Figure 11 shows that much more of the variability in the BBL flux can be explained by using multiple regression. At the two sites the R 2 values improve from about 0.3 and 0.1, to 0.8 and 0.6 respectively ( Figure 12). The intercept at Foinhaven is much lower than in the Ekman-only model but this is not the case at Schiehallion, where the value is similar. There still remains some residual variability despite the improvement.
The shelf-mean R 2 = 0.75[0.66, 0.87] for the multivariate model. The distribution of R 2 in the multivariate case is significantly greater (p < 0.01; 2-sample K − S test) than in either single variable models (Ekman-only R 2 = 0.25, V int -only R 2 = 0.68). Improvements to the Ekman-only model (section 3.4.1) are made at over 98 % of locations along the shelf when V int is included as a variable in the multi-regressive model. The median ratio of mean-V int to mean-Ekman fluxes in the multivariate model is 2.2[1.1, 4.8], demonstrating that the mean-flow contributes about twice as much as the Ekman flux to the across-slope BBL flux.
We additionally repeated our multivariate regression at the 400 m isobath (taking a 40 m BBL layer with a 120 m interior current layer above this). The multivariate model shelf-mean R 2 = 0.50[0.41, 0.60], compared to Ekman-and V int -only models with shelf-mean R 2 = 0.29[0.14, 0.43] and 0.43[0.31, 0.57] respectively. The multi-variate model shelf-mean R 2 is a significant improvement over both univariate models (p < 0.01; 2-sample K − S test), though smaller than the multivariate model at 200 m (R 2 = 0.75).
We have shown that for individual locations at the 200 m and 400 m isobath the BBL fluxes can be better explained as the sum of an Ekman flux plus mean-flow flux (by the improvement of multivariate model R 2 over univariate model R 2 ). Though we have explained much of the observed temporal variability of BBL fluxes, we do not claim that this regression model can be used as an effective predictor model. This is because the coefficients in the multivariate regression model of Eq. 5 have high spatial variability (not shown). The Γ (Ekman) coefficient distribution is nearly symmetrical about a mean 0.46[0.27, 0.63]. The α (V int ) coefficient distribution is positively skewed, with a modal values between 18 m to 19 m, with mean 13.8 [10.6,17.9] m. The C (intercept) distribution has mean (0.15 ± 0.20) m 2 s −1 , which is significantly lower than in the quadratic drag Ekman-only case (0.30 m 2 s −1 ; section 3.4.1), though still not significantly different from zero (p < 0.01, 2-tailed t-test).  The domain mean winds are predominantly westerly and so are downwelling favourable (Figure 13). This may potentially drive a 2-D downwelling circulation affecting the variability of BBL fluxes. We tested whether wind-driven surface Ekman transport, Q ⊥ W , would contribute to the BBL flux by the addition of this term to our existing regression model (section 2.4). For the 200 m isobath we found the shelfmean R 2 = 0.75[0.66, 0.87] with the intercept C = (0.15 ± 0.20) m 2 s −1 . There is no change in the parameters, and the goodness of fit is not significantly different from the no-wind multivariate regression model of Eq. 5. This was also true at the 400 m isobath. Thus we neglect the addition of a 2-D wind-forced downwelling term in our model and reject it as a candidate explanation for variability of BBL fluxes. Volume transports (Sv)

Ratio of V int to Ekman transports for global coefficient model
Fig. 14 (a) Shelf-integrated volume transport of the terms in the multivariate regression model (Eq. 5): the across-slope BBL transport (solid), the Ekman term (dashed) and the V int term (dash-x). Colours differentiate which θ s orientation was used. (b) The ratio of V int to Ekman terms from (a), for both orientations, with a constant reference line of unity plotted.
Estimates of the shelf-integrated volume transports of the terms from the multivariate regression of Eq. 5 (BBL, Ekman and V int ), as a function of the time-averaging of the currents, are all positive/off-shelf (Figure 14a). For the θ LG s orientation when currents are averaged over multiple days, the V int /Ekman transport increase/decrease with averaging but the total BBL transport remains constant. For the θ LG s orientation when currents are averaged daily, the BBL transport is 0.53 Sv and the contribution from the Ekman and V int flux is 0.15 Sv and 0.31 Sv (so V int approximately double). Applying our multivariate regression model (section 2.4) to the shelf-integrated transports calculated after daily current averaging, the model R 2 = 0.91. We also considered the addition of a shelf-integrated wind-driven Ekman term (c.f. section 3.4.5) where we found that the wind was not a contributing factor to the shelf-integrated BBL transport (not shown). The BBL transport for the θ T PT s orientation is 0.45 Sv. The Ekman and V int fluxes are 0.20 Sv and 0.14 Sv respectively (Ekman approximately 1.5 times larger). It may appear contradictory that there is a mean-flow contribution in the TPT case, given that the TPT slope angle is calculated from the time-mean angle (recall θ T PT s = θ T PT ; see Figure 2). This apparent contradiction can be resolved by noting that the magnitude of V int can change without a change in direction. Therefore in the TPT case, when integrating the fluxes over time, there can be a non-zero net V int transport. Only in the limiting case where all the off-shelf flow is equal in magnitude to the on-shelf flow will there be no net V int transport. Our result of non-zero V int transport in the TPT case therefore means that the off-shelf flow is stronger than the on-shelf flow. That being said, using the TPT orientation does cause the net V int and BBL transports to be smaller than in the LG case. This is because in calculating θ T PT s some of the across-shelf transport, that would be present in the LG case, is absorbed.
The ratio of V int to Ekman terms in the multivariate regression is above 2 for daily averages but increases over 10 for longer time averages (Figure 14b). The difference in orientation method is stark: compared to the θ LG s orientation, θ T PT s has nearly even ratios (for current averaging ≥ 4 days the V int term begins to dominate over the Ekman term but not to the same scale as for θ LG s ). In short, for the θ LG s orientation the V int term is dominant over the Ekman term, however for the θ T PT s orientation the Ekman contribution is more pronounced and nearly equal to V int contribution.

Decorrelation length scale
The decorrelation length scales of the transports were determined using two different methods (section 2.5). The e-folding distance of r (the mean correlation coefficient per bin) for Ekman transport is approximately 80 km and near-grid (∼ 5 km) scale for BBL and V int (Figure 15). There is a strong similarity between BBL and V int decorrelation. The standard errors are sufficiently small so that the decorrelation scale of the fluxes (Ekman and V int ) are significantly different. However, to one standard deviation the variability in e-folding distance is large: from ∼ 5 km to over 200 km for the Ekman flux, and from sub-grid scale to over 50 km for V int . The decorrelation length of ε (the residual flux in the regression model) is very small (on the order of grid spacing).
The second method of determining the decorrelation length scale for Ekman transport (section 2.5) shows the normalised RMSD grows rapidly with distance to 1.0 at 90 km and then slowly rises. We therefore take the distance over this rapid growth as the decorrelation length scale, approximately 90 km. As this is the same order of magnitude as the first method, we use it as a confirmation of the result from the previous method, i.e. the Ekman transport has a decorrelation length scale of about 80 km. The grid-scale spatial variability in across-slope BBL transport is due to the mean-flow.

Discussion
The major motivation for this study was to investigate the variability, magnitude and spatial scale of across-slope BBL transport along the Shetland shelf. Concerning short-term variability, a previous hypothesis (Souza et al, 2001;Simpson and Mc-Candliss, 2013) of an Ekman-drain model only partly explains the bottom boundary layer fluxes. We used a regression method to decompose the daily BBL transport into different physically based components. We introduce a multivariate regression model of the BBL transport (Eq. 5) as a linear combination of the Ekman flux and a mean-flow flux (the across-slope component of the interior flow). Adjusted-R 2 in the multivariate model is higher than in either the Ekman-only or V int -only single regression models (section 3.4.4), and this was confirmed for both the 200 m isobath (near the shelf break) and a deeper 400 m isobath. We find that the mean-flow term dominates over the Ekman term in explaining the variability (section 3.4.4 and 3.5).
One additional component of the BBL transport previously considered was alongshore winds driving upwelling and downwelling (e.g. Ekman, 1905;Niebauer et al, 1977). At some locations the BBL and wind-driven Ekman transports have been shown to be in balance (e.g. Perlin et al, 2005) but only when the interior current is weak (e.g. Smith, 1981). This has not been the case at other locations, e.g. Schaeffer et al (2014), who additionally report poor correlations between the BBL and wind-driven Ekman transports. Here we have shown the addition of a wind-driven Ekman transport term, driving a local 2-D downwelling circulation, contributes little to further explaining the daily variability of the bottom fluxes.
The power of our regression approach on model data has enabled us to disentangle the components of the BBL fluxes, for many locations all on the shelf-break and for a few years of data-which is not yet achievable through observations alone given their sparsity. We do not claim that our regression model and coefficients can be used to accurately predict BBL fluxes for observational data that do not extend to the near sea-bed. In practice the coefficients that we have established for our regression coefficients are not spatially fixed, and so it is not clear which to use for a predictive model. What we have shown is that any model of BBL fluxes that does not extend to the seabed should take into consideration both Ekman and mean-flow term, but a 2-D wind-forced downwelling term is not required.
The seasonal variability of the Shetland slope current inflow (i.e. flow parallel to the slope rather than across the slope) was investigated by Gould et al (1985) who reported a sinusoidal seasonal variability of inflow with maximum in winter and minimum in summer. Sherwin et al (2008) find a similar low in summer but consistently high inflow for the majority of the year. The seasonal variability of the Hebridean slope current was studied by Souza et al (2001). They find the acrossslope velocity is offshore and stronger in winter than in summer (over the entire water column and not just the BBL). They link this to changes in the windstress which become more off-shore favourable in the winter than summer. We provide for the first time a seasonal perspective of across-slope BBL transport on the Shetland shelf ( Figure 10; section 3.3) and find it also has a sinusoidal pattern with a winterhigh and summer-low. This is also reflected in the components of BBL transport (Ekman and mean-flow terms). Sherwin et al (2008) provide a time-series of south-westerly wind stress and attempt to connect it to the inflow variability. There is a clear similarity between their wind time-series and our BBL transport time-series. It may be that the wind influences the interior slope current on seasonal time-scales and therefore the across-slope BBL transport (a result of our regression analysis). However, this connection between the wind and the BBL transport via the interior current is implicit rather than explicit: we have shown that an explicit wind-driven 2-D downwelling term is not required in explaining the BBL transport (section 3.4.5 and 3.5).
The magnitude of BBL transport was previously estimated in a model study of the north-west European continental shelf by Holt et al (2009). They used a σ -coordinate model which ran for 44 years at ∼ 12 km resolution. For their section of the Shetland shelf they found, for the same 20 m BBL height, a similar Ekman transport (0.16 Sv vs. 0.15 Sv here). This is encouraging given one limitation with our transport estimates is that we have only simulated four years. However they do report less BBL transport (0.28 Sv vs. 0.53 Sv here). Differences between estimates could be due to inter-annual variability and the lateral extent of the isobath used. The correlation between single Ekman-drain models is also similar (R 2 = 0.24 vs. shelf-mean R 2 = 0.25 here). We show here that a large part of the unaccounted BBL transport is due to the V int component of the interior current which we estimate to be 0.31 Sv. This was previously not attributed.
For both the variability and magnitude, the orientation method used to define the across-and along-slope direction of the interior flow will affect flux estimates. An assumption for the Shetland slope current (e.g. Turrell et al, 1992;Simpson and McCandliss, 2013) is that the interior current direction (here θ T PT ) can be used as proxy for the true bathymetric direction (here θ LG ), and from that estimates of the across-slope transport can be made. If this is not the case (i.e. the mean flow direction is not parallel to the slope), then such estimates will be different to the true acrossslope transport. Souza et al (2001) find the Hebridean slope current is closely parallel to the isobaths. We find the Shetland slope current at the 200 m isobath is directed 5 • off-shelf, so nearly parallel but not exactly (section 3.2). We have shown that using a mean-flow orientation method (θ T PT ) dampens the mean-flow contribution to acrossslope BBL transport, whilst enhancing that of the Ekman transport (section 3.5). This effect must be better considered in future estimates.
The spatial scale of the Ekman drainage is approximately 80 km whereas the V int transport is more local (Figures 15 and 16). Brink and Robinson (2005) have previously stated that the decorrelation length of the along-slope current is larger than for the across-slope current. The large along-shelf coherence of the Shetland slope current (Figures 3 and 9) causes a large spatial coherence of the Ekman transport. For V int , the mean e-folding distance is about 5 km but within errors this is O(10 km) demonstrating a range of scales for this process. Sherwin et al (2006) analyse slope current meanders of O(10 km) and their contribution to shelf-edge exchange. They suggested it was likely these meanders break down to smaller eddies, hinting a range of scales for shelf-edge exchange. Our decorrelation plots for the V int and BBL transport are similar ( Figure 16) and it is apparent that the short scale of V int has consequences for BBL transport. It is precisely the localisation of V int transport that causes localisation of BBL transport. By capturing the mean-flow on smaller scales we capture more of the short-term variability in the BBL transport.
Utilising a multivariate regression approach, we have provided insight into the variability, magnitude and spatial scale of the various BBL transports. We now consider potential sources of local residual behaviour. Residuals lead to variability of the BBL transports around the best fit which decreases local correlation. Firstly, one of the primary assumptions of Ekman theory is that the interior ocean current does not accelerate: it is in steady state. Acceleration of the interior current, V int , may not correspond directly to the BBL transport. Stewart (2004) states acceleration of currents are important for horizontal scale of less than 50 km and for less than a few days. We attempt to remove the effect of acceleration by time averaging our currents but we find that the correlation does not improve with time averaging, making this an unlikely explanation.
One limitation in the multivariate regression is a larger veering layer in our simulation. If this increases into the assumed depth we have taken for the interior current V int (Section 2.3), then the two sides in our regression equation (Eq. 5) are not independent. A veering layer develops because of the turbulent transfer of momentum upwards from the boundary stress, parametrised by the vertical eddy viscosity: increasing this will increase the height of the veering layer. The basic Ekman model (Ekman, 1905) uses a constant eddy viscosity to derive a veering layer height. Cushman- Roisin and Malačič (1997) show in an unstratified boundary layer the eddy viscosity may vary proportionally to the boundary distance. It is also known that the veering layer height varies with vertical mixing and stratification, for both the surface (e.g. Lentz, 2001) and bottom layers (e.g. Perlin et al, 2005). Vertical mixing and stratification are linked: the presence of stratification inhibits vertical (diapycnal) mixing. Conversely mixing of light and denser bottom waters increases the vertical turbulence and decreases stratification by homogenising the waters. Increased vertical turbulence increases the vertical eddy viscosity, thus increasing the veering layer height. Stratification and veering have been considered in Perlin et al (2007) for the Oregon coast who report veering layers of around 20 m, and observed by Hosegood and van Haren (2003) in the FSC with veering up to 50 m above the seabed. Hosegood and van Haren (2003) additionally considered the effect of the slope on the veering height for the Shetland shelf (using Trowbridge and Lentz, 1991), providing a minimum value of 8.3 m with the strongest stratification and increasing inversely with the buoyancy frequency to potentially O(100 m). In our model, stratification profiles are reasonable compared to observations ( Figure 8). Furthermore, we have attempted to mitigate the assumption of a fixed veering layer of 20 m by taking the interior current as the depth-averaged current in the 80 m layer above the proposed BBL. We also tested our regression model at the 400 m isobath where it was also valid.
An additional complication is the difficulty of the MITgcm, a z-coordinate model, in allowing dense water to flow down-slope, as discussed by Legg et al (2006). They show that the amount of dense water overflowing in the MITgcm is dependent on the model resolution-coarse models generate excessive spurious mixing preventing dense water from descending. At intermediate resolutions the model produces less mixing and more dense water can descend, though less than their non-hydrostatic simulations run at the highest resolution. Intermediate resolutions also produced sim-ilar levels of mixing as compared to isopycnal models. Our model was run (section 2.1) between intermediate to high resolution using their benchmarks, and it is therefore likely to be suitable to simulate downslope flows.

Conclusions
We have presented a high-resolution model of the Faroe-Shetland Channel that simulates many features seen in observations (e.g. Sherwin et al, 1999;Berx et al, 2013). We have developed a multivariate regression model of across-slope BBL transport, to understand its daily variability, as a function of the along-and across-slope components of the interior flow. These terms correspond to an Ekman and mean-flow term. The inclusion of the mean-flow term represents an improvement over previous Ekman-only models of the BBL transport Simpson and McCandliss, 2013). We also reject 2-D wind-forced downwelling as a variable in explaining the daily variability of the BBL transport. We also attribute a greater portion of the BBL transport budget to the mean-flow compared to the Ekman transport. A previous study with similar Ekman transport (Holt et al, 2009) had not done this.
For the first time we have presented a seasonal cycle for the across-slope BBL transport and found it has a winter-high and summer-low, similar to the interior flow of the Hebridean slope current (Souza et al, 2001). Sherwin et al (2008) connects the wind on seasonal timescales with the seasonal inflow variability of the slope current. From our study the interior slope current directly explains the BBL transport variability. We have also rejected the explicit role of the wind, in the form of a 2-D downwelling term, in explaining the variability of the BBL transport. Therefore if the wind does play a role in explaining the seasonal BBL transport, this connection is implicit and via a modification of the interior current.
We have estimated the average scale of the Ekman transport to be over 80 km, near-grid scale (< 5 km) for the mean-flow and the BBL transport, and sub-grid scale for the residual behaviour. Local mean-flow is therefore important in determining shelf-edge exchange.
Our regression method is general and may be used diagnose the contribution of Ekman, mean-flow and the wind to across-slope BBL fluxes in models. In the simple Ekman-drain model it is only the change in interior current speed that explains the BBL flux variability. Here we show that for this region the variation in direction of the interior flow is of greater importance in understanding the BBL transport.