Numerical Modelling of Oil Spill Transport in Tide-Dominated Estuaries: A Case Study of Humber Estuary, UK

: Oil spills in estuaries are less studied and less understood than their oceanic counterparts. To address this gap, we present a detailed analysis of estuarine oil spill transport. We develop and analyse a range of simulations for the Humber Estuary, using a coupled hydrodynamic and oil spill model. The models were driven by river discharge at the river boundaries and tidal height data at the offshore boundary. Satisfactory model performance was obtained for both model calibration and validation. Some novel ﬁndings were made: (a) there is a statistically signiﬁcant ( p < 0.05) difference in the inﬂuence of hydrodynamic conditions (tidal range, stage and river discharge) on oil slick transport; and (b) because of seasonal variation in river discharge, winter slicks released at high water did not exhibit any upstream displacement over repeated tidal cycles, while summer slicks travelled upstream into the estuary over repeated tidal cycles. The implications of these ﬁndings for operational oil spill response are: (i) the need to take cognisance of time of oil release within a tidal cycle; and (ii) the need to understand how the interaction of river discharge and tidal range inﬂuences oil slick dynamics, as this will aid responders in assessing the likely oil trajectories.


Introduction
Estuaries can be defined as semi-enclosed bodies of water in which incoming saline ocean water is diluted by freshwater [1]. Estuaries hold ecological, recreational and economic value given their free connection to the open sea and they are rich in biological and ecological resources [2,3]. As a result, estuaries are often centres of human activities, which can increase the risk of environmental problems [4,5]. For example, most estuaries because of the services they can provide, increasingly host shipping vessels of different sizes [6]. These vessels include crude and refined oil carriers, with potential risks of oil spills during routine processes such as engine maintenance, transfer of contents, and tank cleaning [7]. There is also the risk of accidents, such as grounding and collision which can be major sources of oil spills [8][9][10]. These raise significant concerns as oil is considered one of the most detrimental sources of anthropogenic pollution in estuaries [11,12].

Estuary Characteristics Fluid Dynamics Principles Criteria for Model Selection
Salt wedge estuary R/V ≥ 1 Highly stratified estuary R/V~0.1-1.0 Include the vertical dimension in at least two-layer model partially mixed (weakly stratified) estuary R/V~0.005-0.1 Can include the vertical dimension in a multi-layered model Well-mixed estuary R/V < 0.005-0.1 Neglect vertical dimension, unless water quality process dictates vertical resolution where R is runoff volume of freshwater entering the estuary during high tide; and V is runoff volume of saltwater entering the estuary during high tide.
Considering the significant risks that oil spills present to estuaries, it is important to better understand the influences of both tidal flows and seasonal freshwater discharges on oil spill transport in estuaries. Hence, our study aims to avail some novel understanding of (a) the implications of intermingling tidal currents and riverine flows on oil spill transport; (b) the influence of river discharge variability on oil spill transport; and (c) how the time of release affects oil slick transport-in estuaries. Because of the complexities of estuaries, numerical models have become essential tools to understand the processes occurring within them [21,36,37]. Therefore, this study employs numerical models to address the study aim.

Study Area
The Humber Estuary, located on the northeast coast of the UK, is one of the largest estuarine systems in the UK with a mean fluvial flow of~250 m 3 /s (and high flow of 1600 m 3 /s), with Q10 and Q95 (i.e., 10th and 95th percentile flows) of 610 and 58 m 3 /s, respectively, over the period 1980-2015 [26,38,39]. These volumes are delivered from several rivers (including the R. Trent and the R. Ouse), with a combined catchment area of 24,240 km 2 (draining approximately 20% of England's landmass) [40]. Several smaller rivers such as the R. Hull, R. Ancholme and R. Freshney, also discharge into the Humber Estuary. The estuary meets the North Sea at Spurn Head, located approximately 62 km downstream from the Trent-Ouse confluence at Trent Falls [41] (Figure 1). The distance of tidal influence extends from the mouth to 120 km on the River Ouse and 147 km on the River Trent [42,43]. The width of the Humber Estuary is approximately 8 km at the mouth (Spurn Head) and is less than 0.5 km upstream of the Ouse-Trent confluence, making it a funnel-shaped estuary [26,35]. Its mean spring tidal range is 5.7 m, making the Humber a well-mixed macro-tidal estuary [26,43]. km 2 (draining approximately 20% of England's landmass) [40]. Several smaller rivers such as the R. Hull, R. Ancholme and R. Freshney, also discharge into the Humber Estuary. The estuary meets the North Sea at Spurn Head, located approximately 62 km downstream from the Trent-Ouse confluence at Trent Falls [41] (Figure 1). The distance of tidal influence extends from the mouth to 120 km on the River Ouse and 147 km on the River Trent [42,43]. The width of the Humber Estuary is approximately 8 km at the mouth (Spurn Head) and is less than 0.5 km upstream of the Ouse-Trent confluence, making it a funnelshaped estuary [35,26]. Its mean spring tidal range is 5.7 m, making the Humber a wellmixed macro-tidal estuary [26,43]. Figure 1. The Humber Estuary system. Inset shows the location of the Humber Estuary in the UK [44,45]. Tidal limits are indicated as red bars crossing the channel.
The Humber Estuary has one of the UK's biggest constellations of oil refining industries [44,46]. In the outer Humber Estuary is the Tetney Monobuoy from which crude oil is unloaded and transferred through underground pipelines to the ConocoPhillips Ltd.'s Humber refinery [47]. Furthermore, three marine terminal facilities, located at the Southbank of the Humber (South Killingholme Jetty, Immingham gas jetty and Immingham oil terminal) provide services to two oil refineries at Immingham [46]. Consequently, the Humber Estuary plays host to oil tankers which berth at the oil terminals [47]. In 2001, over 40 million tonnes of oil and chemicals were transported in and out of the Humber [48] which makes it the main east coast port for hydrocarbon landing in the UK [47]. The estuary thus is prone to oil spill risk, and several incidents have occurred in the past. The Sivand tanker incident, which occurred in September 1983, spilled 6000 tonnes of Nigerian Figure 1. The Humber Estuary system. Inset shows the location of the Humber Estuary in the UK [44,45]. Tidal limits are indicated as red bars crossing the channel.
The Humber Estuary has one of the UK's biggest constellations of oil refining industries [44,46]. In the outer Humber Estuary is the Tetney Monobuoy from which crude oil is unloaded and transferred through underground pipelines to the ConocoPhillips Ltd.'s Humber refinery [47]. Furthermore, three marine terminal facilities, located at the Southbank of the Humber (South Killingholme Jetty, Immingham gas jetty and Immingham oil terminal) provide services to two oil refineries at Immingham [46]. Consequently, the Humber Estuary plays host to oil tankers which berth at the oil terminals [47]. In 2001, over 40 million tonnes of oil and chemicals were transported in and out of the Humber [48] which makes it the main east coast port for hydrocarbon landing in the UK [47]. The estuary thus is prone to oil spill risk, and several incidents have occurred in the past. The Sivand tanker incident, which occurred in September 1983, spilled 6000 tonnes of Nigerian light crude into the Humber Estuary [49]. Moreover, 46 tonnes of crude oil spills from nine different incidents have been recorded on the Humber Estuary between 1989 and 1997 [47]. These oil spill records and the presence of oil refineries and crude oil transport along the Humber Estuary further highlight the significance of this study.

Hydrodynamic Model
The hydrodynamics of the Humber Estuary were computed with TELEMAC3D, a three-dimensional open-source finite element model [50]. TELEMAC3D solves the Navier-Stokes equations with or without hydrostatic pressure approximation [51,52]. The main originality of TELEMAC, which lies in the flexibility and efficiency of its finite element algorithms, makes it stand out in comparison to other comprehensive modelling systems [50,53,54]. In addition, it provides ready integration with TELEMAC's built-in oil spill model [55,56] (Section 3.2). This study employs the hydrostatic pressure approximation due to the intense computational power demanded by non-hydrostatic pressure models. Furthermore, the impact of using hydrostatic pressure approximation on the reliability of the model is expected to be relatively insignificant, as there is no indication that the processes (i.e., small scale ocean processes) [57,58] emphasised using non-hydrostatic pressure models influence estuarine oil slick transportation. TELEMAC3D uses a sigma transformation to resolve the vertical direction and unstructured triangular grid in the horizontal direction [50,59]. It provides a wide variety of numerical options to solve various processes within the model [52] including numerous turbulence models [56]. Here, horizontal turbulence was resolved using the Smagorinsky model, which is best suited for tidal systems that involve highly non-linear flow [60]. Vertical turbulence was resolved using the Nezu and Nakagawa mixing length model, which offers a good representation of wind drift [60]. While TELEMAC3D offers four vertical turbulence models (Prandtl model; Nezu and Nakagawa model; Quetin model; and Tsanis model), Rahman and Venugopal [60] demonstrate that the differences between the models are almost negligible when used with the Smagorinsky model as the horizontal turbulence model.
The simulated model domain extends from the mouth of the Humber Estuary to upstream Blacktoft on River Ouse and Keadby on the River Trent ( Figure 2). The Association of British Ports (ABP) provided the bathymetry data collected in 2008 with reference to chart datum. All positions are aligned to the UK Ordnance Survey National Grid (OSGB) 1936/British national grid and depth with reference to chart datum. Each dataset is referenced using a local chart datum; a local reference level. This datum varies over the area (and sometimes over time) ( Figure S1). To obtain a consistent reference level, an interpolated grid of values is employed to convert all data to ODN (Ordnance Datum Newlyn), which is the UK's national reference level.
light crude into the Humber Estuary [49]. Moreover, 46 tonnes of crude oil spills from nine different incidents have been recorded on the Humber Estuary between 1989 and 1997 [47]. These oil spill records and the presence of oil refineries and crude oil transport along the Humber Estuary further highlight the significance of this study.

