Next Article in Journal
Understanding the Mechanical Biases of Tipping-Bucket Rain Gauges: A Semi-Analytical Calibration Approach
Previous Article in Journal
Occurrence of Pharmaceuticals and Personal Care Products in the Water Environment of Poland: A Review
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Coupling Rivers and Estuaries with an Ocean Model: An Improved Methodology

1
+Atlantic, Molhe Leste, 2520-620 Peniche, Portugal
2
MARETEC/LARSyS, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisbon, Portugal
*
Author to whom correspondence should be addressed.
Water 2021, 13(16), 2284; https://doi.org/10.3390/w13162284
Submission received: 7 July 2021 / Revised: 18 August 2021 / Accepted: 18 August 2021 / Published: 21 August 2021
(This article belongs to the Section Oceans and Coastal Zones)

Abstract

:
Freshwater sources are essential inputs for regional ocean models covering coastal areas such as the western Iberian Peninsula. The problem is how to include the mixture between fresh and salt water, typically performed by estuaries and in the adjacent areas of river mouths, without unsustainable increases of computational time and human setup errors. This work provides a proof-of-concept solution to both these problems through the use of an offline two-way methodology, where local schematic rivers and estuaries are responsible for mixing river freshwater with salt water of a regional model application. Two different offline upscaling methodologies—which focus on the implementation of tidal fluxes from local domains to regional domains in the context of operational modelling—are implemented in the Portuguese Coast Operational Modelling System (PCOMS) regional model application as well as in a version without rivers. A comparison between results produced by these methodologies, field data, and satellite imagery was performed, which confirmed that the proposed methodology of using schematic rivers and estuaries, combined with the new offline upscaling methodology proposed herein, represents a good solution for operational modelling of coastal areas subject to a high dominance of freshwater inputs.

1. Introduction

Coastal areas have been in the centre of human development due to the availability of water, fertile grounds, and abundance of food—marine and terrestrial. This is especially true for the Iberian coast, with its many rivers and estuaries serving as natural marine life maternities, and the seasonal upwelling phenomena on the western coast, which is responsible for the replenishment of nutrients and subsequent high marine biodiversity [1,2] and abundance of fish for local populations, whose socioeconomic development greatly depended on it [3]. There are also essential areas for carbon storage, with an important role in the fight against climate change, where an increase in human population will undoubtedly require a significant effort towards a sustainable management of coastal ecosystems resources [4]. In this context, coastal areas such as the Iberian Peninsula (IP) have been extensively studied through monitoring and modelling, but often separating between watershed and coastal area. This separation is done mainly at the level of the rivers that feed watersheds. One of the reasons is the difficulty of coupling watershed models with coastal models, which require not only specific programming and software to couple them but also multidisciplinary teams (or research projects) with members specialised in both systems. The study of this whole system, from watershed to coastal area—also known as “watershed-coast continuum” [5]—is especially important in areas with strong variabilities of the fluxes between the two systems [6,7], which is the case of the western Iberian Peninsula, with its tides, seasonal variations of river flows, and its coastal morphology.
Due to their importance as a source of energy, drinking water, for irrigation, and of nutrients and sediments to the estuaries and beaches, rivers were intensively monitored throughout the 20th century and the first decade of the 21st century, which generated a sizeable database on rivers’ flow and properties. However, economic constraints from the 2010s resulted in a reduction of maintenance operations of the hydrological stations, which led to a decline in the number of operational stations, following a global tendency of reduction in hydrometric networks [8]. This decline in the number of stations in Portugal has been steadily offset by the improvement of watershed modelling solutions, an example of which is the work of [9]. In their work, a watershed model application was implemented for the IP, whose results are extremely important as they provide flow and water properties of all the major river systems in the IP to the operational ocean modelling systems (and in forecast mode), with a more natural variability than the more standard method of climatology.
A move towards integrated modelling solutions considering both the watershed and ocean mediums can provide a better representation of the coastal processes due to its better representation of flood and drought events and their variability in time and space, which affect water mixing processes in density driven currents such as in Region Of Freshwater Influence (ROFI) areas [10,11]. In fact, performance of modelling forecasts of salinity in these areas has room to improve, as [12] concluded in their analysis of the Iberia-Biscay-Ireland (IBI) system, who linked it with uncertainties in the river flow data, especially during storm events [13,14].
An accurate representation of river and estuarine fluxes into the coastal area is paramount for modelling of biogeochemical, morphological, and ecological processes in coastal areas due to the amount of nutrients discharged through them. This is important when increased loads of land-based nutrients are increasing eutrophication problems [15,16], which can become worse when combined with fish farming activities [17,18]. As such, correctly estimating the nutrient loads into coastal areas will (1) improve the knowledge available for decision and policy makers at national and transnational levels [19], as currents will transport nutrients and organic matter across maritime borders, and (2) allow scientists to improve their quantifications of nutrient budgets, paths, and influence of river load variability on fish recruitment [2]. Furthermore, improving the coupling between watersheds and oceans is paramount for climate change studies focused on the variability of both the ocean (sea level and temperature rise) and land (lower precipitation and subsequent river flow rate), as well as for coastal morphology studies, which are highly dependent on the inputs from land and ocean. Neither climate change nor morphodynamics studies can be done without the use of ocean and watershed models (preferably coupled), but these are only as good as the inputs given by the modeller [20]. For climate change scenarios, the inputs to these models are what make the scenarios. In this case, scientists would benefit from better represented regional models taking as many freshwater sources as possible into account to create more realistic coastal circulation patterns, essential for studies focused in highly populated areas near river deltas or estuaries [21]. This also includes the study of the climate change impact in biological dynamics in estuaries such as the Great Lakes estuary [22]. As for morphodynamic studies, coastal currents are key actors in sediment transport dynamics, as they will determine the amount of sediment transported along shore and their direction. The other most important factor is the input of sediment from the watersheds. Typically these studies are focused on a particular lagoon, estuary, or river delta, such as in [23,24,25]; however, in some cases, two adjacent estuaries or river deltas can influence one another, as is the case of the Tagus and Sado estuaries in Portugal [26] or the northern Portuguese rivers [27]. Therefore, improving regional ocean models would help scientists to produce better quality information for decision makers acting on the local and regional level, through the study of management policy scenarios taking (or not) into consideration climate changes [20,28] and using regional model inputs that can better represent the combined effect of multiple freshwater sources in coastal currents.
Regarding the implementation of river discharge into an ocean model, [29] presents a good review of the different approaches, where the most common is a point-source volume discharge on the coastline with zero salinity. A more accurate approach is to also include the velocity through an artificial channel added to the coastline [30], as it provides momentum to the discharges and an initial mixing between fresh and salt waters. Naturally, the most accurate solution would be to properly represent the mixing area by including, for example, an estuary in the regional model. This would reproduce the mixing processes and hydrodynamics through proper bathymetric features and total volume of an estuary, as shown in MacCready et al. [31] and Liu et al. [32], who applied the method to the Colombia River. This is not evident in areas such as the Mediterranean Sea, whose very low tides are not the main driver for mixing between fresh and salt water and where a simple discharge of volume without prior mixing can be implemented with reasonable accuracy [33].Vertically, discharges can be either uniform along the water column or implemented in specific layers, the most common of which is a surface discharge due to the lower density of fresh water. However, Herzfeld [29] proposed a dynamic adjustment to obtain more realistic inputs in the coastal area, which tries to better estimate the depth of the model cell where the discharge is implemented and modify the flow profile and tracer properties accordingly.
Recently, the European component of the Global Ocean Observing System of the Intergovernmental Oceanographic Commission of UNESCO (EUROGOOS) group analysed the different methodologies for including river discharges into regional models and offered a recommendation [34] towards a more integrated watershed–ocean approach, as its improved temporal variability of river inputs led to better regional ocean modelling results. Campuzano et al. [10], who compared the use of river climatology with modelled data from a watershed modelling solution, also concluded that using watershed models improved the overall solution of the Portuguese Coast Operational Modelling System (PCOMS).
As typical grid resolutions of regional ocean models are greater than 1 km, a proper representation of all major rivers and estuaries in the western Iberian Peninsula would prove to be very complicated. Another issue is the collection of local bathymetric data for initial setup as well as their update over the years due to sediment dynamics. This study aimed to overcome this barrier by improving land–ocean interactions through the input of river flows into schematic rivers and estuaries.
Schematic rivers and estuaries are run separately from the regional model following an offline downscaling [35,36] followed by an offline upscaling that considers only a momentum discharge in the regional model, the combination of which is hereby defined as an offline two-way system. Two different upscaling approaches are tested in this article and applied to the PCOMS model application. The PCOMS operational version does not yet include freshwater discharges (hereafter no-rivers methodology) which compromises the validation of the model in the nearshore. A first approach to include these discharges was developed in [37], which intakes estuarine discharges from the operational estuarine model applications. However, this approach proved difficult to maintain due to the amount of time required to simulate the bigger estuaries—increased total time of a day’s simulation—and has not been reactivated since. The methodology presented in this work provides a solution to this problem by using schematic representations of the main freshwater sources in the western Iberian Peninsula—Minho, Douro, Mondego, Tagus, Sado, Guadiana, and Guadalquivir—which mix fresh and salt water before being discharged into PCOMS. The first upscaling approach follows the methodology proposed by Campuzano et al. [38] (hereafter detached methodology), and the second is an improved methodology proposed in this article (integrated methodology). Schematic rivers and estuaries have the advantage of a very simple representation of the river channels and estuaries—easy to update over time—but the new upscaling methodology allows for a simpler implementation of the freshwater discharges into the regional model, which only needs to receive an HDF or NetCDF from each separate estuary. With this methodology, the PCOMS system can now receive land boundary conditions from a watershed modelling system covering the entire western Iberian Peninsula, as well as from in situ and climatological data, and can generate an initial mixing with salt, therefore improving the representation of the haline fronts characteristic of this region. This work can have a strong impact on long-term studies, and more so in operational forecast modelling systems, due to the lower computational time, robustness, and ease of land sources implementation in regional model domains, which are important scientific issues faced by the modelling community in coastal areas.

