Validation of the 3D-MOHID Hydrodynamic Model for the Tagus Coastal Area

Abstract

The hydrodynamics of the TagusROFI (Regions of Freshwater Influence) is affected by the coastal upwelling, the estuarine tidal flow, the thermohaline circulation that is modulated by the Tagus freshwater discharge, and by its complex bathymetry. The use of numerical models is the best way to explain the processes that characterize this region. These models are also crucial to answer important scientific and management questions. Nevertheless, the robustness of the products derived from models depend on their accuracy and therefore models must be validated to determine the uncertainty associated. Time and space variability of the driving forces and of bathymetry enhance flow complexity increasing validation difficulties, requiring continuous high-resolution data to describe flow and thermohaline horizontal and vertical variabilities. In the present work, to increase the precision and accuracy of the coastal processes simulations, the sub-systems coastal area and the Tagus estuary were integrated into a single domain, which considers higher resolution grids in both horizontal and vertical directions. The three-dimensiosal (3D)-MOHID Water model was validated for the TagusROFI by comparing statistically modelling results with in situ and satellite L4 data. Validation with a conductivity, temperature, and depth probe (CTD), an acoustic doppler current profiler (ADCP) and satellite data was performed for the first time. Validation against tidal gauges showed that the model is able to simulate tidal propagation inside the estuary with accuracy. A very good agreement between CTD data and surface sea water temperature (SST) and salinity simulations was observed. The validation of current direction and velocity from ADCP data also indicated a high model accuracy for these variables. Comparisons between model and satellite for SST also showed that the model produces realistic SSTs and upwelling events. Overall results showed that MOHID setup and parametrisations are well implemented for the TagusROFI domain. These results are even more important when a 3D model is used in simulations due to its complexity once it considers both horizontal and vertical discretization allowing a better representation of the heat and salinity fluxes in the water column. Moreover, the results achieved indicates that 3D-MOHID is robust enough to run in operational mode, including its forecast ability, fundamental to be used as a management tool.

1. Introduction

Modelling coastal systems particularly in regions of freshwater influence (ROFI), as it is the case of the Tagus estuary and its coastal adjacent area, is paramount to understand the physical and biogeochemical processes that characterize them. The simple measurement of environmental variables in situ alone or time-limited studies is not sufficient to translate and fully comprehended coastal processes since coastal systems are complex, decomposable, and large-scale [1]. Indeed, the combination and interaction of the physical factors (bathymetry; coastline topography; estuarine outflow; heat flux through the sea surface) and the associated physical transport and dispersal processes of the coastal flow field, makes coastal systems complex and unique in its hydrodynamics [2,3]. Therefore, and like other large-scale systems, the coastal system must be partitioned or decomposed into a number of small systems and the best way to gain insights into coastal system structure, organization, and functioning is through the use of numerical models. These models, besides allowing a better understanding of the interactions between the various components of the system, are also indispensable to set scientific scenarios and to answer important scientific and management questions [4]. Thus, only the proper knowledge of the processes that characterize the coastal systems, and their evolution in space and time, allow the design of robust coastal management plans. In this regard, hydrodynamic models assume a preponderant role since the establishment of currents and density fields is fundamental for transport models, whether they are Eulerian or Lagrangian, due to the importance of advection and/or diffusion in the study of issues related to the coastal system.

Among the numerical models that have been used to simulate the hydrodynamics of coastal systems (e.g., MARS3D [5]; MIKE3 [6]; SELFE [7]; SHYFEM [8]; TELEMAC [9]; MIKE 21/3 [10]; CSIRO [11]; Delft3D [12]), MOHID Water (Portuguese acronym of MOdelação HIDrodinâmica) is one of the most widely used. The MOHID Water model was developed by MARETEC (Marine and Environmental Technology Research Center). MOHID Water is open-source and most relevant key strengths are its ability to: Deal with two-dimensional (2D) and 3D simulations, with sigma, cartesian or lagrangian vertical coordinates; deal with eulerian or lagrangian transport references; use the same biogeochemical formulations independently of the number of spatial dimensions or space reference; and incorporate alternative formulations for every process (due to its modular approach). The versatility and flexibility of the MOHID Water eases its implementation in any type of system and to accomplish different modelling requirements depending on the objectives of the work to be performed. Indeed, in the last four years the MOHID Water was applied in multiple studies worldwide, namely in Canada [13], Colombia [14,15], Brazil [16], Argentina [17,18], Uruguay [19], Holland [20], France [21], Spain [22], Croatia [23], Australia [24,25], Malaysia [26] and Korea [27]. In Portugal, this model has been implemented in the entire Portuguese coast [28,29,30] and in all main estuaries [31,32,33,34,35].

Specifically, in the coastal system that comprises the Tagus estuary, most of the studies that used MOHID Water considered the coastal area and the Tagus estuary as sub-systems. In the case of the sub-system estuary, studies focused on sediment dynamics [36,37,38], on chlorophyll and/or nutrients [39,40], on microbiology [41,42], on water level [43], residence time [44] and CO₂ fluxes [45]. Concerning the sub-system coastal area, studies were directed to analyse the effects of the Tagus estuarine plume [45,46,47].

All the above mentioned works regarded the sub-systems estuary and adjacent coastal area as different domains due to the spatial resolution that is required to represent local physical processes [48,49]. However, computational advances allowed to integrate both sub-systems into a single domain, which considers higher resolution grids in both horizontal and vertical directions, essential to increase the precision and accuracy of the coastal processes simulations [50]. For this reason, only recently those two sub-systems were considered and modelled as an estuarine-coastal continuum [47] using a variable grid in the horizontal axis, in order to have a higher horizontal resolution within the estuary and within the adjacent continental shelf.

