Coupled Hydrodynamic and Biogeochemical Modeling in the Galician Rías Baixas (NW Iberian Peninsula) Using Delft3D: Model Validation and Performance

1. Introduction

Coastal zones and estuaries are among the most important ecosystems globally. These transitional environments, where freshwater meets seawater, are highly productive and provide essential food resources and habitats for species of both ecological and economic significance [1,2]. Despite their importance, estuaries worldwide face numerous environmental and anthropogenic perturbations that alter habitats and disrupt biotic communities [3,4]. Key drivers of these changes include rising sea levels, eutrophication, and direct human activities within estuarine systems [5,6]. As these pressures intensify, the stability of aquatic communities and ecosystems is increasingly at risk, leading to deteriorating water quality and changes in ecosystem structure, function, and productivity [7,8,9,10].

The Galician Rías Baixas, located in the northwest Iberian Peninsula, are four estuaries of significant ecological, economic and social importance. These estuaries support rich marine biodiversity, including commercially valuable species, such as mussels and clams, that play a vital role in the economy of the region [11,12]. For example, in 2023, small-scale bivalve fisheries contributed approximately €62 million and sustained the livelihoods of around 4,300 fishers [13]. Moreover, the Rías Baixas support approximately 32% of the world’s production of Mytilus galloprovincialis [14]. However, like other coastal zones, the Rías Baixas are also experiencing increasing environmental pressures. Numerous studies have reported adverse effects and even mortality events in key species related to rising seawater temperatures [15,16,17] and salinity drops [18,19,20,21]. Other stressors, such as the increasing frequency and intensity of harmful algal blooms [22] and increased concentrations of metals and metalloids in sediments [23], also pose significant risks to biodiversity and the local economies [24].

In response to these challenges, numerical models have become essential tools for understanding estuarine dynamics and predicting the impacts of human activities, extreme events and climate change [25,26]. Models provide valuable insights into how changes in one or various factors influence estuarine systems [27,28]. Additionally, they help overcome limitations in field observations, allowing for a comprehensive assessment of hydrodynamic, water quality and sediment transport [29]. Such information is critical for improving ecosystem management, enhancing service delivery, and ensuring ecosystem resilience. These models have been applied to assess pollution and contamination responses in China (e.g., [30,31,32]); to predict algal blooms and eutrophication in the Dutch North Sea [33,34]; and to enhance estuarine knowledge and evaluate climate change impacts on aquaculture in Portugal (e.g., [35,36,37,38]).

In the Rías Baixas, previous research has extensively utilized numerical models to improve the understanding of hydrodynamic processes (e.g., [10,39,40,41,42,43]). Additionally, simplified models have provided valuable insights into biogeochemical cycles and nutrient dynamics in the region [43,45,46,47]. To our knowledge, however, only [48] have addressed the complex task of jointly simulating hydrodynamic and water quality processes in three dimensions.

This study presents the validation of a coupled hydrodynamic and water quality model developed using Delft3D-4 software, designed for a comprehensive spatial and temporal analysis of environmental conditions within the Rías Baixas. The model, which integrates the Delft3D-FLOW and Delft3D-WAQ modules, was calibrated and validated against observed data for the four estuaries of the Rías Baixas. Simulations were conducted from September 1, 2016, to June 1, 2018, with the first 6 months serving as spin-up and the study period spanning from March 2017 to May 2018. A key advancement of this model, compared to other water quality models developed for this area, is the inclusion of the Minho River, which has a significant influence on the Rías Baixas impacting residence times and primary production dynamics through plume intrusions. Additionally, this study advances the validation process by including nutrient dynamics, a crucial factor in processes such as eutrophication, which can severely impact biodiversity and ecosystem health. The insights gained from this research will support future studies and inform management strategies to protect and improve water quality in these ecosystems.

2. Study Area