Hydrodynamic Model
The hydrodynamics of the Humber Estuary were computed with TELEMAC3D, a three-dimensional open-source finite element model [50]. TELEMAC3D solves the Navier-Stokes equations with or without hydrostatic pressure approximation [51,52]. The main originality of TELEMAC, which lies in the flexibility and efficiency of its finite element algorithms, makes it stand out in comparison to other comprehensive modelling systems [50,53,54]. In addition, it provides ready integration with TELEMAC's built-in oil spill model [55,56] (Section 3.2). This study employs the hydrostatic pressure approximation due to the intense computational power demanded by non-hydrostatic pressure models. Furthermore, the impact of using hydrostatic pressure approximation on the reliability of the model is expected to be relatively insignificant, as there is no indication that the processes (i.e., small scale ocean processes) [57,58] emphasised using non-hydrostatic pressure models influence estuarine oil slick transportation. TELEMAC3D uses a sigma transformation to resolve the vertical direction and unstructured triangular grid in the horizontal direction [50,59]. It provides a wide variety of numerical options to solve various processes within the model [52] including numerous turbulence models [56]. Here, horizontal turbulence was resolved using the Smagorinsky model, which is best suited for tidal systems that involve highly non-linear flow [60]. Vertical turbulence was resolved using the Nezu and Nakagawa mixing length model, which offers a good representation of wind drift [60]. While TELEMAC3D offers four vertical turbulence models (Prandtl model; Nezu and Nakagawa model; Quetin model; and Tsanis model), Rahman and Venugopal [60] demonstrate that the differences between the models are almost negligible when used with the Smagorinsky model as the horizontal turbulence model.
The simulated model domain extends from the mouth of the Humber Estuary to upstream Blacktoft on River Ouse and Keadby on the River Trent ( Figure 2). The Association of British Ports (ABP) provided the bathymetry data collected in 2008 with reference to chart datum. All positions are aligned to the UK Ordnance Survey National Grid (OSGB) 1936/British national grid and depth with reference to chart datum. Each dataset is referenced using a local chart datum; a local reference level. This datum varies over the area (and sometimes over time) ( Figure S1). To obtain a consistent reference level, an interpolated grid of values is employed to convert all data to ODN (Ordnance Datum Newlyn), which is the UK's national reference level. The Blue Kenue software tool [61] was employed to design the unstructured grid and set boundary conditions. The model was implemented with five equidistant σ-coordinate layers in the vertical. The application of five σ-coordinate levels is sufficient and efficient for developing an operational oil spill system [62]. In addition, the Humber Estuary is a well-mixed estuary [26], hence, it does not necessitate multi-layers to resolve the vertical direction [35] ( Table 1). The resulting computational domain consists of 92,369 nodes and 183,925 elements, with a mean edge length of 54.57 m varying from 11 m to 803 m ( Figure 2). The choice of mesh size was also influenced by computational power demand as the computational system employed in this study could not cope with the computational demand of edge length less than 50 m. Brown and Davies [63] indicate that a highresolution mesh ranges between 20-100 m as this allows high resolution of the significant bathymetric features. Furthermore, an edge length of 50 m is enough to adequately capture flow propagation [60]. Consequently, the mesh size is not expected to have a significant impact on the results. Chezy's law of bottom friction was employed as it is most suited for TELEMAC3D models applying equidistant vertical layering [60]. In an analysis of over 15 estuaries, Savenije [64,65] indicate that the Chezy friction coefficient typically ranges from 45-70 m 0.5 /s. In our study, the Chezy friction coefficient is a calibration parameter as it is an important and sensitive TELEMAC3D model parameter (Section 3.3.2). A simulation time step of 45 s was employed which satisfies the Courant-Friedrichs-Lewy criterion.

Oil Trajectory Model
The movement of an oil slick on water is mainly governed by advection, i.e., transport due to the water movement [66][67][68][69]. This process is dominated by current, wind and waves, which can act independently (i.e., in the same or opposite direction) and concurrently on oil slick transport [69,70]. Horizontal and vertical displacement (dispersion) of an oil slick is also influenced by turbulent diffusion [71][72][73]. The current generation of oil spill modelling tools employed to support operational oil spill response and impact assessment include GNOME [74]; ADIOS2 [75]; OSCAR [76]; OILMAP [77]; MEDSLIK [78][79][80], amongst others. This generation of oil spill modelling tools mostly uses Lagrangian formulation [81] to compute oil transport (advection and dispersion) and utilises individual formulations to compute crude oil weathering processes [82]. The Lagrangian approach involves representing oil slicks by several constituents (particles) that are transported by advection and dispersion [83]. In contrast, some oil spill models are based on Eulerian approaches [84,85], which either apply the mass and momentum equation to the oil slick layer or apply a convection-diffusion equation [86,87]. In general, Lagrangian models are more popular in oil spill modelling because complexities that arise from irregular and discontinuous oil slick shapes are more easily simulated than with Eulerian models [88]. Furthermore, dispersion phenomena are not directly resolved by Eulerian models, thereby resulting in numerical diffusion [89]. Compared to Eulerian models, the Lagrangian approach is more efficient and cost-effective [90]. Here, we employed the TELEMAC oil spill model because of its integration within the TELEMAC3D hydrodynamic model. Also, the TELEMAC oil spill model combines the Lagrangian and Eulerian approaches [55], which make its choice more suitable in this study.
Integrated into TELEMAC3D is a two-dimensional oil spill model for short-term forecasting of oil spill behaviour in continental waters (lakes, rivers and estuaries) [55,91]. The oil spill model accounts for advection, diffusion, spreading, evaporation and dissolution processes [91]. A set of hydrocarbon particles, considered as a mixture of discrete non-interacting hydrocarbon components (soluble and insoluble components), represent the oil slick [92] while the TELEMAC3D hydrodynamic model (Section 3.1) provides the required hydrodynamic data. TELEMAC's capability for simulating oil spill transport was validated by Goeury [91] using the ERIKA oil spill incident that occurred in the Bay of Biscay in December 1999. A full description of the validation of TELEMAC oil spill processes is detailed in Goeury [91] and Goeury et al. [55]. While TELEMAC accounts for spreading, evaporation and dissolution of the oil spill, it does not model other key short term weathering processes such as dispersion, emulsification, viscosity, density and buoyancy. The TELEMAC oil spill model can, therefore, be considered rudimentary with regards to modelling oil weathering processes. Consequently, the focus of this study is on oil spill transport in a tide-dominated estuarine environment. Advection of an oil slick can be expressed as: where U oil represents the oil slick velocity vector, U w represents wind velocity vector above the water surface, represents the drift of the surface slick due to the wind (β = 0.036; [92]), and U c represents the current velocity vector at the free surface. A stochastic approach is employed to account for turbulent diffusion [92]. The turbulent diffusion equation can be expressed as: where C represents a depth-averaged pollutant concentration, U represents the depthaveraged velocity vector replaced by U oil , h represents water depth computed from TELEMAC3D, and v t represents turbulent viscosity computed from TELEMAC3D.

Calibration and Validation Design
We examined the temporal variability of river discharge in the Humber Estuary to characterise the estuary's seasonality. Hydrodynamic model scenarios were developed for summer and winter conditions in the estuary. In selecting a typical summer and winter month, it is important to pick a representative month for comparison with reallife data. Consequently, stage data available at the river boundary points were analysed to determine the best representative summer and winter months for the hydrodynamic models. Because of limited high-resolution discharge data, the models were driven by constant river discharge at the river boundaries and 15 min tidal height data at the offshore boundary. In the absence of current data, the models were calibrated and validated with water levels during the representative summer and winter months. Consequently, we are unable to assess the reliability of the simulated current. This suggests that there might be some discrepancies between simulated and actual oil spill transport, as current is the main driver of oil spill advection. However, in large estuaries, reasonable reproduction of the depth-averaged currents can be achieved if accurate bathymetry is employed and tidal elevation is accurately simulated [93].

River Discharge
This study examines the temporal variability of discharge from the Rivers Ouse and Trent over an 11-year period (2007-2017). The UK Environment Agency provided highresolution discharge data recorded at 15 min intervals from 2007 to 2017 on Rivers Ouse at Skelton which lies 17 km from Naburn [94] and Trent at North Muskham which lies 0.75 km from Cromwell Weir [95] (Figure 1). At Skelton, the mean monthly discharge from 2007 to 2017 was 56.78 m 3 /s, while at North Muskham it was 89.20 m 3 /s. Analyses of discharge data at Skelton and North Muskham reveal that both Rivers Ouse and Trent exhibit strong seasonal discharge variations, with the highest discharge in the winter and lowest in the summer (Table S1, Figure S2b).
Stage data obtained from the UK Environment Agency from 2007-2017, at Blacktoft on River Ouse and at Keadby on River Trent, indicate that the influence of seasonality reaches the areas under tidal influence ( Figure S3a). We simulated full 28-day semi-diurnal tidal cycles combined with two river flow conditions, representative of summer and winter. The monthly stage averages for each season at Blacktoft and Keadby were compared (using standard deviation) to find the month that best represents the seasonal flow condition ( Figure S3b). For each season, the average monthly stage and stage variability were calculated. Subsequently, we selected from the observed data the month which most resembles the seasonal average using standard deviation. This procedure results in representative seasonal data (Table 2), without losing the fidelity and variability of the observed empirical data, whilst maintaining compatibility and comparability with observed tidal data for the selected periods. Table 2. Simulation period employed to develop hydrodynamic models for summer (characterised by low river discharge) and winter (characterised by high river discharge).

Season
Representative Because of limited high-resolution discharge data at Blacktoft and Keadby (Figure 2), the open river boundary was driven by constant discharge obtained from the literature ( [45,96]; Table 3). Tidal gauge measurements were extracted from the interaction of 34 tidal components in the Finite Element Solution model (FES 2014; [97]). The open offshore boundary was driven by 15 min tidal height data extracted from the FES2014 model (at National Grid References 541,049, 402,468 and 540,841, 409,420) with respect to Mean Sea Level (see Eke et al. [98] for details). Constant temperature (17.4 • C and 6.6 • C for summer and winter respectively) and salinity (29.3 ppt and 27.2 ppt for summer and winter, respectively) were assumed along the open offshore boundary.

Model Calibration and Validation Scenarios
The hydrodynamic model was calibrated and validated against 15 min tidal height data for Immingham station located at National Grid Reference 520,049, 416,473 (Figure 1), provided by the British Oceanographic Data Centre (BODC). In large estuaries (e.g., the Humber Estuary), reasonable reproduction of the depth-averaged currents can be achieved if accurate bathymetry is employed and tidal elevation is accurately simulated [93,100]. To improve the model's ability to reproduce accurate depth-averaged currents, the model performance was compared using three metrics: Root mean square error (RMSE), which is appropriate to determine the agreement between model output and measured data for scalar quantities (e.g., water levels) [101,102]; coefficient of determination, R2; and regression coefficient, b. Bottom friction (i.e., Chezy's friction coefficient) is an important and sensitive TELEMAC3D model parameter [54,60] and formed the focus of the model calibration. Furthermore, Skinner et al.'s [43] hydrodynamic model of the Humber Estuary revealed that the estuary is sensitive to bottom friction. Sensitivity analysis can be undertaken to understand other TELEMAC3D model parameters that are important and sensitive within the Humber Estuary, however, this is beyond the scope of this study. The calibration was done independently for each representative month in each season (Table 4) to allow for seasonal variation in optimal parameter values. Model validation was undertaken by comparing model results with observed data from Immingham station for the second most representative month for each season (Table 1). To investigate the relative impacts of river discharge variation and tidal flow on oil slick dynamics for the Humber Estuary, we simulate several scenarios. Each scenario simulates an instantaneous oil release from one of two arbitrary locations, an upstream and a downstream one (L1 and L2; Figure 2). The oil is released from these locations under eight scenarios, representing spring and neap tidal conditions at Immingham, during both high and low tide, under both summer (low discharge) and winter (high discharge) flows (Table S2). Immingham (Figure 1) was chosen as the reference point because of data availability. As a result, 16 oil spill scenarios are simulated in total. Each scenario simulates an instantaneous release of 10,000 m 3 (arbitrarily chosen) of Brent Crude (883.6 kg/m 3 ), and the resulting simulated oil slick is monitored for a 48 h period (over two semi-diurnal tidal cycles). Due to the lack of a real-case oil spill, this study was unable to validate the oil spill results. Consequently, there might be some discrepancies between the simulated and actual oil slick transport but the study, nevertheless, highlights the importance of understanding the dynamics of oil transport in estuaries as a step towards improving operational response.
To visualise the oil spill transport, a python script was utilised to convert the TELEMAC oil spill displacement output file (tecplot ® format) into an ArcMap readable (xyz) format. The measure tool within ArcMap was employed to measure (a) the oil slick impacted area on the water surface, (b) length of the oil slick over time (distance from one end of the slick to the other), and (c) overall distance travelled (maximum upstream and downstream displacement from the point of release). We are interested in studying the dynamics of the oil slick in the estuarine environment and these parameters are useful to understand the dynamics, which have not been previously studied (Section 1). Furthermore, we employ a one-way ANOVA followed by the Holm-Sidak method (p < 0.05) to understand the relative influence of hydrodynamic conditions river discharge variation, water level and tidal range) on these parameters. We defined oil beaching as occurring when oil interacts with the edge of the computational grid. The length of the polluted estuary bank line was measured within ArcMap.