2. Materials and Methods

2.1. Study Area

The study area extends from Cape São Vicente (37° N) to Cape Finisterre (43° N) (Figure 1), with a predominantly north–south coastal orientation. Circulation in the western Iberian coast is mainly driven by the North Atlantic Oscillation (NAO). It is caused by an atmospheric pressure gradient between a low-pressure system in Iceland and a high-pressure system in the Azores [39,40], and its seasonality and intensity changes with the relative position of these two systems.
Westerly winds dominate the western Iberian Peninsula during the winter, as the Azores High travels southeast and north winds prevail during the summer, when the Azores High moves northward (due to the clockwise rotation associated with a high-pressure system). North winds during summer generate an equatorward flow [41] and, when combined with the Coriolis force and the north–south orientation of this coastline, force surface waters offshore, upwelling deep, cold, and nutrient-rich waters to the surface near the coast [42,43]. During autumn and winter, the region is characterised by the Portuguese Coastal Counter Current [44] and the Iberian Poleward Current (IPC) [45,46], which flows along the continental slope from Cape Carvoeiro northwards during the winter and until April–May [47].
A large part of the Iberian rivers flow into the western coast of the Iberian Peninsula, which include the Tagus, Mondego, Aveiro, Douro, Lima, and Minho. During the winter period, freshwater discharges from these rivers (except the Tagus) produce a strong surface freshwater plume with an upper salinity limit below 35.7–35.8, denominated the West Iberia Buoyant Plume (WIBP). The WIBP causes the stratification of surface waters near the coast, along with strong haline fronts [45,48], and combined with the IPC generates a convergence area in the shelf break [2] that is essential for fish eggs and larvae [49]. Although the Tagus and Sado estuaries are excluded from the WIBP, they have been found to produce their own combined signature in very wet years, producing a buoyant plume denominated the West Iberian Central Plume (WICP) [37]. In these wet years, the size of both the WIBP and WICP is such that they form a single buoyant plume along the western coast of the Iberian Peninsula [37].

2.2. MOHID Water

The MOHID Water model is a finite-differences family model that evolved to a finite-volume model. Variables are placed in space using the Arakawa-C staggered grid. Velocity components are computed over the faces of the scalars control volume and so are the advective and diffusive fluxes. Water-volume fluxes across the faces of the control volume are computed using the last known velocities, and the concentration at the face is obtained by interpolation of the values known inside the finite-volumes using different alternative methods, including upwind, central differences, QUICK, and TVD. New values can be computed using explicit or semi-implicit algorithms. This approach ensures that the flux leaving a cell across a face is the flux entering in the cell sharing that face, guaranteeing that advection is conservative. Diffusivities are also computed over the finite-volume faces, assuring that diffusion is also conservative. Momentum fluxes are computed using the same rationale. The model solves the 3D incompressible primitive equations [50,51,52,53], built and developed using an object-oriented philosophy [54]. Its modular structure includes more than 40 modules, representing hydrodynamic, biogeochemical, wave, and sediment processes. Some of these processes and their correspondent descriptions include: biogeochemical model [55]; wave integration with currents [56]; bivalve modelling [57]; benthic marine systems modelling [58]; sediment dynamics [59]; and oil spill modelling [60].
The model assumes hydrostatic equilibrium and partially assumes the Boussinesq approximation. The density value used in the horizontal momentum equation is time and space dependent, but momentum fluxes are computed per unit of mass. Null velocity divergence is also assumed,
u i x i = 0
where u i are the three velocity components. Density is computed as a function of temperature, salinity, or suspended matter. Integrating the along the water column between the bottom h and a depth z , one gets the vertical velocity at that level.
u 3 x 3 = z = x 1 h z u 1 d x 3 x 2 h z u 2 d x 3 + h z q d x 3
where q are natural or anthropogenic local water sources. Integrating up to the free surface, one gets the free surface level η rate of change:
η t + x 1 h η u 1 d x 3 + x 2 h η u 2 d x 3 = h η q d x 3
In this equation, q includes rain and evaporation. Using the hydrostatic approximation, the vertical momentum equation becomes
p x 3 = ρ g
That integrated over the water column between the free surface and a depth z gives
p x 3 = z = p a t m + z η ρ g d x 3
if one performs the derivative along the horizontal axis, one gets the pressure gradient,
p x i x 3 = z = p a t m x i + g ρ x 3 = η η x i + g z η ρ x i d x 3
and the horizontal momentum equations can be written as:
ρ u 1 t + u 1 u j x j f u 2 = p a t m x 1 g ρ x 3 = η η x 1 g z η ρ x 1 d z + x j μ u 1 x j
ρ u 2 t + u 2 u j x j + f u 1 = p a t m x 2 g ρ x 3 = η η x 2 g z η ρ x 2 d z + x j μ u 2 x j
where f and μ represent the Coriolis parameter, a parallel line to the surface elevation, and the turbulent viscosity, respectively. These are the equations solved by the model version used in this work. The non-hydrostatic pressure component permitted by MOHID would require extra computing power and is not relevant for the purpose of this work. Temporal discretization is done using a semi-implicit scheme—alternating direction implicit (ADI) with two time levels per iteration. The numerical schemes implemented include the Abbot’s four equations system [61], where the water level is computed at Δt/2 and the velocities at Δt, and the Leendertse’s six equations system [62], where both the water level and the velocities are computed at Δt/2. Currently, and unless specifically indicated by a user, the model uses the latter.
Tracer properties such as temperature and salinity are calculated with the computed flow using
P t + P u i x i = x i μ P x i + S P   i = 1 , 2 , 3
where SP stands for sink-sources of the property (P) in question.
Horizontal and vertical advection are computed using the TVD advection scheme—a combination of a high order scheme (in the case of MOHID, a hybrid between first- and third-order upwind schemes) with a flux limiter—Superbee [63]—which is found to reduce spurious oscillations generated from shocks, discontinuities, or sharp changes in bathymetry by imposing monotonicity of the solution.

2.3. Offline Upscaling by Discharge