The Rías Baixas comprise four specific subtypes of estuary, known as rias, located in the NW Iberian Peninsula, at the northern boundary of the North Atlantic Upwelling System (Figure 1b, [49]). From North to South, these estuaries are named Ría de Muros, Ría de Arousa, Ría de Pontevedra, and Ría the Vigo. They have a NE-SW orientation along their main axis, are bordered by steep hills, and were formed by the flooding of river valleys due to rising sea levels [50,51,52]. The average water depths within the rias range from 5-10 m in the inner parts, where the main rivers flow into the estuary, to 40-60 m at the outer (SW) entrance to the sea. The rias, with the exception of the Ría de Muros, are partially protected by small islands at their mouths, which corresponds to structural highs [51,53,54]. More information about the morphologic characteristics of the Rías Baixas are shown in Table 1.

Figure 2. (a) Location of the study area, with the box indicating the modeled region. (b) Close-up view of the study area, showing the locations of the in-situ data stations (dots) used to perform model calibration and validation.

Figure 2. (a) Location of the study area, with the box indicating the modeled region. (b) Close-up view of the study area, showing the locations of the in-situ data stations (dots) used to perform model calibration and validation.

Hydrodynamically, Rías Baixas typically behave as partially mixed estuaries with a two-layered positive circulation [39]. Freshwater flows offshore through the surface while dense, saline ocean water enters through the bottom layers. However, this circulation is strongly influenced by upwelling and downwelling events [50]. During the upwelling season, from April to October, equatorward winds over the shelf enhance the positive circulation by driving Eastern North Atlantic Water into the estuaries along the bottom increasing the net inflow of cold, nutrient-rich water [45,55,56]. This nutrient influx makes the rias highly productive, with elevated remineralization and sedimentation due to intensified physical and biological activity [44,57]. Contrarily, during downwelling events, when poleward winds dominate, the positive circulation weakens or even reverses, trapping water inside the rias. This leads to nutrient depletion, except in the innermost areas, where river outflow becomes the main nutrient source, boosting biological activity [44,58]. The Minho River can also induce a reversal of the positive circulation through freshwater intrusions during downwelling conditions and high river discharge [59,60,61]. Under these conditions, water exchange with the continental shelf is limited, potentially reducing water quality [59]. Additionally, intrusions can also increase nutrient levels, promoting phytoplankton blooms that can enter the rias [62].

Nutrient dynamics within the Rías Baixas are also significantly influenced by upwelling and downwelling events. Upwelling introduces nutrient-rich water into the deeper layers, while thermal stratification limits the upward movement of nutrients, resulting in high concentrations at depth and lower levels near the surface. In contrast, during downwelling, the influx of less saline surface waters from continental runoff—especially in winter—elevates nutrient levels in the upper layers. This spatial variability in nutrient distribution is shaped by the interaction of continental runoff and upwelled water [63]. Primary production follows this variability, with lower production in winter but reaching levels exceeding 10 g C m² d⁻¹ in spring and summer as phytoplankton growth accelerates. Notably, about 55% of this production is available to higher trophic levels [64,65]. Phytoplankton biomass peaks seasonally, reaching around 8 mg m⁻³ in spring and autumn, while winter values remain below 1 mg m⁻³ and summer values average around 5 mg m⁻³ [66,67].

The tidal regime in the Rías Baixas is semi-diurnal and mesotidal, with tides occurring approximately every 12.25 hours. The tidal amplitude ranges from 1.3 m during neap tides to 3.4 m during spring tides. Tides are the most energetic force driving water circulation within the rias.

3. Methodology

3.1. In Situ Data for Model Calibration and Validation

In situ data measured by the Technological Institute for Monitoring the Marine Environment in Galicia [69] at representative points in the inner, middle, and outer sections of each estuary (Figure 1b) were used for model calibration and validation. Samples for nitrate, ammonium, phosphate, and silicate were collected weekly during the study period (March 2017 to May 2018) using a PVC hose to sample three depth intervals: 0-5m, 5-10m, and 10-15m. Additionally, biweekly samples were collected from the inner (V1) and outer stations (M3, A4, P3, and V4) at 1m depth and at the bottom. Nutrient analysis was based on colorimetric methods via continuous segmented flow analysis, as described by [68], using unfiltered samples. Weekly vertical profiles of seawater temperature, salinity, dissolved oxygen, and fluorescence (used to estimate the chlorophyll-a), were collected using a SBE25 CTD [69].

3.2. Numerical Modeling

