Modelling the Influence of Riverine Inputs on the Circulation and Flushing Times of Small Shallow Estuaries

Abstract

Simple flushing time calculations for estuarine systems can be used as proxies for eutrophication susceptibility. However, more complex methods are required to better understand entire systems. Understanding of the hydrodynamics driving circulation and flushing times in small, eutrophic, temperate estuaries is less advanced than larger counterparts due to lack of data and difficulties in accurately modelling small-scale systems. This paper uses the microtidal Christchurch Harbour estuary in Southern UK as a case study to elucidate the physical controls on eutrophication susceptibility in small shallow basins. A depth-averaged hydrodynamic model has been configured of the estuary to investigate the physical processes driving circulation with particular emphasis on understanding the impact of riverine inputs to this system. Results indicate circulation control changes from tidally to fluvially driven as riverine inputs increase. Flushing times, calculated using a particle tracking method, indicate that the system can take as long as 132 h to flush when river flow is low, or as short as 12 h when riverine input is exceptionally high. When total river flow into the estuary is less than 30 m³ s⁻¹, tidal flux is the dominant hydrodynamic control, which results in high flushing times during neap tides. Conversely, when riverine input is greater than 30 m³ s⁻¹, the dominant hydrodynamic control is fluvial flux, and flushing times during spring tides are longer than at neaps. The methodology presented here shows that modelling at small spatial scales is possible but highlights the importance of particle tracking methods to determine flushing time variability across a system.

Introduction

Despite the introduction of European directives to improve ecosystem health, surface water quality across Europe remains poorly controlled with only 53% of European Union water bodies having attained the Water Framework Directive’s (WFD) required ‘good’ status by the appointed deadline (European Commission 2012). Due to attention from the WFD and Marine Strategy Framework Directive (MSFD), there has been increased awareness and research into eutrophication and its consequences; examples include estuaries in Denmark (Nielsen et al. 2003; Krause-Jensen et al. 2005), along the Portuguese coast (Ferreira et al. 2006; Cabrita et al. 2015; Caetano et al. 2015) and along the north coast of Spain (Borja et al. 2004; Ménesguen et al. 2018). However, there remain regions across Europe where, despite frequent data collection of water quality parameters, there has been a lack of published data synthesis. In particular, many of the small estuaries of southern England have not been subject to recent published assessment of water quality. This paper considers the implementation of hydrodynamic and particle tracking models as a means of analysing potential water quality issues in a small south coast shallow estuary, Christchurch Harbour, that has been reported to be subject to eutrophication.

Of critical importance to the development of eutrophic conditions is the ability of the water body to expel the cause of the problem. The flushing time of a system reflects its ability to recover from a period of reduced water quality, providing a unique value for the entire water body (Bárcena et al. 2012). Flushing time is an important physical component of an estuary as the length of time water resides in a basin can influence both chemical and biological processes (Anderson et al. 2003; Lee et al. 2011; Rynne et al. 2016), thus making the water quality dependent on the flushing time of the system (Choi and Lee 2004; Defne and Ganju 2015). In theory, a well-mixed system with a short flushing time is assumed resistant to extended periods of low dissolved oxygen (DO); however, Verity et al. (2006) have shown that even small estuaries characterised by short flushing times have observed gradual decreases in DO and eventual hypoxia when placed under continuous anthropogenic stress.

Due to difficulties in calculating estuarine retention times, numerical models are often used as the primary and most reliable method for estimating flushing durations. Numerous methods have been proposed with some calibrated modelling techniques coming into wide use (Sanford et al. 1992). Many of these complex modelling techniques are unsuitable for small relatively simple estuaries (i.e. a one entrance and exit system with non-complex bathymetries and tides). In such cases, a more simplified modelling approach, in which flushing capabilities are underestimated, would be more appropriate (Sanford et al. 1992). However, there will be cases in which small estuaries require complex modelling due to one or more non-simple factors influencing the system. As such, the method used to calculate flushing time should depend on the complexity of the water body in question (Herman et al. 2007). Traditional simple non-modelling methods include the tidal prism or salt balance concepts (Ketchum 1951; Dyer 1973) with newer mathematical methods such as particle tracking models also being used.

For large estuaries, such as the Baltic Sea (Meier et al. 2011) and the Gulf of Mexico (Fennel et al. 2011), extensive ecological data and complex models are available to provide a comprehensive picture of surface water quality in the past, present and future; however, in most small, shallow estuaries, eutrophication has not been as well documented resulting in limited data collection and a lack of configured numerical models (Nezlin et al. 2009; Evans and Scavia 2013). To improve understanding of eutrophication susceptibility in less well-studied systems, smaller estuaries and water bodies are now under scrutiny with modelling being an ideal method to analyse water quality. To elucidate the processes that underpin eutrophication development, the physical drivers influencing circulation must first be understood. The rationale of this paper is to apply hydrodynamic modelling to predict the susceptibility of a small, shallow estuary to declines in water quality under varying riverine flow conditions.

The overall aim of this paper is thus to assess how physical drivers influence circulation and water flushing times in Christchurch Harbour, a small shallow enclosed microtidal estuary on the south coast of England. The aim will be achieved by addressing the following two objectives:

