1. Introduction

The hydrodynamics and salinity dynamics play a major role in estuarine ecosystems. For instance, they control the exchanges of waterborne material with the adjoining oceans and the residence times, determine the spatial distribution of biota and drive the sediment dynamics. Understanding the water and salinity dynamics is thus a crucial step in the understanding of estuaries.

Water circulation and salinity dynamics are strongly coupled through complex mechanisms ^[1]. First-order mechanisms, such as tidal propagation, are well understood and numerical models have long been able to accurately reproduce them. Yet, the long-term dynamics can be significantly affected, if not controlled, by second-order processes, which play a major role in the generation of residual currents (e.g., ^[2]). Also, vertical mixing and stratification are controlled by small-scale processes. As a result, accurately modeling estuarine dynamics remains a challenge, in particular for stratified conditions ^[3].

Simultaneously, obtaining observational datasets extensive enough to provide a clear understanding of the dynamics of a particular estuary is often too costly due to the space and time variability of estuaries at different time scales. Hence, there is a clear need for detailed and extensive validations of three-dimensional baroclinic estuarine models that highlight their accuracy.

Thus, the goal of this paper is to present a detailed validation and assessment of a three-dimensional baroclinic model of the Tagus estuary (Portugal) using the modeling system SCHISM. The model validation was performed for very distinct conditions in order to evaluate the robustness of the parametrization. Moreover, the influence of the model setup (e.g., vertical grid resolution, advection scheme) and forcings (e.g., waves, atmospheric forcing) in the simulation results, in particular on the ability to reproduce the stratification, was analyzed. The model is then used to provide a preliminary description of salinity distribution in this estuary for extreme river discharges (floods and droughts).

The Tagus estuary (Figure 1), located in the western Portuguese coast, is one of the largest estuaries in Europe with an area of about 320 km² (e.g., ^[4]). The area surrounding the Tagus estuary is intensively used and the estuarine margins support diverse uses and activities, such as urban, industrial/harbors, agriculture and shellfish harvesting. The ecological and natural values of the Tagus estuary, associated with the wetlands, the terrestrial habitats and the biological diversity, are well recognized. The Tagus estuary hosts a natural reserve, the Tagus Estuary Natural Reserve, which is one of Europe's most important sanctuaries for wintering or staging of birds and is also an important nursery area for fish and shellfish juveniles. The western and northern margins of the estuary are densely urbanized, contrasting with the productive agricultural areas along the eastern side ^[5]. Human activities in the intertidal zones and along the estuarine margins are the main driver for the loss of natural areas and estuarine shoreline changes during the last 60 years ^[6]. Many studies showed the importance of the Tagus estuary, in terms of sediments, nutrients, contaminants, plankton and fishery dynamics, and of its interaction with the adjacent coastal area (e.g., ^{[7,8,9,10,11,12]}).

The hydrodynamics and circulation patterns in the Tagus estuary have been addressed in several studies based on data analysis or numerical modelling. In general, much is known about tidal propagation, water levels and general circulation patterns in the Tagus estuary, both from the analysis of a large tide gauge dataset from 1972 ^[13,14] and from the application of several numerical models (e.g., ^{[4,13,15,16,17,18]}). However, the present knowledge of estuarine circulation, salinity dynamics and low-frequency motions in this estuary remains limited. Detailed datasets on salinity during stratification events (e.g., ^[14]) and three-dimensional baroclinic applications in this estuary (e.g., ^[19,20,21]) are still scarce. Developments towards a more accurate three-dimensional baroclinic numerical model of the estuary are still recent ^[22], and are further extended and assessed in this paper.

---

2. Materials and methods

---

2.1. Characterization of the hydrodynamics of the study area

The Tagus estuary is included in the Metropolitan Area of Lisbon and comprises about 1.6 million of total inhabitants along its margins. The detailed land-use cartography of the Tagus estuarine fringe, that covers 130 km², shows that agricultural areas (35%), urban zones (34%), and industrial, port and airport facilities (24%) are the main uses along the margins ^[6]. In the intertidal area, occupancy includes ^[6]: salt marshes (13%), anthropogenic structures (e.g., salt pans, old tide mills or aquaculture installations, 15%) and beaches (1%).