This study uses the open-source Delft3D-4 suite [70] Access to to perform a three-dimensional (3D) hydrodynamic and water quality modeling of the Rías Baixas. The following modules were used:

The Delft3D-FLOW module was used to simulate water circulation, accounting for key physical processes such as tides, currents, and coastal structures. The 3D approach allows for a detailed spatial analysis of hydrodynamic conditions, including vertical stratification and circulation patterns within the estuaries. The grid and parametrization used in this study were presented and validated by [61], ensuring accurate flow predictions across the Rías Baixas.

The main focus of this study is implementing the 3D Delft3D-WAQ module, which simulates water quality dynamics through three-dimensional advection-diffusion equations. The model captures a wide range of ecological and chemical processes, including nutrient cycling, phytoplankton dynamics, eutrophication, and oxygen levels. It also considers the behavior of substances at air-water and water-sediment interfaces. One of the model’s keys is its extensive library of constituents, allowing users to customize the chemical and biological composition of aquatic systems effectively.

3.1.1. Hydrodynamic Module

This section summarizes the setup of the hydrodynamic module, with a detailed description available in [61]. The model uses a structured curvilinear grid that spans from 41.18° to 43.50° N and from 10.00° to 8.33° W (Figure 1a). The horizontal grid resolution varies, starting at 2200 × 800 m at the western boundary and refining to 220 × 140 m in the inner areas of the Rías Baixas. The grid is vertically divided into 16 sigma layers with refined top layers. The computational time step is 0.5 minutes. A spin-up period of 6 months was considered to ensure numerical stability and accuracy.

Bathymetric data were sourced from several providers. Nautical charts from the Hydrographic Institute of the Spanish Navy provided the data for the Muros and Arousa rias, while bathymetric data for the Vigo and Pontevedra rias were obtained from the Spanish General Fishing Secretary. The Hydrographic Institute of the Portuguese Navy provided data for the Minho estuary. To fill data gaps data from the General Bathymetric Chart of the Oceans was incorporated [71].

Daily salinity and water temperature at the open ocean boundary and for initial conditions were sourced from the Atlantic-Iberian Biscay Irish-Ocean Physics Reanalysis (product ID: IBI_MULTIYEAR_PHY_005_002), available through the Copernicus Marine Environment Monitoring Service (CMEMS) website [72]. This dataset has a spatial resolution of 0.083° × 0.083° and 50 vertical levels.

Astronomical forcing at the oceanic boundary was applied using thirteen tidal harmonic constituents (M2, S2, N2, K2, K1, O1, P1, Q1, MgF, MM, M4, MS4, MN4) based on a high-resolution model derived from TOPEX/Poseidon Altimetry data, with a spatial resolution of approximately 25 meters.

Atmospheric surface boundary conditions (air temperature, relative humidity, net solar radiation, surface pressure, and wind components) were imposed using hourly data from the ERA5 dataset [73] with a spatial resolution of 0.25° x 0.25°.

River discharges for the five main tributaries flowing into the Rias Baixas were imposed as time series with a daily resolution, sourced from Meteogalicia [74]. Minho river discharge data was provided by the Conferencia Hidrográfica Miño-Sil [75] and imposed as time series with a daily resolution.

3.1.1. Water Quality module

The Delft3D-WAQ model configuration was optimized to balance computational efficiency, result quality, and data availability. The input variables selected for the simulation were ammonium (NH₄⁺), nitrate (NO₃^-), phosphate (PO₄³⁻), dissolved silica (Si), dissolved oxygen (DO), total inorganic carbonate (TIC), alkalinity, biogenic silica (Opal-Si), diatoms, green algae, particulate organic carbon (POC1, fast decomposing fraction), nitrogen (PON1, fast decomposing fraction) and phosphorus (POP1, fast decomposing fraction). Key processes selected for the simulation of each variable are summarized in Table S1. Detailed descriptions of the formulations for these processes can be found in the Delft3D-WAQ technical reference manual [76].