Wind Data for Oil Spill Scenarios
The average wind speed and direction for the representative months were computed from hourly wind speed and direction data at Donna Nook station near the Humber's mouth (National Grid Reference 542,900; 399,700) obtained from the UK Meteorological Office. Wind speed and direction were assumed constant in time and space in the oil spill simulations (Table S3). However, this assumption is not expected to significantly affect overall simulation results, as the influence of wind speed on oil slick transport is approximately 3.6% of the wind speed [92]. The calibrated results revealed that the best agreement between the measured data and model results was obtained with Chezy C of 70 m 0.5 /s and 75 m 0.5 /s in summer and winter, respectively ( Table 4). The model was able to replicate the observed tidal elevation cycles (Figure 3). This occurred in all seasons. Minor discrepancies in simulated high and low tides might have been influenced by the exclusion of surges in the hydrodynamic model [103], the use of simulated tidal height data rather than measured data to drive the offshore boundary, and the use of constant discharge-driven boundary condition. To improve model performance, the computational domain can be extended upstream to reach a point where high-resolution water discharge measurement is available-a subject of future research. mouth (National Grid Reference 542,900; 399,700) obtained from the UK Meteorological Office. Wind speed and direction were assumed constant in time and space in the oil spill simulations (Table S3). However, this assumption is not expected to significantly affect overall simulation results, as the influence of wind speed on oil slick transport is approximately 3.6% of the wind speed [92].

Calibration
The calibrated results revealed that the best agreement between the measured data and model results was obtained with Chezy C of 70 m 0.5 /s and 75 m 0.5 /s in summer and winter, respectively ( Table 4). The model was able to replicate the observed tidal elevation cycles (Figure 3). This occurred in all seasons. Minor discrepancies in simulated high and low tides might have been influenced by the exclusion of surges in the hydrodynamic model [103], the use of simulated tidal height data rather than measured data to drive the offshore boundary, and the use of constant discharge-driven boundary condition. To improve model performance, the computational domain can be extended upstream to reach a point where high-resolution water discharge measurement is available-a subject of future research.

Validation
The validation was undertaken using different time periods for summer and winter. The calibrated Chezy friction coefficients of 70 m 0.5 /s and 75 m 0.5 /s were employed in secondary summer and winter simulations, respectively ( Table 2). Evans [104] classified the RMSE value that is 15% of the spring tidal range as good, while Brown et al. [105] and Skinner et al. [43] agree that an RMSE value that is less than 20% of the tidal range is considered good. The simulated water level showed good agreement with the recorded data (Table 4; Figure S4) resulting in an RMSE of 9.09% and 12.86% of the spring tidal range in summer and winter, respectively. Furthermore, Guo et al. [106] and Maréchal [107] classify a model performance with R 2 > 0.85 as excellent

Validation
The validation was undertaken using different time periods for summer and winter. The calibrated Chezy friction coefficients of 70 m 0.5 /s and 75 m 0.5 /s were employed in secondary summer and winter simulations, respectively ( Table 2). Evans [104] classified the RMSE value that is 15% of the spring tidal range as good, while Brown et al. [105] and Skinner et al. [43] agree that an RMSE value that is less than 20% of the tidal range is considered good. The simulated water level showed good agreement with the recorded data (Table 4; Figure S4) resulting in an RMSE of 9.09% and 12.86% of the spring tidal range in summer and winter, respectively. Furthermore, Guo et al. [106] and Maréchal [107] classify a model performance with R 2 > 0.85 as excellent

