A tidally driven estuary close to an amphidromy

This paper describes the implementation of a high-resolution three-dimensional model of the estuary “Sundalagið Norður”. This estuary is bound by narrow sills in both ends and has a large semidiurnal tidal variation. The proximity to an amphidromic region, results in a periodically varying difference in sea level height between both ends of the estuary, which generates strong semidiurnal tidal currents across the sills. The available observational data support the model results. The model results verify the dominance of tidal forcing with semidiurnally varying volume fluxes into and out of the estuary. The 5 amplitudes of these fluxes vary in strength with the fortnightly variation between spring and neap tides. More surprisingly, the model also indicates a strong fortnightly variation of net fluxes averaged over 25 hours reducing both diurnal and semidiurnal tidal currents. These variations are caused by fortnightly variations in sea level difference between both ends of the estuary, which are verified by comparison with observed sea level variations. This rather surprising result implies that exchanges within the estuary and with its surroundings vary systematically; typically with one week of net northward flow followed by one week 10 of net southward flow. This variation also appears to affect the mixing processes in the estuary and should be taken into account in planning development or activities. More observational data would be beneficial to validate the model more thoroughly and we recommend that a dedicated experiment with combined observations and numerical modeling is implemented.

2 Material and methods

The model
In this paper we have applied a model setup based on the open-source Regional Ocean Model System (ROMS, http://myroms.org, Shchepetkin and McWilliams, 2005;Haidvogel et al., 2008). This is a state-of-the-art three-dimensional hydrostatic, free-  2013)), the first nesting has a resolution of 800 m in the horizontal and is further described in Erenbjerg et al. (2020).
The second nesting contains a 160 m horizontally resolved grid and is run for five months in 2013. This second nesting is used 75 as forcing for the ultra-high-resolution setup used in our current study.
Atmospheric forcing is provided by the Weather Research and Forecasting (WRF) model on the surface. The WRF-model is set up with a configuration that has a resolution of 9-3-1 km in the horizontal and the area with 1 km ⇥ 1 km resolution covers the entire Faroe Islands. More detail on configuration can be found in Myksvoll et al. (2012).
The freshwater input to the estuary from runoff is assumed to be constant in time and based on two reports: Erenbjerg (2020) 80 and Davidsen et al. (1994) as well as data from the local energy supplier. The two freshwater sources are a hydropower plant on the eastern coast and a river on the western coast (Fig. 1b). Their average annual supplies are 1.7⇥ 10 8 m 3 and 6.3 ⇥ 10 7 m 3 , respectively. The entire annual freshwater input from runoff is averaged and divided into a constant daily input to ROMS. The computationally demanding model was run from the 11th of February until the 12th of March 2013. To avoid potential spin-up effects, results from the first day have been omitted. Model output was stored every hour and includes horizontal 85 velocities, temperature, and salinity for each grid cell (Arakawa C grid). Sea level height was also stored for each surface grid cell.

Observations
During the modeling period, two upward-looking Acoustic Doppler Current Profilers (ADCPs) were deployed on the bottom of the estuary (Fig. 1b). Details are documented in Larsen et al. (2014a, b). From 2012 to 2018 there were a high number 90 of CTD (Conductivity Temperature Depth) observations in the estuary, most of them documented in Simonsen et al. (2018)).
Unfortunately, no hydrographic observations were made during the modeling period. Sea level measurements are available at two sites "Eiði" and "Tórshavn" (Fig. 1a) from the Faroese Office of Public Works (Landsverk). These measurements have many gaps, but they are sampled every ten minutes, which allows generation of the hourly averages that are used here.