The simulations were run with a computational time step of 30 min and a 6-month spin-up period. The 16 vertical layers from the FLOW grid were aggregated vertically into the WAQ module, where each layer from the hydrodynamic module directly corresponds to one layer in the WAQ module (1:1 aggregation approach). Default parameters were used, except for those listed in Table 2. Carbonaceous biological oxygen demand (CBOD5), sediment oxygen demand (SOD), and zooplankton parameters were incorporated into the model as spatially homogeneous time series, following the values from [77,78,79] respectively (see Figure S1). Irradiation was applied uniformly across the water surface using a net solar radiation time series. Julian day numbers were calculated starting from January 1 of the corresponding year.

Daily water quality data for the oceanic open boundary was obtained from the Atlantic-Iberian Biscay Irish Ocean BioGeoChemistry NON-ASSIMILATIVE Hindcast (product ID: IBI_MULTIYEAR_BGC_005_003), provided by CMEMS. This product provides 3D biogeochemical data with a horizontal resolution of 0.083º x 0.083º and 50 vertical layers. Variables such as NH₄⁺, NO₃^-, PO₄³⁻, Si, and DO were directly retrieved from this dataset. TIC, alkalinity, diatoms, green algae, POC1, PON1, POP1 and Opal-Si were derived using other variables from IBI_MULTIYEAR_BGC_005_003 (see Table 3).

Due to the lack of available data for river boundary water properties during the simulation period, climatological and bibliographic data were used instead (see Table 4 for further details).

Finally, the atmospheric condition for surface radiation was forced using hourly net solar radiation as a time series from ERA5 [72].

3.3. Model Calibration

The model was run from September 1, 2016, to June 1, 2018, with the first 6 months serving as spin-up period. Calibration covered the study period from March 2017 to May 2018, focusing on accurate simulation of water temperature, salinity, and water quality variables in the Rías Baixas. Calibrated variables included NH₄⁺, NO₃⁻, PO₄³⁻, Si, DO, and Chl-a. Model outputs were compared with in-situ measurements, with the best fit achieved using the parameterization shown in Table 2.

Model performance was assessed using qualitative and quantitative methods based on [35,86,87]. Evaluation metrics included Root Mean Square Error (RMSE), relative RMSE (rRMSE), percentage Bias (pBias), cost function (CF), Pearson correlation coefficient (r), and reliability (F).

RMSE quantifies the difference between predicted ($P_i$) and observed ($O_i$) values. Lower RMSE values indicate better model performance. It was calculated as:

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

(1)

where N is the total number of observations.

rRMSE normalizes the RMSE by the range between the 1st and 99th percentiles of observed values, O_99 and respectively O_1, reducing the impact of outliers. Acceptable thresholds for model performance are ±20% for water temperature, salinity, and DO; ±50% for nutrients, and ±100% for algae biomass following [88]. It was calculated as:

$$rRMSE(\%)=\left({\displaystyle \frac{RMSE}{O_{99}-O_1}}\right),$$

(2)

Cost Function (CF) is a dimensionless metric that assesses the fit between model outputs and observations. According to the [89], CF < 1 indicates very good model performance, 1–2 good, 2–5 reasonable, and CF > 5 poor:

$$CF={\displaystyle \frac{1}{N}}{\sum }_{n=1}^N{\displaystyle \frac{|O_i-P_i|}{{\sigma }_D}},$$

(3)

were ${\sigma }_D$ is the standard deviation of the data.

pBias measures whether the model systematically underestimates or overestimates observed values. According [90], an absolute value of pBias < 10% is excellent, 10–20% is very good, 20-40% is good, and values above 40% are weak. It was calculated as:

$$pBias\left(\%\right)={\displaystyle \frac{{\sum }_{n=1}^N\left(P_i-O_i\right)}{{\sum }_{n=1}^N{O}_i}}x100,$$

(4)

Pearson correlation coefficient (r) quantifies the correlation between observed and simulated data. It was calculated as:

$$r={\displaystyle \frac{{\sum }_{i=1}^N\left(P_i-{\underset{\_}{P}}_i\right)(O_i-{\underset{\_}{O}}_i)}{\sqrt{{\sum }_{i=1}^N{\left(P_i-{\underset{\_}{P}}_i\right)}^2}\sqrt{{\sum }_{i=1}^N{\left(O_i-{\underset{\_}{O}}_i\right)}^2}}},$$

(5)