In order to use model outputs, models should be scientifically sound and robust, and therefore must be validated to inform users how much confidence can be placed in the products. Model validation is a process that aims to demonstrate that a given site-specific model is capable of making sufficiently accurate simulations [51]. Nevertheless, despite the importance of the validation process and although several statistical indexes have been developed and used for model evaluation [51,52,53], there is a lack of a common validation procedure [53]. Notwithstanding, the robustness of the models must be demonstrated which is only possible by comparing the model simulated results with in situ or remote observations such as satellite or HF radar products. Although the MOHID model is periodically validated, none of the works published until now focus exclusively on the validation of the 3D-MOHID for the TagusROFI. Only Vaz et al. [47] validated the coastal area circulation model through the comparison of the superficial velocity predictions with HF Radar data for this domain in the vicinity of the estuary mouth for 2012. These authors found that model predictions do not fit perfectly with the HF radar data probably due to the high spatial variability of velocity that changes in spatial scales shorter than the model cell size. Nonetheless, in 3D hydrodynamic models, validation requires the comparison of the model results with at least sea level, current speed, and current direction [54,55], and since the model application is baroclinic, validation should also consider seawater temperature and salinity. Janeiro et al. [30] validated a 3D application of the MOHID Water model for southern Portugal by comparing in situ measurements and HF Radar with predictions for all the five variables above. Although the TagusROFI domain has been applied in several research studies and as a tool for its management, this study will provide a higher insight concerning the validation of the results.

The present study aims to validate a 3D baroclinic application of MOHID Water model for the TagusROFI domain by using the above variables from a spatial-temporal perspective. With this purpose, model results will be compared with in situ measurements from several types of sensors (tidal gauges, CTD, and ADCP) and satellite L4 products (OSTIA, ODYSSEA, and MUR). The present work uses CTD, ADCP and satellite data to validate the TagusROFI domain for the first time. The main results are discussed, and conclusions are drawn concerning the model implementation in representing the coastal physical processes, its robustness to run in operational mode, including its forecast capability and consequently its ability to be a useful and powerful tool for coastal environmental management.

2. Materials and Methods

2.1. Study Area

The study area is located in the Portuguese west coast (30.6° N–9.5° W) and comprises a ROFI zone (Region of Freshwater Influence) (Figure 1) that encompasses the Tagus estuary, which is a highly populated region, subject to both natural and human pressures. This is the most extensive estuary of the Iberian Peninsula that occupies a volume of 1.9 × 10⁹ m³ and a surface area of approximately 320 km² [7,56]. This estuary is composed of a central body 30 km long and 10 km wide, connected to the Atlantic Ocean through a channel 12 km long and 2 km wide [57] forming a NE-SW oriented talweg [58]. Its width varies from 400 m at its head to 15 km in the central bay, with a depth range of 0 to 20 m (average depth of 5.1 m) [59], and presents extensive intertidal zones that occupy 40% of the estuary [60]. The coastal zone adjacent to the open ocean of the study area comprises 70 km of continental shelf that extends along the 39° N parallel, limited by the Cape Raso in the north and the Cape Espichel in the south, with the shelfbreak at 140 m depth [61]. It is a narrow platform that extends, on average, until 180 m depth and with a width ranging from 3 km (near the Cape Espichel) to 30 km [61]. This area is dominated by two important geomorphological features, the Cascais and Lisbon submarine canyons (the Cascais canyon is located NW of the Lisbon canyon and separated from the latter by the 1000 m high ridge of the Afonso Albuquerque Plateau [62]).

The Tagus estuary is a semi-diurnal mesotidal estuary, presenting an average amplitude of 2.4 m at the river mouth with the tidal range varying from 0.9 m to around 4.1 m during neap tides and spring tides, respectively [63]. According to Guerreiro et al. [7], river discharge may significantly influence water levels until 40 km upstream of the mouth, whereas downstream the levels are mainly controlled by tide and storm surges. The annual average river flow is 258 m³ s⁻¹, with monthly averages ranging from 6 to 2090 m³ s⁻¹ (2006–2018; National Water Resources Information System, SNIRH), and although being usually a well-mixed estuary, stratification may occur at high flow rates [64]. Flow is related to the amplitude of the tide which is also responsible for the mixture processes within the estuary. Indeed, the primary driver of the circulation in the Tagus estuary are tides, but other factors also influence it, such as river flow, atmospheric pressure, and the wind [65]. In addition, the estuary is strongly ebb-dominated due to the large extension of tidal flats [59].

In the adjacent area of the estuary, the configuration and orientation of the coastline protect the estuarine outflow from the waves, being only the southern part of the estuarine channel exposed to the swells [66]. Coastline geometry also influences plume dynamics, leading to a trajectory near the shore which is displaced offshore during northern winds events [46]. On the other hand, according to Jouanneau et al. [61], the continental shelf may also affect coastal currents, which in turn influences the mixture processes of the plume. Similarly, Fernández-Nóvoa et al. [58], stated that submarine canyons will modify the circulation on the shelf, whilst wind patterns in the ROFI are somehow controlled by the Sintra Mountains. In this area, predominant wind during winter is from northwest, north, and southeast (ranging from 1.5 m s⁻¹, 8.5% of occurrence, to 14 m s⁻¹, 7.5% of occurrence), whilst in summer winds from the north or from northwest prevail (ranging from 1 m s⁻¹, 5% of occurrence, to 15 m s⁻¹, 20% of occurrence) [47,67].