Oil Spill Scenarios
The calibrated and validated Chezy friction coefficients are used for all our oil spill scenarios. It is the first time that these resulting dynamics have been illustrated and confirmed ( Figure 4). As these patterns of oil behaviour in estuaries have not been empirically and formally demonstrated until now, this is a novel contribution to the field. Furthermore, the capability of the simulations to confirm some of these dynamics lends additional credibility to the more surprising findings, e.g., because of seasonal variation in river discharge, winter slicks released at high water did not exhibit any upstream displacement over repeated tidal cycles, while summer slicks travelled upstream into the estuary over repeated tidal cycles (Section 5.1). Oil slicks are seen moving back and forth in the estuary as tides flood and ebb ( Figures S5 and 4), which is not surprising. However, the exceptions are oil spills released at L1 at high water, which leave the computational domain during the first tidal cycle (Figure 4(1a,1c,2a,2c)). Overall, spills travel downstream over time. Nevertheless, spills released at low water first travel upstream as the tides come in (Figure 4(columns b and d)) and the oil slicks spread over time. After 48 h, i.e., after two semi-diurnal tidal cycles, the simulated oil slicks cover an area between 13.7 and 67.8 km 2 on the water surface, with an average of 37.8 km 2 (Figures 4 and S6; Table 5). Within these overall results, there are systematic patterns and trends that relate to the scenario configurations. These nuances are discussed in detail in the next section.

Oil Spill Scenarios
The calibrated and validated Chezy friction coefficients are used for all our oil spill scenarios. It is the first time that these resulting dynamics have been illustrated and confirmed ( Figure 4). As these patterns of oil behaviour in estuaries have not been empirically and formally demonstrated until now, this is a novel contribution to the field. Furthermore, the capability of the simulations to confirm some of these dynamics lends additional credibility to the more surprising findings, e.g., because of seasonal variation in river discharge, winter slicks released at high water did not exhibit any upstream displacement over repeated tidal cycles, while summer slicks travelled upstream into the estuary over repeated tidal cycles (Section 5.1). Oil slicks are seen moving back and forth in the estuary as tides flood and ebb ( Figures S5 and 4), which is not surprising. However, the exceptions are oil spills released at L1 at high water, which leave the computational domain during the first tidal cycle (Figure 4(1a,1c,2a,2c)). Overall, spills travel downstream over time. Nevertheless, spills released at low water first travel upstream as the tides come in ( Figure  4(columns b and d)) and the oil slicks spread over time. After 48 h, i.e., after two semidiurnal tidal cycles, the simulated oil slicks cover an area between 13.7 and 67.8 km 2 on the water surface, with an average of 37.8 km 2 (Figures 4 and S6; Table 5). Within these overall results, there are systematic patterns and trends that relate to the scenario configurations. These nuances are discussed in detail in the next section.  Figure 4c is the same as shown in detail in Figure S5.  Figure 4c is the same as shown in detail in Figure S5.  The results show that an oil slick is likely to remain in the estuary within the first 48 h (Figure 4). This is dependent on the point of release, as an oil slick released at L1 is observed to leave the computational domain (Figure 4(1a,2a,1d,2d)). An oil slick released at L1 could re-enter the estuary on the flood tide, depending on the North Sea current magnitudes and directions. However, this scenario is beyond the scope of this paper, because of the computational domain under consideration.
Using the oil slick released from L2 as the reference case, oil slick lengths were observed to increase over time, albeit not at a uniform rate ( Figure S7). It was observed that the maximum oil slick lengths in summer were longer than the winter slicks for the two HW scenarios ( Figure S7), while maximum oil slick lengths in winter were longer than the summer slicks for the two LW scenarios ( Figure S7). This indicates that the influence of seasonal discharge on oil slick spreading is dependent on the time of release within a tidal cycle.
The oil slick impacted area depends on the release point, time of oil release, neap/spring tide and high/low water scenario and the season (Table 5; Figure S6). After 48 h, the oil slick impacted area is predominantly greater in the summer (by an average of 4%; Table S4). This contrasts with the findings that winter slicks were displaced (i.e., travelled) farther (by an average of 12% than summer slicks) after 48 h (Table S4). Winter slicks thus are longer but narrower than summer slicks. This shows that river discharge variation significantly influences oil slick spreading in the tide-dominated Humber Estuary.
Overall, the distances travelled by winter slicks released at low water are 20.4% and 36.9% greater than those of summer slicks (Figures S9 and S10; Table S4). The differences in distance travelled under high water are less pronounced, as the distance covered by winter slicks is only 1.6% greater than summer slicks when released at high water (HW) neap tide (NT) while the distance travelled by summer slicks is only 0.7% greater than winter slicks when released at HW spring tide (ST) (Figures S9 and S10; Table S4). However, these distances obscure an important qualitative difference in the oil slick dynamics. Oil slicks released under winter conditions (high river discharge) were observed to travel further downstream compared to oil slicks released under summer conditions (low river discharge). After 48 h, the downstream displacement of winter slicks was, on average, 4.19 km further than summer slicks (Figures S9 and S10 ; Table S5). However, summer slicks reached an average of 2.51 km further upstream (Figures S9 and S10; Table S5). While winter slicks released at high water do not experience any upstream displacement, summer slicks were able to travel upstream over repeated tidal cycles (Figures S9 and S10; Table S5). This indicates that river discharge variability has a key influence the upstream and downstream displacement of the oil slicks.
Thus, there is a distinct difference in oil slick dynamics within the estuary, whereby summer spills (released during low river discharge) tend to travel further upstream and winter spills (released during high river discharge) tend to migrate downstream ( Figures S8-S10). As a consequence, oil slicks released under high river discharge are likely to leave the estuary more quickly than oil slicks released under relatively lower river discharge. This dynamic can be attributed to the relatively higher ebb velocities as a result of higher freshwater flow [45]. It also supports Townend and Whitehead [38], who point out that in the Humber Estuary, flood asymmetry dominates in summer, while flow in winter months becomes ebb dominant as gravity flow dominates.