Reliability (F) indicates the confidence in the correlation between model and in situ data. Reliability categories are 99%, 95%, and 75%, with higher F values indicating greater confidence in the model’s correlation. It is derived from the p-value and was calculated as:

$$F=100*(1-P_{value})$$

(6)

4. Results

From the complete set of model results, only outputs for Ría de Muros and Ría de Arousa are presented in the main text. Ría de Arousa, with extensive mussel rafts and significant aquaculture, contrasts with the minimally impacted Ría de Muros, allowing for a comparison between areas with and without major aquaculture influence. Detailed results for the other two rias, Ría de Pontevedra and Ría de Vigo, as well as additional depth layers, are provided in the Supplementary Material.

4.1. Model Calibration and Performance

4.1.1. Hydrodynamic

Near surface (averaged from 0 to 5 m depth) model outputs showed good agreement with in-situ observations across all stations. RMSE values for both variables (temperature and salinity) tended to be higher at inner stations and decreased toward middle and outer stations, reflecting spatial variations within the estuaries (Figure 2 and Figure 3). Specifically, the RMSE for water temperature ranged from 0.43 to 0.63 ºC in the Ría de Muros (Figure 2) and from 0.57 to 0.78 ºC in the Ría de Arousa (Figure 3). For salinity, RMSE ranged from 0.32 to 1.69 in Ría de Muros and from 0.57 to 1.73 in Ría de Arousa. The best fit was observed at middle and outer stations. Similarly, the model accurately captured the seasonal water temperature and salinity patterns observed in the field data, with higher values during summer (June, July, and August) and lower values in winter (December to February). Spatial gradients were also well reproduced, with higher water temperatures recorded at inner stations (M1, A1, and A2) and higher salinity levels observed at outer stations (M3 and A4). Additionally, the model effectively simulated the vertical distribution of water temperature and salinity providing a robust representation of surface and subsurface dynamics (Figures S2 to S9). These results suggest that the model is capable of replicating both horizontal and vertical hydrodynamic variability in the Rías Baixas, essential for further water quality simulations.

Pearson’s correlation coefficient (r) for salinity followed a similar pattern to RMSE, with higher values at the middle and outer stations compared to the inner ones (Figure 2 and Figure 3). In Ría de Muros, r for salinity ranged from 0.80 to 0.87, while in Ría de Arousa it ranged from 0.71 to 0.94. Water temperature, however, showed consistently high correlations across most stations, with r values of 0.9 or higher, except for slightly lower correlations at M1 (0.88) and V3 (0.8; Figure S15).

The rRMSE for both salinity and water temperature was below the threshold of 20% (Figure 4 and Figure S20), demonstrating the model’s strong predictive capabilities. However, for water temperature, rRMSE was slightly higher in deeper layers compared to near-surface (0–5 meters), while for salinity, rRMSE was lower at greater depths. At specific stations, such as A1 for salinity (Figure 4) and V3 (Figure S20) for water temperature, rRMSE exceeded 20% in the bottom and surface layers. Despite these discrepancies, both variables showed 99% reliability across all stations.

The cost function (CF) values for salinity and water temperature ranged from 0.15 to 0.6 (Figure 5 and Figure S21), while the pBIAS remained below 6% (Figure 6 and Figure S22), indicating very good model performance for both metrics. Notably, the highest CF values were observed for salinity at inner stations of the Ría de Arousa (0.85 at A1) and Ría de Pontevedra (0.78 at P1; Figure S21), as well as the bottom layer in Ría de Muros (0.66 at M3). Importantly, despite these elevated CF values, the pBIAS for these stations remained below 1% reflecting a strong model accuracy.

4.2. Water Quality