Current velocity validation
The two ADCPs were located on a transect crossing the estuary (Fig. 2a), but quality data were only obtained for the deep parts of the velocity profiles (Larsen et al., 2014a, b). Modelled and observed average cross-estuary profiles are similar and close to zero ( Fig. 2b and 2c). For mooring site BW (Fig. 2b), the modelled and observed along-estuary profiles are also fairly similar for the depths reached by the ADCP. For mooring site BE (Fig. 2c), the discrepancy is larger, but Fig. 2a indicates that 100 this site is close to a transition zone between northward and southward flow at depth, which makes the comparison sensitive to modeling details.
Instantaneous velocities in the estuary will depend highly on the phase of the tidal signal. In the model, this is determined by the forcing along the outer boundaries. Thus, a model-observation comparison of instantaneous velocities is not very meaningful. Instead, the top panels in Fig. 3 compare observed and simulated along-estuary velocities averaged over consecutive 105 25-hour periods (to average out the tides). The observed and simulated velocities are certainly not identical, but they are generally of the same magnitudes. All of the correlation coefficients are also positive and three out of four are significantly higher than zero at the 95 % level (p < 0.05).
On daily time scales, the model is able to reproduce velocity variations to a significant extent. A similar conclusion can be reached for shorter time scales including the tidal velocity variations, which is documented in the bottom panels of Fig. 3. For 110 both the top and the bottom panels, the largest discrepancy between observations and model is for 18 m depth at mooring site BW. From Fig. 2a it appears that this is close to the average interface between northward flow in the upper layers and southward flow at depth. Thus, details of model topography, setup, and forcing may affect the comparison for this case disproportionally.

Salinity validation
As for other estuaries, a key feature of our estuary is the interplay (mixing, transport) between the freshwater added from land 115 or the atmosphere and the saline seawater entering the estuary over the two sills. This interplay is especially reflected in the salinity distribution. To assess the model performance for this parameter, we collected the data from all the CTD observations in the estuary 2012 -2018 (details are documented in Simonsen et al. (2018)).
In the shallow regions on either side of the estuary, detailed bottom topography and proximity to a river outlet may affect the salinity disproportionately. We therefore considered only CTD casts with bottom depth at least 50 m. To exclude situations with 120 a stagnant bottom layer, we only used observations from winter (November -April). This gave a total of 92 salinity profiles, which are shown as the thin black lines in Fig. 4.
The semi-transparent red area in Fig. 4 shows the range of daily averaged salinity profiles from model grid cells that fulfill the same criterion for bottom depth. One clear difference is that the observed salinities are lower than the range from the model even at depth. This is not surprising since the salinity of the shelf water decreased considerably from 2013 to 2018 (Erenbjerg 125 et al., 2020). The decreasing salinity cannot explain the large differences in the upper layers where the observations show that the water often was much fresher than indicated by the model. Still keeping in mind, that the two data sets have different conditions. The model was run over one month with nearly constant freshwater supply, whereas the observations are over a six-year period with variable precipitation, river runoff, hydropower outflow.

Model results
There are no indications of serious spin-up effects, but for precautionary reasons we omit the first day of the simulations. Thus, the model output comprises 672 hourly values (from 12 February 01:00 to 12 March 00:00) of velocity and hydrographic parameters for each grid cell as well as sea level. This period will in the following be referred to as the "simulation period".
We define our "estuary" to be the region between the two sills. The sills are defined by minimum cross-sectional (east-west)   S3) give an indication of the estuarine characteristics, but may be rather misleading due to the large and systematic changes through the period.

Hourly variations 140
On hourly time scales, the flow through the estuary is clearly dominated by the tides. This is evident in Fig. 5a, which shows the hourly variations of volume fluxes across both the northern, q N (t), and the southern, q S (t), sill, during the 4-week simulation in 2013, where positive values indicate northward volume flux and negative values indicate southward flux. The semidiurnal variation is clear in the hourly output, as is a fortnightly variation in the amplitude of the flux. The amplitude of the flux across the northern sill is much higher than the amplitude in flux across the southern sill. From a regression analysis, the highest 145 correlations are found when q S (t) lags one hour after q N (t), and the amplitude ratio is then 0.46 (Table 1).
Typically, the volume flux is in one direction for six hours before changing sign. Adding up all the water flowing into or out of the estuary across the northern sill during one of these tidal phases, we find that this typically is around 10 % of the volume, but occasionally the phases may last longer and transport more water, such as the period by the end of the first week of simulation where an episode of excess southwards flow flushes 30 % of the volume out ( Fig. 5b).