The general circulation on the western coast of the Iberian Peninsula presents a marked seasonal variability. Between May and September, the Portuguese west coast presents a generally equatorward flow (Portuguese Current; [68]) which is associated with the dominant regime of winds. During autumn and early winter, dominant winds are reversed becoming predominantly directed to E–N and a north poleward flow is observed (Portugal Coastal Counter Current; [69]). This current weakens when the northerly winds prevail developing conditions for the occurrence of seasonal upwelling (Portugal Coastal Current [70,71]), promoting an offshore advection of surface waters and vertical advection of rich nutrient upwelled waters toward the surface. To underline that the coastal area of the present study is characterized by significant and frequent upwelling events originated by jet-like flow extending more than 20 km seaward [72].

2.2. Model

The MOHID Water model assumes a hydrostatic equilibrium and Boussinesq approximation. To discretize the equations, it uses the finite volume approximation, so the discrete form of the equations that govern the flow can be applied macroscopically to the control volume (cell), permitting the equations to be independent of the cell geometry, allowing the use of a generic vertical coordinate [71]. According to Blazek [73], the main advantage of the finite-volume method over the finite-difference method is that the spatial discretization is carried out directly in the physical space, which avoids issues related to the transformation between the physical and the computational coordinate system. Moreover, the finite-volume method is very flexible and therefore can be easily implemented on both structured and unstructured grids [73].

The momentum balance equations for horizontal velocities are, in differential form and cartesian coordinates:

$$\frac{\partial u}{\partial t}+\frac{\partial \left(uu\right)}{\partial x}+\frac{\partial \left(uv\right)}{\partial y}+\frac{\partial \left(uw\right)}{\partial z}-fv=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial x}+\frac{\partial }{\partial x}\left(v_H\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(v_H\frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial z}\left(v_v\frac{\partial u}{\partial z}\right)$$

(1)

$$\frac{\partial v}{\partial t}+\frac{\partial \left(uv\right)}{\partial x}+\frac{\partial \left(vv\right)}{\partial y}+\frac{\partial \left(vw\right)}{\partial z}-fu=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial x}+\frac{\partial }{\partial x}\left(v_H\frac{\partial v}{\partial x}\right)+\frac{\partial }{\partial y}\left(v_H\frac{\partial v}{\partial y}\right)+\frac{\partial }{\partial z}\left(v_v\frac{\partial v}{\partial z}\right)$$

(2)

where, u, v and w are the components of the velocity vector in the x, y and z directions, f the Coriolis parameter, vH, and vv are the turbulent viscosities in the horizontal and vertical directions, p is the pressure. Assuming hydrostatic pressure the vertical momentum equation becomes an equation for pressure:

$$\frac{\partial p}{\partial z}+\rho g=0$$

(3)

And vertical velocity is computed using the continuity equation (assuming constant density, according to the Boussinesq approach):

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$

(4)

MOHID Water solves the seawater density using a non-linear state equation, depending on pressure, salinity and potential temperature using the algorithm of Millero and Poisson [74]. For the 3D momentum (zonal and meridional velocities), heat and salt balance equations in the vertical direction are computed implicitly while the horizontal directions are calculated explicitly for enhanced stability purposes [75].

The implementation of the study domain uses a one-way downscaling strategy of nested domains: the child high-resolution domain, TagusROFI, receives its open boundary conditions with 900 s temporal resolution from the lower resolution parent domain, the Portuguese Coast Operational Model System (PCOMS), which open boundary condition combines low-frequency Copernicus Marine Environment Monitoring Service (CMEMS) results (levels, salinity and temperature) with tidal levels computed by a 2D model with the same horizontal resolution (WestIberia domain). The downscaling strategy followed in this application was developed in the framework of the projects “European coastal-shelf sea operational observing and forecasting system” (ECOOP) and “Development and pre-operational validation of upgraded GMES Marine Core Services and capabilities” (MyOcean) and is subjacent to the CMEMS forecast modelling service rationale. This service provides large scale ocean forecasts of hydrodynamic and biogeochemical parameters which are refined locally by higher resolution models that can add the interaction with other forcing functions such as waves [9]. In our application, tide is the new forcing. The downscaling strategy using an intermediate regional model application between CMEMS and a local application has the advantage of running many local applications in parallel. In particular, PCOMS provides boundary conditions for all the major Portuguese estuaries such as the TagusROFI domain. In operational modelling this allows forecasts to be available in a smaller time period, as PCOMS can run separately from all its receiving child applications. Tide is imposed at the open boundary of the WestIberia domain using the version of the Finite Element Solution tide model (FES2004) global tide model [76,77]. Feedback into the parent domains is not considered. At the open boundary between nested domains, a Flather radiation method for the barotropic flow combined with a flow relaxation scheme to temperature, salinity, and total velocities was applied to the first 10 grid cells [78]. Turbulent diffusion coefficients are computed in MOHID using its embedded version of the General Ocean Turbulence Model (GOTM) [79].

The PCOMS and the TagusROFI domains applications are a fully 3D baroclinic hydrodynamic and ecological model, capable of simulating a wide range of processes (e.g., hydrodynamics, transport, water quality, oil spills) in surface water bodies (oceans, coastal areas, estuaries, and reservoirs). The PCOMS domain covers the Continental Iberian Atlantic coast and its contiguous ocean. This domain presents a constant horizontal resolution of 0.06° (≈5.2 km) populated with bathymetric information derived from the EMODnet Hydrography portal and covers the area comprised by the latitudes 45.00° and 34.38° N and the longitudes 12.60° and 5.10° W, resulting in a grid of 177 × 125 cells and maximum depths around 5300 m. The TagusROFI domain presents a variable horizontal resolution (200 m–2 km), which is an area that covers the interior of the estuary and the coastal zone (coast W–E) resulting in a grid of 121 × 146 cells (Figure 2) comprised by the range of latitudes 39.21° to 38.16° N and of longitudes 10.02° to 8.90° W. In the present study, it was used the best bathymetric information currently available, which comes from bathymetric data collected by the Hydrographic Institute (IH) between 1964 and 2009 for the Tagus estuary, and from the General Bathymetric Chart of the Oceans (GEBCO) for adjacent coastal area. The bathymetry was generated through a Delaunay triangulation using both data sets which are very consistent in the overlapping area.