This mesotidal coastal plain estuary has a complex morphology with a deep, long and narrow inlet connecting the Atlantic Ocean to a broad and shallow inner basin with extensive tidal flats and marshes (Figure 1). About 40 km upstream, the estuary markedly narrows at the bay head. The intertidal area represents about 43% of the total estuarine surface ^[23].

The circulation in the Tagus estuary is primarily driven by tides but it is also influenced by the river flow, wind, atmospheric pressure and surface waves. Tides are predominantly semi-diurnal, as in the whole Portuguese coast. The tidal form number (ratio between the amplitudes of the two major diurnal constituents and the two major semi-diurnal tidal constituents) is close to 0.1. Tidal ranges at the nearest coastal gauge (Cascais) vary between 0.55 m and 3.86 m ^[4] but resonance significantly amplifies this range within the estuary ^[13,15]. With a resonance period between 8 to 9 hours, the Tagus estuary amplifies the 4^th-diurnal and semi-diurnal constituents, while leaving the amplitude of diurnal constituents almost constant ^[4]. As a result, the tidal form number decreases along the estuary, enhancing the diurnal inequality of the tides. The estuary is strongly ebb-dominated due to the large extent of the tidal flats ^[13]. In its central section, ebbs are shorter than floods by over one hour. In contrast, the upper estuary is slightly flood-dominated.

River discharge may significantly influence water levels, but only farther than 40 km upstream of the mouth ^[18]. Downstream, the water levels are mainly controlled by tides and storm surges. The Tagus River is the main tributary of the Tagus estuary with a mean river discharge of 336 m³/s ^[24]. Two other rivers contribute to the water inflow to the estuary: the Sorraia River (about 5% of the Tagus river discharge ^[25]) and the Trancão River.

Wind, atmospheric pressure and waves affect the Tagus estuary dynamics less than tides, and are therefore less well studied. Atmospheric pressure is more important than wind in determining the surge within the estuary ^[26]. Oliveira ^[27] applied a wave model to the estuary mouth and found that short waves were dissipated in the ebb-delta sand bars at the entrance of the estuary. At the estuary mouth, the currents increase the significant wave height and the mean wave period on ebb and decrease them on flood ^[28]. The wave-induced setup can reach tens of centimeters even in the upper reaches of the estuary, with a significant tidal modulation ^[29]. Inside the estuary, locally-generated waves can reach significant wave heights of 0.8 m ^[30,31]. Wind waves can affect the salinity dynamics of estuaries. For instance, mixing created by waves can weaken estuarine circulation ^[32], and longshore currents can transport the estuarine plume away from the inlet, thereby promoting the penetration of ocean water in the estuary on flood. To our knowledge, the effect of wind and waves on the Tagus salinity dynamics has not been studied.

The estuary is often considered well-mixed, even though stratification has been observed at high flow rates ^[14]. To estimate the conditions for which the estuary could be stratified, the empirical model of ^[33] was used, following ^[34]. According to ^[33], an estuary is well-mixed for Ur/Ut < 0.01, stratified for Ur/Ut > 0.1, and partially mixed otherwise. Here, Ur is the river velocity defined as the river discharge divided by the cross-section and Ut is the tidal velocity defined as the mean velocity due to tides. Both Ur and Ut were estimated along the estuary using the results of the model of ^[13] for neap, mean and spring tides conditions, and for a range of river flows. Results indicate that the estuary is well-mixed for average and lower river flows, but that it can be partially-mixed for higher river flows, and even stratified for extreme river flows, particularly under neap tide conditions (Figure 2). For river flows lower than 1000 m³/s, ^[35] and ^[36] also estimated that the Tagus estuary is well-mixed.

---

2.2. Numerical model

---

2.2.1. SCHISM modeling system