• 1. To assess the influence of river flow and changing tidal conditions on circulation patterns in the Christchurch Harbour estuary; and

• 2. To investigate the impact of river flow on flushing times in Christchurch Harbour under a range of tidal conditions.

The microtidal Christchurch Harbour estuary is principally fed by two rivers, the Hampshire Avon and the Stour, and meets the English Channel through a 47-m-wide narrow opening known as the Run (Fig. 1). As a tidally driven estuary, salinity within the harbour depends greatly on the state of the tide, resulting in fresh conditions at low tides in the winter and near fully saline conditions at high tides (ABP 2009). Christchurch Harbour is a shallow and mostly well-mixed system (Gao and Collins 1994; ABP 2009) with a history of oxygen under-saturation in warmer summer months when river fluxes are low and average flushing times long (Murray 1966). Although not well studied in the past, the estuary has recently been part of an extensive monitoring programme by the Christchurch Harbour Macronutrients Project (British Oceanographic Data Centre 2017), established to improve understanding of macronutrient behaviour in the Hampshire Avon, River Stour and Christchurch Harbour over a range of spatial and temporal scales. This monitoring programme detected summer low oxygen saturations and high chlorophyll concentrations in the upper estuary. The observational data sets obtained from both a discrete sampling programme (Panton et al. 2018; Panton et al. 2020) and automated continuous monitoring (Beaton et al. 2017) conducted during the project (between 2012 and 2017) will be most beneficial for use in model calibration, validation and parameterisation. Of particular interest is the interaction between tidal and fluvial dominance in small shallow estuaries and as such is something we look into in this system. Due to its comparatively high river flow, it is possible that under certain conditions, the small tidal amplitude of Christchurch Harbour may be overwhelmed by the riverine input, leading to the appearance of a constantly ebbing tide. We would define this as the point at which the usually tidally driven estuary becomes fluvially driven. Furthermore, we propose that interaction between riverine inputs and the volume of water within the system (i.e. spring/neap tidal cycle) will have a significant impact on circulation and can provide information on controls within a small system due to the flushing effectiveness of the tide compared to the fluvial flows.

Fig. 1

Christchurch Harbour and the surrounding areas. Located in the centre of the UK’s south coast, Christchurch Harbour exits into Christchurch Bay in the English Channel as shown in a. b The sample location at the ferry pontoon where a RBR duo pressure sensor was deployed (50°43′11.69″N 001°44′37.16″W)

This paper is structured as follows. “Model Configuration and Validation” describes the hydrodynamic and particle tracking model configuration and hydrodynamic model validation. The simulations undertaken and details of flushing time methods used are described in “Methodology”. Results are presented in “Results” and discussed in “Discussion”. Finally, conclusions are given in “Conclusions”.

Model Configuration and Validation

Two different models are used in this study: a hydrodynamic model and a particle tracking model. The configuration and validation of these are described in the next two sections.

Model Configuration

Hydrodynamic Model

The depth-averaged hydrodynamic model was configured using the MIKE 21 Flexible Mesh (FM) suite of modelling tools. MIKE 21, developed by the Danish Hydraulic Institute (DHI), is a two-dimensional modelling tool for use in estuarine, coastal water and sea environments (DHI 2017b). MIKE 21 Hydrodynamic FM is the basic module for use in the field of free-surface flows (DHI 2017b; Fitri et al. 2017), simulating water level variations and flows in the computational area in response to pre-defined forcing functions. Flows are calculated in the x and y direction based on vertically integrated equations of conservations of volume and momentum (Warren and Bach 1992; Fitri et al. 2017). The MIKE 21 modelling system is based on a cell-centred finite volume method on an unstructured mesh (DHI 2017b) allowing for graded mesh resolution making it an ideal modelling framework for shallow estuarine environments (Robins et al. 2011). For simplicity, we consider a two-dimensional approach, justifiable for this case study of Christchurch Harbour due to it being mostly well-mixed, small, shallow and microtidal. This is similar to methods applied by Umgiesser et al. (2004) and Cucco and Umgiesser (2006).

The model mesh, configured using the mesh generator tool in MIKE Zero, includes the estuary, the River Avon and River Stour up to their respective tidal limits and a large coastal area outside of the estuary (Fig. 2). The unstructured grid, based on linear triangular elements, allows for variable resolution with regions such as rivers and narrow passages to be well resolved whilst off-shore regions utilise larger grid-spacing to maximise computational efficiency (Robins and Davies 2011; Robins et al. 2011). The grid was created such that the resolution increases from 1000 m in the off-shore regions to less than 5 m inside the estuary.

Fig. 2

Model grid and bathymetry. The model area and entire grid are shown in a. The thick line represents the open tidal boundary. b The model grid within the estuary. The bathymetry of the entire domain and bathymetry of the estuary are shown in c and d, respectively. The location of the 507 particle sources (filled circles) used for the particle tracking model (described in "Methodology") are also shown in d