Both PCOMS and TagusROFI domains were forced in the surface by meteorological models which provide hourly fields of surface wind, temperature, relative humidity, pressure, and solar radiation. PCOMS is forced with MM5 (PSU/NCAR mesoscale model) meteorological model with 9 km of spatial resolution, while the TagusROFI domain was forced by the WRF (Weather Research & Forecasting) meteorological model with 3 km of spatial resolution [80].

Both domains share a common vertical discretization consisting on a mixed vertical geometry composed of a sigma domain with 7 layers from the surface until 8.68 m depth, with variable thickness decreasing up to 1 m at the surface, on top of a cartesian domain of 43 layers with thickness increasing towards the bottom [77]. Intertidal areas computation is subjected to land masks which take the value 0 (the model does not compute forces for that cell) when the water column is below 20 cm. For the bottom friction, a rugosity of 0.0025 m² s⁻¹ was used.

River discharges in the TagusROFI application were imposed as volume discharges (without momentum). The Tagus river discharge using the Almourol hydrological station (39.22° N, 8.67° W) data, located 70 km off the head of the estuary, and obtained from the National Water Resources Information System (SNIRH) service with a frequency of 1 h. Other minor rivers (Trancão and Sorraia) were imposed using climatological values, with river flow ranging between 3 and 60 m³ s⁻¹ for Sorraia river and between 1 and 9 m³ s⁻¹ for Trancão river. Monthly climatological values of freshwater temperature and salinity for the three rivers and for water discharges from 14 wastewater treatment plants (WWTPs) were also imposed. Model configurations for the TagusROFI are summarized in

Table 1. Model configuration for the TagusROFI.

Settings	Level 1 –WestIberia	Level2–PCOMS	Level 3-TagusROFI
Model characterization	2D-Barotropic	3D–Baroclinic	3D-Baroclinic
Grid corners	33.50° N–49.90° N 1.00° W–13.50° W	34.38° N–45.00° N 12.60° W–5.50° W	38.16° N–39.21° N 10.02° W–8.90° W
Cells dimension	208 × 156	177 × 125	121 × 146
Bathymetry	EMODneta Hydrography portal	EMODneta Hydrography portal	Delaunay triangulation with IH Data and GEBCO
Horizontal Grid	Regular: (≈5.7 km)	Regular: (≈5.7 km)	Irregular: 200 m to 2 km
Vertical Grid	1 layer	7 Sigma Layer (0 m–8.68 m) 43 Cartesian layers	7 Sigma Layer (0 m–8.68 m) 43 Cartesian layers
Δt	60 seconds	60 seconds	6 seconds
Tides	FES2004b & FES2012c	From Level1	From Level2
OBC Water		From MercatorOcéan PSY2V4 (Releases 1–4)d	From Level2
Assimilation		Flow relaxation scheme of 10 cells with a time decay of 1 week at the open boundary and 0 inside the domain	Flow relaxation scheme of 10 cells with a time decay of 1 week at the open boundary and 0 inside the domain
OBC Atmosphere	MM5e (9 km)	MM5e (9 km)	WRFf (3 km)
Discharges	No	No	Tagus (hourly), Sorraia, Trancão and WWTP (monthly)
Turbulence		GOTMg	GOTMg
Bottom	Rugosity of 0.0025 m2 s−1	Rugosity of 0. 0025 m2 s−1	Rugosity of 0. 0025 m2 s−1

a—EMODnet Hydrography portal; b—[76], c—[81]; d—[82]; e—MM5: PSU/NCAR mesoscale model, [83]; f—WRF: Weather Research & Forecasting [80]; g—General Ocean Turbulence Model [79].

.

The actual implementation has been running operationally without interruption since 2011 providing a 3-day forecast once per day, for water velocity, water level, surface seawater temperature, salinity, and biogeochemical properties from a nutrient-phytoplankton-zooplankton- detritus (NPZD) model. The daily operational cycle consists of a 5-day simulation, starting the day before, in order to get the latest meteorological modelling results. Discharges are also updated every day with field data provided by the Almourol station thus providing the best flow value available. Forecasts are simulated using the last available river flow. The Level 3 configuration takes approximately 1 h and 30 min for each day of simulation under the configuration presented in Table 1.

2.3. Available Observations Data

Modelling performance was assessed through comparison with observations. All in-situ data and remote observations have advantages and disadvantages: The main advantage of tide gauges is that continuously record the height of the surrounding water level with a fine resolution, whereas the main disadvantage is its coarse spatial resolution in complex coastal areas. Similarly, the main advantage of CTD (salinity, temperature and depth data) is the frequency of data acquisition, whilst its disadvantage is the lack of spatial information. The ADCP (current velocity data) has the same advantages and disadvantages of a CTD but has as additional advantage the obtention of depth profiles along the entire water column. Satellite images has the advantage of providing daily surface maps (2D) of a region, but with very low frequency. This could be an issue for an area such as the TagusROFI that is dominated by small temporal scales. The data sets used for this purpose were the following:

2.3.1. Water Level

The water level was obtained from two tidal gauges, one located within the Tagus estuary (Lisbon harbor-LTG; 38.71° N and –9.13° W) and the other placed in the adjacent coastal area (Cascais-CTG; 38.69° N and −9.42° W) (Figure 1). Data from LTG was provided by the Hydrographic Institute, whereas data from CTG is available from the General-Directorate of the Territory web site. Distance between these tidal gauges was approximately 25 km (in a straight line). Modelling results were compared with data obtained in October 2012 from both tidal gauges with a time interval of regular 10-minute, comprising one spring and one neap tide.

2.3.2. Seawater Temperature and Salinity

Surface seawater temperature (SST) and salinity were recorded at 1 m depth every 10 min using a CTD (YSI 6600 V2) placed at the entrance channel of the Tagus estuary (38.69° N and –9.23° W) (Figure 1). Both variables were compared with the results obtained from MOHID, in the period between November 2012 and April 2013.

2.3.3. Current Velocity

Current velocity was registered using an ADCP (Teledyne Marine WorkHorse Sentinel, 600 kHz, WHS600 with four beams) placed at surface, off Cascais (38.67° N and –9.46° W) (Figure 1). This equipment recorded eastward and northward current, at 15 min sampling intervals, during the period from 2 to 17 July 2009.

The velocity vector was measured along the entire water column, every 1 m. In the present study, after analysing the profile of the current velocity, the water column was split into two layers from 2.5 and 15 m and from 15 to 30 m. In the top layer, current velocity is higher and dominated by the action of the wind, whilst in the bottom layer the effect of wind very reduced. For each layer, the mean value was obtained for the eastward and northward components as well as for the intensity. Modelling results were interpolated using the ADCP grid disregarding data obtained in the first 2.5 m and in the last 5 m water depth, to avoid bias due to the Doppler noise.

2.3.4. Seawater Temperature Based on Satellite Images

In order to compare the results obtained from MOHID with satellite images, three level 4 gridded RS data products were used, namely, OSTIA ([84], https://doi.org/10.5067/GHOST-4FK01), ODYSSEA (http://marine.copernicus.eu/; [85]), and MUR (http://doi.org/10.5067/GHGMR-4FJ01).

These products have a spatial resolution of 5, 2 and 1 km, respectively. Level 4 products are more accurate than other products since they combine complementary information from several satellite sensors and in situ instruments, being generated through statistical interpolation and temporal averaging [86]. According to Parkinson et al. [87], these products have also the advantage of providing gap-free gridded outputs. Each of these products, establish a reference depth used to correct bias and to integrate data into a single product. Daily SST satellite images were compared with the daily average of SST derived from the model (hourly outputs). The comparisons were performed for the period 2014–2016 and the statistical metrics correspond to the annual average obtained from the daily averages of SST.

2.4. Statistics

Various statistical indexes exist to evaluate the quality of the model’s calibration and validation results. Nonetheless, there is a lack of standard validation procedure. In the present study, observed and modelled results were compared using the correlation coefficient (r), root mean square error (RMSE), and BIAS index (mean error). The formulas that express these statistics followed the ones presented in Williams and Esteves [55] where n represents the number of model/observation data pairs, whilst M and O stand for model and observation, respectively:

$$r=\frac{{{\displaystyle \sum }}_{i=1}^n\left(M_i-\overline{M_i}\right)\left(O_i-\overline{O_i}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(M_i-\overline{M_i}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-\overline{O_i}\right)}^2}}$$

(5)

The correlation coefficient (also referred as Pearson product-moment correlation) measures the strength and direction of the relationship between two variables and can be either positive or negative depending on the directionality of the estimate. The correlation coefficient ranges from −1 to 1. A value of 1 indicates that the linear equation describes the relationship between M and O perfectly, whereas a value of 0 indicates that there is no linear correlation between the variables. Negative values indicate that M decreases as O increases.

$$BIAS = {\displaystyle \sum }_{i=1}^n\frac{1}{n_i}\left(M_i-O_i\right)$$

(6)

The model bias is a simple statistic, which measures the mean deviation between model-predicted values (M) and the observation data (O) and can be either positive or negative reflecting an over- or underestimation of observations by the model, respectively. Similar to RMSE, the smaller the absolute values of bias the better the agreement between model-predicted values and observation data.

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(M_i-O_i\right)}^2}$$

(7)

The RMSE is a generalized form of standard deviation that measures the deviation between model and observations in a least-squares sense. The square of the deviations ensures that negative as well as positive contributions are added, whereas the square-root restores the units to that of the original variable. RMSE has the benefit of penalizing large errors. Small absolute values of RMSE indicates a good agreement between model and observations.

3. Results and Discussion

3.1. Water Level–Time Series Data

Figure 3 shows the comparison between the results of the model with the data obtained for water level through the tide gauges for Cascais (CTG) and Lisbon (LTG) stations. The comparison interval refers to a full month (October 2012) covering one ebb and one flood tide period. The mean difference between modelling results and observations for CTG varied between −0.28 m and 0.26 m with an average of 0.003 m, whereas for LTG ranged from −0.24 m to 0.27 m with a mean average of 0.004 m (Figure 3). In both stations, the higher differences were registered during slack water. During these periods, the effect of the tide is less pronounced and thus the effect of wind on water level is evident. Model results showed an excellent agreement with the data measured in amplitude and tidal phase for both stations (CTG and LTG), reproducing accurately the measurements obtained in situ and for all periods of the tide. This is corroborated by the results of the statistics that were performed to evaluate the coherence between simulated and observed data for water level (Figure 3;

Table 2. Water level. Data summary and results of statistics used to assess the level of agreement between measured data and modelling results. n—number of observations; r—correlation coefficient; BIAS—Average bias; RMSE—Root mean square error.

Tide Gauge	Average (min–max) MOHID	Average (min-max) Tide Gauge	n	Pearson (r)	BIAS	RMSE
CTG	2.02 (0.39–3.73)	2.02 (0.31–3.69)	4445	0.995	0.003	0.17
LTG	2.01 (0.32–3.70)	2.00 (0-17–3.77)	4341	0.994	0.005	0.45

).