Offline upscaling in this work is characterized by the upscaling of a local (or child) model application (hereafter CD) into a regional (or parent) model application (hereafter PD) but where the two applications are run separately, like the offline downscaling [35] often applied in ocean modelling. In MOHID, the online upscaling was developed by Sobrinho et al. [64] following a similar methodology as in [65,66]. However, the online upscaling algorithm was made to upscale the entire overlapped area between two nested model domains, which means that for a PD to benefit from CDs simulating rivers or estuaries, it would need to upscale directly from full-scale bathymetric implementations of these coastal features. CDs would also need to cover a large enough area of the PD to improve the transfer of information and avoid issues at the open boundaries (for a more extensive review of the issues in online two-way couplings and upscaling, readers are referred to [67,68]). The problem is that online upscaling from large CDs increases the total simulation time of the PD, which is unsustainable for time-sensitive operational forecast modelling. This work aimed to provide a proof-of-concept that could solve this issue. To accomplish it, we propose the upscaling of only the information on the few numerical cells of a schematic version of the CD that simulates the exchange of water in tidal inlets, such as in estuaries. To further improve the methodology, an offline upscaling approach is proposed, in which horizontal velocities and tracer properties computed by the schematic CDs are used to produce a lateral momentum discharge in one or many numerical cells of the PD.
Offline upscaling through momentum discharges in MOHID was first developed in Campuzano et al. [38], where a full description of the flow and tracer properties computation is provided. The difference to the offline upscaling methodology proposed herein is procedural only, and the reader is referred to [38] for additional details. In both methodologies, the flow rate, velocity, and tracer properties are computed from the output of the schematic rivers and estuaries domains (Figure 2) through the definition of a cross-section perpendicular to the flow direction of the local domain in a suitable location, such as the mouth of the river or estuary. However, while in the detached methodology, the discharge properties are obtained through an external tool [38], where the user must define the location of the cross section, and one discharge file (time series containing flow, velocity, and tracer properties) is produced for each vertical layer of the estuary model application; in the integrated methodology the discharge properties are obtained directly in run-time by the MOHID model, where the user only needs to provide the HDF (or netCDF) outputs of the schematic rivers and estuaries along with the coordinates of the regional model numerical cell where the discharge is to be made (Figure 2—blue cell). In this case, MOHID will automatically detect the schematic river or estuary model cells that are required for the computation of the flow, velocity, and tracer properties and apply them to the user-specified regional model cell(s). Computation also takes into account the vertical structure by making a correspondence between vertical layers of both model applications. This means that different vertical discretizations can be used in the regional and local model applications, as the model will identify to which PD vertical layer each CD vertical layer belongs to. The flow is computed using the equation:
F = C D ( u i * · A k h )   i = 1 , 2
which is a simple summation of all local domain fluxes (ex: u C D * multiplied by the vertical area of its correspondent cell face—in the right-most green dashed line square of (Figure 2) crossing the cell faces of the regional domain cell identified by the user; in this case, it is imposed in the blue filled PD cell of Figure 2. F is computed for every vertical layer of the PD, creating one discharge point per vertical layer. In this equation, u i * represents the normal velocity of the cross section (in Figure 2 the dashed green line, where u C D * is computed. A k h stands for the vertical area of the CD cell associated with u C D * , where k is the vertical layer and h is the depth of the cell at the face. As these fluxes are obtained by output files of the CD, a linear temporal interpolation is performed using the nearest two time instants of the output file, to be imposed in the current time instant of the PD.
By automating this entire procedure, the implementation of several schematic rivers and estuaries becomes easier for the user and reduces the chance of human error in the implementation of multiple land discharges in a regional model application, such as the PCOMS application used in this article.

2.4. Grid Communication

Both methodologies tested in this article can be considered as offline two-way. However, in the implementation that follows the detached methodology, the CDs never receive the updated open boundary condition from the PD, while in the new integrated methodology they do. Adaptations to the detached methodology could be performed to enable a full two-way coupling, but it was not in the scope of this work. As such, in the detached methodology, a simulation of the entire period is performed by the PD, followed by the simulation of the same period by the CDs using the newly created PD results. Then, a second simulation (also of the entire period) is performed by the PD, but now with upscaling of the just made available CD results. This means that the grid coupling in the detached methodology is very limited (there is now a daily update of open boundary conditions from the PD or tidal fluxes from the CD). In the integrated methodology, results produced by the parent domain (PD) are downscaled into its nested child domains (CD) in an offline manner and for the entire simulation period. However, and in contrast with the detached methodology, the grid communication is performed on a daily basis. This means that after a one-day simulation of the PD without upscaling (because the CD results are typically behind in time), the CDs use the PD results as open boundary conditions to run the same day. After the CD simulation is finished, the PD runs a second time the same day, but it now upscales the CD results. After this second PD run, the updated PD results are again downscaled into the CDs (Figure 3).
In both cases, the information is transferred on a daily basis, where each model simulates one day at a time and then transfers the information. This procedure can be repeated any number of times, increasing the grid communication of the coupling. However, in this work, only one repetition is performed, which means upscaling is only performed one time.

2.5. Automatic Running Tool

In order to run all model applications (PCOMS and schematic rivers and estuaries), the automatic running tool (ART) [69,70] developed in MARETEC was improved. The tool is responsible for all file transfers (including meteorological model input files and open boundary conditions), invocation of external tools, and the MOHID model itself, as well as for storing all results. In the case of the detached methodology, the ART is responsible for calling a specific tool designed to extract the flow, velocity, and tracer properties from an HDF file produced by a schematic channel. This tool takes as user inputs an HDF file, a cross section that the user must identify, and outputs n time series containing date, flow, velocity, and tracer properties, where n is the number of vertical layers in the HDF. This is not necessary in the new integrated methodology proposed herein, as all these computations are done directly by MOHID, but the ART tool, in this case, is responsible for copying the HDF files produced by the schematic rivers and estuaries into the designated input folders of PCOMS implementation. The main improvement made to the ART tool refers to the launch of multiple runs per simulation (each day must be run more than once), combined with a trigger system which forces the model domains to wait for their parent (PCOMS) or son (schematic rivers or estuaries) domain (Figure 3) results.

2.6. Model Setup

In this study, the (3D)-MOHID water model was applied to the PCOMS system for the period between 1 October and 1 May 2018, using two different methodologies whose results were then compared between them. The PCOMS system is a nested grid configuration comprised of two domains: (i) 2D barotropic regional domain with 5.7 km constant grid resolution for the Portuguese Coast (WestIberia-Level1) (33.5° N–49.9° N, 1.0° W—13.5° W) running only with tidal forcing from FES2004 [71,72]; (ii) 3D full baroclinic regional domain for the Portuguese Coast (hereafter PCOMS) (34.4° N–45.0° N and 12.6° W–5.5° W) with a grid resolution of 5.7 km and 50 vertical layers (7 sigma at the surface and 43 Z coordinate layers below). A full description of the PCOMS is available in [36,37].
A new initialization was performed with a spin-up period of 3 months for both the PCOMS and the schematic rivers and estuaries. The PCOMS domain receives its lateral open boundary conditions from the 2D barotropic domain and the Mercator Ocean Psy2V4, linearly superimposed on the level barotropic velocities, as proposed in Leitão et al. [50]. Regarding PCOMS, a Flather radiation scheme [73] was applied at the open boundary to radiate water level, followed by a flow relaxation scheme (FRS) applied to the baroclinic velocities and tracer properties, as proposed by Martinsen and Engedahl [74], provided by the Mercator solution, as described in Leitão et al. [50]. The flow relaxation scheme applied decreases exponentially, a relaxation coefficient from 105 s in the outmost open boundary cell to 109 s in the 10th outmost cell (perpendicular to the open boundary). The rest of the domain was nudged to the Mercator solution with a coefficient of 109 s. Additionally, a biharmonic filter of 5.5 × 109 m4s−1 was applied to reduce high frequency noise inside the domain. As for the schematic rivers and estuaries, they receive their open boundary conditions from PCOMS HDF outputs holding 3D data every 900 s. The same radiation and relaxation schemes are applied, but due to the small number of numerical cells located inside the first water point of PCOMS, only three or four cells are using the FRS (depending on the channel). Their bathymetries were created taking into consideration the bathymetric values of the PCOMS domain near the discharge location. Thus, they have the same values in the overlapped area of water points. They were also based on the work of Campuzano et al. [38]. At the atmospheric boundary, both PCOMS and the schematic rivers and estuaries were forced by the meteorological model MM5 [75,76], with 9 km resolution for the Portuguese coast. Time steps used in the two levels of the PCOMS system were 120 s and 60 s, respectively, and the time step of the several channels was 30 s. A summary of the implementation is presented Table 1 and Table 2.
Discharge flow, temperature, and salinity were obtained from different sources according the river in question, which are described in Table 3. Sources include in situ data, climatological data, and a watershed model, which was implemented and validated in Campuzano et al. [10]. A main event can be observed in March (Table 3), which is visible in all rivers with the exception of Sado, whose watershed was not fed during this event. Most of these rivers are highly artificialized, which, combined with the absence of dams in the watershed models, can lead to unrealistic flow values. However, this is not in the scope of this work.

2.7. Validation Procedure

To validate the schematic rivers and estuaries methodology and the two methods of offline upscaling, a comparison was made against in situ sea surface temperature (SST) and salinity data for the periods between 1 January and 1 May and 18 March and 1 April 2018. Satellite SST data regarding the peak flow period (20 March) was used to validate SST fields produced by the integrated methodology and the no-rivers implementation of PCOMS. A separate salinity comparison was also made between the two offline upscaling methods (detached and integrated) against the no-rivers version of PCOMS in order to assess the advantages of including the freshwater inputs and to compare results from the two offline upscaling methods.
Regarding in situ data, the comparison was made against the moored buoys of Sillero and Estaca de Bares, located near the Galician rias and Cape Finisterre, respectively (see Figure 1), and for the period of 1 January to 1 May 2018. Statistical parameters of bias, Pearson correlation and root-mean-square error (RMSE) were computed and are shown in the Results section. Satellite images were retrieved from the ODYSSEA product, a 2 km gridded RS data product (http://marine.copernicus.eu/) [78] and were compared with the PCOMS solution under the integrated and the no-rivers methodologies. In order to obtain a better comparison between methodologies, satellite validation focused on the peak flow (20 March).

3. Results and Discussion

3.1. Comparison with In Situ Data

Temperature and salinity obtained by the three PCOMS implementations (without rivers, detached, and integrated offline upscaling) were compared against the Silleiro and Estaca de Bares moored stations during the period between 1 January and 1 May 2018 (Figure 4 and Figure 5).
As the moored station of Silleiro is closest to the last Portuguese river included in PCOMS (Minho), the riverine plume signal is easier to detect. Figure 4 shows the signal of the rivers especially during the period between 18 March and 1 April, as a result of the higher flow rates of the northern rivers (from Mondego upwards—Table 3). Furthermore, the separation between waters is best seen in salinity, as the gradient between ocean waters (≈36) and freshwater (≈0) is much higher than the temperature gradient (usually less than 10 °C). Regarding salinity, PCOMS without river inputs failed to reproduce the variability verified during the event, which is clearly dominated by land discharges, as it correlates well with results shown in Table 3. The PCOMS solution was improved when rivers were introduced in the detached methodology, which is sufficient to produce a signal similar to that observed by the in situ data. However, although the variability and timing of the events are well represented, the model produced an excess of salinity in the first event (18 March–1 April) and underestimated it in the second event (23 April–1 May). The solution improves again with the integrated methodology, which produces lower salinities during the first event and higher salinities during the second, following the in situ data. This improvement is mainly because in the new methodology the schematic rivers and estuaries run a second time, with the updated PCOMS solution, a step that was not present in the detached methodology. As such, the new integrated methodology can produce lower salinities simply because its open boundary condition values of salinity decrease with the river flow rates. During the second event, and due to the stronger influence of the northern river plumes during the entire simulation period (the schematic rivers and estuaries received feedback from the updated PCOMS model, lowering the open boundary salinity values), it led to a more intense combined plume being obtained by the integrated methodology, which remains closer to shore in the particular latitude of the moored buoy. With the exception of the no-rivers PCOMS methodology, which fails to represent the variability measured by the moored buoy, temperature results are more consistent between methodologies. In the two previously identified events, a shift in temperature can be observed when both offline upscaling methodologies are compared with field data and the no-rivers PCOMS methodology. While in the first, the freshwater inputs are responsible for a decrease in temperature, in the second—due to the increase of freshwater temperature during the spring—the river inputs contributed to the increase of temperature in the coastal area.
The statistical analyses (Table 4 and Table 5) focused on the global simulation period (1 January and 1 May 2018) and on the main peak flow period (18 March–1 April). Results produced by the new integrated methodology show smaller average bias than the detached methodology for temperature regarding the entire simulation period and salinity. RMSD values improved in the integrated methodology for all periods and parameters, with the exception of temperature during the peak flow. Nevertheless, differences between bias and RMSD values obtained by these two methodologies are quite small and indicate a similar accuracy. However, the correlation coefficients obtained for the different periods and parameters improved substantially with the use of the integrated methodology (r2 always above 0.7). As expected, the no-rivers PCOMS methodology provided the lowest accuracies and correlations with the exception of the temperature bias and RMSD for the peak flow period, which is associated with low variability of the time series, combined with the casual similarity of temperature during this particular time period.
Further north, near Cape Finisterre, the freshwater influence from the most northern Portuguese rivers is hardly felt (Figure 5), not only due to the distance (which leads to a higher diffusion of the freshwater plumes) but also due to the absence of the northern Spanish river inputs which were not included in this work. However, the influence of these freshwater discharges is noticed in the in situ data which show a salinity drop at the same time as the temperature increases sharply. As such, for an accurate representation of density, currents near the north-western coast of the IP the main rivers outflowing in the Spanish coast are required. An effort is already being made towards solving this issue with the recently finished Lambda project, where a data base on river flow inputs has been compiled using watershed models to cope with the lack of hydrographic stations [79].
The statistical analysis of the results obtained by the different methodologies for the two selected periods confirms what was observed in Figure 5, where the absence of the Spanish rivers reduces the reach of the freshwater plume along the north-western coast of the IP.

3.2. Surface Salinity Maps

A comparison between salinity maps (Figure 6) produced by the three methodologies is important for a broader analysis of their impact on the evolution of river plumes in the western coast of the IP. As such, time-average surface maps of salinity were computed from the model outputs for the period of 18 March to 1 April, which corresponds to the strongest precipitation event, as observed in Table 3. Results show the obvious impact of adding river inputs in the PCOMS system, which were responsible for the formation of a buoyant plume (WIBP) produced by the three Portuguese rivers north of the Tagus mouth (Mondego, Douro, and Minho). Figure 7 (bottom-right) shows this impact, which is most visible in the coastal areas adjacent to the river mouths, where salinity differences can reach 10 salinity units (although the legend starts at −4) near the river and estuaries mouths. It also shows the importance of including freshwater inputs in regional model applications, which produce a significant change in the local currents, essential for the transport of nutrients, plankton, and sediments required for biological and geomorphological studies alike.
Very small differences were observed between results obtained by the two offline upscaling methodologies (detached and integrated—Figure 7). The first is in the vicinity of the Tagus estuary, where the new methodology produced higher salinity values (1 to 1.5) instead of the expected lower values (due to the higher level of grid communication). The explanation for this difference is more complex and is related to the schematic domains volume available, combined with the difference in the methodology. Because in the detached methodology the rivers and estuaries open boundary conditions were always fed by the no-rivers PCOMS version, their open boundary values of salinity are constantly high (above 35). When combined with the lower volume of the channel in comparison with the current operational version of the Tagus estuary model domain (described and validated in de Pablo et al. [80]), lower salinity values will be produced near the estuary mouth (from where the fluxes are retrieved and added to PCOMS). As such, there will be less vertical mixing between ocean and river water, which leads to a higher surface gradient of salinity when high river flow rates are present. In the northern rivers, this effect does not take place, most likely due to a better volume correlation between the schematic rivers and estuaries and reality and due to their higher flow rate in comparison with their tidal prism. Therefore, in the northwest coast of the IP, the integrated methodology produced slightly lower values (around 0.5) of surface salinity—as expected when the schematic rivers and estuaries receive updated boundary conditions from PCOMS. In the adjacent area of the Guadiana and Guadalquivir mouths (Figure 6 and Figure 7), the impact of the two offline upscaling methodologies is quite similar, which can be caused by the local hydrodynamics—through higher mixing in the area—which can increase diffusion and therefore reduce the differences between methodologies. A zoom of the area of interest, which covers a part of the west coast of the IP, is shown in both Figure 6 and Figure 7 for a better visualization of this analysis in the ROFI areas.

3.3. SATELLITE SST Analysis

An analysis of surface temperature map was done near the peak period (20 March) against satellite data with an aim not only to obtain validation of the schematic rivers and estuaries methodology and the integrated offline upscaling but also to obtain a broader view of the surface temperature patterns over the coastal area of the western coast of the IP when river inputs are considered. The comparison (Figure 8) was made only for the no-rivers and the integrated methodologies for the period of peak flow (18 March–1 April) due to the high similarity between the latter and the detached methodology. At first glance, the results (Figure 8) show the typical Northern Hemisphere north–south temperature gradient, with increasing temperatures towards the equator in both methodologies. However, major differences were obtained in the coastal area of the IP and near its ROFI areas (Figure 9), where the influence of the freshwater inputs (present only in the PCOMS version where the integrated offline upscaling methodology was applied). River inputs played a key role in surface temperature near the coast and up to 80 km offshore in the latitude of the Douro and Minho Rivers mouths, which have the highest flow rates (Table 3). A zoom of the coastal area of the IP is provided in both Figure 8 and Figure 9 for an easier analysis of these results.
In fact, the comparison with satellite data demonstrates the accuracy of the model and the offline upscaling methodology in representing the north–south surface temperature gradient. However, a more detailed analysis of the SST pattern (Figure 8 and Figure 9) obtained by model and satellite data for the coastal areas suggests that the no-rivers methodology compared better against satellite than the integrated methodology and suggests an overestimation of SST by the satellite, as confirmed by the comparison against the Silleiro moored buoy (Figure 4). In this location, the in situ data show a value of around 12.2 °C, while the satellite suggests a value above 13 °C. Although this difference is small, Figure 4 shows a decrease in temperature of 1.3 °C associated with the plumes around 18 March, whereas the satellite observed a constant temperature during these days for that location. In these coastal areas, a possible issue is related to the accuracy of satellite SST near the coastal area. This was observed by Brewin et al. [81], where in situ SST collected from autonomous buoys showed a bias on the order of 0.45–0.51 °C with satellite SST values, which was addressed in the validation procedure of [10]. Another reason is the solar radiation effect on surface temperature. When a discharge of fresh water takes place and its temperature is lower than the ocean, the solar radiation will warm the surface water, mainly in the first millimetres during the day, and the signal will soon disappear from the satellite (but not from in situ stations measuring below this depth). As satellite images are taken only once per day, punctual river signals can be lost from one day to the next. Differences in the entire domain can be explained by the fact that satellite SST measures the temperature at the “skin” of the ocean, which represents approximately the top 0.01 mm, while the model SST is the average temperature in the top cell of the model domain, 0.95 m plus tidal range. This is also a possible explanation for the detection, by the satellite, of the Douro and Minho plumes traveling along the coastline to the Galician rias in the north, while the model results suggest this influence is felt further away from the coast.

4. Conclusions

A good model representation of the coastal hydrodynamic processes in the IP—namely, those associated with freshwater discharges, such as the WIBP and the WICP—requires the input of all major river sources. These processes are important to subsequent ecological and morphodynamic processes which will then affect all human coastal activities and respective coastal management policies, especially relevant for climate change studies focusing on the impact of extreme precipitation events, sea level rise, and storm surges. However, to properly add freshwater discharges in regional ocean models, horizontal and vertical mixing between freshwater and ocean water due to tidal prism is necessary so as to correctly compute the haline fronts. In the case of the western IP, the Portuguese Coast Operational Modelling System (PCOMS)—a regional mode application—has been running since 2009 without freshwater sources (no-rivers methodology), which has partially compromised its results with regard to surface temperature and salinity near the coast. A first approach to solve this issue was presented in Campuzano et al. [38], who added freshwater discharges through the coupling of the ocean model with real case estuaries and rivers (detached methodology). However, this improvement of the regional model came at a great cost in computational time, as results produced by PCOMS now needed to wait for the slower estuarine model application. Another complication was related to the implementation of the several discharges files (one per river per vertical layer) by the user, as well as the need to run an external tool also implemented by a user. All these steps increase the chance for human error, making the entire system more demanding in terms of maintenance. A new solution is proposed in this work, whereby the heavier, real river and estuary applications are replaced by simpler and less time-consuming schematic versions of these applications. As such, the total amount of time required for PCOMS to simulate one day increased by only around two minutes, making this methodology much more attractive for operational modelling purposes. Results obtained by the three aforementioned methodologies were compared against in situ sea surface salinity and temperature data and against satellite observations of sea surface temperature. They showed an improvement from the no-rivers methodology to the detached and integrated methodologies, particularly in the coastal area of the IP. It was also found that the satellite had some difficulty representing the surface temperature signal from the northern Portuguese rivers in comparison with in situ and model data, reinforcing the need for land discharges in coastal modelling simulations, which can fill in the gaps of satellite data. Additionally, the two methodologies where river inputs were added to PCOMS produced very similar results, suggesting that this new integrated methodology can be used for operational modelling of coastal areas where freshwater sources play an important role. Nevertheless, this work consists of a proof-of-concept in which the estuaries volumes have not been adjusted to their real volume, which means that the tidal prisms are not properly represented. This is another research field and should be addressed in future work considering this methodology. Longer simulation periods should also be performed in order to study the impact of the methodology over time, especially for use in climate change and morphodynamic studies which typically need to simulate several years.

Author Contributions

Conceptualization, J.S., F.C. and R.N.; methodology J.S. and F.C.; software, J.S. and F.C.; validation, J.S. and H.d.P.; formal analysis, J.S., H.d.P. and F.C.; investigation, J.S., H.d.P. and F.C.; writing—original draft preparation, J.S. and H.d.P.; writing—review and editing, J.S., H.d.P., F.C. and R.N.; visualization, J.S. and H.d.P.; project administration, R.N. All authors have read and agreed to the published version of the manuscript.

Funding

The present work was supported by FCT/MCTES (PIDDAC) through the project LARSyS—FCT Pluriannual funding 2020–2023 (UIDB/50009/2020).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data regarding surface temperature and salinity from in situ data are available in https://www.emodnet-physics.eu/map/platinfo/piroosplot.aspx?platformid=7340—Silleiro moored buoy (accessed on 1 March 2021) and in https://www.emodnet-physics.eu/map/platinfo/PIROOSDownload.aspx?PlatformID=7339—Estaca de Bares (accessed on 1 March 2021). The ODYSSEA satellite data were obtained from http://marine.copernicus.eu/ (accessed on February 2019).

Acknowledgments

The authors would like to acknowledge the E.U. Copernicus Marine Service Information and Satellite Monitoring of the Northwest Shelf Sea for providing the ODYSSEA SST product. The authors would also like to acknowledge Manuela Juliano from the University of Azores for her invaluable help in the creation of the validation scripts.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Santos, A.M.P.; Chícharo, A.; Dos Santos, A.; Moita, T.; Oliveira, P.B.; Peliz, Á.; Ré, P. Physical-biological interactions in the life history of small pelagic fish in the Western Iberia Upwelling Ecosystem. Prog. Oceanogr. 2007, 74, 192–209. [Google Scholar] [CrossRef]
  2. Santos, A.M.P.; Peliz, A.; Dubert, J.; Oliveira, P.B.; Angélico, M.M.; Ré, P. Impact of a winter upwelling event on the distribution and transport of sardine (Sardina pilchardus) eggs and larvae off western Iberia: A retention mechanism. Cont. Shelf Res. 2004, 24, 149–165. [Google Scholar] [CrossRef]
  3. Canuel, E.A.; Cammer, S.S.; McIntosh, H.A.; Pondell, C.R. Climate change impacts on the organic carbon cycle at the land-ocean interface. Annu. Rev. Earth Planet. Sci. 2012, 40, 685–711. [Google Scholar] [CrossRef]
  4. Ward, N.D.; Megonigal, J.P.; Bond-Lamberty, B.; Bailey, V.L.; Butman, D.; Canuel, E.A.; Diefenderfer, H.; Ganju, N.K.; Goñi, M.A.; Graham, E.B.; et al. Representing the function and sensitivity of coastal interfaces in Earth system models. Nat. Commun. 2020, 11, 1–14. [Google Scholar] [CrossRef]
  5. Salomons, W. Sediments in the catchment-coast continuum. J. Soils Sediments 2005, 5, 2–8. [Google Scholar] [CrossRef]
  6. Malara, G.; Zema, D.A.; Arena, F.; Bombino, G.; Zimbone, S.M. Coupling watershed—coast systems to study evolutionary trends: A review. Earth-Sci. Rev. 2020, 201, 103040. [Google Scholar] [CrossRef]
  7. Chen, Y.; Cheng, W.; Zhang, H.; Qiao, J.; Liu, J.; Shi, Z.; Gong, W. Evaluation of the total maximum allocated load of dissolved inorganic nitrogen using a watershed–coastal ocean coupled model. Sci. Total Environ. 2019, 673, 734–749. [Google Scholar] [CrossRef]
  8. Mishra, A.K.; Coulibaly, P. Developments in hydrometric network design: A review. Rev. Geophys. 2009, 1–24. [Google Scholar] [CrossRef]
  9. Brito, D.; Campuzano, F.J.; Sobrinho, J.; Fernandes, R.; Neves, R. Integrating operational watershed and coastal models for the Iberian Coast: Watershed model implementation—A first approach. Estuar. Coast. Shelf Sci. 2015, 167, 138–146. [Google Scholar] [CrossRef]
  10. Campuzano, F.; Brito, D.; Juliano, M.; Fernandes, R.; de Pablo, H.; Neves, R. Coupling watersheds, estuaries and regional ocean through numerical modelling for Western Iberia: A novel methodology. Ocean. Dyn. 2016, 66, 1745–1756. [Google Scholar] [CrossRef]
  11. Sotillo, M.G.; Campuzano, F.; Guihou, K.; Lorente, P.; Olmedo, E.; Matulka, A.; Santos, F.; Amo-Baladrón, M.A.; Novellino, A. River freshwater contribution in operational ocean models along the european atlantic façade: Impact of a new river discharge forcing data on the cmems ibi regional model solution. J. Mar. Sci. Eng. 2021, 9, 401. [Google Scholar] [CrossRef]
  12. Aznar, R.; Sotillo, M.G.; Cailleau, S.; Lorente, P.; Levier, B.; Amo-Baladrón, A.; Reffray, G.; Álvarez-Fanjul, E. Strengths and weaknesses of the CMEMS forecasted and reanalyzed solutions for the Iberia-Biscay-Ireland (IBI) waters. J. Mar. Syst. 2016, 159, 1–14. [Google Scholar] [CrossRef]
  13. Ruiz-Parrado, I.; Genua-Olmedo, A.; Reyes, E.; Mourre, B.; Rotllán, P.; Lorente, P.; García-Sotillo, M.; Tintoré, J. Coastal ocean variability related to the most extreme Ebro River discharge over the last 15 years. Section in Copernicus Marine Service Ocean State Report, Issue 4. J. Oper. Ocean. 2020, 13 (Suppl. 1), s160–s165. [Google Scholar]
  14. Sotillo, M.G.; Mourre, B.; Mestres, M.; Lorente, P.; Aznar, R.; García-León, M.; Liste, M.; Santana, A.; Espino, M.; Álvarez, E. Evaluation of the Operational CMEMS and Coastal Downstream Ocean Forecasting Services During the Storm Gloria (January 2020). Front. Mar. Sci. 2021, 8, 1–27. [Google Scholar] [CrossRef]
  15. Andersen, J.H.; Carstensen, J.; Conley, D.J.; Dromph, K.; Fleming-lehtinen, V.; Gustafsson, B.G.; Josefson, A.B.; Norkko, A.; Villn, A. Long-term temporal and spatial trends in eutrophication status of the Baltic Sea. Biol. Rev. 2017, 92, 135–149. [Google Scholar] [CrossRef] [Green Version]
  16. Le Moal, M.; Gascuel-Odoux, C.; Ménesguen, A.; Souchon, Y.; Étrillard, C.; Levain, A.; Moatar, F.; Pannard, A.; Souchu, P.; Lefebvre, A.; et al. Eutrophication: A new wine in an old bottle? Sci. Total Environ. 2019, 651, 1–11. [Google Scholar] [CrossRef] [Green Version]
  17. Smith, V.H.; Schindler, D.W. Eutrophication science: Where do we go from here? Trends Ecol. Evol. 2009, 24, 201–207. [Google Scholar] [CrossRef]
  18. Du, J.; Shen, J.; Park, K.; Wang, Y.P.; Yu, X. Worsened physical condition due to climate change contributes to the increasing hypoxia in Chesapeake Bay. Sci. Total Environ. 2018, 630, 707–717. [Google Scholar] [CrossRef]
  19. Desmit, X.; Thieu, V.; Billen, G.; Campuzano, F.; Dulière, V.; Garnier, J.; Lassaletta, L.; Ménesguen, A.; Neves, R.; Pinto, L.; et al. Reducing marine eutrophication may require a paradigmatic change. Sci. Total Environ. 2018, 635, 1444–1466. [Google Scholar] [CrossRef]
  20. Samaras, A.G.; Koutitas, C.G. An integrated approach to quantify the impact of watershed management on coastal morphology. Ocean. Coast. Manag. 2012, 69, 68–77. [Google Scholar] [CrossRef]
  21. Duong, T.M.; Ranasinghe, R.; Walstra, D.; Roelvink, D. Assessing climate change impacts on the stability of small tidal inlet systems: Why and how? Earth-Sci. Rev. 2016, 154, 369–380. [Google Scholar] [CrossRef]
  22. Liu, Q.; Anderson, E.J.; Zhang, Y.; Weinke, A.D.; Knapp, K.L.; Biddanda, B.A. Modeling reveals the role of coastal upwelling and hydrologic inputs on biologically distinct water exchanges in a Great Lakes estuary. Estuar. Coast. Shelf Sci. 2018, 209, 41–55. [Google Scholar] [CrossRef]
  23. Yin, Y.; Karunarathna, H.; Reeve, D.E. Numerical modelling of hydrodynamic and morphodynamic response of a meso-tidal estuary inlet to the impacts of global climate variabilities. Mar. Geol. 2019, 407, 229–247. [Google Scholar] [CrossRef] [Green Version]
  24. Franz, G.; Delpey, T.M.; Brito, D.; Neves, R.; Leitão, P.; Pinto, L. Modelling of sediment transport and morphological evolution under the combined action of waves and currents. Ocean. Sci. 2017, 13, 673–690. [Google Scholar] [CrossRef] [Green Version]
  25. Bruneau, N.; Fortunato, A.B.; Dodet, G.; Freire, P.; Oliveira, A.; Bertin, X. Future evolution of a tidal inlet due to changes in wave climate, Sea level and lagoon morphology (Óbidos lagoon, Portugal). Cont. Shelf Res. 2011, 31, 1915–1930. [Google Scholar] [CrossRef]
  26. Ribeiro, A.S.; Sousa, M.C.; Lencart, E.; Silva, J.D.; Dias, J.M. David and goliath revisited: Joint modelling of the Tagus and Sado Estuaries. J. Coast. Res. 2016, 1, 123–127. [Google Scholar] [CrossRef]
  27. Oliveira, V.H.; Sousa, M.C.; Picado, A.; Mendes, R.; Ribeiro, A.S.; Morgado, F.; Dias, J.M. Coupled modelling of the interaction between dissolved substances emitted by Minho and Lima estuarine outflows (Portugal). J. Mar. Syst. 2021, 222, 103601. [Google Scholar] [CrossRef]
  28. Samaras, A.G.; Koutitas, C.G. Comparison of three longshore sediment transport rate formulae in shoreline evolution modeling near stream mouths. Ocean. Eng. 2014, 92, 255–266. [Google Scholar] [CrossRef]
  29. Herzfeld, M. Methods for freshwater riverine input into regional ocean models. Ocean. Model. 2015, 90, 1–15. [Google Scholar] [CrossRef]
  30. Lacroix, G.; Ruddick, K.; Ozer, J.; Lancelot, C. Modelling the impact of the Scheldt and Rhine/Meuse plumes on the salinity distribution in Belgian waters (southern North Sea). J. Sea Res. 2004, 52, 149–163. [Google Scholar] [CrossRef]
  31. MacCready, P.; Banas, N.S.; Hickey, B.M.; Dever, E.P.; Liu, Y. A model study of tide- and wind-induced mixing in the Columbia River Estuary and plume. Cont. Shelf Res. 2009, 29, 278–291. [Google Scholar] [CrossRef]
  32. Liu, Y.; MacCready, P.; Hickey, B.M.; Dever, E.P.; Kosro, P.M.; Banas, N.S. Evaluation of a coastal ocean circulation model for the Columbia river plume in summer 2004. J. Geophys. Res. Ocean. 2009, 114, 1–23. [Google Scholar] [CrossRef] [Green Version]
  33. Estournel, C.; Kondrachoff, V.; Marsaleix, P.; Vehil, R. The plume of the Rhone: Numerical simulation and remote sensing. Cont. Shelf Res. 1997, 17, 899–924. [Google Scholar] [CrossRef]
  34. Capet, A.; Fernández, V.; She, J.; Dabrowski, T.; Umgiesser, G.; Staneva, J.; Mészáros, L.; Campuzano, F.; Ursella, L.; Nolan, G.; et al. Operational Modeling Capacity in European Seas—An EuroGOOS Perspective and Recommendations for Improvement. Front. Mar. Sci. 2020, 7, 1–19. [Google Scholar] [CrossRef] [Green Version]
  35. Mateus, M.; Riflet, G.; Chambel, P.; Fernandes, L.; Fernandes, R.; Juliano, M.; Campuzano, F.; de Pablo, H.; Neves, R. An operational model for the West Iberian coast: Products and services. Ocean. Sci. 2012, 8, 713–732. [Google Scholar] [CrossRef] [Green Version]
  36. Campuzano, F.J.; Fernandes, R.; Leitão, P.C.; Viegas, C.; de Pablo, H.; Neves, R. Implementing local operational models based on an offline downscaling technique: The Tagus estuary case. In Proceedings of the 2nd Jornadas Engenharia Hidrográfica, Lisbon, Portugal, 20–22 June 2012; pp. 105–108. [Google Scholar]
  37. Campuzano, F. Coupling Watersheds, Estuaries and Regional Seas through Numerical Modelling for Western Iberia. Ph.D. Thesis, University of Lisbon, Lisbon, Portugal, 2018. [Google Scholar]
  38. Campuzano, F.J.; Juliano, M.; Sobrinho, J.; de Pablo, H.; Brito, D.; Fernandes, R.; Neves, R. Coupling Watersheds, Estuaries and Regional Oceanography through Numerical Modelling in the Western Iberia: Thermohaline Flux Variability at the Ocean-Estuary Interface. Estuary 2018. [Google Scholar] [CrossRef] [Green Version]
  39. Ottersen, G.; Planque, B.; Belgrano, A.; Post, E.; Reid, P.C.; Stenseth, N.C. Ecological effects of the North Atlantic Oscillation. Oecologia 2001, 128, 1–14. [Google Scholar] [CrossRef]
  40. Hurrell, J.W.; Dickson, R.R. Climate Variability Over the North Atlantic. In Marine Ecosystems and Climate Variation; Oxford University Press: Oxford, UK, 2005; ISBN 9780198507499. [Google Scholar]
  41. Fiúza, A.F. de G.; de Macedo, M.E.; Guerreiro, M.R. Climatological space and time variation of the Portuguese coastal upwelling Upwelling Portugal Coastal winds Sea surface temperature Sardine Upwelling Portugal Vents côtiers Température de surface de la mer Sardine. Oceanol. Acta 1982, 5, 31–40. [Google Scholar]
  42. Fiuza, A.F.G. Upwelling Patterns Off Portugal. NATO Conf. Ser. 4 Mar. Sci. 1983, 10 A, 85–98. [Google Scholar] [CrossRef]
  43. Ruiz-Villarreal, M. Hydrodynamic Model Study of the Ria de Pontevedra Under Estuarine Conditions. Estuar. Coast. Shelf Sci. 2002, 54, 101–113. [Google Scholar] [CrossRef]
  44. Ambar, I.; Fiúza, A.F.G. Some features of the Portugal current system: A poleward slope undercurrent, an upwelling-related summer southward flow and an autumn-winter poleward coastal surface current. Proc. Second Int. Conf. Air-Sea Interact. Meteorol. Oceanogr. Coast. Zo. 1994, 286–287. [Google Scholar]
  45. Peliz, Á.; Dubert, J.; Santos, A.M.P.; Oliveira, P.B.; Le Cann, B. Winter upper ocean circulation in the Western Iberian Basin—Fronts, Eddies and Poleward Flows: An overview. Deep. Res. Part. I Oceanogr. Res. Pap. 2005, 52, 621–646. [Google Scholar] [CrossRef]
  46. Frouin, R.; Fiúza, A.F.G.; Ambar, I.; Boyd, T.J. Observations of a poleward surface current off the coasts of Portugal and Spain during winter. J. Geophys. Res. Ocean. 1990, 95, 679–691. [Google Scholar] [CrossRef]
  47. Relvas, P.; Barton, E.D.; Dubert, J.; Oliveira, P.B.; Peliz, Á.; da Silva, J.C.B.; Santos, A.M.P. Physical oceanography of the western Iberia ecosystem: Latest views and challenges. Prog. Oceanogr. 2007, 74, 149–173. [Google Scholar] [CrossRef] [Green Version]
  48. Peliz, Á.; Rosa, T.L.; Santos, A.M.P.; Pissarra, J.L. Fronts, jets, and counter-flows in the Western Iberian upwelling system. J. Mar. Syst. 2002, 35, 61–77. [Google Scholar] [CrossRef]
  49. Santos, A.M.P.; Re, P.; Dos Santos, A.; Peliz, Á. Vertical distribution of the European sardine (Sardina pilchardus) larvae and its implications for their survival. J. Plankton Res. 2006, 28, 523–532. [Google Scholar] [CrossRef] [Green Version]
  50. Leitão, P.; Coelho, H.; Santos, A.; Neves, R. Modelling the main features of the Algarve coastal circulation during July 2004: A downscaling approach. J. Atmos. Ocean. Sci. 2005, 10, 421–462. [Google Scholar] [CrossRef]
  51. Leitão, P. Integracção de Escalas e de Processos na Modelação do Ambiente Marinho; Universidade Técnica de Lisboa, Instituto Superior Técnico: Lisbon, Portugal, 2002. [Google Scholar]
  52. Martins, F.; Leitão, P.C.; Silva, A.; Neves, R. 3D modelling in the Sado estuary using a new generic vertical discretization approach. Oceanol. Acta 2001, 24, 51–62. [Google Scholar] [CrossRef] [Green Version]
  53. Martins, F.A.; Neves, R.J.; Leitão, P.C. A three-dimensional hydrodynamic model with generic vertical coordinate. Hydroinformatics 1998, 1403–1410. [Google Scholar]
  54. Braunschweig, F.; Leitao, P.C.; Fernandes, L.; Pina, P.; Neves, R.J.J. The object-oriented design of the integrated water modelling system MOHID. Dev. Water Sci. 2004, 55, 1079–1090. [Google Scholar] [CrossRef]
  55. Mateus, M. A Process-oriented Biogeochemical Model for Marine Ecosystems: Development, Numerical Study, and Application; Universidade Técnica de Lisboa—Instituto Superior Técnico: Lisbon, Portugal, 2006. [Google Scholar]
  56. Delpey, M. Etude de la Dispersion Horizontale en Zone Littorale sous l’Effet de la Circulation Tridimensionnelle Forcée par les Vagues—Application à la Baie de Saint Jean de Luz—Ciboure et au Littoral de Guéthary-Bidart: HAL Id. Ph.D. Thesis, Université de Bretagne Occidentale-Brest, Brest, France, 2013. [Google Scholar]
  57. Saraiva, S. Modelling Bivalves in Estuaries and Coastal Areas. Ph.D. Thesis, University of Lisbon, Lisbon, Portugal, 2014. [Google Scholar]
  58. Ascione, I. Development and Application of a Process-oriented Model for Benthic Marine Systems. Ph.D. Thesis, Lisbon University, Lisbon, Portugal, 2014; 214p. [Google Scholar]
  59. Franz, G. Numerical Modelling of Hydrodynamics and Sediment Transport in Coastal Systems. Ph.D. Thesis, Instituto Superior Técnico, Lisbon, Portugal, 2017; pp. 2–3. [Google Scholar] [CrossRef]
  60. Fernandes, R. Risk Management of Coastal Pollution from Oil Spills Supported by Operational Numerical Modelling. Ph.D. Thesis, University of Lisbon, Lisbon, Portugal, 2018. [Google Scholar]
  61. Abbott, M.B.; Damsgaard, A.; Rodenhuis, G.S. System 21, “Jupiter” (A design system for two-dimensional nearly-horizonal flows). J. Hydraul. Res. 1973, 11, 1–28. [Google Scholar] [CrossRef]
  62. Leendertse, J. Aspects of a Computational Model for Long Water Wave Propagation; The RAND Corporation: Santa Monica, CA, USA, 1967. [Google Scholar]
  63. Roe, P.L. Characteristic-Based Schemes for the Euler Equations. Annu. Rev. Fluid Mech. 2003, 18, 337–365. [Google Scholar] [CrossRef]
  64. Sobrinho, J.; de Pablo, H.; Pinto, L.; Neves, R. Improving 3D-MOHID water model with an upscaling algorithm. Environ. Model. Softw. 2021, 135, 104920. [Google Scholar] [CrossRef]
  65. Oey, L.Y.; Chen, P. A nested-grid ocean model: With application to the simulation of meanders and eddies in the Norwegian Coastal Current. J. Geophys. Res. 1992, 97, 20063–20086. [Google Scholar] [CrossRef]
  66. Spall, M.; Holland, W.R. A Nested Primitive Equation Model for Oceanic Applications. J. Phys. Oceanogr. 1991, 21, 205–220. [Google Scholar] [CrossRef] [Green Version]
  67. Debreu, L.; Blayo, E. Two-way embedding algorithms: A review. Ocean. Dyn. 2008, 58, 415–428. [Google Scholar] [CrossRef]
  68. Debreu, L.; Marchesiello, P.; Penven, P.; Cambon, G. Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. Ocean. Model. 2012, 49–50, 1–21. [Google Scholar] [CrossRef]
  69. Kenov, I.A.; Campuzano, F.; Franz, G. Advances in Modeling of Water Quality in Estuaries Chapter 10 Advances in Modeling of Water Quality in Estuaries. In Remote Sensing and Modeling; Springer: Cham, Switzerland, 2014; ISBN 9783319063263. [Google Scholar]
  70. Franz, G.; Femandes, R.; de Pablo, H.J.; Viegas, C.; Pinto, L.; Campuzano, F.; Ascione, I.; Leitão, P.; Neves, R. Tagus Estuary Hydro-biogeochemical Model: Inter-Annual Validation and Operational Model Update. 3.as Jornadas Engenharia Hidrográfica; Extended Abstracts: Lisbon, Portugal, 2014. [Google Scholar]
  71. Lyard, F.; Lefevre, F.; Letellier, T.; Francis, O. Modelling the global ocean tides: Modern insights from FES2004. Ocean. Dyn. 2006, 56, 394–415. [Google Scholar] [CrossRef]
  72. Lefèvre, F.; Lyard, F.H.; Le Provost, C.; Schrama, E.J.O. FES99: A global tide finite element solution assimilating tide gauge and altimetric information. J. Atmos. Ocean. Technol. 2002, 19, 1345–1356. [Google Scholar] [CrossRef]
  73. Flather, R.A. A tidal model of the northwest European continental shelf. Mem. Soc. R. Sei. Liège 1976, 10, 141–164. [Google Scholar]
  74. Martinsen, A.E.; Engedahl, H. Implementation and testing of a lateral boundary scheme as an open boundary condition in a barotropic ocean model. Coast. Eng. 1987, 11, 603–627. [Google Scholar] [CrossRef]
  75. Grell, G.A.; Dudhia, J.; Stauffer, S.R. A Description of the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5); NCAR Technical Note NCAR/TN-398+STR; University Corporation for Atmospheric Research: Boulder, CO, USA, 1994; p. 121. [Google Scholar] [CrossRef]
  76. Trancoso, A.R. Operational Modelling as a Tool in Wind Power Forecasts and Meteorological Warnings. Ph.D. Thesis, University of Lisbon, Lisbon, Portugal, 2012. [Google Scholar]
  77. Ruiz Villarreal, M.; Bolding, K.; Burchard, H.; Demirov, E. Coupling of the GOTM turbulence module to some three-dimensional ocean. In Marine Turbulence: Theories, Observations, and Models—Results of the CARTUM Project; Baumert, H.Z., Simpson, S., Sündermann, J., Eds.; Cambridge University Press: Cambridge, UK, 2005; pp. 225–237. [Google Scholar]
  78. Autret, E.; Paul, F.; Tandéo, P.; Hackett, B. For. Level 3 and 4 ODYSSEA SST Products over the Global Ocean and North Western Shelves—SST _ GLO _ SST _ L3 _ NRT _ OBSERVATIONS _ 010 _ 010—SST _ NWS _ SST _ L4 _ NRT _ OBSERVATIONS _ 010 _ 003; EU Copernicus Marine Service–Public: Brussels, Belgium, 2017. [Google Scholar]
  79. Lambda Project. CMEMS Service Evolution Project. 2020. Available online: http//www.cmems-lambda.eu/ (accessed on 6 April 2021).
  80. de Pablo, H.; Sobrinho, J.; Garcia, M.; Campuzano, F.; Juliano, M.; Neves, R. Validation of the 3D-MOHID hydrodynamic model for the Tagus coastal area. Water 2019, 11, 1713. [Google Scholar] [CrossRef] [Green Version]
  81. Brewin, R.J.W.; de Mora, L.; Billson, O.; Jackson, T.; Russell, P.; Brewin, T.G.; Shutler, J.D.; Miller, P.I.; Taylor, B.H.; Smyth, T.J.; et al. Evaluating operational AVHRR sea surface temperature data at the coastline using surfers. Estuar. Coast. Shelf Sci. 2017, 196, 276–289. [Google Scholar] [CrossRef]
Figure 1. Representation of the study area as well as the most important bathymetric features and location of the time series used in the validation of the results. The area corresponds to the PCOMS grid domain (Level 2) and its bathymetry (on the map). Level 1 is not shown in the figure, as it is equal to that of Level 2 but with 2 more numerical cells in each direction. On the right side, a representation of the schematic river and estuary domains is shown. In this case the dimensions of these estuary and river domains are shortened in order to fit in the image. Their real dimensions are presented in Table 1.
Figure 1. Representation of the study area as well as the most important bathymetric features and location of the time series used in the validation of the results. The area corresponds to the PCOMS grid domain (Level 2) and its bathymetry (on the map). Level 1 is not shown in the figure, as it is equal to that of Level 2 but with 2 more numerical cells in each direction. On the right side, a representation of the schematic river and estuary domains is shown. In this case the dimensions of these estuary and river domains are shortened in order to fit in the image. Their real dimensions are presented in Table 1.
Water 13 02284 g001
Figure 2. Schematic domain (dashed green squares) nested in the PCOMS domain (white, blue, and grey boxes). The lateral discharge is computed using the schematic domain cell velocities u C D * located inside the land cells of the regional domain (grey), and the discharge is implemented in the first regional domain cell adjacent to the land cell (blue square). The momentum produced by the discharge velocity is added to regional domain in the face of the cell opposed to the land cell u P D * , which is located in the centre of the u velocity cell (red dashed line). For a more detailed description of these computations, the reader is referred to [38].
Figure 2. Schematic domain (dashed green squares) nested in the PCOMS domain (white, blue, and grey boxes). The lateral discharge is computed using the schematic domain cell velocities u C D * located inside the land cells of the regional domain (grey), and the discharge is implemented in the first regional domain cell adjacent to the land cell (blue square). The momentum produced by the discharge velocity is added to regional domain in the face of the cell opposed to the land cell u P D * , which is located in the centre of the u velocity cell (red dashed line). For a more detailed description of these computations, the reader is referred to [38].
Water 13 02284 g002
Figure 3. Time scheme and grid communication of the offline two-way methodology proposed in this article. Each day is run twice (or more) in order to receive updated open boundary conditions (schematic channels) and discharges (PCOMS).
Figure 3. Time scheme and grid communication of the offline two-way methodology proposed in this article. Each day is run twice (or more) in order to receive updated open boundary conditions (schematic channels) and discharges (PCOMS).
Water 13 02284 g003
Figure 4. Salinity and temperature time series produced by the three PCOMS implementations and the Silleiro in situ moored buoy during the period between 1 January and 1 May 2018.
Figure 4. Salinity and temperature time series produced by the three PCOMS implementations and the Silleiro in situ moored buoy during the period between 1 January and 1 May 2018.
Water 13 02284 g004
Figure 5. Salinity and temperature time series produced by the three PCOMS implementations and the Estaca de Bares in situ moored buoy during the period between 1 January and 1 May 2018.
Figure 5. Salinity and temperature time series produced by the three PCOMS implementations and the Estaca de Bares in situ moored buoy during the period between 1 January and 1 May 2018.
Water 13 02284 g005
Figure 6. Average surface salinity maps for the entire domain (top) and a zoom over the area of interest on the west coast of the IP (bottom) produced by the three implementations: no-rivers (left), detached (middle), and integrated (right) methodologies over the peak flow period (18 March to 1 April 2018).
Figure 6. Average surface salinity maps for the entire domain (top) and a zoom over the area of interest on the west coast of the IP (bottom) produced by the three implementations: no-rivers (left), detached (middle), and integrated (right) methodologies over the peak flow period (18 March to 1 April 2018).
Water 13 02284 g006
Figure 7. Difference between average surface salinity maps for the entire domain (top) and a zoom over the area of interest in the west coast of the IP (bottom) produced by: integrated minus detached (left) and integrated minus no-rivers (right) over the peak flow period (18 March to 1 April 2018).
Figure 7. Difference between average surface salinity maps for the entire domain (top) and a zoom over the area of interest in the west coast of the IP (bottom) produced by: integrated minus detached (left) and integrated minus no-rivers (right) over the peak flow period (18 March to 1 April 2018).
Water 13 02284 g007
Figure 8. Average surface temperature maps for the entire domain (top) and a zoom over the area of interest on the west coast of the IP (bottom) produced by the two implementations: Satellite SST (left), no-rivers (middle) and integrated (right) methodologies and at the peak flow (20 March 2018).
Figure 8. Average surface temperature maps for the entire domain (top) and a zoom over the area of interest on the west coast of the IP (bottom) produced by the two implementations: Satellite SST (left), no-rivers (middle) and integrated (right) methodologies and at the peak flow (20 March 2018).
Water 13 02284 g008
Figure 9. Difference between surface temperature maps for the entire domain (top) and a zoom over the area of interest on the west coast of the IP (bottom) produced by: no-rivers minus satellite (left) and integrated minus satellite (right) at the peak flow period (20 March 2018).
Figure 9. Difference between surface temperature maps for the entire domain (top) and a zoom over the area of interest on the west coast of the IP (bottom) produced by: no-rivers minus satellite (left) and integrated minus satellite (right) at the peak flow period (20 March 2018).
Water 13 02284 g009
Table 1. Dimensions of the schematic river and estuary domains.
Table 1. Dimensions of the schematic river and estuary domains.
Model DomainDimensionsHorizontal Grid
Minho14 × 33.4 km × 630 m
Douro15 × 31.9 km × 300 m
Mondego14 × 33.8 km × 170 m
Tagus14 × 34.3 km × 5.3 km
Sado14 × 36.8 km × 2.2 km
Guadiana3 × 14300 m × 3.9 km
Guadalquivir14 × 39.4 km × 520 m
Table 2. Model setup configuration for the TagusROFI area.
Table 2. Model setup configuration for the TagusROFI area.
SettingsLevel 1—West IberiaLevel 2—PCOMSSchematic Rivers and Estuaries
Model
characterization
2D—Barotropic3D—Baroclinic3D—Baroclinic
Grid corners33.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 dimension208 × 156177 × 125a
BathymetryEMODnet b
Hydrography portal
EMODnet b
Hydrography portal
[38]
Horizontal GridRegular: (≈5.7 km)Regular: (≈5.7 km)a
Vertical Grid1 layer7 Sigma Layer (0–8.68 m)
43 Cartesian layers
7 Sigma Layer (0–8.68 m)
43 Cartesian layers
Δt60 s60 s30 s
TidesFES2004From Level1From PCOMS
OBC WaterFrom MercatorOcéan PSY2V4 (Releases 1–4)From PCOMS
AssimilationFlow relaxation scheme of 10 cells with a time decay of 1 day at the open boundary and 109 s inside the domainFlow relaxation scheme of 3 or 4 cells with a time decay of 900 s at the open boundary and 109 s inside the domain
OBC Atmosphere MM5 (9 km)
DischargesNoFrom schematic rivers and estuariesMinho, Douro, Mondego, Tagus, Sado, Guadiana, Guadalquivir
TurbulenceGOTM cGOTM c
BottomRugosity: 0.0025 m2 s−1Rugosity: 0.0025 m2 s−1Rugosity: 0.0025 m2 s−1
a dimensions are presented in Table 1. b https://www.emodnet-bathymetry.eu/ (accessed on 1 March 2017). c [77].
Table 3. River flow between January and May of 2018 and respective source of data used to feed the schematic channels.
Table 3. River flow between January and May of 2018 and respective source of data used to feed the schematic channels.
Flow Time SeriesSource
Water 13 02284 i001Watershed model
Water 13 02284 i002In situ (Crestuma dam)
(https://snirh.apambiente.pt/)
Water 13 02284 i003In situ (Pte Coimbra)
(https://snirh.apambiente.pt/)
Water 13 02284 i004In situ (Almourol)
(https://snirh.apambiente.pt/)
Water 13 02284 i005Watershed model
Water 13 02284 i006Watershed model
Water 13 02284 i007In situ (Alcala)
(https://www.emodnet-physics.eu/map/platinfo/piroosplot.aspx?platformid=974539)
Table 4. Statistical parameters obtained for the validation of the three methodologies against field data (Silleiro moored buoy).
Table 4. Statistical parameters obtained for the validation of the three methodologies against field data (Silleiro moored buoy).
ComparisonParameter1 January–1 May 2018 (n = 2876)18 March–1 May 2018 (n = 1053)
rr2RMSDBIASrr2RMSDBIAS
Data-DetachedSal0.7520.5650.9820.0130.8630.7450.454−0.522
Temp0.8370.7010.6710.4150.8420.7090.3920.647
Data-IntegratedSal0.8650.7480.9630.0310.9200.8470.333−0.357
Temp0.8450.7140.6530.3940.8560.7330.4040.660
Data-No-RiversSal0.2270.0521.047−0.111−0.7120.5070.693−0.833
Temp0.7530.5670.6840.420−0.1780.0320.2960.392
Table 5. Statistical parameters obtained for the validation of the three methodologies against field data (Estaca de Bares moored buoy).
Table 5. Statistical parameters obtained for the validation of the three methodologies against field data (Estaca de Bares moored buoy).
ComparisonParameter1 January–1 May 2018 (n = 2876)18 March–1 May 2018 (n = 1053)
rr2RMSDBIASrr2RMSDBIAS
Data-DetachedSal0.4870.2370.1620.025−0.267−0.071−0.267−0.015
Temp0.8260.6830.4400.0850.1590.0250.4990.432
Data-IntegratedSal0.4490.2020.1640.025−0.1430.0210.265−0.014
Temp0.8270.6840.4170.0660.1920.0370.4640.393
Data-No-RiversSal0.3710.1380.1660.020−0.4540.2060.268−0.023
Temp0.8280.6860.4330.0810.1820.0330.4940.428
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sobrinho, J.; de Pablo, H.; Campuzano, F.; Neves, R. Coupling Rivers and Estuaries with an Ocean Model: An Improved Methodology. Water 2021, 13, 2284. https://doi.org/10.3390/w13162284

AMA Style

Sobrinho J, de Pablo H, Campuzano F, Neves R. Coupling Rivers and Estuaries with an Ocean Model: An Improved Methodology. Water. 2021; 13(16):2284. https://doi.org/10.3390/w13162284

Chicago/Turabian Style

Sobrinho, João, Hilda de Pablo, Francisco Campuzano, and Ramiro Neves. 2021. "Coupling Rivers and Estuaries with an Ocean Model: An Improved Methodology" Water 13, no. 16: 2284. https://doi.org/10.3390/w13162284

APA Style

Sobrinho, J., de Pablo, H., Campuzano, F., & Neves, R. (2021). Coupling Rivers and Estuaries with an Ocean Model: An Improved Methodology. Water, 13(16), 2284. https://doi.org/10.3390/w13162284

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop