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

: 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.


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 areaalso 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 estuarieseasy 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.

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 highpressure system in the Azores [39,40], and its seasonality and intensity changes with the relative position of these two systems. 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. 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].

MOHID Water
The MOHID Water model is a finite-differences family model that evolved to a finitevolume 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, where 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 (−ℎ) and a depth , one gets the vertical velocity at that level.
where q are natural or anthropogenic local water sources. Integrating up to the free surface, one gets the free surface level rate of change: In this equation, includes rain and evaporation. Using the hydrostatic approximation, the vertical momentum equation becomes That integrated over the water column between the free surface and a depth gives if one performs the derivative along the horizontal axis, one gets the pressure gradient, and the horizontal momentum equations can be written as: 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 where SP stands for sink-sources of the property (P) in question.
Horizontal and vertical advection are computed using the TVD advection schemea 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.

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 crosssection 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: which is a simple summation of all local domain fluxes (ex: * 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, * represents the normal velocity of the cross section (in Figure 2 the dashed green line, where * is computed. stands for the vertical area of the CD cell associated with * , 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. . The lateral discharge is computed using the schematic domain cell velocities * 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 * , 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].
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.

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.

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.

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 1 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 10 5 s in the outmost open boundary cell to 10 9 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 10 9 s. Additionally, a biharmonic filter of 5.5 × 10 9 m 4 s −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 Tables 1 and 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.

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).

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 (Figures 4 and 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 (March 18-April 1) 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 (Tables 4 and 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 (r 2 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. Table 4. Statistical parameters obtained for the validation of the three methodologies against field data (Silleiro moored buoy).  Table 5. Statistical parameters obtained for the validation of the three methodologies against field data (Estaca de Bares moored buoy). 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].

Comparison
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.

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.  [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 (Figures 6 and 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 Figures 6 and 7 for a better visualization of this analysis in the ROFI areas.

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 norivers 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 Figures 8 and 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 (Figures 8 and 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.

Conclusions
A good model representation of the coastal hydrodynamic processes in the IPnamely, those associated with freshwater discharges, such as the WIBP and the WICPrequires 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 applicationhas 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.