A strong correlation was observed for both stations between modelling results and observed data (CTG: r = 0.995; LTG: r = 0.994) which is in concordance with the RMSE error determined for both stations (Table 2). The BIAS determined was close to zero indicating that modelling results and observed values do not differ significantly. Altogether, these results highlight the capacity of the MOHID Water model to simulate the tidal propagation in a highly dynamic estuary mouth.

The tidal wave is amplified when it moves into the estuary as a consequence of the reduction of the depth which, along with the high velocities of the flow, causes the refraction of the tidal wave and consequently the increase of its amplitude. This phenomenon can be observed in Figure 4, which corresponds to a zoom of 7.5 h of Figure 3 that encompasses an ebb period. Data from CTG and LTG revealed a wave time-lag of approximately 20 min and an increase of the amplitude of about 15–20 cm between both stations. Both tidal wave time-lag and differences on the tidal wave amplitude between stations were well reproduced by the model (Figure 4).

3.2. Seawater Temperature and Salinity–Time Series Data

Between November 2012 and April 2013, seawater temperature and salinity were obtained at surface (1 m depth) with a CTD located in the area further downstream of the outlet channel of the estuary. These data were compared with a time series extracted from the model for the same local. During the study period, the average river flow was 522 m³ s⁻¹, ranging from 44 to 8711 m³ s⁻¹ (Figure 5). From the analyses of Figure 5, it can be seen that the observed and simulated results present a similar pattern for both variables. Daily fluctuations observed in SST and salinity were related to different river flows during flood and ebb periods. The high Pearson correlations coefficients obtained for SST (r = 0.91) and salinity (r = 0.86) reveal the good agreement between observations and simulated results (Figure 5;

Table 3. Surface seawater temperature and salinity. Data summary and results of statistics used to assess the level of agreement between measured data and modelling results. n—number of observations; r—correlation coefficient; BIAS—Average bias; RMSE—Root mean square error.

Parameter	Average (min-max) MOHID	Average (min-max) CTD	n	Pearson (r)	BIAS	RMSE
Temperature	14.2 (11.4–17.4)	14.3 (11.9–17.3)	21436	0.91	0.1	0.4
Salinity	28.0 (2.1–35.6)	27.4 (1.3–35.3)	21436	0.86	–0.9	2.9

).

The low values obtained for RMSE error and BIAS (Table 3) indicates that modelling results and observed values for SST are in good agreement. Regarding salinity, although the BIAS obtained was extremely low reflecting an overall low difference between observed and simulated results, the RMSE error determined, despite being considered low, was higher than 2 (Table 3). Model results were obtained with the operational forecasting model in which the river flow is daily updated. However, sometimes it occurs that the flow is not updated on the database in due time and therefore the model uses the last data available. Whenever this happens, a time-lapse between the real river flow and the one used in simulations occur. Time lapses are not uncommon but usually do not exceed 3 days. This may constitute an issue if during this period the river flow is atypical which will affect negatively forecasting. Indeed, this was verified during the period under study which may explain the low value of RMSE determined.

Although the salinity within the estuary depends on the entrance of freshwater, the main driver that controls the mixing processes is the tide. Considering that the tidal prism inside the Tagus estuary is 640 × 10⁶ m³, the river flow is only relevant to the non-linear mixing processes during high river flow events. At the end of the comparison period, an extreme event occurred where the mean daily river flow over three consecutive days attained 7500 m³ s⁻¹ which represents 25% of the volume of the tidal prism. This has led to the abrupt decrease of salinity (Figure 5).

Data from CTD put in evidence this extreme event, registering very low salinity values outside the estuary’s mouth (Figure 5), indicating that MOHID is well implemented for the estuary regarding transport and mixing processes. Furthermore, also reveals that the boundary conditions obtained from the Portuguese Coastal Operational Modelling System are appropriate.

3.3. Currents–ADCP Analysis

During the period of comparison, the prevailing wind was from the north quadrant with an average intensity of 3.4 m s⁻¹ (ranging from 0.1 to 9.2 m s⁻¹; data recorded in MS, Figure 1), and the mean flow rate of the river was 104 m³ s⁻¹ (varying between 27 and 451 m³ s⁻¹; data recorded in HS, Figure 1). Figure 6 compares raw data from ADCP with model results. There is a good qualitative agreement between the present model and ADCP measurements since predictions follow the same pattern of the observation data for current intensity and direction, as well as for both components of velocity (u and v). Periodical oscillations on the direction of the current due to tide are also perceptible in Figure 6.

ADCP data were obtained in an area that is protected from the North winds by the Sintra Mountains and weakly influenced by the Tagus plume due to the morphological characteristics of the estuary that leads it southwest which in turn is forced to move north due to Coriolis force. Nevertheless, measurements and modelling results showed that higher superficial current velocities had a southern direction, mainly due to the prevailing north winds during the study period. North winds, and consequently upwelling conditions, intensifies the surface velocity´s v component and forces it to move southwards almost during the entire period. The distribution pattern of the current direction results fundamentally from the tide that confers a harmonic behaviour to the flow. There is a mismatch in the direction of the prevailing current between the surface and the bottom, due to the variation of current velocity with the increase of depth which makes the fluid inertia to be different. Indeed, surface layers are subject to a stronger wind stress than the middle and bottom layers, since the influence of the wind decreases as depth increases. Thus, the water surface usually presents higher current velocities, which are directly related to wind direction and intensity. This can also be observed in Figure 7, where velocity profiles at each moment are presented for both observations and modelling results. The model reflects with high accuracy the individual and mean profiles obtained through ADCP.

