Hydrological 0D and 2D Modelling of the Navio-Quebrado Coastal Lagoon (La Guajira, Colombia): A Challenging Exercise

Abstract

In a previous paper, we presented a very simple monitoring system. That system aimed at providing the basic information needed to set up a hydrological simulation model that would be used to understand the behavior of the system and to explore possible changes linked to future climates. After more than 1.5 years of observations, the gathered information was utilizable to test a 0D hydrological model. This exercise is presented here, together with a comparison with a more refined, and more demanding, 2D hydrodynamic model, enlightening strengths and weaknesses. Several details are revealed in this paper showing how even a simple case needs some type of art in order to overcome the unavoidable difficulties. Surprisingly nice answers have been obtained that illuminate the behavior of the lagoon and shed light on the key points to be improved both in monitoring and modelling. This exercise can be of help in several similar situations.

1. Introduction

This paper can be considered a second chapter of the exercise already described in a previous paper [1].

Coastal lagoons are delicate environments governed by energy, nutrients, and, most of all, water exchanges amongst, mainly, the atmosphere (precipitation and evaporation), the river basin (fluvial freshwater inputs), and the sea (tidal water exchange). The hydrodynamics of these systems plays a key role in regulating water quality, sediment transport, and ecological processes [2]. Understanding this dynamic is essential for effective coastal management, especially regarding climate change and sea level rise.

Hydrodynamic modelling is an essential tool for simulating the behavior of lagoons under different environmental and anthropogenic conditions [3]. Recent studies have shown that numerical modelling can capture complex interactions within coastal lagoons, including water circulation, interactions with adjacent marine environments, and responses to external forces [4]. For example, a two-dimensional model was successfully used to simulate the hydrodynamic process of Mediterranean lagoons and provided insight into the impact of artificial regulation on water level and flow patterns [5,6].

Advances allowed also for improving the calibration and validation of numerical hydrodynamic models, which has improved the representation of the flow and exchange processes in lagoons. For example, studies of the Cabras lagoon in Italy showed that factors such as bathymetry, tides, and wind stress were important factors to improve model accuracy [7]. These models help not only understand current hydrodynamics but also provide predictability of future climate change and anthropogenic changes in coastal lagoons [8].

A hydrological monitoring exercise has been conducted in the coastal lagoon Navio-Quebrado (Colombia) [1], and the obtained data allowed us to set up a 0D simulation model useful to understand its dynamics by firstly replicating quantitatively its behavior during the recorded period. Such a model can then be used to investigate future scenarios under climate change, which, particularly in Colombia, is considered to be a very serious threat [9,10]; indeed, a key question is what will happen to lagoons with a modified freshwater input from the feeding river, a different (increased) evaporation rate, a different pattern of precipitation, a changed morphology of the lagoon (possibly undergoing a sedimentation process), and a risen sea level.

Another lower-level, more technical issue concerns the assessment of strengths and weaknesses of 0D hydrological modelling (dynamic water balance) with respect to a hydrodynamic 2D model: these two tools are here compared in critical terms to provide hopefully useful hints. In the cited previous paper [1] on the monitoring system, the different components have been presented, which will be recalled in what follows.

Although based on preliminary results, the investigation described in this paper may contribute significantly to the management practice of coastal lagoons that, particularly in Colombia, constitute a peculiar and extremely important environmental asset. Indeed, the experience can help understand how to set up and improve monitoring systems; what can be expected from a simulation model like this; what are its delicate points; and how to choose between seemingly more attractive but more complex and computationally heavier 2D models rather than a simpler 0D model.

This paper starts with a position on the system schematization, where the relevant variables are introduced. We hence define all the factors relevant for a hydrological model of the lagoon (basically a water balance). Then, some important aspects of how missing information has been dealt with to set up a complete and consistent dataset to feed the model are provided. Some structural key differences between the 0D and the 2D models adopted are then pointed out. Results are finally presented and discussed for both the simplified and the full model and conclusions derived.

2. System Approach: Two Views of the Same System

As shown in the system sketch of Figure 1a, the hydrological dynamic system “Navio-Quebrado (alias Camarones) coastal lagoon” is subject to a number of exogenous, uncontrollable (or “natural”) input variables, namely: the freshwater inflow from the main river basin Tomarrazón-Camarones; the precipitation, p, determining the direct inflow input on the water surface S(t) and the runoff q(t) from the local watershed; the climatic variables determining the evaporation rate e(t) (and the evaporation flow E(t) = e(t)S(t)); and the sea level (y_S), governing the water exchange with the sea. There is also a human-induced variable, treated here as an exogenous variable: the occasional manual modification of the lagoon mouth. This intervention is sometimes carried out to lower the lagoon’s water level when the mouth has not naturally opened. This latter input affects the dynamic block of the “boca” (i.e., the mouth) that periodically opens, grows, then shrinks and again closes, so governing the mouth cross-section area A(t). This complex dynamic behavior is dependent on the elevation y_L of the lagoon water surface and that (y_S) of the sea. These two variables, together with the cross-section A(t) of the mouth, and possibly other factors (wind, salinity, …), determine the water exchange Q_B(t) between the sea and the lagoon (positive, when outgoing, or negative, when entering the lagoon).

Hence, the system has two state variables: the storage volume V(t) and the cross-section area A(t) of the boca. The dynamic behavior of the storage is relatively simple, given all the flows affecting it, as it is a straightforward dynamic water balance. Clearly, a key component is the morphological characteristic of the lagoon shelf; that is, its volume–elevation relationship V = V(y_L), as well as its area–elevation relationship S = S(y_L).

The dynamic of the boca cross-section is much more complicated: In a dry period, it closes completely because of a sand bar set by the sea waves in front of it; then, it comes the time when it opens again, usually after a significant river flood that rises the level over a given threshold (but sometimes locals intervene manually to “help” the process, or sometimes it is the sea that rises anomalously and opens from the outside); then, it remains opened for several weeks or months, allowing for an alternate lagoon–sea water exchange depending on the tide; and then it shrinks and closes again when the freshwater input drops. This process is often, but not always, repeated twice a year according to the bimodal climatic pattern.

A much more simplified version of the Camarones lagoon system, with just one state variable (the storage), is presented in Figure 1b. This version simply ignores the dynamics of the boca that hence must be specified as an exogenous variable. Clearly, this is not suited for the simulation of future scenarios, but it is very well suited to calibrate the core hydrological block (the storage dynamics) by using historical data; i.e., a record of opening and closing sequences of the boca.

Notice that, up to now, we referred to the “mono-dimensional” and “bi-dimensional” attributes, thinking of the system concept of “number of state variables” (one for the simplified version, two for the full version). But both versions can also be classified with regard to their spatial dimensions: the scheme can be 0D, when no spatial direction of flow is actually considered (i.e., the lagoon is just a tank, filling and emptying), or 2D, when movement of water in the horizontal x, y directions is considered. The 0D case is managed here as a discrete time dynamic system, and the 2D as a continuous time dynamic system discretized in both time and space.

3. Elements of the Modelled System

3.1. Morphometrics

The purpose of this analysis was to obtain an area–elevation relationship (curve A(y)) and associated storage-elevation relationship (V(y)) for the lagoon. This includes a bathymetric part (under water) and a topographic part (under the open air), varying of course with the particular moment considered. We relied on three types of data: (i) a set of available large-format aerial images captured with a Vexcel UltraCam large digital camera; (ii) sporadic data from our survey in situ with an RTK along the limits of the water-shore and on the bottom (Topcon HiPer V Dual Frequency GNSS Receiver, n = 93 along the shore, n = 111 on the bottom); (iii) a sonar survey of the lagoon shelf (Garmin STRIKER™) with 2 transects and a total distance of 12.6 km, while lagoon width is about 3.3 km and its length 3.9 km, hence a point density of 41 points per km². In addition, an already installed hydrometer recorded the elevation of the water surface, referring to the average sea level reference adopted. The latter has been used as a validation tool, as explained further ahead. Figure 2 shows the three-step process adopted to generate the topo-bathymetric Digital Elevation Model (DEM) of the whole lagoon, including its floodable areas (dry at the moment of the survey).

Basically, altimetric data of the floodable area obtained from a crewed aircraft via photogrammetric reconstruction were integrated with topo-bathymetric data obtained from a survey in situ.

3.1.1. Process

Firstly, we generated the bathymetric surface by interpolation of (orthometric) data obtained with the GNSS-RTK system (n = 111 points, Figure 2a) from a boat; these data, linked to the local geodesic network, were taken on 19 and 20 September 2022 (Datum MAGNA-SIRGAS central and geoidal model GeoCol2004). Secondly, these data were processed in GIS (Geostatistical Analyst-ArcGIS Pro version 3.0 packed) via kriging interpolation. Thirdly, smoothed elevation curves were extracted that allowed the construction of a Terrain Digital Model in the form of a Triangular Irregular Network (TIN) as adherent to the flat lagoon bottom as possible (Figure 2b). Lastly, this TIN was converted into a bathymetric raster and integrated with the points cloud from photogrammetry, so obtaining the desired topo-bathymetric DEM (Figure 2c) with a 1 m spatial resolution. The photogrammetric points cloud is the output of the processing of georeferenced images obtained with a pixel size expressed in ground or Ground Sample Distance (GSD) image of 20 cm/pixel, taken by the governmental agency Fondo de Adaptación during March 2017.

Finally, the surface-elevation curve A(y) and the volume-elevation curve V(y) were built from the topo-bathymetric DEM by using the r.lake.xy module of GRASS_GIS version 7.8.7 [11].

3.1.2. Validation

A checking of our DEM elevations in the dry floodable areas around the lagoon, via a survey carried out with GNSS-RTK (n = 93) during a posterior dry season (13 March 2023), reveals an altimetric average error of 3 cm and a RMSE of 32 cm. As this metric is smaller than 2 times the GSD, we concluded that the absolute precision of our DEM in those areas is satisfactory.

Independently, during the bathymetric campaign, we obtained data of the lagoon depth through a sonar (Garmin Striker) in order to validate the DEM generated by comparing the values in n = 532 points. The raster map of depths from the topo-bathymetric DEM has been calculated again with the r.lake.xy tool of GRASS GIS, taking as a reference the water mask obtained from a Landsat 8 image the same day of the campaign. The average error of 11 cm and the RMSE of 21.2 cm speak of an acceptable consistency.