150
Sea-level variations also follow the tidal cycle. We have defined four hourly sea level time series: h N (t) is sea level north of the northern sill. h S (t) is sea level south of the southern sill. h IS (t) is sea level just north of the southern sill. h I (t) is sea level averaged over all surface grid points between the two sills. In addition, h S (t) is the sea level difference across the southern Supplementary Fig. S4).
Correlation coefficients (R) between some of these parameters are listed in Table 1. Average sea level between the two sills, 155 h I (t), is very highly correlated with sea level north of h N (t), with zero lag (less than one hour) and a regression coefficient (↵) close to one. Even h IS (t), just north of the southern sill ( Supplementary Fig. S4), follows h N (t) almost identically. Thus, sea level within the estuary responds more or less instantaneously to the sea level north of the estuary. Sea level south of the estuary, h S (t), is also highly correlated with h N (t), but with a smaller regression coefficient (Table 1).
Hourly values for volume flux across the two sills are fairly well correlated with the difference in sea level between both 160 ends of the estuary, h N (t) h S (t), consistent with the idea that this difference drives the flow through the estuary. The volume flux across the southern sill, q S (t), is better correlated (R=-0.95) with the sea level change across the southern sill, h S (t) ( Supplementary Fig. S4), but an even better fit is obtained by a quadratic fit, rather than a linear one (red curve in Fig. 6).   Table 1. Lagged correlation and regression analysis of relationships between hourly values of various model parameters. R0 is the correlation coefficient for zero lag. Lagmax is the lag (in hours) that gives the numerically highest correlation coefficient. Rmax, ↵max, and max are the correlation coefficient and the regression coefficients for that lag.

Regression equation
The quadratic fit (red equation in top right corner of Fig. 6) is motivated by a simple model for the flow across the southern sill, in which the speed at any instant, v S (t), is assumed to be homogeneous on the section. If friction is ignored, the Bernoulli  locations for sampling h S (t) (Supplementary Fig. S4) were chosen arbitrarily. The value for was derived by fitting the model data into this simple framework. When friction is included, the speed and volume transport decrease and the value for 170 in the fit changes ( Fig. 6) this indicates that a considerable part of the potential energy is lost to friction .