Impact of Tide (Spring Tide vs. Neap Tide)
After 48 h, oil slick-impacted areas on the water surface are between 86% and 144% larger under spring tide conditions than under neap tides, with an average of 125% (Table S6). This is due to the larger magnitudes of the currents driving the oil slick ( Figure S11). Similarly, after 48 h, the lengths of oil slicks released under spring tide were on average 100% longer than under neap tides ( Figure S7; Table S6), whilst the distance travelled by spring tide oil slicks was on average 78% greater than for neap tide oil slicks (Figures S9 and S10; Table S6). After 48 h, the upstream displacement of oil slicks released under spring tide was on average 3.53 km further than under neap tide (Table S5). Similarly, after 48 h, the downstream displacement of oil slicks released under spring tide was on average 6.08 km further than under neap tide (Table S5). Mendes et al. [108] made similar findings in the Ria De Aveiro Lagoon, pointing out that compared to neap tides, spring tides have a larger influence on oil travel distance.

Impact of Tidal Stage (High Water vs. Low Water)
After 48 h, oil slick-impacted areas are between 23% and 74% larger when released at high water than at low water (Table S7). Similarly, after 48 h, the lengths of oil slicks released at high water were on average 52% longer than at low water ( Figure S7; Table S7). In agreement with Yan et al. [30], this suggests that to efficiently deal with oil spills, responders will have to be made aware of the tidal stage at the time of oil release. Simulated oil slicks released at L2 high water were observed to travel further downstream compared to oil slicks released at low water, while oil slicks released at low water were observed to travel further upstream (Figures 4, S9 and S10). After 48 h, the upstream migration of oil slicks released at low water was on average 7.63 km further than at high water. However, downstream displacement of oil slicks released at high water had migrated on average 13.49 km further than at low water (Table S5). This is likely due to the initial direction of displacement at the time of oil spill release caused by flood/ebb tides. The results suggest that oil slick upstream displacement at low water, and downstream displacement at high water, is enhanced during spring tides. The overall distances travelled by oil slicks released at high water were on average 44% greater than at low water (Figures S9 and S10 ;  Table S7). However, oil spills released at low water will have a longer slick residence within the estuary since the oil slick will first move upstream with the incoming tide.

Relative Impacts of River Discharge vs. Tide vs. Stage
Until now, no empirical comparison has been undertaken on the relative influence of hydrodynamic conditions (river discharge variation, water level and tidal range) on oil slick impacted area, length and distance travelled. However, our results indicate that the influence of time of release in a tidal cycle on oil slick dynamics is greater than river discharge. The above analyses (Sections 5. 1-5.3) indicate that these differences are statistically significant (p < 0.05 in pairwise ANOVA using the Holm-Sidak method). This thus provides the first clear evidence that oil slick impacted area, length and distance travelled are predominantly affected by the tidal range (i.e., spring tide or neap tide) at the time of oil release, and that tidal stage and river discharge are only secondary influencing factors (Tables S5, S7 and S8; Figure 5).
influence of time of release in a tidal cycle on oil slick dynamics is greater than river discharge. The above analyses (Sections 5.1, 5.2 and 5.3) indicate that these differences are statistically significant (p < 0.05 in pairwise ANOVA using the Holm-Sidak method). This thus provides the first clear evidence that oil slick impacted area, length and distance travelled are predominantly affected by the tidal range (i.e., spring tide or neap tide) at the time of oil release, and that tidal stage and river discharge are only secondary influencing factors (Tables S5, S7 and S8; Figure 5).

Influence of Oil Spill Release Location
Our results show that the proximity of the oil spill release point to the estuary mouth influences the oil slick impacted area on the water surface. Oil slicks released near the estuary mouth (L1; Figure 4) covered a greater area than oil slicks released further upstream (L2) (Figure 4; Table S8). This can be explained by the stronger currents observed at the estuary mouth ( Figure S11). This also suggests that the oil slick area decreases as the point of release moves further upstream, due to decreasing current velocities. Since the amount of oil released is the same in all scenarios, it also suggests that the oil slicks remain denser (more concentrated) as the point of release is further upstream.

Oil Beaching
The oil slick model shows the risk of an environmental disaster as oil may spread over the shallow intertidal/saltmarshes on the Southbank (as observed in the oil slick released from L1)- Figure 4, or on the densely populated north bank of the Humber Estuary (as observed in oil slick released from L2) ( Figure 4). The length of bank line affected and time of first contact (Table 4) suggests that although the oil slick reaches the estuary bank quicker in winter compared to summer, oil slick beaching covers a greater distance in summer. This is due to the differences in seasonal discharge. Oil slick beaching occurs in summer when released at L2 under spring tide conditions, affecting 13.21 km and 19.45 km of the estuary bank line under low and high water, respectively ( Table 4). The absence

Influence of Oil Spill Release Location
Our results show that the proximity of the oil spill release point to the estuary mouth influences the oil slick impacted area on the water surface. Oil slicks released near the estuary mouth (L1; Figure 4) covered a greater area than oil slicks released further upstream (L2) (Figure 4; Table S8). This can be explained by the stronger currents observed at the estuary mouth ( Figure S11). This also suggests that the oil slick area decreases as the point of release moves further upstream, due to decreasing current velocities. Since the amount of oil released is the same in all scenarios, it also suggests that the oil slicks remain denser (more concentrated) as the point of release is further upstream.