A final, key validation check was performed by using the reverse-way the A(y) curve and comparing the elevation obtained from it for a given area A, with the elevation provided by our hydrometer [1]. We considered four different situations of the lagoon; for each one, a water mask of the wet area was extracted from Sentinel 2 satellite images (EOS Land Viewer Web Interface for Land/Water band combination: https://eos.com/landviewer (accessed on 5 May 2024): 2 May 2023; 15 March 2023; 31 May 2022 and 20 September 2022). We obtained a satisfactory RMSE of 8.4 cm [1], which supports us in assuming that the obtained model can be reliably used for modelling purposes.

As a final step, we obtained two further relationships needed for modelling purposes: elevation–volume y = y(V) (Figure 3a, Equation (1)); and area–volume A = A(V) (Figure 3b, Equation (2)).

y = 1.32 × 10⁻³⁶V⁵ − 8.37 × 10⁻²⁹V⁴ + 1.98 × 10⁻²¹V³ − 2.22 × 10⁻¹⁴V² + 1.79 × 10⁻⁷V − 7.60 × 10⁻¹(1)

A = 6.11 × 10⁻²⁹V⁵ − 3.82 × 10⁻²¹V⁴ + 8.95 × 10⁻¹⁴V³ − 9.81 × 10⁻⁷V² + 5.46 × 10⁰V + 2.31 × 10⁶(2)

3.2. Mechanism of the Lagoon–Sea Exchange Flows

The idea is assuming that the exchange flow be governed by the height difference (y_L − y_S) between the lagoon and the sea (Figure 4) and the size of the cross-section area A of the boca (given as exogenous data), while other possible factors (e.g., shape of the cross-section, wind, salt concentration gradient, …) are ignored. Of course, the flow would be outgoing (lagoon into the sea) when the difference is positive and vice versa incoming (from the sea into the lagoon) when negative.

When looking for a relationship of the type Q_B = Q_B(y_L, y_S), we encountered an inconsistency as already mentioned in [1]; apparently, it always resulted y_L << y_S, as shown in Table 1, which is physically impossible and most probably due to a mismatching of the reference elevation adopted. Therefore, we looked for a constant φ that, once subtracted from the P.Brisa sea-gauging station data, provided a consistent height difference. We proved that with φ = 0.618 m, all cases, except one, became consistent (see Table 1). Plotting the Q data with the height difference reveals a meaningful relationship, particularly when data are separated into two groups according to the size of the cross-section (not shown for simplicity); this is indeed quite reasonable because the flow must be larger for a larger cross-section.

A simple multiple linear regression linking the flowrate Q_B to the height difference (y_L − y_S*) and to the cross-section area A (Equations (3) and (4)) provides a nice performance (other non-linear relationships can certainly improve this result):

Q_B = Q_B ((y_L − y_S), A)(3)

Q_B = [k₁ (y_L − y_S*) + k₂ A + k₃]⁺(4)

where: Q_B denotes the exchange water flow [m³/s]; the symbol y_S* denotes the corrected sea level datum, that is: y_S* = y_S − φ, with y_S the elevation measured at P. Brisa tide-gauge station, and φ = 0.618 m; [x] + denotes the absolute value of x; and: k₁ = 94.6528; k₂ = 0.6150; k₃ = −6.7008; φ = 0.618 m are calibrated parameters.

Figure 5a shows the matching between the measured flows and the estimated ones through Equation (4), which we consider highly satisfactory, given the rude methods adopted and other uncontrollable factors affecting the outcome (e.g., distance of P. Brisa from the lagoon mouth; presence of wind; etc.). Of course, the relationship is fully valid only within the range of values utilized in its calibration.

To strengthen the quite surprising results obtained, we performed a completely independent verification of the obtained relationship by applying classic hydraulic concepts. We schematized the boca as a short channel where the flow would be governed by the typical Chezy–Manning equation [12], with a length of 100 m and a slope given by the height difference of Table 1 over such length (Figure 4b). With a calibrated Manning coefficient n = 0.040, the consistency of the found values is, once more, surprisingly nice (Figure 5b) (we are aware that, for instance, the flow is not uniform and the actual gradient is lower because of the initial acceleration).

As a conclusion, the relationship (Equation (4)) found by multiple linear regression seems suitable to represent the flow exchange. It adopts as an exogenous variable the sea elevation, while the lagoon elevation is a state variable of the dynamic system, but it also needs the cross-section area A so that a rule to describe its evolution as a function of the other variables should be defined, as discussed in the following paragraph.

3.3. Mechanisms: Dynamics of the Mouth

The challenge here is to achieve an explanation of the dynamic of the boca in terms, first of all, of its process of opening and closing, and then in terms of cross-sectional area. Without this understanding, it is virtually impossible to carry out long-term simulations of future scenarios. Developing a credible theory for this dynamic behavior allows us to switch from scheme (a) to scheme (b) of Figure 1.

From the few observations we obtained during our monitoring period (10 December 2021–30 June 2023), it seems legal to infer (Figure 6) that the boca opens when either the lagoon water surface elevation or the sea elevation passes a threshold ϑ₁ (approximately 0.35 m asl). Then, it keeps open for at least a minimum interval τ_min (about 1.5 months) unless one of the following two conditions occur:

• –. a time superior to the maximum openness duration τ_max elapses (about 4 months);

• –. the lagoon level drops below a minimum threshold ϑ₂ (about −0.24 m asl) and the maximum height difference between the lagoon and the sea is below the threshold σ (about 12 cm), i.e., there is little enough energy for the sea to re-deposit the mouth bar that closes the lagoon.

Implementing this rule provides a quite nice reproduction of the pattern historically observed during our monitoring, as visible in Figure 6a, although the duration does not perfectly coincide with the (supposedly) observed one. No guarantee, however, can be given that this “model” actually behaves reliably well for other periods. Further data will be gathered to explore its validity and possibly to suitably modify it.

The dynamic of the cross-section area A is assumed to adhere to the following rules:

• –. from when it opens, it increases linearly, but quickly (τ_full = 5 days) from 0 to the maximum A_max (this in principle could depend on the strength of the opening process, but this is influenced in turn by human intervention, not recorded up to now, so at the moment it is not possible to establish a relationship); the maximum value recorded is here adopted (A_max = 92 m²);

• –. then it keeps more or less constant for a certain interval τ_keep (about 20 days);

• –. then it diminishes exponentially until the next closure instant.

The results of applying this “theory”, together with the pattern reconstructed (extrapolated) from the few available measurements, is plotted in Figure 6b. It is evident that some additional mechanisms need to be incorporated in order to modulate the maximum area values, but this needs to wait for further observations and data.

3.4. Inflow from Runoff and Direct Precipitation

We assume that the precipitation over the area usually occupied by waters directly flows into the lagoon, even when it shrinks, but not all of it translates into an actual surface water flow because (i) part of it evaporates before it touches the ground; (ii) the available data refer to Camarones station (latitude 11°25′43.75″, longitude 73°3′9.83″, see [1]), and, given the type of meteorological events occurring in the area, certainly cannot be assumed to spread uniformly above the lagoon area. Therefore, a reduction factor μ₁ is needed, and the equation is (Equation (5)):

Q_p [m³/day] = μ₁ p S/1000 [m³/day](5)

where S is the lagoon max surface [m²] (directly observed from aerial photos) and p the daily precipitation at the Camarones gauging station [mm/day].

To determine the flow generated by direct runoff from the local water basin (see [1]), a kind of hydrological rainfall–runoff model is needed. As a first, basic step, we assume that there is no significant delay nor accumulation so that a constant fraction μ₂ of the precipitation over the surface S_R of such a basin (excluding the lagoon itself already accounted for) enters the lagoon; hence, the flow is (Equation (6)):

Q_R [m³/day] = μ₂ p S_R/1000 [m³/day](6)

Clearly, these parameters have to be calibrated, but it is not possible to obtain direct measurements useful to this purpose; the only possibility is to calibrate them through the overall answer of the lagoon model. But this task, as is made evident in what follows, can be undertaken only once the gaps regarding key forcing factors, like the river inflow, are filled.

3.5. Evaporation

The evaporation flow from the free surface of the lagoon is determined as expressed by Equations (7) and (8):

e(t) [mm/day] = e (constant)(7)

E(t) [m³/h] = S(t) e(t)/(24 × 1000)(8)

where E(t) is the evaporation flow from the lagoon in the interval (t, t + 1); S(t) [m²] is the lagoon surface at time t, as evaporation is directly linked to the actual water surface present; and e(t) is the daily evaporation rate [mm/day] depending in principle on the specific climatic conditions of that day t. To determine the evaporation rate e(t), several formulas can be applied; however, their use is limited by the meteo-climatic data available. An official station is present in the area of our lagoon but has no data for our period. The temperature excursion is however quite limited in the area, with about 3 °C degrees only between the “cold” and “hot” seasons (see Figure 7); therefore, evaporation does not vary markedly. This pattern can be considered representative of our lagoon, although the station is quite far away (20 km) and so the absolute values of temperature (and evaporation) can be different, owing to local conditions (mainly winds, for which we have no data). For this reason, they have not been used to determine the actual evaporation in the model. As a preliminary approximation, we adopted instead a constant value of 10 mm/day, consistent with previous works [13].

4. Modelling of the Lagoon

4.1. Characteristics

4.1.1. 0D Spatial Modelling

The “0D” denotes here, as already said, that no spatial dimension is considered: the lagoon is seen as a tank. Additionally, thanks to the simplifying position of Figure 1b, the system has only one state variable: the storage V(t). Therefore, modelling the dynamics of the boca is bypassed.

Such a model of course has limitations. The main one is that it cannot see the distribution of flows and depths across the water body. However, the almost omnipresent wind or breeze and the very low depth (mostly less than 1 m) ensure that the complete mixed schematization is acceptable, at least vertically. Possibly, the freshwater input from the river and the incoming exchange flow from the sea are more complex spatially. But, concerning the latter, we observed only once a situation where the mouth was split into two parts with both an incoming and an outgoing flow present in the two separate sections, so that situation certainly is more an exception than the rule.

Although the physical system lives in the continuous time t, a discrete time schematization with a quite short time step Δ (an hour) is very likely to be acceptable, and this is what we adopted (Equation (9)):

V(k + 1) = [V(k) + Q_R(k) + Q_p(k) + Q_L(k) − Q_B(k) − E(k)](9)

where k denotes the time step, and all the terms have already been introduced, but are recalled here (all expressed as m³/hour]:

V(k): lagoon water volume; Q_R(k): freshwater input from the Tomarrazón-Camarones River (exogenous input); Q_p(k): freshwater input (exogenous input) from direct precipitation on the maximum lagoon surface S; Q_L(k): freshwater input (exogenous input) from the runoff into the lagoon from its local catchment; Q_B(k) = Q_B(y_L(V(k)), y_S(k), A(k)): the lagoon–sea exchange flow (positive when outgoing), completely determined after specifying the lagoon surface elevation y_L, which is a function of the volume V(k) at that time through the corresponding inverse morphometric relationship so that y_L = y_L(V); the sea elevation y_S(k) and the cross-section area A(k) of the boca are once again exogenous inputs; E(k) = S(V(k)) e(k): evaporation loss from the lagoon depending on the area S(k) of its surface given for any volume V(k) by the corresponding indirect morphometric relationship S = S(V); and e(k): again an exogenous input.

In this schematization, all the exogenous inputs, by definition, are completely known given the time k.

The approximation of this formulation does not lie in the fact of it being discrete time; rather, it lies in the way the terms are computed. Indeed, the approach adopted in Equation (9) is known as the Euler method [14]; that is, given a continuous time system of the type given by Equation (10) (i.e., that describing the lagoon physical system):

dx(t)/dt = f(t,x(t))(10)

the discrete time version is approximated as (Equation (11)):

x(k + 1) = x(k) + f(t,x(k)) Δ(11)

where dx/dt is the total derivative of time variable x with respect to the continuous time t; f(t,x(t)) is the transition function; and x(t) is the state vector of the system. Clearly, this is in essence a rude approximation, as the continuous state may indeed change during the interval (t; t + 1), while in the discrete version, the new discrete state vector x(k + 1) is computed by disregarding such further changes as if x(t) remained fixed at time t.

Many other approximations exist still within explicit schemes, like that of the Runge–Kutta method [14], or jumping to implicit schemes [14], which are more precise but heavier computationally. The Euler method seems, however, sufficiently precise in our case.

4.1.2. 2D Spatial Modelling

For water bodies such as coastal lagoons, where the vertical dimension is much smaller than the horizontal dimensions and the depth is relatively uniform, a 2D model, rather than a full 3D model is a good cost-effective solution [15]. This kind of model represents hydrodynamic behavior well by simulating the distribution of flows of water across the water body without the need to describe complex vertical stratifications [16,17]. Two-dimensional models have proven adequate in shallow water flow situations, significantly reducing computational time with respect to the full 3D equation [18].

In principle, they are able to provide extensive information: time series of relevant variables (velocity, depth, water quality), velocity fields and current patterns, and sediment transport patterns, including the representation of changes in topo-bathymetry due to erosion and sedimentation processes [19,20,21,22]. Such models can potentially simulate also the opening and closing of the mouth of the lagoon, like in the case of our Navío-Quebrado Lagoon. And, finally, they allow the simulation and analysis of the impact of extreme climate events related to climate change.

Therefore, a 2D model conceptually appears more suitable and attractive than a 0D model. However, the accuracy of numerical models is heavily reliant on the presence and quality of the input data. Indeed, 2D hydrodynamic models for coastal lagoons require a series of detailed and specific input data to configure, implement, and produce accurate results. Among these inputs, topography and bathymetry (topo-bathymetry) are fundamental since detailed information is required on the conformation of the lagoon bed, the marine area, and the surrounding terrain.

Taking a step back, a number of issues hence become relevant referring to our exercise. In particular, a first set of issues refers to the ability to reproduce the historical behavior (observed period):

Is it technically possible to model a time-varying mouth according to a measured time series ? (i.e., a mouth geometry treated as a known exogenous input; actually, an interpolated time series obtained from few, irregularly scattered measures of the mouth width):

• ▪ Models could dynamically adjust the size, shape, cross-sectional area, and other characteristics of the lagoon mouth based on input time series data as handled by the XBeach model [20,23]. But this exercise involves significant challenges and complexities [23,24].

• ▪ One simple way to deal with this is to adopt a (topo-bathymetric) DEM including the largest mouth and then run a number of simulations, each starting with the final state of the previous one and including a border condition with a progressively larger width on top of the DEM mouth while the mouth is widening and vice versa when it is narrowing. A more elegant solution is to embed directly in the model code (usually not accessible) a time-varying border condition.

The answer is hence “positive”, although not straightforward.

Is the 2D model able to reproduce the observed historical behavior (lagoon water level time series) while adopting the same input information utilized in the 0D exercise ? (that is: river inflow, sea level, precipitation time series, mouth width and depth, and constant evaporation rate).

• ▪ This implies that the model determines automatically the exchange flows across the mouth, given the sea level and the current lagoon water level. This seems quite feasible according to the very nature of the 2D model approach: at any time, in any cell of the discretization modelling grid, a gradient is computed and a corresponding flow determined (in x,y directions).

• ▪ Winds are crucial in the phenomena of coastal barrier overflow and favor water exchange between the lagoon and the sea (ebb and flow), so for the case at hand, where winds are scarcely known, the 2D model output can be negatively affected.

Therefore, the answer is in principle positive, but assessing the actual output is another story and deserves a check.

Is it possible to incorporate within the 2D model a formalized mechanism that determines mouth opening and closure ? This is by far a more difficult requirement:

• ▪ In the 0D model, we solved this by including a very basic rule stating, essentially, that on top of a certain lagoon level (threshold), the mouth would open, while it would close below another threshold.

• ▪ In the 2D model, the mechanisms of bed morphology change can be embedded in the “physics” of the model itself, so that it autonomously and dynamically adjusts the mouth geometry in response to changing conditions, simulating sediment transport and morphological changes [24]. Several studies on lagoon barrier overtopping events have employed 2D hydrodynamic models to reconstruct the entire process and variability of water flows, including dam formation, breach flow hydrographs, and downstream flooding. This demonstrates the models’ abilities to capture complex morpho-dynamic and hydro-sedimentological processes [19,25].

• ▪ However, obtaining the data required by these sophisticated models (e.g., XBeach) is expensive and logistically challenging, limiting the model’s applicability in practical situations where data are scarce and generally unreliable [26,27]. In addition, sophisticated algorithms and robust numerical solutions [28] are required, as well as specialized knowledge and experience to configure and run the model, interpret its results [23], and occasionally access source code to perform the necessary adjustments and configurations to the computational algorithm. This level of expertise may not always be available, especially in resource-limited regions, further complicating the practical application of these models.

• ▪ To model in particular the closure of the mouth, essential input data detailing exogenous forces inside the lagoon as well as in the marine domain (such as tide levels, waves, solids concentration, and sediment transport by marine currents) are fundamental [20,25]. Without detailed, extensive, high-quality data input, hardly obtainable in general, as well as sufficient computational resources, the feasibility and accuracy of these types of models are significantly limited.

We can therefore conclude that, in particular, for the opening stage, the mechanism can possibly work, but it hardly could incorporate the fact that usually opening is supported by human intervention (the 0D model does not see this explicitly, but it does by appropriately setting the threshold that already incorporates this human behavioral mechanism). For the closing stage, then, the exercise is by far more difficult because the model should incorporate the dynamic mass balance of the sea shore, including the mass transfer along the coast, as it is this input that, when waves prevail on the lagoon emptying process, build up the mouth bar.

Hence, the answer may again be positive in principle, but in practice, it becomes quite challenging.

Finally, another big question deserves attention, this time concerning the future:

Is it possible to use the 2D model to simulate the lagoon behavior according to (several) long-time horizons (and associated hydrological and sea level inputs) ?

• ▪ Here, the key issue concerns the computing and processing speed.

• ▪ Dimensionality reduction can be a way out (as explored, for instance, by [29]). However, especially concerning the simulation of coastal barrier breaches, significant computational resources are required due to the need for high-resolution numerical meshes and the complexity of the physical processes involved [19,20,21,22]. Mesh refinement in the breach area is generally required, implying an increase in computational cost, as it necessitates decreasing the time step to maintain numerical stability. This can be a limiting factor, especially in the case of long-term simulations.

• ▪ Most probably, the answer becomes negative when dealing with a personal computer, while it may be positive when thinking of cloud computing [30].

To explore the practical feasibility of adopting a 2D model, we proceeded with the Environmental Fluid Dynamics Code (EFDC) (more precisely, the EFDC+ Explorer software package version 11.3) because it is a coupled model that solves hydrodynamics, contaminant transport, and water quality kinetics in surface water systems, all in one integrated code, eliminating the need for external coupling between the hydrodynamics and transport modules, with a combination of finite volume and finite difference techniques on a Staggered Cell type grid to solve fluid momentum equations [31,32].

The EFDC has emerged as one of the most extensively applied numerical models for hydrodynamic and water quality simulations. The EFDC+ Explorer model, used in this work, stands out for its efficiency in data pre-processing, high-performance computing capabilities, numerical stability, and advanced post-processing functionalities. Furthermore, its adaptability has been validated through successful implementations across various case studies worldwide, demonstrating its effectiveness in modelling complex lagoon and estuarine systems [33,34,35].

In this study, EFDC+ was used to simulate hydrodynamics coupled with the temperature module to consider the atmospheric conditions, mainly evaporation and precipitation. The EFDC+ hydrodynamic module solves the three-dimensional, vertically hydrostatic primitive equations of motion for turbulent flow in a curvilinear orthogonal coordinate system. It incorporates a second-order turbulence closure model, which relates turbulence viscosity and diffusivity to turbulence intensity and length scale. The model employs a state equation that correlates density with pressure, salinity, temperature, and suspended sediment concentration, enabling comprehensive simulation of complex hydrodynamic processes [35,36].

4.2. Replicating the Historical Dataset

The simulation period was established from 20 January 2022 to 30 June 2023, covering several wet and dry periods, as well as incorporating the temporal variability that characterizes the dynamics of the lagoon due to the influence of multiple factors such as tides, winds, river flows, and mouth opening/closing.

4.2.1. The 0D Case (Excel^® Model)

In order to proceed with a simulation, once the model is implemented in a software (our 0D was made in Excel^®), the harshest difficulty found is that the model needs a “continuous” time series of the river freshwater input; more precisely, it needs it on a discrete hourly time basis, but with no missing data. The reality, unfortunately, is that we can hardly count with few measures (actually, just estimations) of flowrate (gauging exercises). For the major part of the remaining time, instead, we only have one level measurement a day at an irregular time moment (usually between 6 am and 9 am, but sometimes even in the afternoon) at the chosen gauging station of P. Troncal (and that, closer to the lagoon, of P. Viejo; see map in [1]). But, in both cases, there is no information about what happened the remaining 23 h of that day. Finally, this time series has “holes” because some days (indeed, not many), no measurement was taken at all. A first easy step adopted was to transform the level measurements into flowrate data thanks to the stage–discharge relationship obtained (see [1]).

Then, to “fill the rest of the day”, we assumed that the flowrate would remain more or less the same by repeating the same measured or estimated value multiplied by a factor λ to be set by trial and error; this factor would in principle be smaller than unity in case the recorded value was a flood (because it is unlikely that the flood lasts several hours, as the river basin is just about 500 km² and events are rapid), while it would be greater than unity in the case the measured value was a low flow case. In the application, only one value of λ was adopted for the whole time series. For days with no measurement, we then adopted the following reconstruction criterion: if the no-data period is one or two days, the reconstructed value is the average of the previous and following measured values; when the period is longer, the flow is null. This is consistent with the policy agreed with our partners that each time a flow was visible in the river, a measurement had to be taken, which implies that in all periods of no flow, no measurement is taken. It is clear, however, that the reverse statement is not necessarily true; this means that structurally it is expected that the river inflow be underestimated.

The sea level is measured on an hourly basis so that no problem arises except the few days where a datum is missing, in which case it is interpolated between the closest data. The level of the lagoon was measured usually just twice a day (in the morning and in the afternoon), but this is not a problem because we use it only to compare the model answer with reality and assess its reliability, while the lagoon–sea water exchange module uses the state of the system (storage) calculated by the model itself.

In the simplified version of the system scheme (Figure 1b), there is a need to assign a time series of the status of the mouth (opening/closure and cross-section area); in the periods of “boca cerrada” (closed mouth), it is easy because A = 0 and no exchange can occur; in the other periods, we interpolated linearly the few measurements available.

For the direct precipitation inflow and runoff from local catchment and the evaporation rate, we adopted the recorded time series of relevant meteo-climatic variables (temperature and so on) from a nearby station (see [1]).

4.2.2. The 2D Case (EFDC+ Model)

We discretized the domain area with curvilinear orthogonal horizontal coordinates. The grid was generated using the Curvilinear Grid Generator for EFDC (Grid+). The model has 203 elements in the X direction and 186 elements in the Y direction, with ΔX varying from 29.2 to 44.9 m and ΔY varying from 13.4 to 38.5 m, for a total of 25,564 active cells (Figure 8). The topo-bathymetry of the study area obtained as described in Section 3.1 was processed into this high-resolution grid.

Input conditions employed for hydrodynamic simulation are shown in Table 2.

• Atmospheric Forcing Data

A constant evaporation rate of 0.01 m/day was applied throughout the simulation time, consistently with the 0D exercise. Wind forcing was based on some measured wind speed and direction data available from 1 February to 7 September 2022 and kept constant throughout the period.

• Initial Conditions

The water surface was set at the level reported from the field measurements recorded. This implied the assumption of a horizontal water surface and no movement (v = 0) everywhere.

• Boundary Conditions

The two main boundary conditions used were, consistent with what had been performed for the 0D, a Dirichlet-type boundary for the discharge of the Tomarrazón River and a Neumann-type boundary for the connection with the sea, where sea level variation was imposed. The river’s boundary condition included the freshwater inflow from the main river basin Tomarrazón-Camarones and the runoff from the local watershed. Figure 8 shows the location of the flow river boundary condition. The boundary condition at the sea interface needs further clarification. To represent the closing and opening of the lagoon mouth, we use the time-varying field for bathymetry option available in EFDC+. This option allows us to specify the time (or space) varying forcing factors for each cell over the whole model domain. By using this capability, EFDC+ can receive time-varying topo-bathymetric updates within its computational domain from external sources. The area of the sand bar modified to reproduce the opening and closing of the mouth is demarcated with a red polygon in Figure 8. First, the topographic variation in the mouth zone was determined using the field-measured mouth surface area (m²). To do this, a fixed width was established to obtain an approximated depth; then, with the measured lagoon free surface level data, this depth was converted into elevation above sea level. All the cells in that zone were modified by assigning a new bottom elevation each hour. Thus, a file with 527 different time-varying fields for bathymetry was created to reproduce the behavior and variation of the mouth throughout the simulation.

5. Results

5.1. Results for 0D Model

Suitable values of the parameters involved (reduction factor for direct precipitation μ₁ = 0.5 and local runoff μ₂ = 0.2; modulation factor of the river inflow λ = 1.2) were assigned (not through a formal calibration, but just via an approximate trial and error procedure), and the main results obtained are presented in Figure 9 and commented here below.

The general behavior of the model is coherent with the system physics: when there is a significant river inflow (flood: dots of the Figure 9a) or precipitation (histogram “Q rain”, Figure 9a), the level of the lagoon water surface (both measured, “y Lag” dots, and simulated, “y Lag simulated”, in part b of the figure) suddenly rises and then soon starts dropping with a first order behavior (exponential decay) because of the evaporation rate, possibly drying out the water body. The same occurs when the mouth is closed (code C in the horizontal line above the time line) and there is no flood or precipitation input the lagoon: the level drops. When the mouth is open (O) or semi-open (S), the level fluctuates according basically to the 24 h moving average (“Sea mov average”, bottom figure b) of the sea level (“Sea”, bottom figure).

The model represents the measured elevations (pink dots) quite well, although with a general sub-estimation, clearly due to missing freshwater inputs. In some occasions, particularly in the beginning period (t ≤ 3000 h, when the hydrometer and piezometer still suffered from some structural defect), the model overestimates the values (but they still lie between the hydrometer and piezometer measurements, not reported in the figure to keep it readable: as explained in [1], we finally adopted as a measure of the lagoon level y_L, a combination of both according to a specified criterion).

The tide is very “nervous”, with an amplitude of oscillations significantly larger than that of the lagoon; this is more linked to the 24 h moving average. In addition, it is evident that between times 7000 and 9300 approximately, the oscillations of the lagoon elevation are much wider, but so is the tide (actually, its 24 h moving average value). In general, it must be considered that the lagoon elevation measurements only occur twice a day; hence, with such a low frequency sampling, it is very unlikely that the maxima and minima are indeed captured (as shown in the inset, Figure 9b).

The model represents measurements quite well, although there is a general sub-estimation due to lack of river inflow inputs (particularly between times 8300 and 11,000, inset in Figure 9b); there also is a sporadic small overestimation (this is still lying between the hydrometer and piezometer data boundaries not reported in the figure to keep it readable; notice that in the initial period until t = 3000, the two instruments needed some adjustments and were therefore less reliable).

Further analysis was conducted to explore the behavior, namely:

• ▪ The system appears to be very sensitive to evaporation, as shown in Figure 10. This means that more attention should be paid to the determination of the evaporation rate.

• ▪ Finally, the cross-section area of the boca determines the amplitude of the water surface oscillations (and the associated flowrates) but not the general pattern of the elevation, as shown in Figure 11, for an increased value of 50% in the period 7000–9000 (as an example).

In Table 3, the water balance is reported for different time intervals. It is apparent that, by far, the major water input comes from the sea exchange; the river input covers the second majority, while evaporation is approximately double the precipitation input.

It is worth noting that the mass balance shows a small error. This may be due to (i) mutual inconsistency amongst the morphological relationships (the numerical V(y) is not really the integral of the S(y) as it should be) and numerical approximations adopted for their functional forms; (ii) errors introduced by the Euler integration method; (iii) approximations in the numerical reckonings (as it works at hourly time steps, in a considered period of about one and a half years, there are more than 13,000 steps). In a further stage of this research, we may investigate this further because if the model is to be used for longer horizon simulations, errors might sum up and affect the results.

When the dynamics of the boca are incorporated (Section 3.3), results are a bit worse but still acceptable, as shown in Figure 12. It is noteworthy that the main difference is an underestimation after the closure of the boca, particularly around time 1000–2500; this is due to the too-early closure that impedes a recharge of the lagoon by sea water. Additionally, the oscillations when the boca is open are in general a bit exaggerated. This latter problem can be solved by modulating the A_max parameter as already noted (see Figure 10); the former problem may improve by adopting a less gradual reduction of the section after reaching the maximum opening (i.e., not the linear pattern adopted, but for instance an exponential one), but the main point is to capture the right closure, and with our dataset—where the initial data are lacking (particularly the river inflow and the sea level)—this is not possible. This can be addressed once the observation period will be further extended. On the contrary, around time 9500–10,500, the full model (incorporating the dynamics of the boca) seems better than the partial one (adopting the extrapolated area based on measurements); this is due to the opposite fact: the closure of the boca occurs later, and so an increased volume from the sea can fill the lagoon.

A final observation is that the full model answer perfectly coincides with the partial one (without the dynamics of the boca) in the periods right after the dry out (time 3000 and 11,000): this is where the system memory is nullified and a new cycle begins.

5.2. Results for 2D Model

The comparative plot of the time series of water levels in the Navio-Quebrado lagoon (Figure 13), which includes measured data and results from 0D and 2D models, highlights the key differences between both modelling approaches. The observed fluctuations in water levels are clearly influenced by the temporary opening and closing of the mouth that connects the lagoon to the sea, as well as by the seasonal contributions of the Tomarrazón River. These factors constitute the most important external drivers affecting the dynamics of the lagoon, as they regulate both the water level peaks during periods of highest flow and the decreasing levels when the lagoon is isolated, and evaporation becomes the main mechanism of water loss.

Curiously, the 2D model shows an adequate fit at moments when levels in the lagoon suffer gradual changes or when there are no large fluctuations induced by intense external forcing; that is, during periods without sudden changes. However, it tends to overestimate some peaks and has less detailed responses with smaller variations compared to measured data and the 0D model.

Again, as expected, there is an underestimation of the water level towards the end, most probably due to “escaped floods”.

Four fitness indicators were adopted to assess the goodness of fit of models by comparing simulated values and in situ measurements: root mean squared error (RMSE), where RMSE yields only positive values (a zero value (rarely obtained) indicates a perfect fit to the data, values close to zero represent an ideal fit, and values far from zero indicate a deficient fit to the data [37]):

(12)$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}\sum _{\mathrm{n}}^1(\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}-\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}{)}^2}$

MAE is the average of the absolute values of all differences between simulated and observed values:

(13)$\mathrm{MAE}={\displaystyle \frac{1}{\mathrm{n}}}{\sum }_{\mathrm{n}}^1\left|\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}-\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}\right|$

The Kling–Gupta Efficiency (KGE) [38] is an evaluation criterion that can be decomposed into the contribution of mean, variance, and correlation to model performance (Equation (5)). KGE ranges in (–∞, +1]. Essentially, the closer to +1, the more accurate the model is:

(14)$\mathrm{K}\mathrm{G}\mathrm{E}=1-\sqrt{{\left(r-1\right)}^2+{\left({\displaystyle \frac{\overline{P}}{\overline{O}}}-1\right)}^2+{\left({\displaystyle \frac{\mathrm{S}\mathrm{D}\mathrm{P}}{\mathrm{S}\mathrm{D}\mathrm{O}}}-1\right)}^2}$

Finally, the Nash–Sutcliffe coefficient of efficiency (NSE) [39] varies in (–∞, +1], with higher values indicating better agreement, and NSE = 1 is the optimum value. The NSE ≤ 0 indicates unsatisfactory performance, and NSE closer to 1 is considered a better model performance:

(15)$\mathrm{N}\mathrm{S}\mathrm{E}=1-{\displaystyle \frac{\sum _{i=1}^N{\left(P_i-O_i\right)}^2}{\sum _{i=1}^N{\left(O_i-\overline{O}\right)}^2}}$

where: $O_i$ and $P_i$ represent the measured and calculated data, respectively, $\overline{O}$ is the mean of observed data, $\overline{P}$ is the mean of simulated data, SDP and SDO are the simulated and observed standard deviations, r is the linear correlation coefficient between $P_i$ and $O_i$, and $N$ is the number of the analyzed values.

This study did not explicitly evaluate different grid resolutions as part of the modelling workflow. However, we prioritized ensuring that the numerical representation of the lagoon accurately captured its key hydrodynamic characteristics. Specifically, we focused on precisely representing the dimensions of the lagoon’s connection to the sea, a critical factor controlling tidal exchange and salinity transport. We also paid close attention to the Tomarrazón River delta’s branched structure, where distributary channels influence the spatial distribution of flows. To balance model accuracy and computational efficiency, we used the curvilinear grid with refined zones in these critical areas. This refinement allowed for a more detailed resolution of hydrodynamic processes in the “La Boca” region and the Tomarrazón delta, while maintaining a coarser resolution in the lagoon’s less dynamic areas to optimize computational performance. The time step was selected by assessing numerical stability using the Courant–Friedrichs–Lewy (CFL) number, ensuring model stability throughout the computational domain.

6. Discussion and Conclusions

6.1. Discussion

An interpretation of the statistical results in the context of hydrodynamic modelling of the lagoon is presented in Table 4.

Evidently, the 0D model shows lower values in the error metrics (MAE and RMSE), indicating that, in absolute terms, it has a better fit to the measured data than the EFDC model. This apparent improvement should be interpreted with caution, since 0D models usually present less variability in their results by smoothing out extreme fluctuations, which can artificially reduce these errors without necessarily better capturing the real dynamics of the system. However, the NSE metric—which evaluates the capacity of the models to reproduce the variability of the observed data—is, like the previous indicators, higher in the 0D model, indicating a good overall fit to the measured levels. Only the KGE value is slightly higher in the EFDC model, suggesting that this model has a better balance in terms of the combination of correlation, bias, and variability, which is key in hydrodynamic models influenced by external forcing.

From the point of view of lagoon dynamics, a more critical point of view may be taken, as these results may hide the limitations and strengths of each model. The better performance of the 0D model in terms of the three indicators (MAE, RMSE, and NSE) clearly says that, for the case at hand and the specific dataset available, the 0D model outperforms the 2D; quite a surprising result, particularly considering that the dataset includes quite abrupt oscillations, in particular of both the sea tide and the river input. In general, indeed, 2D models have a better representation of the spatial and temporal variability of levels, which may be crucial in coastal lagoons where tide, wind, and river contributions generate complex dynamics. But this is not the case for our case study.

From a more general point of view, both models replicate quite satisfactorily the patterns of the measured variables when the dynamics of the boca are considered as an exogenous input (i.e., known from observations).

This result may appear obvious as the system and the model are extremely simple; however, reconstructing the inputs is not a trivial exercise. In particular, the lagoon–sea exchange process is quite delicate and, a priori, one might think that too many factors come into play (say: wind, waves, changing bathymetric morphology, salinity). In spite of this reality, as a matter of fact, the simplified relationship found well described the observed behavior. The water balance points out, in turn, the role of the different processes in the water balance.

A slightly worse result is obtained when the mouth dynamics are explicitly modelled in the 0D model, but it can certainly be improved and can be considered sufficient for exploring future scenarios.

The comparison between 0D and 2D hydrodynamic models unveils the strength of the 0D model: it performs slightly better than the 2D and is evidently much simpler. Possible reasons for this are: (i) although the area of the mouth between the lagoon and sea has been represented as varying, following the actual system information (exogenous input), its width has been assumed constant, while it should have been modulated as well; and (ii) the wind was assumed constant throughout the whole period because of lack of detailed information, which it is certainly not, and this may increment the actual evaporation.

Moreover, the 0D is able to replicate the dynamics of the mouth, although not perfectly, while for the 2D model, this predictive exercise has not been undertaken and would certainly prove much more difficult (notice that we included the mouth dynamic—as an exogenous input—to replicate the historical dataset, but we did not conduct the very prediction of the mouth dynamic; that instead was performed for the 0D, as explained in Section 3.3). Finally, even if this additional element could (theoretically) be incorporated, simulating very long time series with the 2D model would be a challenge owing to the computational burden.

Importantly, the 0D model enabled us to ascertain that our comprehension of the physical behavior was correct and actually added several pieces of information. For instance, we ascertained that the tide, although determining the major water exchange, does not affect the general time pattern of the surface elevation; a variation of the mouth cross-section area influences the amplitude of the oscillations, while it is its moving average that introduces some real forcing effect, but of course only when the mouth is open. The evaporation flow, although by far smaller than the sea exchange flow and the river inflow, is the governing factor when the mouth is closed and outside of the rainy season. Another finding, although to be ascertained with further data, is that the effective input from the precipitation (directly to the lagoon or indirectly through runoff from its catchment) is much lower than what could be expected from the recorded data, contrary to the initial intuition; the main reason for that probably is the nonspatial representativeness of the data. Finally, or most of all, it is essential to determine the full river hydrogram, otherwise an underestimation of the lagoon water surface elevation is unavoidable.

However, to face this latter problem, an indirect approach can be adopted: set up a hydrological rainfall–runoff model of the river basin based on satellite precipitation data (now available) and reconstruct the flowrate based on precipitation data. Then, the problem would be how to calibrate such a model, as there exists no official systematic river flow gauging station in our basin. In principle, the low granularity flowrate data available from our own gauging station can serve this purpose, but adding more data up to the whole hourly time series of river flow data would offer a better resource.

But an even more empirical approach can be explored by observing that the quite primitive monitoring system set up, together with the 0D hydrological model developed, can be used to implement a reconstructor of the river flowrate hydrogram: the idea is that, from the variations of the lagoon level not explained by the causal variables that can be easily measured (precipitation and evaporation and sea exchange), one can infer what river flow must have occurred for the model to fit the discrete measurements of the lagoon water surface elevation. In other words, the lagoon monitoring system plus the model constitutes a gauging system for the river (Figure 14).

The most basic approach would be to write Equation (9) recursively, starting from a measured value y_L*(k) of the elevation at discrete time k, until the following measurement y_L*(k + Δ) after Δ steps, so obtaining an articulated non-linear function of the type:

y_L*(k + Δ) = g(k,y_L*(k); u(k), u(k + 1),…, u(k + Δ − 1))(16)

where the function g(•,…,•) is a combination of the recursive use of Equation (1) and the inverse morphological relationship V = V(y), while [u(k), u(k + 1),…, u(k + Δ − 1)) ] is the set of the unknown values sought for, u(k) being the generic river inflow, while all other exogenous inputs are supposedly known and hence embedded in the g(•,…,•) function. This is a problem with one non-linear equation and Δ unknowns, which by construction has no unique solution (if any); this is not at all surprising, as we are here trying to estimate a high granularity input (the river hydrogram at a hourly level) from low granularity observations (the lagoon elevation measured discretely just twice a day), taking advantage of the fact that these variables are interlinked by the continuity equation. Adding the information related to the following measurement (or to all available measurements) does not improve the situation because a new set of unknowns is also added.

However, if the fitting is assessed not on a single value y_L*(k) at each real-world time t, but on the whole set of measurements taken from the last good fitting model time k*(t) until current time t (e.g., with k*(t) located at the end of a significant period of no river inflows), and by introducing some reasonable, practical hypothesis about the shape of the hydrogram portion, there may be a chance that a reasonable solution be found. In any case, that would be a much better estimation than…nothing. The discontinuous measurements at the gauging station may greatly help, although they are not essential in principle. Rigid schemes like the Kalman filter [40] (originally conceived to reconstruct the state, but likely to be extended to estimate model parameters—e.g., [41]—and hence even its inputs, which can be seen as special time-varying parameters) cannot be applied here owing to the irregular and discontinuous measurement protocol adopted for the lagoon (and even more for the river).

This is a new open challenge to explore further.

6.2. Conclusions

The exercise we conducted allowed us to set up a simulation model of the hydrological behavior of the Camarones coastal lagoon in two versions—0D and 2D—which proved to respond in a quite sensible and reliable manner (the MAE, RMSE, NSE, KGE indicators are indeed quite satisfactory). More importantly, the mere model development helped us to reveal some peculiar behaviors of the system. This exercise can therefore be considered successful.

It is perhaps worth summarizing the main concepts explored to set up the different components:

• Hydrology: the inputs due to direct precipitation on the lagoon surface and to direct runoff from the local basin are estimated in a rough manner by assuming an instantaneous relationship and by assigning by trial and error a reasonable value to the parameters introduced so that the hydrological inputs be as far as possible consistent with the general behavior of the system when represented by the 0D model. The main input due to the river inflow is measured indirectly only once a day and transformed into a flowrate by a stage–discharge relationship calibrated on the few actual measurements of discharge carried out (which are indeed based on a rough transit velocity measurement). The main validation of this whole component is found in the observation of the consistency of the model behavior with all measurements, as illustrated in Figure 10.

• DEM: (i) Comparison of DEM elevation and field (GNSS-RTK) measurements along the shores in a particular situation when a field survey was carried out; (ii) comparison of bathymetric DEM elevations with the sonar data obtained in another field (boat) survey (532 points); (iii) comparison of the area provided by the relationship A(y) for the elevation y measured at our hydrometer in each of four given situations where a Sentinel satellite image was available, from which the wet area was determined.

• Lagoon–sea exchange: The flowrate is determined by a bi-dimensional relationship linking it to the height difference of the lagoon–sea and the mouth cross-section area. This relationship has been calibrated on a set of few measurements of the area and velocity at the mouth. The validity of this approach is corroborated by the (surprisingly) nice fitting of data to a smooth interpolating curve. In addition, the so-found relationship has been independently validated by adopting a completely different criterion based on the standard hydraulic relationship (Manning’s equation) and observing that very consistent results are obtained.

• Mouth dynamics: This process has been split into the process of opening-closing and the process of changing, accordingly, the area of the mouth. The “algorithm” found is just a formalization of what can be observed from our data. A validation on additional, independent data should be carried out.

The 0D model has the potential to be used with future climatic scenarios to assess, for instance, the frequency of events of dry out and their duration, the frequency of opening and closure events, or the duration of the lagoon–sea exchange phase. All these statistics can be very useful to predict ecological consequences hampering both the natural value of the system and its provisioning ecosystem service.

On the other hand, it is evident that the 2D model can be key when simulating water quality issues, as horizontal mixing in a shallow lagoon is a complex process that the 0D model cannot reproduce at all. Indeed, 0D models present significant limitations because they treat the system as if it were homogeneous, leading to sometimes underestimating the magnitude of response to extreme events and being temporally delayed under changing conditions, as reported by [42]. These characteristics make 0D models unsuitable for simulating phenomena that require high spatial resolution, especially when interactions with the environment vary significantly. Additionally, the flexibility offered by 2D models to handle complex geometries, particularly using curvilinear or nested grids, allows for higher resolution in areas of interest, improving the accuracy in the simulation of water accumulation and redistribution zones [43].

In this context, the integration of 0D and 2D models has been proposed as a hybrid approach to take advantage of the benefits of both methods. According to [44], this combination can provide a more balanced representation of coastal systems, allowing the use of 2D models to accurately simulate extreme events, while 0D models can be employed for more computationally efficient long-term simulations. This strategy not only maximizes the accuracy in resolving the most complex hydrodynamic processes but also optimizes computation time, which is especially relevant in complex environments such as coastal lagoons.

Before our model can actually be used to explore future climate change scenarios, even if disregarding the mentioned (very important) effects of sediment transport, several improvements have to be made: Evaporation figures deserve special attention, as the model is particularly sensitive to them, and actual measurements would be particularly helpful; direct and indirect precipitation inputs need to be better estimated, taking into account the spatial patterns of the climatic variables and actually calibrating the key parameters μ₁ and μ₂ and λ; the Euler computational integration scheme should be improved; the dynamics of the mouth needs to be better represented and consolidated (parameters τ_min,τ_max, ϑ₁, ϑ₂; A_max, τ_keep), including the dependence on anthropogenic interventions. This latter point is particularly important because a relatively simple artificial intervention to keep the mouth opened during the no-flood periods might greatly increase the productivity of the lagoon as the dry down process, and the associated over-salinization with death of all fishery inside, would be avoided. Clearly, this is again an exercise that would see a beneficial coupled use of the 0D model to assess the frequency and duration of events and the 2D model to specifically assess the model answer in a typical cycle in terms of water quality change. The findings would then have to be evaluated in terms of ecological behavior and impacts on the broader coastal lagoon system.

In terms of data acquisition, most importantly, a better method to measure the river input has to be put in place. Certainly, a more automatized gauging station [45,46] would solve our problem, but as pointed out in [1], unfortunately, the local socio-cultural context ensures that such a system would have a short life. Satellite telemetry methods are not applicable here owing to the small width of the river and the dense riparian vegetation. For the purpose of modelling the lagoon, an alternative, promising approach yet to be explored is the reconstruction scheme presented above.

Author Contributions

Conceptualization, A.G.C.N.; methodology, A.G.C.N., F.T.-B. and J.I.P.-M.; topo-bathymetry: J.R.E.V., J.M.F.A. and J.I.P.-M.; software (specific Excel sheets to analyze and plot data): A.G.C.N., J.I.P.-M. and F.T.-B. for the EDFC+ software; field surveys and measurements: J.I.P.-M., A.G.C.N., J.R.E.V. and J.M.F.A.; writing—original draft preparation: A.G.C.N.; writing—review and editing: J.I.P.-M., F.T.-B., R.-A.R.F. and J.R.E.V.; graphics, J.I.P.-M., A.G.C.N. and J.R.E.V.; supervision: A.G.C.N.; project administration: J.I.P.-M. and F.T.-B.; funding acquisition: J.I.P.-M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data generated within this paper are available under direct request to the corresponding authors.

Acknowledgments

This research is grateful to the Universidad de La Guajira for the contributions and grant to J.I.P.-M; to the CREACUA Foundation (Riohacha, Colombia) for the impulse given to this initiative; to the Fondo Adaptación of Colombia for the images from the 2017 aerial flight of Riohacha as inputs to obtain the DTM; and to the Dirección General Marítima y Portuaria (DIMAR) for the mean sea level information from the Puerto Brisa tide gauge. A special greeting goes to Parques Nacionales Naturales de Colombia and particularly to all the public servants of the Santuario de Fauna y Flora los Flamencos (SFFF) and their illuminated director Nianza Angulo Paredes. Last, but not least, a warm thank you to Cristian Diaz Moscote for her help in the preparation of the figures.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. System approach: our lagoon is a dynamic system with several natural exogenous inputs (river freshwater inflow; direct precipitation inflow; climatic variables determining the evaporation rate; sea level with its astronomic and meteorological pattern)—all affected by climate change—and possibly a human-driven exogenous input when interventions are applied to open or close the lagoon mouth. (a) Full scheme with two state variables: the storage volume V(t) and the area A(t) (and shape) of the mouth exchanging water with the sea; (b) simplified scheme (“core model”) useful to test and calibrate the model, with just one state variable (storage), and the mouth area is considered a known exogenous variable, greatly simplifying the analysis.

Figure 2. Process adopted to generate the topo-bathymetric DEM: (a) topo-bathymetric GNSS-RTK campaign; (b) interpolation of points to generate the bathymetric surface y (c) integration with topographic DEM for the dry areas, obtained from a points cloud from photogrammetry of images obtained by a crewed aircraft.

Figure 3. Relationships: (a) elevation–volume y = y(V), (R2 = 0.999); (b) area–volume A = A(V) (R2 = 0.986).

Figure 4. Water exchange with the sea: (a) the boca (open) of the lagoon during an “open period” (ancho: width; longitude: length); (b) hydraulic schematization adopted for the alternative method to estimate flows; the slope is approximated as Δy/L (with reference to (a)).

Figure 5. Correspondence between measured and estimated flowrates: (a) Matching between measured exchange flows (QB) and estimated flow (q) through Equations (3) and (4); that is, as a function of the height difference and the cross-section area A (R2 = 0.9103). Data utilized from 8 June 2022 until 16 December 2022; (b) matching from the simplified hydraulic method (R2 = 0.9279).

Figure 6. Temporal pattern of lagoon elevation and mouth area: (a) time pattern of lagoon elevation (gray graph) and boca status (1: closed; 2: semi open; 3: open) during our monitoring period (10 December 2021–30 June 2023) estimated from real-world data (orange graph), and modelled according to the logic explained above (blue graph); (b) time pattern of the boca cross-sectional area A extrapolated from the few measurements available (orange graph) and determined by the model implementing the theory presented here by assuming the same starting point (blue graph).

Figure 7. Multiannual monthly temperature pattern (data from IDEAM, Almirante Padilla airport station of the town of Riohacha, about 20 km away from the lagoon; period 1 October 1972–1 March 2020).

Figure 8. Forcing factors and boundary conditions for model 2D configuration.

Figure 9. 0D model performance: (a) River inflow and precipitation measured; (b) General behavior of the model coherent with the system physics when there is a significant river inflow or precipitation. See text for detailed explanations. Notice that to allow discerning the lagoon elevation from that of the sea, the latter is reported without the shift correction (φ = 0.618 m) that the model sees; so, the sea graphs appear in the air, but have to be seen exactly shifted downwards of that quantity. For the sea, the light dense line on the background is the actual sea level measured; the darker shorter amplitude line on top is its 24 h moving average. The bottom line denotes the status of the boca: O: Open; C: Closed; S: Semiopened.