Simulations were performed with the modeling system SCHISM (^[37], http://ccrm.vims.edu/schism/), version 5.3.1. This system evolved from SELFE (^[38], version 3.1dc), with many upgrades including a new extension to large-scale eddying regime and a seamless cross-scale capability from river to ocean ^[37]. SCHISM (Semi-implicit Cross-scale Hydroscience Integrated System Model) is a community modeling system based on unstructured grids in the horizontal dimension and hybrid SZ coordinates or a new LSC² (Localized Sigma Coordinates with Shaved Cell) in the vertical dimension. This modelling system is fully parallelized and uses semi-implicit finite elements and finite volume methods, combined with Eulerian-Lagrangian methods, to solve the shallow water equations. The numerical scheme also naturally incorporates wetting and drying of tidal flats. SCHISM includes modules for surface waves ^[39], morphodynamics ^[40,41] and water quality ^[42,43,44]. Further details about SCHISM's physical and numerical formulations can be found in ^[37].

SCHISM was applied in three-dimensional baroclinic mode to the Tagus estuary. SCHISM and its predecessor SELFE have been extensively used to simulate the three-dimensional estuarine circulation (e.g., ^{[20,42,45,46,47]}). Detailed assessments of the ability of SELFE to model stratified tidal flows are described in ^[3] and ^[48].

---

2.2.2. Data for model calibration and validation

Water levels measured at 13 tidal gauges installed along the estuary in 1972 were used to calibrate the model (Figure 3). These time series were harmonically analyzed to extract the amplitudes and phases of the tidal constituents ^[13]. Time series were obtained by harmonically synthesizing the tidal elevations using 11 tidal constituents (Z0, MSF, O1, K1, N2, M2, S2, M4, MS4, M6, 2MS6). Water levels measured at the Cascais tidal gauge (Figure 3) in 1983 and 1988 were also used to validate the model.

Salinity and water temperature were validated with data from previous field campaigns in the estuary ^[14,49]. Two simulation periods were used to validate the model, covering distinct environmental conditions. Data included salinity and water temperature observations obtained in February 1988, a period during which the observed maximum daily-mean discharge in the Tagus river was about 2000 m³/s. Vertical profiles of salinity and water temperature were measured along three longitudinal profiles (BC: Barra-Corredor, CN: Cala do Norte, CS: Cala de Samora; Figure 3) during high-tide and low-tide. Measurements were made from the surface to the bottom at 1 m intervals. The 1988 simulations were performed between January 15 and February 15. Monthly observations of salinity acquired in 1983 were also used to validate the model. These data, which cover a broader area of the estuary, were measured at six sampling stations (Figure 3). The mean Tagus river discharge during these surveys was 150 m³/s (the maximum discharge was 345 m³/s). The 1983 simulations were performed between January 1 and December 31.

---

2.2.3. Model setup

An horizontal grid composed of triangular elements was generated with the softwares xmgredit ^[50] and nicegrid ^[51], based on grids from previous applications (e.g., ^[4,29]). The domain is over 110 km long and has a surface area of 1442 km². It extends from the ocean, 27 km away from the estuary mouth, to the river (Figure 3). An analysis of the estuary plume shown in satellite images ^[52] indicates that the plume is normally contained within the domain. An a posterior verification of the extent of the domain showed that for typical river flows the estuary plume does not reach the ocean boundary. For the largest river flow simulated herein, which corresponded to a return period of 50–100 years ^[53], the minimum depth-averaged salinity reaching the ocean boundary is about 30. Under those circumstances, the model may tend to slightly overestimate the inflow of salinity through the boundaries. Further details regarding the grid generation can be found in ^[29,54]. The resulting grid has 14,000 nodes and a typical resolution of 15–25 m, which is lower in the channels to allow the proper propagation of the tide. For simulations in which the inundation of areas protected by dykes is not expected, a smaller grid was generated to reduce the computational costs. The smaller grid has 83,000 nodes.

The vertical grid was setup based on preliminary simulations in which three hybrid vertical grids were tested: 20 SZ levels grid (15 S levels and 5 Z levels), 39 SZ levels grid (30 S levels and 9 Z levels) and 54 SZ levels grid (45 S levels and 9 Z levels). Taking into account both preliminary data-model comparisons and the computational times, in most runs the vertical domain is discretized with a hybrid grid with 39 SZ levels: 30 S levels in the upper 100 m, and 9 Z levels between 100 m and the maximum depth. Further results and discussion regarding the vertical grids tested are presented in section 3.

Three open boundaries were considered in the simulations. The model is forced by surface water elevations at the ocean boundary and river flow at the river boundaries (Tagus and Sorraia rivers). Salinity and water temperature are also imposed at the open boundaries. Tidal elevations at the ocean boundary are imposed from an application of SCHISM to the NE Atlantic Ocean ^[55]. The Tagus river flow was established based on the data available at SNIRH (http://snirh.pt), namely at the Ómnias and Almourol stations. For lack of data, the river flow in the Sorraia River was taken as 5% of the river flow in the Tagus River, based on the ratios between annual averages in the two rivers.

The bottom stress was parameterized using a drag coefficient. Two approaches were evaluated to establish the drag coefficient: (ⅰ) drag coefficient based on a previous applications of SELFE 3D in the Tagus estuary ^[19,20,22]; (ⅱ) and drag coefficient determined from the Manning coefficient adopted for a 2D depth-averaged model application ^[54]. To evaluate the influence of the drag coefficient in the model results, preliminary simulations were performed for a period of 20 days. These simulations were forced by 11 tidal constituents (section 2.2.2.) at the ocean boundary and by river flows of 517 m³/s and 20 m³/s at the Tagus and Sorraia boundaries, respectively. The time step was set to 30 s and the Generic Length Scale KKL turbulence closure scheme, with the Kantha and Clayson's stability function, was used. Data and model results were both harmonically analyzed and synthesized. The comparison of the model results with the 1972 water levels data showed that the second approach led to the smallest RMSE, and this approach was retained hereafter.

In all validation simulations the oceanic boundary was forced with 23 tidal constituents (Z0, SSA, MM, MF, MSF, O1, K1, P1, Q1, N2, M2, S2, K2, 2N2, Mu2, Nu2, L2, M3, MN4, M4, MS4, M6, 2MS6). Daily-mean discharges were used at the riverine boundaries. Initial conditions of salinity were set to decrease from 36 in the coastal area to 0 at the river boundaries. Salinity was set as constant at the open boundaries (36 at the oceanic boundary and 0 at the riverine boundaries). Initial conditions for temperature were set as constant. For the 1988 simulation, water temperature at the boundaries was set at a constant 15 ℃ at the oceanic boundary, based on daily data measured at the Cascais tidal gauge, and 13 ℃ at the riverine boundaries, estimated from monthly data measured at the SNIRH stations. For the 1983 simulations, since no data were available for the periods simulated, water temperature was forced with climatological time series. The water temperature was estimated with the Atlantic–Iberian Biscay Irish–Ocean Physics Reanalysis (^[56], http://marine.copernicus.eu/) and with SNIRH data at, respectively, the oceanic boundary and the riverine boundaries.

The heat exchange with the atmosphere is computed using the model of ^[57]. The 1983 and 1988 simulations were forced with two different atmospheric datasets, NCEP-NCAR Reanalysis (~50 km/6h resolution, provided by the NOAA/OAR/ESRL PSD, http://www.esrl.noaa.gov/) and BINGO project dataset (12 km/3h resolution, ^[58]). For the simulated periods, local atmospheric observations were not available. The comparison between the two datasets aimed at evaluating their influence on the salinity and water temperature results.

Two numerical schemes for the simulation of the advection of salinity and temperature, upwind and TVD-Total Variation Diminishing (both available in SCHISM), were evaluated. The influence of these numerical schemes was analyzed for the 1988 period.

In order to assess the potential importance of wave effects on salinity in the Tagus estuary, the simulation 1988 was repeated considering wave effects. Waves were simulated with WWM fully coupled with SCHISM ^[39]. The wave spectra were discretized with 24 directions and 25 frequencies and boundary conditions were provided by an application of WaveWatch Ⅲ ^[59] to the North Atlantic and the Portuguese shelf ^[53]. Further details on the application and validation of the coupled SCHISM-WWM are given in ^[29].

A summary of the simulations performed is presented in Table 1.

---

2.2.4. Model skill assessment

Different skill indices were computed to evaluate the model performance in the simulation of water levels, salinity and temperature. These metrics included the mean absolute error (MAE), the root mean square error (RMSE) and the Wilmott Skill Index of agreement (SKILL ^[60]):

$MAE =\frac{1}{N}\mathop \sum \limits_{n =1}^N \left| {{P_n} - {O_n}} \right|$

(1)

$RMSE ={\left[ {\frac{1}{N}\mathop \sum \limits_{n =1}^N {{\left( {{P_n} - {O_n}} \right)}^2}} \right]^{0.5}}$

(2)

$SKILL =1 - \frac{{\mathop \sum \nolimits_{n =1}^N {{\left| {{P_n} - {O_n}} \right|}^2}}}{{\mathop \sum \nolimits_{n =1}^N {{\left( {\left| {{P_n} - \bar O} \right| + \left| {{O_n} - \bar O} \right|} \right)}^2}}}$

(3)

where P is a given model result, O is a given observation, N is the total number of points used, and Ō is the average of the observations.

Following ^[61] and ^[62] target diagrams were also used to assess the model behavior. Target diagrams provide summary information about pattern statistics using unbiased RMSE (on the x-axis) and bias (on the y-axis), both normalized by the standard deviation of the observations (U-RMSE* and BIAS*). The RMSE is also multiplied by the sign of the difference between the standard deviations of the model and the observations. These diagrams thus provide information about whether the model standard deviation is larger (X > 0) or smaller (X < 0) than the observations' standard deviation, and if the bias is positive (Y > 0) or negative (Y < 0). The closer the points are to the origin, the better the model performs. Since the information in these diagrams is normalized, they also allow the comparison of the model performance between different variables.

---

2.3. Estuarine classification for extreme river discharge events using the Venice system

The Venice system (after ^[63]) is often used to classify waters with respect to the salinity (e.g., ^[34,64,65]). This system considers six classes: limnetic (freshwater, < 0.5), oligohaline (0.5–5), mesohaline (5–18), polyhaline (middle area, 18–25; lower area, 25–30), and euhaline (30–35).

To evaluate the role of the freshwater discharge in the salinity distribution in the Tagus estuary, namely during extreme situations regarding the river flow, two scenarios were established. These scenarios were based on two distinct periods, which were chosen due to their extreme characteristics regarding river discharge: July 2005 (low river discharge scenario) and February 1979 (high river discharge scenario). The estimated average Tagus river flow was 22 m³/s in July 2005 (estimation based on ^[66]), although some uncertainty remains since no observations are available. In contrast, in February 1979 the average river flow was 5464 m³/s at Ómnias (near the Tagus riverine boundary of the model) ^[67]. The estimated return period of the 1979 flood in the lower Tagus river (Vila Velha de Ródão) is 50–100 years ^[53]. Daily-mean river discharges were used in the simulations, with maximum values reaching about 40 m³/s in July 2005 and about 13,900 m³/s in February 1979. Simulations were performed for 40 days. Horizontal salinity classification maps were performed using model average results for the last 30 days of simulation.

---

3. Results and discussion

---

3.1. Tidal propagation

Table 2 presents a comparison between the RMSE obtained for the present application and the previous three-dimensional baroclinic application of SELFE in the Tagus estuary by ^[19]. For the 1983 and 1988 simulations the water levels errors at the Cascais tidal gauge are of the same order of magnitude as those obtained with the 1972 data. Results show that the model represents adequately the tidal propagation. At most stations, particularly the ones located upstream, the errors obtained are lower than the ones from previous applications (Table 2). The improvements obtained may result from the higher horizontal grid resolution, namely upstream, and the use of more detailed bathymetric datasets ^[29], both contributing for a better propagation of the tide, and from the model itself.

---

3.2. Salinity and water temperature: horizontal gradients and stratification

The comparison between the observations and the model results for the 1988 simulations suggests that the model is able to represent the main stratification patterns during high river flow periods (Figures 4–7). The model, however, tends to overestimate the vertical mixing of salt. For the same vertical grid, the vertical mixing increases slightly in the simulation forced by waves. The underestimation of stratification had already been observed in a previous application ^[19], but significant improvements were achieved herein. A satisfactory agreement was also reached for temperature, with significant improvements compared to the application by ^[19], in particular regarding the vertical stratification. Evaluation simulations showed that the setup of the vertical grid and the numerical method used to solve the advection of the tracers were particularly important for the accurate representation of the vertical patterns. The increase of the vertical grid resolution, from 20 SZ to 39 SZ levels, contributed to reduce the simulated mixing in the water column (Figure 8.). The use of the TVD numerical scheme, which is less diffusive than the upwind, significantly improved the simulated stratification along the estuary, in particular during low-tide (Figures 4–7).

For the 1988 period, depth-averaged salinity and depth-averaged temperature at each station also provide information about the model's ability to represent the horizontal gradients in the downstream and central areas of the estuary. Thus, skill indices and target diagrams between the observations and the model results were also computed for the depth-averaged salinity and depth-averaged water temperature using the data from all stations (Figure 9 and Table 3). Salinity errors are about 2.0, except for the simulation with the lowest vertical grid resolution for which RMSE is about 3.0, showing the need for a high vertical grid refinement. Regarding the water temperature, errors ranged from 0.3 ℃ to 0.6 ℃. The largest error was obtained with the BINGO's atmospheric forcing, which leads to an underestimation of the water temperature at the surface, and a lower performance of the model (Figure 9). These results highlight the role of the atmospheric forcing in simulating the heat exchange between the water surface and the atmosphere and, consequently, in the simulated water temperature. Overall, target diagrams suggest a good performance of the model. In particular, the use of the TVD method leads to a more accurate representation of the salinity fields. However, the TVD method increases the run time twofold relative to the use of the upwind method. This computational cost prevents the use of the TVD method for very long simulations or for the simulation of a large number of scenarios. The combined use of the upwind method with a higher resolution vertical grid could improve the computational cost; however, the model still significantly overestimates the vertical mixing (Figure 8) comparatively to the use of the TVD method (Figure 5).

The 1983 simulations allow the assessment of the model over 1-year in a broader area of the estuary, in particular in the upper estuary, and for periods of lower river discharge. Results indicate that the model is able to represent the main spatial and temporal patterns observed (Figure 10). The global RMSE and MAE (i.e., considering the entire set of observations) are, respectively, about 3 and 2 for salinity and about 1.5 ℃ and 1 ℃ for temperature. The skill index is 0.98 for salinity and 0.95–0.97 for temperature. Differences between the model results using the NCEP-NCAR and BINGO atmospheric forcings are small, although the model tends to represent the salinity better with the BINGO atmospheric forcing and the temperature with the NCEP-NCAR atmospheric forcing (Figure 11, Table 4). A local analysis at each station suggests that the differences between the model results and the observations tend to be larger in the upstream stations, in particular regarding the mean errors (Table 4). These differences may be due to the large uncertainty regarding the bathymetric data in these areas, which are shallower and thus more affected by variations in the bathymetry. There is also some uncertainty relative to the time and location of the samplings, which may contribute to the differences observed. The model skill at the downstream station (station 8.0) is lower, although mean errors are about 1 and indicate a good performance of the model (the mean salinity in this station is about 35). The lower skill in this station is due to the small variance in the observations and may be related with a bias in the oceanic boundary condition. Moreover, there is also some uncertainty in the observations themselves, related to the methods and equipment used at the time.

Overall, results indicate a good performance of the model, within the range obtained in recent similar modelling applications (e.g., ^[47,62]).

---

3.3. Salinity classification for extreme river discharge events

Results show a pronounced progression of the saltwater upstream for low river discharges and a contrasting situation for high river discharges (Figure 12). For the low river discharge scenario most of the estuary exhibits polyhaline or euhaline characteristics (salinity > 18), while for the high river discharge scenario the mesohaline area extends downstream until the narrow area of the estuary (near Lisboa, Figure 3).

For the high river discharge scenario, the simulated mean salinity at Cascais was 31, with a range between 14 and 36, the minimum value observed during the peak of the Tagus river discharge (13900 m³/s). The large volume of freshwater entering the estuary significantly reduces the salinity in the downstream area, typically with salinities between 30 and 35, and a pronounced stratification is observed near the inlet (Figure 13). Note that in the riverine area the marginal inundation is not simulated due to both the lack of topographic data to extend the model domain and very demanding computational costs. The marginal inundation in the upstream area of the domain can contribute to reduce the peak of the freshwater flow reaching the downstream area of the estuary and influence the estimated salinity. However, salinities lower than 5 (unpublished data) were observed between Pedrouços and Lisboa (Figure 3) during a period of high river discharge in the beginning of April 2013 (daily-mean river discharge ranging from 3500 to 7500 m³/s), which supports the present results.

Comparatively, for the low river discharge scenario, the simulated mean salinity at Cascais was 36 (range: 35–36) and the estuary is well-mixed (Figure 13). In the upstream area, in particular near the agricultural lands where a major water uptake for irrigation is located (Conchoso), simulated mean and minimum salinities are about 7 and 3, respectively. Although there is a significant lack of salinity data during drought periods in the estuary, particularly upstream, a brief comparison between model results and a few observations near a water uptake for water supply (Valada) in July and August 2005, suggest a reasonable agreement between observations and model results with differences lower than 1. The upstream progression of the salinity during low river discharges may affect the ecosystem dynamics and also impact negatively some of the economic activities in the upstream area of the estuary. In particular, it may prevent the water uptake for irrigation of agricultural land (which requires salinity < 1). Some constraints have already occurred in the past, namely during 2005 when Portugal was under severe drought conditions. Sea level rise and more severe droughts may exacerbate this problem in the future.

---

4. Conclusions

Understanding of how changes in the hydrodynamics and salinity dynamics affect the estuarine ecosystems dynamics, uses and activities is fundamental to support the management of these areas. However, water circulation and salinity are driven by several processes (e.g., tides, wind, freshwater discharge) and are strongly coupled through complex mechanisms, which makes their accurate modelling a challenge. In this study, a three-dimensional hydrodynamic baroclinic model was implemented and assessed in the Tagus estuary in order to improve the representation of such dynamics in this estuary.

The model was implemented using the modeling system SCHISM and assessed for contrasting environmental conditions (i.e., periods of high and low river discharges) in order to evaluate the robustness of the parametrization. Overall, results show the ability of the model to represent the main horizontal and vertical salinity and temperature gradients in the estuary over different periods.

The analysis of the influence of the model setup and forcings in its ability to reproduce stratification suggests that the vertical grid resolution and the advection scheme are of major relevance. When the river discharge is high enough to lead to stratified conditions, using the TVD numerical scheme improves the model performance, in particular in the representation of vertical gradients. The horizontal gradients and the mean variations along the estuary are adequately represented by both the upwind and the TVD methods. Since the use of TVD method is computationally more demanding, which may prevent its use for several applications (e.g., very long simulations, simulation of a large number of scenarios or in operational models providing forecasts), these findings suggest that the upwind method can be used in well-mixed conditions for computational efficiency and the TVD method should only be used when required, in particular, when stratification is expected. To support this approach further investigations are required to better understand the stratification conditions in the Tagus estuary.

Under extreme river discharges (floods or droughts) the salinity distribution along the estuary shows very contrasting characteristics. During low river discharge, in particular, salt water progresses significantly upstream, which may negatively affect agricultural activities in the upper estuary's margins. Moreover, these conditions may be aggravated by climate change, namely sea level rise, as suggested by preliminary results ^[22] and this requires further research.

The model developed herein will thus be used in further research studies regarding the evaluation of the effects of climate change in the Tagus estuary dynamics. The model will also be further extended to study the biogeochemical dynamics in this estuary and integrated in a real-time observatory, providing information that can support both the daily and the long-term management of the Tagus estuary.

---

Acknowledgements

This work was funded by the European Commission through the H2020 project BINGO (Grant Agreement Number 641739) and by the FCT—Fundação para a Ciência e Tecnologia project UBEST—Understanding the biogeochemical buffering capacity of estuaries relative to climate change and anthropogenic inputs (PTDC/AAG-MAA/6899/2014). The first author was partly funded by a post-doctoral grant from FCT—Fundação para a Ciência e a Tecnologia (SFRH/BPD/87512/2012). This work took advantage of the Portuguese National Distributed Computing Infrastructure (INCD). The authors would like to thank Dr. Paula Freire for the pre-processing of topographic data and useful discussions throughout the implementation and validation of the numerical model, Dr. Elsa Alves and Dr. Rui Rodrigues for useful discussions about the river discharges, Prof. Isabel Ambar and Prof. Vanda Brotas for making the 1988 data available, and Prof. Joseph Zhang for making the numerical model SCHISM available. The climate simulations for the BINGO dataset were performed at the German Climate Computing Centre (DKRZ). NCEP Reanalysis data was provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at http://www.esrl.noaa.gov/psd/. Two anonymous reviewers are also acknowledged for their valuable suggestions and comments.

---

Conflict of interest

The authors declare there is no conflict of interest.

---