Oil Beaching
The oil slick model shows the risk of an environmental disaster as oil may spread over the shallow intertidal/saltmarshes on the Southbank (as observed in the oil slick released from L1)- Figure 4, or on the densely populated north bank of the Humber Estuary (as observed in oil slick released from L2) (Figure 4). The length of bank line affected and time of first contact (Table 4) suggests that although the oil slick reaches the estuary bank quicker in winter compared to summer, oil slick beaching covers a greater distance in summer. This is due to the differences in seasonal discharge. Oil slick beaching occurs in summer when released at L2 under spring tide conditions, affecting 13.21 km and 19.45 km of the estuary bank line under low and high water, respectively ( Table 4). The absence of beaching in winter under the same scenario is likely due to the narrower slick formed in winter (Section 5.1).
However, the estuary bank affected by the oil slick is likely to depend on the crosssectional location of the release point, which for this study was arbitrarily chosen. Further studies are needed to fully assess the impact of cross-sectional release point location.

Conclusions
We present the first detailed analysis of oil spill dynamics in a large well-mixed macrotidal estuary with some significant contributions to the field. Using the Humber Estuary as a case study, we assessed the influence of tidal cycles and freshwater discharge variability on the dynamics of hypothetical oil spills. This study examined the temporal variability of discharge from Rivers Ouse and Trent over an 11-year period (2007-2017).
TELEMAC3D was employed to develop hydrodynamic models that represent the behaviour of the Humber Estuary in summer (low river discharge) and winter (high river discharge). Satisfactory model performance (R 2 = 0.883 and 0.852 for summer and winter respectively) was attained using constant discharge at the river boundary. The models were validated against measured tidal height data at Immingham station for summer and winter, also with a satisfactory agreement (R 2 = 0.912 and 0.848 for summer and winter respectively). The two-dimensional oil spill model for the Humber Estuary was developed with the TELEMAC oil spill module.
Our results show new findings that, for this large well-mixed macro-tidal estuary: • because of variation in river discharge, slicks released under high river discharge at high water did not exhibit any upstream displacement over repeated tidal cycles, while slicks released under low river discharge travelled upstream into the estuary over repeated tidal cycles; • there is a statistically significant (p < 0.05) difference in the influence of hydrodynamic conditions (river discharge variation, water level and tidal range) on oil slick impacted area, length and distance travelled; • the tidal range has a key influence on oil slick impacted area, with spring tide slicks being 125% bigger than neap tide slicks, on average; • oil slick impacted area, length and distance travelled is predominantly affected by the tidal range (i.e., spring or neap) at the time of oil release, and only then by the stage or river discharge; • the influence of river discharge on oil slick spreading is dependent on the time of release within a tidal cycle; and • the possibility of oil beaching on the banks of the estuary exposes environmental risks, with up to 24.6 km of shoreline affected in our simulations.
Some of these findings conform to intuitively expected dynamics, at least qualitatively. Others are arguably more surprising insights. However, all findings represent novel contributions in the sense that none had been formally and explicitly demonstrated or quantified, until now.
It should be noted that this study investigated instantaneous oil slick releases in a large well-mixed macro-tidal estuary. To reliably generalize the findings, further studies should be undertaken, using continuous oil slick release and extending the analyses to other estuarine systems.
The implications of these findings for operational oil spill response are: (a) the need to take cognisance of time of oil release within a tidal cycle; and (b) the need to understand how the interaction of river discharge and tidal range influences oil slick dynamics, as this will aid responders in assessing the likely oil trajectories and cost-effective deployment of response tools.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/jmse9091034/s1, Figure S1: Outline of the Humber Estuary showing metres Above Ordnance Datum (m AOD) at various locations, Table S1: Seasonal mean discharges of freshwater (m 3 s −1 ) into the Humber Estuary in the period 2007-2017 (UK Environment Agency data), Figure Table S2: Summary of the oil slick release scenarios (tidal stage is with reference to Immingham station). Scenarios are repeated for both L1 and L2, Table S3: Wind speed and direction applied to the oil spill simulation, Figure S4: Observed (line) and simulated (points) free surface elevations at Immingham for the validation summer period (August 2016). Inset: Point correlation between observed and simulated tidal elevations (R 2 = 0.912), Figure S5: Example simulated oil spill trajectory, released from L2 at high water spring in winter, Figure S6 Table S5: Summary of maximum oil slick displacement (km) from point of release, Table S6: Relative properties of comparable oil slicks released at spring and neap tides (ratio = ST value/NT value), as measured 48 hours after the spill, Figure Table S7: Relative properties of comparable oil slicks released at high and low water (ratio = HW value/LW value), as measured 48 hours after the spill, Table S8: Ratios between impacted areas of oil slick released at L1 and L2. Acknowledgments: FES2014 was produced by Noveltis, Legos and CLS Space Oceanography Division and distributed by Aviso, with support from Cnes (http://www.aviso.altimetry.fr/-assessed: assessed on: 10/01/2018). The authors are grateful to Dr Hailiang Shen, P.Eng. for developing the script used to visualise the TELEMAC oil spill output in ArcMap. The authors are also grateful to the Association of British Ports (ABP) for the provision of Humber Estuary bathymetry data; the British Oceanographic Data Centre for provision of tidal data at Immingham station as well as salinity data; the UK Meteorological Office for wind data at Donna Nook; and the UK Environmental Agency for the provision of stage and water quality (temperature and salinity) data. The data was provided on a non-commercial license.

Conflicts of Interest:
The authors declare no conflict of interest.