Figure 10. The system, represented by the model, appears to be very sensitive to the evaporation rate (same parametrization): the adopted rate e = 1.0 cm/day; e = 1.5 cm/day; and e = 0.5 cm/day. The bottom line denotes the status of the boca: O: Open; C: Closed; S: Semiopened.

Figure 11. The amplitude of oscillations of the water surface elevation (during periods with open mouth) depends very much on the cross-section area of the boca, but not so the general pattern; for example, for the period 7000–9000: the pattern is almost identical for the reconstructed (adopted) cross-section areasand for the same area pattern increased by 50% (inset top-right). The bottom line denotes the status of the boca: O: Open; C: Closed; S: Semiopened.

Figure 12. Comparison of the model answer based on the measured area of the boca (light blue line, “Lag simulated”) and that based on the theory and model adopted for the boca dynamics (orange line, “y simulated mouth dynamics”). The bottom line denotes the status of the boca: O: Open; C: Closed; S: Semiopened.

Figure 13. Results of the 0D and 2D models compared to measured water level data.

Figure 14. The reconstructor scheme of the river flowrate Q: (a) overall view; (b) detailed view. It needs the input time series virtually from a far instant before (t0) until current time t to reconstruct the flowrate time series Q(t0,t − 1) until the last instant t − 1.

Table 1

Data obtained in the measurement campaign at “la boca” where the flow rate is the flow exchanged with the sea (colors distinguish dates; grey lines had to be discarded because of internal inconsistency or too low, and hence unreliable, height difference values). The value in the blue cell is the correction constant ϕ found and applied to obtain the corrected sea elevation y*, and the “delta y*” is (y_lagoon – y_{sea corrected}) when the flow is outgoing (“OUT” from the lagoon into the sea) and vice versa when it is incoming (“IN”), therefore it must be always positive (only September 30, does not match). The (arbitrary) classification of the cross-section area (Large, above 40 m²; Small, below) partitions the data in two groups with meaningful graphical relationships (while one global is more confused).

Date	Time	Flow Direction		COD	Mouth Area	Small (S) or Large (L)	ExchangeFlowrate QB	y Sea P. Brisa	y* Sea P. Brisa	y* Lagoon Piezometer	Delta y* Piezometer
		OUT	IN		(m2)		(m3/s)	(m.a.s.l)	(m.a.s.l)	(m.a.s.l)	(m.a.s.l)
									0.618
10 December 2021	08:46	x		1	14.53	S	6.9	-	-	-	-
10 December 2021	14:28		x	0	17.98	S	7.9	-	-	-	-
8 June 2022	09:20	x		1	38.04	S	19.4	0.400	−0.218	−0.076	0.142
8 June 2022	11:00	x		1	35.31	S	27.2	0.350	−0.268	−0.092	0.176
7 July 2022	05:35		x	0	0.04	S	0.0	0.480	−0.138	−0.166	0.028
30 September 2022	09:30		x	0	41.31	L	27.0	0.630	0.012	0.034	-0.022
30 September 2022	15:30		x	0	56.04	L	36.0	-	-	-	-
30 September 2022	15:30		x	0	56.04	L	44.0	0.845	0.227	0.115	0.112
28 October 2022	09:37	x		1	92.23	L	49.1	0.640	0.022	0.053	0.031
28 October 2022	09:37	x		1	92.23	L	47.9	-	-	-	-
28 October 2022	12:10		x	0	97.92	L	30.6	0.750	0.132	0.120	0.012
25 November 2022	10:30		x	0	68.26	L	61.0	0.70	0.082	−0.140	0.222
25 November 2022	10:30		x	0	68.26	L	61.1	-	-	-	-
25 November 2022	14:30		x	0	63.18	L	50.4	0.675	0.057	−0.088	0.145
16 December 2022	10:30	x		1	19.40	S	5.9	0.455	−0.163	−0.159	0.004

Table 2

Input data and key parameters for hydrodynamic module.

Parameter	Value
Time step, ∆t (s)	2.0
Manning Coefficient	0.01
Background horizontal eddy viscosity (m2/s)	0.3
Background vertical eddy viscosity (m2/s)	1 × 10−5
Vertical molecular diffusivity (m2/s)	1 × 10−6
Total computation time (h)	18

Table 3

The water balance of the lagoon obtained through the simulation model (period refers to 20 December 2021 until 30 June 2023). Notice that the figures “Yearly” and “Monthly” are derived afterwards as a simple proportion to give an idea of the amounts referred to usual comparison periods.

		Period	Yearly	Monthly
IN	river	94.00	61.54	5.13
	precipitation and runoff	30.78	20.15	1.68
	incoming sea	195.10	127.72	10.64
	Total	319.88	209.41	17.45
OUT	outgoing to sea	265.09	173.55	14.46
	evaporation	58.46	38.27	3.19
	storage loss	−1.90	−1.24	−0.10
	Total	321.65	210.58	17.55
BALANCE	IN-OUT	−1.78	−1.17	−0.10
Error	as a % of IN	0.56	0.55	0.55

Table 4

Performance of the two models.

Metric	0D Model	EFDC Model
MAE	7.784	9.718
RMSE	9.897	12.584
NSE	0.846	0.751
KGE	0.812	0.822

References