Since the influence of the wind on the water column is higher in the top meters, the statistical analysis was performed for two depth layers (between 2.5 m–15 m and 15 m–30 m) (

Table 4. Current velocity and direction. Data summary and results of statistics used to assess the level of agreement between measured data and modelling results. n—number of observations; r—correlation coefficient; BIAS—Average bias; RMSE—Root mean square error.

	Average (min-max) MOHID	Average (min-max) ADCP	n	Pearson (r)	BIAS	RMSE
Depth (2.5 m–15 m)
Vel. Modulus (m s−1)	0.16 (0.0087–0.41)	0.17 (0.0054–0.43)	1440	0.68	0.014	0.10
Direction (rad)	151 (0.2–360)	134 (2.4–359)		0.63	–9.8	45
Velocity u (m s−1)	0.097 (−0.069–0.30)	0.099 (−0.14–0.34)		0.63	0.0019	0.092
Velocity v (m s−1)	−0.089 (−0.38–0.11)	−0.11 (−0.43–0.11)		0.73	−0.025	0.11
Depth (15 m–30 m)
Vel. Modulus (m s−1)	0.10 (0.0039–0.30)	0.12 (0.0025–0.34)	1440	0.62	0.018	0.065
Direction (rad)	151 (0.2–360)	164 (2–358)		0.66	–14	50
Velocity u (m s−1)	0.084 (−0.046–0.30)	0.092 (−0.086–0.29)		0.63	0.0082	0.066
Velocity v (m s−1)	−0.052 (−0.15–0.098)	−0.035 (−0.26–0.13)		0.71	–0.029	0.071

). The coefficient of correlation was relatively high and varied between 0.63 and 0.73 and between 0.62 and 0.71 for the first and second water column layer, respectively. For both layers, the RMSE and BIAS obtained for velocity modulus, and for the current components of the velocity (u and v) was always lower than 0.2 which indicates a statistically significant fit [55]. Regarding the current direction, the RMSE and BIAS determined indicate the occurrence of slight differences between observed and simulated data. According to Williams and Esteves [55], the minimum-level model performance should be ±15° for BIAS which is in accordance with the results obtained.

The results obtained clearly shows the reproducibility of the model in representing the 3D water column in the coastal zone, where nonlinear processes (interaction between wind/ tides/ oceanic circulation/ bathymetry interaction) become more important. Since vertical stratification of the velocity was well captured by this MOHID Water application, it can be inferred that bottom rugosity and vertical diffusion are correctly implemented in the model.

3.4. Seawater Temperature–Validation with Satellite

In the last decade, remote sensing data have been used to characterize coastal systems because they provide information on several variables with higher comprehensive spatial scale that any in situ equipment could provide. In the present study, an interannual comparison of the MOHID predictions with the three satellite products (OSTIA, ODYSSEA, and MUR) was performed, aiming at identifying SST patterns in the TagusROFI domain. The annual modelling results for the period 2014–2016 correspond to a daily average (of an hourly output) for the surface layer (1 m) which is represented in Figure 8 along with the three L4 products used. From the simple visual analysis of this figure it can be concluded that the best correspondences were obtained between ODYSEEA and MOHID, and between OSTIA and MUR.

To evaluate interannual differences in the upwelling process, a time series of the upwelling index was analysed, derived from two meteorological models. Upwelling index time series was provided by the Spanish Institute of Oceanography (IEO) and was calculated using the sea level pressure obtained from the Spanish Meteorology Agency (Meteogalicia) WRF atmospheric model and from the Fleet Numerical Meteorology and Oceanography Center (FNMOC) GFS atmospheric model. Data were obtained for a site located off Roca Cape. The upwelling index was considerably lower in 2014 (<100 m³ s⁻¹ km⁻¹) compared with the following years, which presented very similar annual average values (≈300 m³ s⁻¹ km⁻¹) (Figure 9). The surface seawater temperature drop in the upwelling region can be detected by remote sensing since the upwelled water is characterized by lower temperature comparing to the surrounding seawater. Although satellite products and MOHID results allow identifying the part of the coastal area where the upwelling process occurred, this phenomenon is more evident in ODYSEEA and MOHID, where a dark blue strip parallel to the N-S shoreline is clearly visible, which corresponds to upwelled water (Figure 8). Furthermore, the visual analysis of this figure also allows concluding that upwelling was less intense in 2014 than in the subsequent years analysed.

The Pearson correlation coefficients (r) obtained were very high and always superior to 0.86, indicating a very strong relationship between MOHID and the remote sensing products (

Table 5. Surface seawater temperature from satellite data. Data summary and results of statistics used to assess the level of agreement between remote sensing data and modelling results. n–number of observations; r–correlation coefficient; BIAS–Average bias; RMSE–Root mean square error.

Year	Satellite	Average MOHID	Average L4 products	n (per day)	Pearson (r)	BIAS	RMSE
2014	OSTIA	17.11	17.18	355	0.937	−0.064	0.846
	ODYSSEA	17.11	17.17	2095	0.948	−0.059	0.773
	MUR	17.12	17.20	8356	0.934	−0.078	0.894
2015	OSTIA	16.50	16.91	255	0.919	−0.407	0.946
	ODYSSEA	16.51	16.83	2095	0.924	−0.320	0.889
	MUR	16.56	16.87	8356	0.912	−0.359	0.992
2016	OSTIA	16.74	16.19	355	0.930	−0.176	0.866
	ODYSSEA	16.71	16.83	2095	0.864	−0.127	0.978
	MUR	16.73	16.83	8356	0.916	−0.102	0.914

). Regarding BIAS, the differences between modelling results and products data were always low and varied between −0.407 and −0.102 (Table 5). For all comparison, the RMSE error was always lower than 1 and higher than 0.7, which also indicates high accuracy of MOHID results. These results show that MOHID captures well the dominant pattern of the horizontal superficial distribution of the surface seawater temperature.