The model’s performance in predicting water quality variables exhibited a similar trend to that of salinity, with RMSE values generally higher at inner stations and decreasing towards the middle stations and the estuary mouth. NO₃^- and Si had the highest average RMSE, ranging from 2.46 to 4.02 µmol L^-1 and from 2.27 to 6.61 µmol L^-1 in Ría de Muros (Figure 7), and from 2.96 to 3.18 µmol L^-1 and from 3.97 to 6.43 µmol L^-1 in Ría de Arousa (Figure 8). Underestimations of nutrient concentrations were observed during autumn and winter, particularly in inner areas and upper layers of the rias of Arousa and Vigo, as well as for NH₄⁺ and Si in the Ría de Pontevedra (see stations P1, P2, V1 and V2; Figures S10 and S15). Deeper layer concentrations at the mouth of the estuaries were also slightly underestimated during summer (Figures S5 and S9). Peak concentrations of NH₄⁺ and Si were recorded in the rias of Muros and Pontevedra in October 2017, and in the rias of Muros, Arousa and Vigo in spring 2018. RMSE for Chl-a ranges from ~3 mg m-3 in outer stations to ~4.4 mg m-3 in inner stations (Figure 7 and Figure 8), with a tendency to underestimate winter concentrations in the Ría de Muros and overestimate spring peaks in the rias of Arousa and Pontevedra (Figure S10). Nevertheless, the model was able to satisfactorily reproduce inorganic nutrient concentrations (for the upper 15 meters) in spring and summer, and Chl-a concentrations in summer and autumn. Similarly, spatial gradients were also well reproduced with the highest (lowest) concentrations occurring at internal (external) stations.

Regarding rRMSE, all water quality variables exhibited average values below or equal to acceptable thresholds across all stations (~21% for NO₃^-, ~23% for PO₄^3-, ~26% for NH₄⁺, Si, and DO, and 30% for Chl-a; Figure 9 and Figure S20). The cost function (CF) indicated strong model performance, with values below 1 across the region (~0.49 for NO₃^-, ~0.62 for PO₄^3-, ~0.64 for Si, ~0.73 for NH₄⁺, ~0.83 for DO, and ~0.99 for Chl-a; Figure 10 and Figure S21). CF exceeded 1 only at inner stations (V1, A1, and M1), particularly for NH₄⁺, Si, DO, and Chl-a. NO₃^-, PO₄^3-, and DO exhibit the best overall performance, with pBIAS values of ~21% for NO₃^- and PO₄^3-, and ~6% for DO (Figure 11 and Figure S22). Correlation coefficients (r) were above 0.60 for NO₃^- and 0.50 for both PO₄^3- and DO, with reliability levels averaging ~97% for PO₄^3-, ~98% for DO, and ~99% for NO₃^-. In contrast, Si presented higher pBIAS values, exceeding 40% in most cases, while NH₄⁺ exhibited a similar pattern in deeper layers (10-15 meters and bottom) in the rias of Vigo, Pontevedra, and Arousa. Both nutrients showed positive correlations, with r values of 0.24 for NH₄⁺ and 0.47 for Si, while their reliability was lower, around 80% for NH₄⁺ and 86% for Si. Chl-a had the lower reliability (~77%; Figure 11) and r value (0.26), with an average pBIAS of ~37% across the region (Figure 11 and Figure S22).

5. Discussion

The coupled Delft3D-FLOW and Delft3D-WAQ modules exhibited higher accuracy at the estuary mouths compared to inner stations. This trend aligns with modeling errors reported for other estuarine systems, such as in Mondego [91] and Sado [38,87] estuaries. Discrepancies between model outputs and observations can be due to uncertainties in riverine boundary conditions, the exclusion of minor tributaries due to data limitations, and inadequate bathymetric data in shallower areas [92]. Moreover, in this study, water quality inputs in the riverine boundary were based on existing datasets from the region and, due to the limited availability of data, often represent different periods or climatological conditions and may not fully capture specific conditions of the simulation period. For instance, October 2017 was marked by the second-highest wildfire outbreak in Galicia [93], which likely influenced nutrient dynamics by contributing NH₄⁺, NO₃^-, and PO₄^3- through atmospheric deposition and post-fire runoff [94,95,96].

Challenges in reproducing nutrient dynamics in deeper layers likely derive from the simplified parameterization of processes such as phytoplankton production and remineralization. This study’s approach, which considered only two phytoplankton groups (diatoms and green algae) with most rates and parameters fixed may not fully capture the intra- and inter-annual variability within the phytoplankton community. Multiple phytoplankton groups differ in their carbon-specific photosynthetic rates and their responses to light and nutrient uptake [98,99,100,101,102]. Consequently, the simplified model configuration could lead to overestimated production and nutrient uptake, leading to depletion in the water column [103,104,105]. The model’s constant rates for remineralization processes, like nitrification and ammonification, also omit seasonal variations, potentially explaining nutrient peaks, particularly NH₄⁺, observed in October and November in upper layers [66,106].