1. Nardini, A.G.C.; Escobar Villanueva, J.R.; Pérez-Montiel, J.I. Hydrological Monitoring System of the Navío-Quebrado Coastal Lagoon (Colombia): A Very Low-Cost, High-Value, Replicable, Semi-Participatory Solution with Preliminary Results. Water; 2024; 16, 2248. [DOI: https://dx.doi.org/10.3390/w16162248]

2. Zoidou, M.; Kokkos, N.; Sylaios, G. Dynamics of Water, Salt, and Nutrients Exchange at the Inlets of Three Coastal Lagoons. J. Mar. Sci. Eng.; 2022; 10, 205. [DOI: https://dx.doi.org/10.3390/jmse10020205]

3. Michielotto, A.; Tognin, D.; Matticchio, B.; Carniello, L.; D’Alpaos, A. Investigating Natural and Anthropogenic Impacts on the Hydrodynamics of the Venice Lagoon (Italy), a Numerical Approach. Proceedings of the EGU General Assembly Conference Abstracts; Vienna, Austria, 23–28 April 2023; EGU-8582.

4. Alba, J.G.; Gómez, A.G.; del Barrio Fernández, P.; Gómez, A.G.; Álvarez Díaz, C. Hydrodynamic Modelling of a Regulated Mediterranean Coastal Lagoon, the Albufera of Valencia (Spain). J. Hydroinformatics; 2014; 16, pp. 1062-1076. [DOI: https://dx.doi.org/10.2166/hydro.2014.071]

5. Sousa, M.C.; Dias, J.M. Hydrodynamic Model Calibration for a Mesotidal Lagoon: The Case of Ria de Aveiro (Portugal). J. Coast. Res.; 2007; 50, 1075. [DOI: https://dx.doi.org/10.2112/JCR-SI50-194.1]

6. Simonetti, I.; Cappietti, L. Influence of Inlets Morphology and Forcing Mechanisms on Water Exchange between Coastal Basins and the Sea: A Hindcast Study for a Mediterranean Lagoon. J. Mar. Sci. Eng.; 2022; 10, 1929. [DOI: https://dx.doi.org/10.3390/jmse10121929]

7. Ferrarin, C.; Umgiesser, G. Hydrodynamic Modeling of a Coastal Lagoon: The Cabras Lagoon in Sardinia, Italy. Ecol. Modell.; 2005; 188, pp. 340-357. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2005.01.061]

8. Huang, P.; Hennig, K.; Kala, J.; Andrys, J.; Hipsey, M.R. Climate Change Overtakes Coastal Engineering as the Dominant Driver of Hydrological Change in a Large Shallow Lagoon. Hydrol. Earth Syst. Sci.; 2020; 24, pp. 5673-5697. [DOI: https://dx.doi.org/10.5194/hess-24-5673-2020]

9. IPCC. Climate Change 2022: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; Pörtner, H.-O.; Roberts, D.C.; Tignor, M.; Poloczanska, E.S.; Mintenbeck, K.; Alegr, A. Cambridge University Press: Cambridge, UK, New York, NY, USA, 2022.

10. Murillo, L.G.; Florián Buitrago, M.; Pabón Restrepo, G.A.; Pérez Álvarez, P.A.; Rojas Laserna, M.; Suárez Castaño, R. Política Nacional de Cambio Climático: Documento Para Tomadores de Decisiones; Ministerio de Ambiente y Desarrollo Sostenible: Bogotá, Colombia, 2017.

11. Nartiss, M. R.Lake.Xy Module. Available online: https://grass.osgeo.org/grass-stable/manuals/r.lake.html (accessed on 5 May 2024).

12. Chow, V.T. Open-Channel Hydraulic; McGraw-Hill: New York, NY, USA, 1959; ISBN 978-0-07-085906-7

13. Corporación Autónoma Regional de La Guajira [CORPOGUAJIRA]. Ajuste y/o Actualización Del POMCA Del Río Camarones y Otros Directos Al Caribe—La Guajira; Corporación Autónoma Regional de La Guajira [CORPOGUAJIRA]: Riohacha, Colombia, 2022.

14. Young, D.M.; Gregory, R.T. A Survey of Numerical Mathematics; Addison-Wesley Pub. Co.: Reading, MA, USA, 1972; Volume 1.

15. Moraru, A.; Rüther, N.; Bruland, O. Investigating Optimal 2D Hydrodynamic Modeling of a Recent Flash Flood in a Steep Norwegian River Using High-Performance Computing. J. Hydroinform.; 2023; 25, pp. 1690-1712. [DOI: https://dx.doi.org/10.2166/hydro.2023.012]

16. Mao, M.; Xia, M. Seasonal Dynamics of Water Circulation and Exchange Flows in a Shallow Lagoon-Inlet-Coastal Ocean System. Ocean Model.; 2023; 186, 102276. [DOI: https://dx.doi.org/10.1016/j.ocemod.2023.102276]

17. Torres-Bejarano, F.M.; Torregroza-Espinosa, A.C.; Martinez-Mera, E.; Castañeda-Valbuena, D.; Tejera-Gonzalez, M.P. Hydrodynamics and Water Quality Assessment of a Coastal Lagoon Using Environmental Fluid Dynamics Code Explorer Modeling System. Glob. J. Environ. Sci. Manag.; 2020; 6, pp. 289-308. [DOI: https://dx.doi.org/10.22034/gjesm.2020.03.02]

18. Tsanis, I.K.; Wu, J.; Shen, H.; Valeo, C. Environmental Hydraulics: Hydrodynamic and Pollutant Transport Modelling of Lakes and Coastal Waters. Environmental Hydraulics; Tsanis, I.K.; Wu, J.; Shen, H.; Valeo, C. 2549; Developments in Water Science; Elsevier: Amsterdam, The Netherlands, 2006.

19. Gunasinghe, G.P.; Ruhunage, L.; Ratnayake, N.P.; Ratnayake, A.S.; Samaradivakara, G.V.I.; Jayaratne, R. Influence of Manmade Effects on Geomorphology, Bathymetry and Coastal Dynamics in a Monsoon-Affected River Outlet in Southwest Coast of Sri Lanka. Environ. Earth Sci.; 2021; 80, 238. [DOI: https://dx.doi.org/10.1007/s12665-021-09555-0]

20. Herrling, G.; Winter, C. Spatiotemporal Variability of Sedimentology and Morphology in the East Frisian Barrier Island System. Geo-Mar. Lett.; 2017; 37, pp. 137-149. [DOI: https://dx.doi.org/10.1007/s00367-016-0462-6]

21. Nienhuis, J.H.; Törnqvist, T.E.; Esposito, C.R. Crevasse Splays Versus Avulsions: A Recipe for Land Building With Levee Breaches. Geophys. Res. Lett.; 2018; 45, pp. 4058-4067. [DOI: https://dx.doi.org/10.1029/2018GL077933]

22. Li, X.; Leonardi, N.; Plater, A.J. Impact of Barrier Breaching on Wetland Ecosystems under the Influence of Storm Surge, Sea-Level Rise and Freshwater Discharge. Wetlands; 2020; 40, pp. 771-785. [DOI: https://dx.doi.org/10.1007/s13157-019-01219-x]

23. Elsayed, S.M.; Oumeraci, H. Combined Modelling of Coastal Barrier Breaching and Induced Flood Propagation Using XBeach. Hydrology; 2016; 3, 32. [DOI: https://dx.doi.org/10.3390/hydrology3040032]

24. Maicu, F.; Abdellaoui, B.; Bajo, M.; Chair, A.; Hilmi, K.; Umgiesser, G. Modelling the Water Dynamics of a Tidal Lagoon: The Impact of Human Intervention in the Nador Lagoon (Morocco). Cont. Shelf Res.; 2021; 228, 104535. [DOI: https://dx.doi.org/10.1016/j.csr.2021.104535]

25. Nienhuis, J.H.; Heijkers, L.G.H.; Ruessink, G. Barrier Breaching Versus Overwash Deposition: Predicting the Morphologic Impact of Storms on Coastal Barriers. J. Geophys. Res. Earth Surf.; 2021; 126, e2021JF006066. [DOI: https://dx.doi.org/10.1029/2021JF006066]

26. Wright, K.A.; Goodman, D.H.; Som, N.A.; Alvarez, J.; Martin, A.; Hardy, T.B. Improving Hydrodynamic Modelling: An Analytical Framework for Assessment of Two-Dimensional Hydrodynamic Models. River Res. Appl.; 2017; 33, pp. 170-181. [DOI: https://dx.doi.org/10.1002/rra.3067]

27. Yang, Q.; Wu, W.; Wang, Q.J.; Vaze, J. A 2D Hydrodynamic Model-Based Method for Efficient Flood Inundation Modelling. J. Hydroinform.; 2022; 24, pp. 1004-1019. [DOI: https://dx.doi.org/10.2166/hydro.2022.133]

28. Kamoun, M.; Zaïbi, C.; Langer, M.R.; Carbonel, P.; Ben Youssef, M. Coastal Dynamics and the Evolution of the Acholla Lagoon (Gulf of Gabes, Tunisia): A Multiproxy Approach. Arab. J. Geosci.; 2021; 14, 2021. [DOI: https://dx.doi.org/10.1007/s12517-021-07590-6]

29. Rohmer, J.; Sire, C.; Lecacheux, S.; Idier, D.; Pedreros, R. Improved Metamodels for Predicting High-Dimensional Outputs by Accounting for the Dependence Structure of the Latent Variables: Application to Marine Flooding. Stoch. Environ. Res. Risk Assess.; 2023; 37, pp. 2919-2941. [DOI: https://dx.doi.org/10.1007/s00477-023-02426-z]

30. Liu, R.; Wei, J.; Ren, Y.; Liu, Q.; Wang, G.; Shao, S.; Tang, S. HydroMP—A Computing Platform for Hydrodynamic Simulation Based on Cloud Computing. J. Hydroinform.; 2017; 19, pp. 953-972. [DOI: https://dx.doi.org/10.2166/hydro.2017.140]

31. Jeong, S.; Yeon, K.; Hur, Y.; Oh, K. Salinity Intrusion Characteristics Analysis Using EFDC Model in the Downstream of Geum River. J. Environ. Sci.; 2010; 22, pp. 934-939. [DOI: https://dx.doi.org/10.1016/S1001-0742(09)60201-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20923109]

32. Mathis, T.J.; Lee, S.-T.; Craig, P.M.; Lam, N.T.; Scandrett, K.; Mishra, A.; Arifin, R.R.; Jung, J.Y. Estuarine Salinity Intrusion and Implications for Aquatic Habitat: A Case Study of the Lower St. Johns River Estuary, Florida. Proceedings of the World Environmental and Water Resources Congress 2019; Pittsburgh, PA, USA, 19–23 May 2019; pp. 13-25.

33. Feng, L.U.O.; Ruijie, L.I. 3D Water Environment Simulation for North Jiangsu Offshore Sea Based on EFDC. J. Water Resour. Prot.; 2009; 1, pp. 41-47. [DOI: https://dx.doi.org/10.4236/jwarp.2009.11007]

34. Shin, S.; Her, Y.; Song, J.-H.; Kang, M.-S. Integrated Sediment Transport Process Modeling by Coupling Soil and Water Assessment Tool and Environmental Fluid Dynamics Code. Environ. Model. Softw.; 2019; 116, pp. 26-39. [DOI: https://dx.doi.org/10.1016/j.envsoft.2019.02.002]

35. Torres-Bejarano, F.; González-Martínez, J.; Rodríguez-Pérez, J.; Rodríguez-Cuevas, C.; Mathis, T.J.; Tran, D.K. Characterization of Salt Wedge Intrusion Process in a Geographically Complex Microtidal Deltaic Estuarine System. J. S. Am. Earth Sci.; 2023; 131, 104646. [DOI: https://dx.doi.org/10.1016/j.jsames.2023.104646]

36. Hamrick, J.M. A Three-Dimensional Environmental Fluid Dynamics Computer Code: Theoretical and Computational Aspects; Special Report in Applied Marine Science and Ocean Engineering; no. 317 Virginia Institute of Marine Science, William & Mary: Gloucester Point, VA, USA, 1992; 64p. [DOI: https://dx.doi.org/10.21220/V5TT6C]

37. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. Trans. ASABE; 2007; 50, pp. 885-900. [DOI: https://dx.doi.org/10.13031/2013.23153]

38. Gupta, H.V.; Kling, H.; Yilmaz, K.K.; Martinez, G.F. Decomposition of the Mean Squared Error and NSE Performance Criteria: Implications for Improving Hydrological Modelling. J. Hydrol.; 2009; 377, pp. 80-91. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2009.08.003]

39. Nash, J.E.; Sutcliffe, J. V River Flow Forecasting through Conceptual Models Part I—A Discussion of Principles. J. Hydrol.; 1970; 10, pp. 282-290. [DOI: https://dx.doi.org/10.1016/0022-1694(70)90255-6]

40. Kalman, R.E. A New Approach to Linear Filtering and Prediction Problems. J. Basic Eng.; 1960; 82, pp. 35-45. [DOI: https://dx.doi.org/10.1115/1.3662552]

41. Bavdekar, V.A.; Deshpande, A.P.; Patwardhan, S.C. Identification of Process and Measurement Noise Covariance for State and Parameter Estimation Using Extended Kalman Filter. J. Process Control; 2011; 21, pp. 585-601. [DOI: https://dx.doi.org/10.1016/j.jprocont.2011.01.001]

42. Protsenko, S.V.; Protsenko, E.A.; Kharchenko, A. V Comparison of Hydrodynamic Processes Modeling Results in Shal-Low Water Bodies Based on 3D Model and 2D Model Averaged by Depth. Comput. Math. Inf. Technol.; 2023; 6, pp. 49-63.

43. Hernandez, I.; Castro-Rosero, L.M.; Espino, M.; Alsina Torrent, J.M. LOCATE v1.0: Numerical Modelling of Floating Marine Debris Dispersion in Coastal Regions Using Parcels v2.4.2. Geosci. Model Dev.; 2024; 17, pp. 2221-2245. [DOI: https://dx.doi.org/10.5194/gmd-17-2221-2024]

44. Hu, D.; Chen, Z.; Li, Z.; Zhu, Y. An Implicit 1D-2D Deeply Coupled Hydrodynamic Model for Shallow Water Flows. J. Hydrol.; 2024; 631, 130833. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2024.130833]

45. Kabi, J.N.; wa Maina, C.; Mharakurwa, E.T.; Mathenge, S.W. Low Cost, LoRa Based River Water Level Data Acquisition System. HardwareX; 2023; 14, e00414. [DOI: https://dx.doi.org/10.1016/j.ohx.2023.e00414]

46. Notter, B.; MacMillan, L.; Viviroli, D.; Weingartner, R.; Liniger, H.P. Impacts of Environmental Change on Water Resources in the Mt. Kenya Region. J. Hydrol.; 2007; 343, pp. 266-278. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2007.06.022]