Fortnightly variations
The fortnightly variation of tidal flux amplitude seen in Fig. 5a is a normal phenomenon in Faroese waters (Hansen, 1978), reflecting the variation between neap and spring tides. However, the fortnightly variations in the 25-hour average fluxes (thick lines in Fig. 5a) are unexpected. The simulation period includes two full periods with net southward flow and one period with 175 net northward flow, all of them lasting about one week (Fig. 7a). Both the first days and the last days during the simulation have a net northward flow. This figure also has a curve (black) that shows the variation of the standard deviation of q N (t) over consecutive 25-hour intervals. This curve should reflect the strength of the tidal amplitude and it is worth noting that the 25-hour average flux does not reach its extremes during spring tide, but rather a few days out of phase (Fig. 7a).
To help understand the fortnightly flux variations, Figure 7b shows 25-hour averages and standard deviations of sea level 180 height. During the periods with average southward flux, average sea level is higher north of the estuary. When the average flux is northward, the average sea level is higher south of the estuary. Consistent with Table 1, the volume flux through the estuary may be seen as forced by the sea level difference between both ends. In this paradigm, the reason for the fortnightly variations in volume flux is the variation in this sea level difference. Figure 8 shows the selected two full periods of northward and southward flow, respectively. Period 1 is the first full period 185 with average northward flow, whereas Period 2 is the following period with average southward flow. To illustrate the differences between these two periods, northward velocity and salinity are averaged over each period and across (east-west) the estuary and then plotted against northward grid number and depth in Fig. 8.
During period 1 (upper panels in Fig. 8 From the average salinity during period 1 (upper panels in Fig. 8), it seems that the entering seawater is highly diluted on the way into the basin. Only a very small part of the entering seawater in Period 1 is undiluted once reaching the bottom of the 195 basin. For the upper most layers in the basin the salinity is low, and it appears that freshwater is added from the south.
During Period 2 (lower panels in Fig. 8), the average velocity is southwards throughout most of the estuary. This implies, that the southern part of the sound receives a net inflow of water from our estuary. Except for some minor areas close to the  bottom, the only northward flow during Period 2 is near the surface over and just south of the northern sill. This appears to be freshwater from the river and hydropower plant in the northern part of the estuary (Fig. 1b).

200
The average salinity distribution during Period 2 (lower panels in Fig. 8), indicates that seawater enters the estuary over the northern sill and fills it up from below. Thus, almost undiluted seawater fills the basin below 30 m depth. Almost no freshwater seems to be mixed below 20 meters depth during this period in contrast to Period 1 where a large amount of the freshwater is mixed into the water column.

Observational validation of fortnightly variations
The fortnightly variations in average fluxes and sea level differences (Fig. 7) were, as mentioned, unexpected. To check that these variations are not purely an artifact generated by the model, we have compared observed sea level at "Eiði", which is just 210 north of the estuary with the observed sea level in Tórshavn, which is south of it. When averaged over 25 hours, the sea level at both sites is dominated by atmospheric pressure variations ( Supplementary Fig. S5). This effect is strongly reduced when the average difference in sea level between both sites is considered, which instead varies consistently with the neap-spring variation of the amplitude of the semidiurnal and diurnal tide (Fig. 10).
For zero lag, the daily (25-hour) averaged sea level difference between the two sites is positively correlated with the tidal 215 amplitude of sea level variations in Eiði with a correlation coefficient R = 0.50, which is significantly different from zero at the 98 % (p <0.02) level. The correlation increases slightly for a lag of 1 day (Fig. 10a), which implies that the 25-hour average sea level difference (north minus south) is highest shortly after spring tides. This is consistent with Fig. 7b, although the lag there seems somewhat larger. Remarkably, the (numerically) highest correlation, R = -0.55, is found when the sea level difference lags after the tidal amplitude by 8 days (Fig. 10a), i.e., shortly after neap tide. This supports our conclusion that the observed 220 fortnightly sea level variations are real and not an artifact of the averaging.

Density inversions downstream from the northern sill
From the salinity distributions in Fig. 8 right, it seems that the seawater entering the estuary across the northern sill is flowing downwards and is then mixed with the water inside the estuary just downstream of the sill. To study this in more detail, the density structure was plotted along a track (Fig. 11a). The track was chosen based on the average maximum velocity close 225 to the bottom, beginning on the shallow part of the sill (track number 0) down-slope to 55 m depth (track number 35). The average velocity during the entire simulation period is more than 10 cm s 1 southwards in the bottom layer of the steep slope of the track (Fig. 11d).
As long as turbulent mixing is weak, the density structure is normally stable with density increasing downwards. Density inversions with density decreasing downwards may therefore be used as a sign of mixing. To utilize this, we define a "density 230 inversion" as the density difference between the third lowest layer and the bottom layer (upper minus deeper) when this value is positive the density is inversed. Under stable conditions, when this value is negative, it is set to zero. The spatial and temporal variations of density inversions along the selected track are illustrated on the Hovmøller diagram (Fig. 11b). If the density inversions are averaged over the track and temporally smoothed, they exhibit systematic variations (red curve in Fig. 11c). Most pronounced is the rather frequent occurrence of inversions in the first part of the simulation from day two to nine. From day nine to 15, the density inversions are less frequent for most part of the track except for track points 20-25.
These variations show some similarity to the variations in Fig. 7a and verified by the blue curve in Fig. 11c. Strong southward flow across the northern sill seems to be a necessary -but not sufficient -condition for density inversion.
When averaged over time, the density inversions mainly occur on the down-hill slope, and are maximal right before the bottom slope levels off at track number 20, Fig. 11e. This could indicate hydraulic jump conditions, although we use a hydrostatic 240 model and cannot conclude on any vertical advection.

Model performance
The comparison between average velocities in the model and in observations (Sect. 3), did not show identical values, but taking into account the spatial variation (Fig. 2a) the corrrespondence is as good as may be expected. For temporal velocity variations, 245 most of the corrrelations between model and observations (Fig. 3) were significant at the 95 % (p <0.05) level, which again is encouraging. No hydrographic observations were made in the estuary during the simulation period.

Hourly variations
On short time scales, the flow through the estuary is clearly of tidal character with semidiurnal dominance (Fig. 5a) and we expect sea level differences between both ends to be the main driving force. This is supported by the high corrrelations between 250 volume fluxes across the sills, q N (t) and q S (t), and sea level difference h N (t) h S (t) ( Table 1).
As the tidal wave enters the estuary from the north, it should propagate southwards as a barotropic wave with sufficient speed to pass through the estuary in less than an hour. This is verified by the high zero-lag corrrelations between h N (t), h I (t), and h IS (t) ( Table 1). Most of the sea level change occurs over a fairly short distance over the southern sill and this is where the highest speeds are observed ( Supplementary Fig. S6).

255
Speeds exceeding 2.5 m s 1 might seem excessive, but according to local sailors, the speeds over the southern sill may occasionally be considerably higher. From the arguments supporting the quadratic fit in Fig. 6, it appears that a considerable part of the potential energy in the sea level is lost to friction over the southern sill. Thus, the flow across the southern sill probably controls how much water the tides can push through the estuary, but this will also make the results sensitive to turbulence and bottom friction parameterization in the model.

260
To reach the high values for speed ( Supplementary Fig. S6) and volume transport (Fig. 5a), strong forcing is required.
This forcing is achieved through the periodically varying large differences in sea level between both ends of the estuary and they may be linked to the amphidromic character of the region. The "textbook" amphidromy is a point in the middle of an ocean basin, around which the tidal wave turns. Here, the islands complicate matters, but Hansen (1978) suggested that an amphidromic region for the semidiurnal tides was located somewhere in the proximity of Tórshavn based on tidal analysis of An alternative explanation may be linked to the tidally rectified residual currrents that circulate the Faroes over the shelf (Larsen et al., 2008). These currrents vary in strength over the spring-neap cycle and would be expected to generate sea level slopes through geostrophy that should have a fortnightly variation. Conceivably, these variations may generate the fortnightly variations in sea level difference between both ends of our estuary that are needed to explain Fig. 7.
On the available evidence, none of these explanations can be validated more quantitatively. However, the effect is important, 305 especially for fish farming activities, and may affect other Faroese sounds. A larger model domain and more comprehensive sea level and currrent observations would be beneficial to increase our knowledge. This is planned to be performed in a follow-up study, but will be too voluminous for this manuscript. Here, we note only that observations verify that there are fortnightly variations between sea level north of the estuary and south of it similar to those in the model. Possibly, the model exaggerates the effect or shifts its timing, but the reality of the fortnightly variations seems well justified.

310
Since southward flow across the southern sill occurs during flood, it has been suggested that the cross-sectional area over the sill -and therefore also volume flux -should be higher during southward than northward flow, which would lead to a net southward volume flux varying between 50 and 175 m 3 s 1 in phase with the strength of the tidal amplitude (VandKvalitetsInstitutet, 1983). This is not supported by the model, which has a net flux that varies between northward and southward direction on a fortnightly time scale (Fig. 7a). To check the effect suggested by (VandKvalitetsInstitutet, 1983), volume transport across 315 the southern sill was calculated with and without varying sea level. On average, the difference was 14 m 3 s 1 . Thus, the effect is real, but much smaller than expected and swamped by the fortnightly variations.

Exchange rates/ flushing rates
For an estuary that is affected by human activity, one of the most important parameters is the flushing rate, i.e., how fast the waters in the estuary (and dissolved contaminants or planktonic organisms) are flushed out of the estuary. An often used 320 measure of this is the "flushing time", defined as the volume of the estuary (or parts of it) divided by the volume flux into or out of the estuary. Combining CTD observations from our estuary with estimated freshwater supply, Hansen (1990) estimated a typical flushing time of 5 days for this estuary, but noted the uncertainty for this value.
From the present model results, there are several ways of obtaining alternative estimates. One way is to use Fig. 5b, which implies that between 5 and 10 % of the estuarine volume is typically flushed in and out over the northern sill every 12 hours. If 325 one assumes no mixing between in-flowing and out-flowing waters, this method gives an average flushing time of around one week, ranging between less than four days and more than 11 days ( Supplementary Fig. S7) over the simulation period.
Alternatively, one could use the average salinity distribution ( Supplementary Fig. S3) to estimate the average total freshwater content in the estuary and combine that with the (almost constant) freshwater supply to calculate a flushing time. This method is, however, very sensitive to the choice of seawater salinity that is diluted by freshwater and also assumes stationary conditions.

330
An extra challenge when trying to estimate a "typical" flushing time is the fortnightly variation in the net flow through the estuary (Sect. 5.3). Considering the two periods defined in Fig. 7, Figure 12 illustrates the two different exchange regimes, that the estuary regularly shifts between.
During Period 2, the seawater from the north seems to experience stronger mixing with ambient waters (Fig. 11c) and to be more spread out through the water column (Fig. 8). The deep arrow in Fig. 12b showing 629 m 3 s 1 represents net southward flow through the whole water column just south of the northern sill, but most of that is below sill level of the northern sill. This is evident from Fig. 9, which also illustrates the difference in deep flow and flushing during the two periods.
If the arrows in Fig. 12 are combined to make budgets, we see that they do not balance. During both periods, Fig. 12 indicates 355 more water leaving the estuary than entering it and the associated sea level changes would be unrealistically large. It has to be kept in mind, however, that the values in Fig. 12 are based on averages for each period with sea level kept constant. The volume flux estimates therefore have an uncertainty, which over the southern sill at least are on the order of 10 %. Within that uncertainty, there is balance. It might also be emphasized once again that the arrows in the figure are net volume fluxes, upon which are superimposed the semidiurnal variations with amplitudes that are considerably higher.

Recommendations
The fortnightly variations in net flow through our estuary revealed by the model may have large impacts on any development or activity that affects the estuary and its surroundings, such as hydropower development or fish farming. They should therefore be considered when planning future operations. We also recommend a more thorough validation of the model. Although the model results show a satisfactory agreement with the available observational data, these data were not targeted at some of 365 the key processes that have appeared. We therefore recommend that a more targeted field experiment is implemented that can reveal more of the strengths and weaknesses of the model.
We also recommend investigating the mixing across the northern sill, to see if there is actually hydraulic jump conditions, with a model that is not using the hydrostatic assumption, thus beeing able to make conclusions on vertical advection.
Code availability. The source code of ROMS is available from ROMS, http://myroms.org.

370
Data availability. Landsverk data of sea level are available upon request (www.landsverk.fo). Data from DMI surface air pressure and precipitation are available upon request (www.dmi.dk). Current measurements and CTD observations are available upon request from SVE.
The data from local energy supplier (SEV) are available upon request (www.sev.fo). The model output from the FC32m model are available from SVE upon request.