Despite these constraints, this study has successfully implemented the coupled Delft3D-FLOW and Delft3D-WAQ model for the Rías Baixas, demonstrating proficiency in reproducing both physical variables (water temperature and salinity) and water quality parameters (inorganic nutrients, DO, and Chl-a). The performance of the hydrodynamic module is consistent with [61], whose grid and parameterization were adopted and aligns with studies by [48,107,108], which used the MOHID, ROMS, and Delft3D models for the same region. The performance of the water quality model lags behind that of hydrodynamic models. It was expected due to the complexity of modeling the natural and biological processes within aquatic environments led to the consensus that water quality models often perform worse than those for hydrodynamics [35,86,87]. Nevertheless, the results of this study are consistent with previous research [35,36,38,82,87,91,109]. For example, Chl-a RMSE values in the Breton Sound Estuary reported by [109] exceeded 10 mg m⁻³, while [82] reported values between 1 and 12 mg m⁻³ for Chl-a, 0.64 to 14.19 µmol L⁻¹ for NO₃^-, and 0.126 to 0.59 µmol L⁻¹ for PO₄³⁻ in the North Sea, all exceeding values obtained for the Rías Baixas. Notably, [38,87] reported similar pBIAS values for inorganic nutrients to those found in this study, indicating that the model demonstrates comparable or even superior performance in nutrient representation compared to other regions. Furthermore, DO results align closely with those obtained by [110] for the same region; however, the set-up used in the present study shows higher performance in reproducing transport conditions.

Usually, the most favorable results for environmental or aquatic models are achieved for water temperature and dissolved oxygen, followed by inorganic nutrients, with algal biomass exhibiting higher errors [111,112]. This study similarly observed the most accurate predictions for water temperature and DO, followed by inorganic nutrients (in the order of NO₃^-, PO₄^3-, NH₄^+, and Si), with the highest errors noted in the predictions for Chl-a across the region.

Finally, the coupled model captured the hydrodynamic and biogeochemical processes accurately, properly representing seasonal and spatial variations in the Rias Baixas. However, accounting for mussel farming dynamics, particularly in the rias of Arousa and Vigo, would likely enhance model performance. Mussels, as filter feeders, remove organic particles (including phytoplankton) from the water, reducing DO levels and depositing particulate organic matter on the seabed, which enriches the nutrient pool [87,113,114]. By incorporating mussel dynamics, the model could better capture phytoplankton fluctuations and nutrient recycling, thereby improving Chl-a and DO predictions in surface layers and helping to correct overestimated nutrients observed during certain periods.

6. Conclusions

This study successfully implemented a coupled Delft3D-FLOW and Delft3D-WAQ model for the Rías Baixas, effectively reproducing key hydrodynamic and biogeochemical processes. By accurately simulating water temperature, salinity, and essential water quality variables such as inorganic nutrients, dissolved oxygen (DO), and chlorophyll-a (Chl-a), the model provides a reliable representation of seasonal and spatial variations within the estuary system.

The model’s higher accuracy at estuary mouths compared to inner regions highlights the importance of refining boundary conditions and model resolution in complex estuarine environments. Discrepancies in nutrient dynamics and Chl-a levels underscore the intricacies of biogeochemical processes in estuarine settings, indicating a need for enhanced parameterization of nutrient cycling and phytoplankton dynamics.

Incorporating factors like phytoplankton diversity, seasonal variation in remineralization rates, and the influence of agricultural runoff and wildfire events could further improve model precision. Additionally, including dynamics related to mussel farming may be particularly beneficial for regions like Arousa and Vigo, where high mussel densities significantly impact nutrient and phytoplankton levels.

In summary, the coupled Delft3D model demonstrates valuable potential for predicting estuarine dynamics. Its application in the Rías Baixas serves as a model for similar studies, supporting better-informed decisions for ecosystem management, climate adaptation, and water quality enhancement in complex coastal systems.