Bathymetry data was interpolated onto the mesh. The bathymetric data implemented was a compilation of the following five datasets: (1) data inside the estuary was taken from a 1980s Admiralty Chart (UK Hydrographic Office 1982); (2) a small bathymetric survey completed in 2000 around the ferry pontoon with a resolution of < 1 m (from Stuart McVey, personal communication, May 5, 2017); (3) LiDAR data, covering shallower regions of the estuary; (4) hydrographic data in the Run and surrounding coastal regions at a resolution of 1 m, courtesy of Channel Coastal Observatory (CCO; downloaded from www.channelcoast.org); and (5) bathymetric data from GEBCO (General Bathymetric Chart of the Oceans; downloaded from the British Oceanographic Data Centre, www.bodc.ac.uk) at a resolution of 20 m for outside of the estuary. As bathymetric data within the Hampshire Avon and River Stour was not available, a constant depth of 2 m was applied throughout the length of both rivers. A uniform bed resistance of 0.032 m^1/3 s⁻¹ using Manning’s formula was applied to the entire domain.

The open boundary was driven using water levels derived from the OTIS (Oregon State University Tidal Inversion Software) Regional European Shelf 1/30° tidal model (Egbert and Erofeeva 2002). The European model provides amplitudes and phases for the eight main diurnal and semi-diurnal tidal constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁ and Q₁) and three non-linear constituents (M₄, MS₄ and MN₄). The harmonic constituents (obtained from http://volkov.oce.orst.edu/tides/ES.html) were implemented to predict the tide using the Tidal Model Driver Matlab routines (Padman 2004; Haigh et al. 2011). Twenty-one points along the tidal boundary were selected and forced with OTIS tidal level data.

Two riverine water sources, with tracers, were added to the model at the tidal limits of the River Avon and River Stour. These water sources inject freshwater into the model, replicating interaction between the saline water outside of the estuary and the fresh river water. The riverine point sources were forced with both real data (shown in Fig. 3) and constant values. River flow data was obtained from the National River Flow Archive (downloaded from https://nrfa.ceh.ac.uk/data). Gauging stations on the River Stour, located at Throop, and the River Avon, located at Knapp Mill (Fig. 1a), have been continuously monitoring freshwater flows since 1973 and 1975, respectively, providing data at 15 min intervals. The River Mude (seen in Fig. 1) also flows into the estuary; however, these flows are negligible compared to the two main rivers (< 1% of annual flow in Stour) and thus were not included in the model.

Fig. 3

River flows in the Hampshire Avon and River Stour. Time series of river flows into Christchurch Harbour from the Hampshire Avon at Knapp Mill (since 1975) and the River Stour at Throop (since 1973) in a and b, respectively. The cumulative percentage flows for the Avon and Stour are presented in c and d, respectively. The thin lines are annual values and the bold line is the average since records began

These simulations were carried out on an individual PC with a 3.6 GHz processor and 16 GB of RAM. Runs took 36 h to complete 30 days of simulation. To compromise between resolution and computational expense, results were saved every 15 min of simulation time.

Particle Tracking Model

A particle tracking model was configured using the MIKE 21 Particle Tracking (PT) module. This module is a Lagrangian solver enabling the simulation of transport of dissolved and suspended substances in the water environment through a discrete-parcel method (DHI 2017a). The PT module uses the hydrodynamic simulation results on a regular grid with a spacing of 10 m (converted from the triangular grid in MIKE Zero) upon which particles, either conservative or non-conservative, can be applied at any time or location. The particle path can then be tracked across the domain. Each particle applied was released on the surface as a point source. As conservative particles were being used, the decay, settling and erosion settings were switched off. Horizontal dispersion sensitivity tests were carried out. Constant horizontal dispersion values were applied (ranging from 0.001 to 0.04 m² s⁻¹). Results indicated that the model was not sensitive to horizontal dispersion, and as such, the default value of 0.01 m² s⁻¹ was thus chosen for all simulations discussed in this paper.

Validation

The model was validated against water level data collected using a bottom mounted RBR duo pressure sensor deployed at the ferry pontoon, near the entrance of the estuary (see Fig. 1b), between June and September 2017, which sampled every 60 s.

Each simulation was run for 30 days from 1st July 2017 00:00:00 to the 31st July 2017 00:00:00. July 2017 was the chosen modelling period as this coincided with the hydrodynamic validation data. The first 2 days of the simulation results, during the model warm up period, were removed and excluded from validation. Statistical tests were used to assess model performance in accurately reproducing water level at the pontoon. To evaluate the model skill, the root mean square error (RMSE), mean absolute error (MAE) and standard deviation of the absolute errors (STD AE) were calculated. The correlation coefficient was also derived.

Comparisons of the predicted and observed water level within the estuary are shown in Fig. 4. The results show good agreement within the estuary with a RMSE of 0.09 m. The model accurately captures the complex tidal regime of the area with double high waters and young flood stands clearly predicted. The MAE is low (0.08 m) as is the STD AE (0.06 m). A strong correlation was observed with a coefficient of 0.96. Overall, the model does a good job of reproducing tidal levels across the domain and we can therefore confidently use the model to analyse physical controls within the estuary.

Fig. 4

Model validation. Model output and measured data for water level at the ferry pontoon for July 2017. Included is the root mean square error (RMSE), mean absolute error (MAE), standard deviation of the absolute error (STD AE) and the correlation coefficient (CC)

Methodology

Hydrodynamics

The first objective is to assess the influence of river flow and changing tidal conditions on circulation patterns in Christchurch Harbour. Using the configured hydrodynamic model, sensitivity tests were undertaken with constant river flows applied to represent observed ranges in the River Avon and River Stour (Fig. 3). In each case, the same value of river flow was applied to both rivers (i.e. 10 m³ s⁻¹ in each river, amounting to a total input into the estuary of 20 m³ s⁻¹). The values chosen represent average flows (10 m³ s⁻¹), intermediate flows (25 m³ s⁻¹), high flows which may be associated with regular storms and high precipitation (50 m³ s⁻¹ and 75 m³ s⁻¹), extreme flows such as those observed during large storm events (100 m³ s⁻¹ and 125 m³ s⁻¹) and a flow that corresponds to the uppermost values observed since recordings began (150 m³ s⁻¹). In addition to these values, simulations of no flow (0 m³ s⁻¹) were carried out as a control. To allow for comparison between the real river flow (validation) simulations and the constant river flow simulations, the same time period and tidal forcings were used.

Wind sensitivity tests were undertaken to examine the effect of wind speed and direction on current velocities within the estuary. Average wind magnitudes for the Dorset region (ranging from 1 to 10 ms⁻¹ which encompasses seasonal and annual averages as well as maximum gust averages; data from Met Office 2019) were applied over each of the main eight wind directions (0, 45, 90, 135, 180, 225, 270 and 315). In each case, water levels across the estuary for wind simulations were compared to water levels for no wind simulations. In addition, circulation patterns between wind and no wind simulations were compared. These tests indicated that wind had minimal effect on water levels and current velocities and circulation and was thus excluded from further simulations.

Flushing Time

The second objective is to investigate the impact of river flow on flushing times in Christchurch Harbour under a range of tidal conditions. To achieve this, the MIKE 21 Particle Tracking model was coupled to the hydrodynamic model and simulations with constant river flow, as described in “Hydrodynamics”, were carried out. For each river flow scenario, 507 point sources within the estuary (shown in Fig. 2d) were applied. Ten particles were released per source, amounting to a total of 5070 particles per simulation. The particles were released on a spring and a neap tide at four different stages of the tidal cycle: (1) flooding tide, (2) high tide, (3) ebbing tide and (4) low tide (exact times are listed in Table 1). Spring and neap tides were assumed to be the 7 days with the highest and lowest tides, with particles released on the maximum and minimum tides for these periods, respectively. For flooding and ebbing tide, the release time occurred when water levels were halfway between the preceding and subsequent high and low tides. These simulations were carried out independently from one another so that the influence of tidal phase on flushing time could be determined. The particles were tracked out of the estuary and flushing time calculated as the amount of time between release and the point at which only 37% (the e-folding flushing time) of the particles remain within the estuary. To gain an understanding of flushing time variability across the estuary (i.e. residence time), 507 boxes (the centre of each corresponding to the 507 point source locations) were considered within the estuary using a 55 m grid. The average amount of time taken for the ten particles in each box to leave the estuary was taken as the flushing time for that box. This provides local flushing times at each box throughout the estuary, similar to methods implemented by Jouon et al. (2006) and Iriarte et al. (2015). The use of the e-folding flushing time is appropriate if the number of particles remaining in the system over time decays exponentially. To confirm if this was an appropriate method to use for this case study, the number of particles in the estuary was observed overtime under a number of varying scenarios and found to decrease exponentially.

Table 1 Particle release times. Implemented in the particle tracking model, 5070 conservative particles were released at 507 point sources across the estuary during the spring and neap tides. In each case, particles were released at high, low flood and ebb tides

We propose the following method of identifying the controls on flushing times within the system. If the neap flushing time is taken away from the spring flushing time, values below zero will indicate that the system is controlled by the tide as neap tides should produce longer flushing times (Shaha et al. 2012; Sridevi et al. 2015; Haddout et al. 2020). Conversely, if the value is above zero, this will indicate that the system is controlled by riverine inputs. This method will be implemented to identify when control on circulation in the system changes from tidal to fluvial.

Although the main focus of this work is the use of a particle tracking model to calculate flushing time, a simple tidal prism calculation method has also been applied as a secondary method of estimating the flushing time of the estuary to provide a simple comparison with the more sophisticated model. Similarity between particle tracking and tidal prism methodologies would indicate that the tidal prism method can be used as an accurate first-order estimate of average basin flushing time. The tidal prism method was applied as follows:

$$ {T}_{\mathrm{fl}}=\frac{V}{\left({V}_{\mathrm{O}}+{V}_{\mathrm{R}}\right)\left(1-b\right)}t $$

where V is the volume of the estuary (assumed constant), V_R is the volume of river water, V_O is the volume of coastal water that enters the system during the flooding tide, b is the return flow factor, a fraction between 0.0 and 1.0 (Sanford et al. 1992; Monsen et al. 2002) and t is the period of the semi-diurnal or diurnal tide. This method was applied under the river flow scenarios described in “Hydrodynamics”. The assumption that water leaving on the ebbing tide does not return on the successive flooding tide can be bypassed with the inclusion of b (Officer 1976; Sanford et al. 1992; Monsen et al. 2002). We justify the use of this comparison method as it is derived from exponential decay and the particle tracking method found particles leave the system in an exponential fashion over time.

Results

Hydrodynamics

To understand physical controls on circulation, numerous sensitivity tests with different river flow values were undertaken. Riverine flows of 0 m³ s⁻¹ indicated very low and consistent current speeds across much of the estuary with values of 0.2–0.4 ms⁻¹ during both flood and ebb tides (Fig. 5). Within the Run, values are much greater, reaching 1.5 ms⁻¹ on the flooding tide and 0.8 ms⁻¹ during the ebb. Downstream flow rates at the narrow section of the estuary increase from values of 0.2 ms⁻¹ with river flows of 10 m³ s⁻¹ to 1.4 ms⁻¹ when flows were set to 150 m³ s⁻¹. At low riverine inputs, flow into the estuary through the Run travels up the deeper channel at values of 0.5 ms⁻¹. As fluvial inputs increase, current speeds in the channel during the flooding tide decrease where the freshwater and seawater meet in the middle of the estuary. The downstream flowing river water does not appear to travel within the deeper channel, instead spreading out over the entire estuary.

Fig. 5

Current speed and current speed difference during flooding and ebbing tides. Estuary current speeds under river flows in each river of 0 m³ s⁻¹ (a, b), 10 m³ s⁻¹ (c, d), 50 m³ s⁻¹ (e, f), 100 m³ s⁻¹ (g, h) and 150 m³ s⁻¹ (i, j). The arrows indicate the direction of current flow. Current speed difference between the higher river flow simulations and the no flow simulation (i.e. current speed when flows were 0 m³ s⁻¹ subtracted from current speeds when river flows were greater than 0 m³ s⁻¹) for river flows in each river of 10 m³ s⁻¹ (k, l), 50 m³ s⁻¹ (m, n), 100 m³ s⁻¹ (o, p) and 150 m³ s⁻¹ (q, r). The red arrows indicate the direction of current flow for the no flow scenario and the blue arrows indicate the direction of the current flow for the higher flow simulation. In each case, the plots on the left side are results for the flooding tide and on the right for the ebbing tide

Current speed difference plots (Fig. 5k–r) show the change in current speeds between the simulations with river flows and zero river flow simulation. Comparison of current speeds between no flow and flow of 10 m³ s⁻¹ shows that, during the flooding tide, currents into the estuary are typically greater when there is no riverine input. Maximum current values occur in the Run and at the narrow sections leading to the estuary from the rivers. When riverine inputs are increased to 50 m³ s⁻¹, positive values are observed at the narrow passage from the rivers and indicate that water is flowing downstream. Values in the Run remain negative and reach − 0.5 ms⁻¹. For the simulations with 100 and 150 m³ s⁻¹, the same trend is observed but the extent of downstream flow into the estuary is greater as river flow increases with current speed differences reaching 1.4 ms⁻¹. During the ebb tide, the difference in current speed increases as the riverine flow value increases. This occurs throughout most of the estuary with the greatest increases observed at the narrow sections of the estuary, where the riverine water enters the basin, and at the Run. Values increase to over 1.4 ms⁻¹ when river flows are 150 m³s⁻¹.

Flushing Time

To constrain flushing times, we used two approaches, the particle tracking and the tidal prism methods. The results of these two methods are now compared. For the particle tracking method, when the river flow is 0 m³ s⁻¹, the estuary is unable to completely flush during the simulation period, resulting in very long flushing times that cannot be completely constrained. Thus, ignoring these results, the average flushing times of the estuary across a tidal cycle range from 13.5 to 124.7 h at neap tides and from 13.1 to 71.9 h at spring tides between river flows of 10 and 150 m³ s⁻¹ per river (Fig. 6a). Maximum flushing times for both neaps (132 h) and springs (88.5 h) occurred at low tide with river flows of 10 m³ s⁻¹ (Fig. 6b). Minimum flushing times occurred at neaps when flows were 125 m³ s⁻¹ during ebbing tide and for springs at flows of 100 m³ s⁻¹ during high tide (Fig. 6b). When river flow is 10 m³ s⁻¹, the neap flushing times far exceed those at spring tides; however, between 25 and 75 m³ s⁻¹, spring flushing times are longer than those at neaps. At river flows over 100 m³ s⁻¹, the neap tides again see longer flushing times. These two switches in the particle tracking approach, visualised in Fig. 6a, occur at river flows of ~ 14 m³ s⁻¹ and ~ 75 m³ s⁻¹. Of note is the apparent rise in flushing time between 100 and 150 m³ s⁻¹ with average times increasing by 3.2 h between the two river flow scenarios and not, as might be expected, decreasing. As such, shortest flushing times are observed when river flow is 100 m³ s⁻¹ and not at river flows of 150 m³ s⁻¹. Although there is some variation in flushing time at high riverine inputs, after flows reach 75 m³ s⁻¹, the flushing time levels out at around 13 h and does not deviate much from this value.

Fig. 6

Flushing time results. The average particle tracking flushing times across a tidal cycle under different river flows applied at spring and neap tides are presented in a. Tidal prism results using a return flow factor of 0.3 (the average return flow factor value expected for Christchurch Harbour as determined by the particle tracking method) are also presented in a. Particle tracking flushing times at different stages of the tidal cycle for springs and neaps under different river flow scenarios (values included on lines) in b. Flushing time for the tidal prism method at spring and neap tides under different river flow scenarios with a range of return flow factors considered in c. To observe any changes in control on circulation, neap flushing times were subtracted from spring flushing times for each river flow value, the results of which can be seen for both the particle tracking and tidal prism methods in d. As with a, the tidal prism results shown in d are those taken using a return flow factor of 0.3. The dashed line indicates the point at which the circulation control mechanism of the estuary changes from tidal (below the line) to fluvial (above the line)

Study of the whole estuary (Figs. 7 and 8) indicates large spatial variability in flushing time (i.e. residence time) across the system for the PT method. The northernmost section of the estuary consistently exhibits long flushing times (over 5 days), regardless of the river flow or tidal state considered. This coincides with an area of limited circulation (seen in Fig. 5). The centre of the estuary flushes the fastest at typically less than half a day with the southernmost sections and channel leading to the rivers taking between 1 and 2 days to flush at flows greater than 25 m³ s⁻¹. From fluxes of 50 m³ s⁻¹, the flushing time does not change significantly, even under river inputs of 150 m³ s⁻¹. River flows of 10 m³ s⁻¹ produce flushing times in the estuary ranging from less than 1 day near the Run to 5 days in the river channel and the southern areas of the system.

Fig. 7

Local flushing times at spring high tide using the particle tracking method. Estuary flushing times under river flows in each river of 0 m³ s⁻¹ (a), 10 m³ s⁻¹ (b), 25 m³ s⁻¹ (c), 50 m³ s⁻¹ (d), 75 m³ s⁻¹ (e), 100 m³ s⁻¹ (f), 125 m³ s^-1 (g) and 150 m³ s⁻¹ (h). For the no flow results, the estuary was not able to completely flush in the simulation period so the maximum time recorded was 20 days. It can be assumed that flushing times were much greater than this value

Fig. 8

Local flushing times at neap high tide using the particle tracking method. Estuary flushing times under river flows in each river of 0 m³ s⁻¹ (a), 10 m³ s⁻¹ (b), 25 m³ s⁻¹ (c), 50 m³ s⁻¹ (d), 75 m³ s⁻¹ (e), 100 m³ s⁻¹ (f), 125 m³ s^-1 (g) and 150 m³ s⁻¹ (h). For the no flow results, the estuary was not able to completely flush in the simulation period so the maximum time recorded was 20 days. It can be assumed that flushing times were much greater than this value

The point in the tidal cycle at which the particles are released greatly influences the flushing time as expected and in agreement with previous works (Oliveira and Baptista 1997). Comparisons at different river flows suggest no overall trend though during spring tides, the shortest flushing times typically occur at high tide and the longest flushing times at low tide (Fig. 6b). The neap tide results appear more unpredictable with both longest and shortest flushing times changing between tidal state. For river flows greater than 50 m³ s⁻¹, flushing times are greater during the spring tide when the particles are released at low tide; however, when particles are released at high tide, the neap flushing times are greater.

Comparison of the two flushing time approaches shows that the overall patterns are similar with both methods exhibiting exponential-type decays in flushing time as river flow changes (Figs. 6 a and c). For both approaches, when river flow is less than 25 m³ s⁻¹, the spring tide flushing times are shorter than those at neap tides; however, at river flows of 25 m³ s⁻¹ and above, a switch is observed between spring and neap tides, resulting in spring tides producing the longer flushing times. Figure 6 d shows the difference between the spring and neap tide flushing times, with the cross-over point for both methods falling at river flows of ~ 15 m³ s⁻¹ per river. There are some discrepancies between the methods; estimated flushing times between spring and neap tides have a much greater range with the particle tracking method than the tidal prism method at low river flows (Figs. 6 a and c) and the tidal prism method produces shorter flushing times than that of the particle tracking method.

Return flow factors were determined using the particle tracking method (results not presented in this work) and varied from 0.04 to 0.534 between river flows of 10–150 m³ s⁻¹ with greatest values occurring under low river flow conditions. It would be expected that very low river flows (less than the 10 m³ s⁻¹ studied here) would yield return percentages of greater than the maximum value mentioned above.

Discussion

Hydrodynamics

The aim of this study was to examine the physical processes driving circulation and flushing time in the shallow semi-enclosed microtidal estuarine basin of Christchurch Harbour. River flow was found to have a major impact on circulation within the estuary. During flooding tides, the estuary remains tidally dominant in the upper reaches until flows exceed 50 m³ s⁻¹ in each river, typical values expected during winter months. At 50 m³ s⁻¹, there is a build-up of water at the estuary mouth in the narrow Run at Mudeford where the incoming tide is slightly stronger than the out flowing river water, preventing the freshwater from leaving the embayment and thus backing up into the estuary. This phenomenon has been observed when the estuary has flooded in the past during periods of high rainfall and river inputs (ABP 2009). Tidal flow into the estuary appears to cease once river flows exceed 100 m³ s⁻¹, at which point circulation is primarily driven by riverine inputs and the estuary acts as a continuously ebbing system. Whilst flows do not usually exceed 100 m³ s⁻¹, they have been observed on numerous occasions within the last decade (twice in 42 years in the River Avon and 13 times in 44 years in the River Stour); thus, it is highly possible that this scenario (a fresh harbour with no apparent flooding tide into the estuary) has been observed in the past.

Most flow within the estuary occurs through the deeper channel with little flow on the much shallower intertidal banks. The northern part of the estuary sees very little in respect to water movement with weak currents that are unchanged by increased riverine inputs. This is observed in Figs. 7 and 8 as this region of the estuary experiences much greater flushing times than the rest of the basin. The only location where current speeds appear to change drastically is in the Run. However, even when river flows are as large as 150 m³ s⁻¹, the modelled current speeds in the Run do not exceed 1.6 ms⁻¹, contrasting with other observations which reported observed values as high as 2.5 ms⁻¹ (ABP 2009). This could be attributed to a lack of up-to-date accurate bathymetry data both within the Run and the estuary. The depth either side of the Run is difficult to accurately constrain due to the transient sandbanks present and as such averages must be used. It is possible values applied in the model are too deep and not fully representative of the true values, potentially culminating in water not being forced through the narrow opening as quickly as it should be.

Flushing Time

The second objective of this study was to investigate the impact of river flow on flushing times in Christchurch Harbour, using a particle tracking method to identify variability across the microtidal system. The results of this suggested changes in the dominant control of circulation in the estuary. At neap tide, when the tidal range is small, the ability to remove pollutants would be expected to be reduced. However, results indicated that this was not always the case, and on occasion, spring flushing times can be longer. In most estuarine systems, when riverine inputs are small, the basin will likely be controlled by the much greater tidal forcing. During spring tides, the tidal input is greater than at neaps, and thus, shorter flushing times will be observed as a greater tidal input will flush the system more effectively. A system with small tidal inputs (e.g. microtidal) studied during a period of high river flow will be controlled by the riverine inputs. When the river water interacts with the pre-existing estuarine water, the additional volume in the estuary at the spring tide will require more time to be removed, promoting longer flushing times than at neap tides. It can be assumed then that this switch indicates the point at which freshwater inputs begin to drive the circulation of the estuary.

Interestingly, the particle tracking method indicates a second switch between springs and neaps. The difference between spring and neap flushing times when river flows exceed 75 m³ s⁻¹ is minimal (typically between 1 and 3 h), unlike during periods of low river flow. The cause of this second switch is unknown. It is possible that the second switch occurs at the point in which the volume of river water entering the system is so great that the change in volume between spring and neap tide is in comparison negligible to the riverine inputs and the tidal inputs become irrelevant. There may be a point just after the switch at 75 m³ s⁻¹ but before tidal irrelevance occurs where the mechanism driving flushing time is the volume of water within the main channel in the estuary. The particles within the harbour where circulation is limited drive the flushing time up considerably so it is flow through the central channel that will bring down the average flushing time of the system. Greater flow at spring tides will reduce flushing times at the margins of the estuary as flow is less restricted to the central channel; however, lower flow at neap tides will restrict flow to the channel, driving flushing times up at the estuary margins. In addition to this second switch at high river flows, flushing times did not conform with the trend as they increase with increasing river inputs. When flows pass a threshold (in this example, 75 m³ s⁻¹), the large freshwater input into the estuary may not be able to leave the Run fast enough, potentially resulting in water backing up into the basin as found by ABP (2009) under high river flow scenarios. This is evidenced in Fig. 5 g and o which show water velocities through the Run slowing under high river inputs due to the freshwater interacting with the incoming tide and appearing to prevent both in- and egressing of water. This would prevent particles from leaving the system. Greater river flows will result in more water backing up in the estuary which would explain the slight rise in flushing time observed between flows of 100 and 150 m³ s⁻¹.

The main question we are addressing here is the behaviour of small shallow systems without the assumption they are the same as larger systems. In addition, we question the practicability of modelling at such proportions. Due to its particularly complex tidal structure, Christchurch Harbour is an ideal location to test the effectiveness of a hydrodynamic model at such small scales and it can be assumed that other similar systems with less complex tidal patterns would also be modelled effectively in the same manner. We show it is possible to model small estuaries at very high temporal and spatial resolution under reasonable time frames using standard computing power. We have established that Christchurch Harbour is controlled by tidal forcings during the summer but can accurately identify when the riverine inputs will drive circulation within the system. This is consistent with studies of larger systems that have shown similar shifts in control in response to changes in riverine inputs (e.g. the Patos Lagoon in Brazil; Fernandes 2001).

Although the work presented here is case study based, the methodologies and results could be used to provide estimates of flushing times and changes in circulation control for other estuaries with similar characteristics such as shallow bathymetry, high riverine input, small size and microtidal. Christchurch Harbour is not alone in its physical characteristics; there are numerous small, shallow semi-enclosed microtidal systems throughout the world. Much of the Mediterranean is microtidal so an understanding of the changes in circulation control in small estuaries along the coast can be gained from the work presented here. Using our results, we can predict how flushing times in small microtidal systems such as Christchurch Harbour will respond to changes in freshwater inputs and thus use this as a proxy for eutrophication susceptibility. With climate change predictions indicating significant declines in summer river flows (von Christierson et al. 2012; Robins et al. 2016) in temperate systems, we can gain a better understanding of how the small estuaries that make up most of our coastlines will respond by projecting the results from Christchurch Harbour to similar estuaries. This will allow for us to be better prepared should a system’s influencing factors appear to be in line with those conducive to eutrophication development. Examples across the world include the Hartenbos and Klein Brak estuaries in South Africa, both of which are temporary open/closed estuaries that have been studied to examine eutrophic condition indicators (Lemley et al. 2015). These systems could be successfully modelled in a similar way to Christchurch Harbour to aid understanding of flushing times and thus eutrophication susceptibility.

To provide a quick and simple comparison with the particle tracking method, we implemented a tidal prism approach to Christchurch Harbour. Although differences in the absolute values between the two methods (shown in Fig. 6a) were observed, this can be attributed to the tidal prism method yielding faster flushing times than are expected in reality (Choi and Lee 2004). The results did demonstrate similar overall trends with a switch in flushing time observed between spring and neap tides. Despite these similarities, this comparison has highlighted the simplicity of non-modelling flushing time calculation methods and emphasises the benefits of using a complex model to study small, shallow systems and their water retention times. The particle tracking flushing times can be used more accurately than the tidal prism results to predict areas susceptible to declines in water quality, particularly since the former method provides spatiotemporal information as opposed to one value for the whole system. If a complex hydrodynamic model could not be applied to these systems, it would be possible to quickly and easily apply the tidal prism method used in this study to obtain reasonable estimates of flushing times and a better understanding of expected patterns across spring-neap cycles. Though the results would lack detailed spatial and temporal information, an inherent problem of using such a method, our results indicate the tidal prism method can provide a sensible estimate of estuarine flushing time patterns in comparison to more complex particle tracking models. If the differences observed between the particle tracking and tidal prism methods presented in this study were applied to tidal prism results for a new system, a better estimate of flushing time than one attained through use of the tidal prism method alone could be obtained.

As Christchurch Harbour is known to be a mostly well-mixed and very shallow system, we consider that use of the 2D model is justified. Numerous other works (e.g. Lee et al. 2013; Nasermoaddeli et al. 2018; Ye et al. 2018) have shown that it is possible to apply a 3D model to shallow systems and as such future work could consider application of a 3D model to Christchurch Harbour to see if better understanding of circulation patterns could be gained. In addition, future work will consider coupling a biogeochemical model to the hydrodynamic model to identify regions of reduced oxygen concentration. We intend to compare these results with those of the particle tracking model to establish similarities and examine the use of a particle tracking approach as a proxy for susceptibility to declines in water quality. A simplification with this study is the omission of wind within the model. Further examination of Christchurch Harbour would include a more in-depth analysis of the influence of wind on circulation patterns and flushing times.

A better understanding of the influence of estuary size on flushing time would be gained from the use of a finite volume box model with one source and one exit. In future, this approach could be used. In addition, to provide a potentially more representative estimate of the flushing of estuarine volume, an interesting addition to the 2D particle tracking model would be to proportionally distribute particles according to depth.

Conclusions

The overall aim of this paper was to assess which physical drivers influence circulation and water flushing time in a small, shallow estuary. This paper has focused on two objectives using Christchurch Harbour as a case study. The first objective was to develop a 2D hydrodynamic model in MIKE 21 to assess the influence of river flow and changing tidal conditions on circulation in Christchurch Harbour. The model accurately replicates the hydrodynamics of Christchurch Harbour proving that it is possible to model the physics of small, shallow, well mixed estuarine systems, even one with a complicated tidal cycle. The second objective was to investigate the impact of river flow on flushing times in the estuary under a range of tidal conditions. The particle tracking method indicated localised areas of the estuary where long flushing times were observed under all river flow scenarios studied and may therefore become susceptible to eutrophic conditions. These results can thus be used as a proxy for reduced water quality in this area of the estuary potentially over extended periods of time, and as such, this should be monitored closely. On a local scale, the results from this work can be used to investigate the physical controls on the biogeochemistry of the system through the use of the ECO Lab module in MIKE 21. This would allow for the study of the impact of macronutrients on the eutrophication status of Christchurch Harbour, the results from which could be used to justify use of flushing time as a proxy for susceptibility to declines in water quality if low oxygen and high chlorophyll concentrations are observed in the long flushing time regions. This would ultimately provide insight to the effect of future climate change scenarios on water quality in small, shallow estuaries. On a larger scale, we can use our findings to predict variation in flushing times for similar semi-enclosed small, shallow basins and use this as a proxy for water quality and system health.