Values in Table 5 represent the mean of the statistical parameters for the entire study domain, but by plotting the spatial distribution of errors it is possible to observe where the differences are found.

Figure 10 shows that these differences (BIAS and RMSE) are more noticeable within the estuary and nearshore, in particular for OSTIA and MUR products that cannot capture the SST gradient. This can also be observed in Figure 8. MOHID and ODYSSEA present a similar distribution pattern of the differences, with BIAS always below 1 °C, including in the estuarine area.

Overall results indicate that the study area is characterized by the recurrent occurrence of upwelling event, which is well reproduced by the model. Therefore, it can be concluded that even in an extremely complex coastal area, as it is a TagusROFI domain, the MOHID Water model is adequately implemented and able to reproduce vertical stratification and horizontal patterns of SST. Moreover, of the three satellite products analysed in the present study, only ODYSSEA reproduces with accuracy upwelling events. Nevertheless, satellite data cannot be used to predict upwelling events as final filtered images are usually only available the next day but can be used to confirm an upwelling event represented by the model. On the contrary, 3D-MOHID Water is running operationally in forecasting mode, forecasting 3 days in every daily simulation.

4. Conclusions

Regions of freshwater influence (ROFIs) are interface areas between terrestrial freshwaters sources and the ocean, where physical processes characteristic of both estuaries and of coastal areas often occur and overlap, making these areas very complex systems. The spatiotemporal scales of these systems are so large and complex that the collection of in situ and/or spatial data alone cannot explain the processes that characterize them. This can be overcome using numerical models that despite being an approximation of the reality, provide a continuous representation in time and space of the system variables. More than giving a better understanding of the interactions between the various components of the system, these models are essential for the monitoring and management of the costal systems and associated resources. Nevertheless, the robustness of model-derived products to address specific scientific and management issues depends on its accuracy. Models assume properties for which no measurements are available and have known inherent uncertainties (e.g., model simplifications, mesh generation, roughness definition, etc.) [88]. Therefore, to minimize these uncertainties models are systematically calibrated by comparing modelling results with monitoring data. Apart from calibration, models should also be validated so as to determine the uncertainty associated with their solutions.

In the present work, the 3D-MOHID Water model was validated for the TagusROFI domain by comparing modelling results with in situ data (water level, surface seawater temperature (SST), salinity, and current direction and velocity) and SST remote sensing gridded L4 products (OSTIA, ODYSSEA, and MUR). These variables were chosen in the validation process because they are the most important hydrodynamics indicators with matching measured data. Indeed, the divergence of the horizontal flux (measurable through the propagation of the tide), in situ measurements of the velocity in water column, and the transport of conservative tracers (e.g., salinity) or of tracers with temporal variability lower than the tide (e.g., seawater temperature), are paramount parameters used to validate hydrodynamics. Seawater temperature and salinity are responsible for water density variability, which plays a central role in flow conditions in ROFI areas. Flow, however, can only be directly characterised with velocity fields and these can only be described through hydrodynamic models.

The statistical analyses used in validation include the calculation of the Pearson’s coefficient of correlation, BIAS index and RMES. Validation against two tidal gauges located at Cascais (coastal area) and at Lisbon (estuarine area) showed that the model is able to simulate observed water levels with accuracy, reproducing the amplitude and tidal phase, as well as its distortion during the propagation of the tidal wave within the estuary. Since the modelling results are accurate, it can be assumed that the model reproduces satisfactorily the horizontal flux divergence. A very good agreement between predictions of SST and salinity, and data from a CTD was also observed. The validation of current direction and velocity results using ADCP data indicated a high model accuracy for these variables. Finally, comparisons between model and satellite L4 products for SST showed that the model produces realistic SSTs and upwelling events. Although the use of more in-situ data would be desirable to increase spatial coverage, in the present paper it was used the most accurate data available. Notwithstanding, overall results showed that the 3D-MOHID Water setup and parametrisations were well implemented for the TagusROFI domain leading to highly accurate simulations. These results are even more important when a 3D model is used in simulations due to its complexity once it considers both horizontal and vertical discretization permitting a better representation of the heat and salinity fluxes in the water column. The validation of the 3D-MOHID Water is important because hydrodynamics is the basis for the transport of any biogeochemical variable through advection and/or diffusion processes, from the eulerian and langrangian point of view, and also indicates that it is robust enough to run in operational mode, including its forecast ability. Therefore, the 3D-MOHID can be considered an important management tool. Future research should be directed to the validation of the model’s biogeochemical component, and to the study of the biogeochemical processes in the TagusROFI domain. Moreover, it would be interesting to couple process-based-models (such as MOHID) with data-driven models [89] that employs machine learning techniques. This approach consists in feeding data-driven models with data from process-based-models that were already validated. Although this hybrid approach is still in its infancy in what concerns hydrodynamics due to the complexity and variability of the inherent processes, it could be a useful approach in the coming future. The results of the present work constitute an important step towards the implementation of this innovative approach.