Natural and Human-Induced Flow and Sediment Transport within Tidal Creek Networks Influenced by Ocean-Bay Tides

Abstract

Improving current understanding of hydrodynamics and sediment dynamics in complex tidal embayments is of major importance to face future challenges derived from climate change and increasing human pressure. This work deepens the knowledge of the hydro-morphodynamics of complex creek networks that connect basins with different characteristics, identifying their morphodynamic trends and the potential impacts of channel deepening. We selected two tidal creeks which flow through salt marshes and tidal flats of the Cádiz Bay (SW Spain) in a singular network due to their double connection to the Atlantic Ocean and the inner bay. We study the interactions between tidal waves that penetrate into the creeks from these two different bodies of water, analyzing the tidal asymmetry and the morphodynamic tendencies of the system. For the analysis, we set up a hydro-morphodynamic model specifically developed for areas with very shallow and complex channels. Results show that the tidal wave penetrates within the tidal network both from the inner Bay and the open ocean with different amplitudes, phases and flow velocities. There is also an asymmetric pattern for the tidal flows caused by the deformation of the dominant astronomical tidal constituents, M2 and M4, due to the non-linear interaction of tidal currents with the irregular creek geometry and bottom topography. Tidal asymmetry promotes the progressive infilling of the area where the tidal waves meet closing the connection between the open ocean and the inner bay, such an infilling trend being accelerated by human interventions.

1. Introduction

Tidal embayments such as bays and lagoons are usually dissected by networks of tidal channels that connect these embayments with the open sea. Furthermore, they provide preferential pathways for the exchange of water, sediments and nutrients among the typical ecogeomorphic features of tidal systems such as salt marshes and tidal flats [1,2,3,4]. Tidal systems are highly valuable for geomorphological, ecological and socioeconomic reasons, and are of great concern for managers and scientists being these systems vulnerable to variations in the environmental forcing as a consequence of climate changes and extensive anthropogenic activities (e.g., [5]). Some of the climate change components, whose outcomes are a great threat to the tidal systems, are sea-level rise, storminess or extreme precipitation, among others. The results of these events change the shoreline development facing coastal hazards such as erosion, flooding or storm waves [6,7]. The deterioration of coastal water quality and reduction of biodiversity, or the impacts on coastal habitats are some of the consequences of these modifications on the tidal systems [8]. These changes are increased when human interventions are carried out in these areas [9,10] since these anthropogenic disturbances and their interplays modify the shape and evolution of these coastal areas [3].

During recent years, many studies have improved the current understanding of the hydro- and morphodynamic behavior of coastal systems such as estuaries, lagoons or bays, and how they are affected by climate change and human interventions [11,12,13,14,15,16,17,18]. The methods employed for such analyses range between field observations (e.g., [11,19,20,21,22]), numerical models (e.g., [18,23,24,25,26,27,28,29,30]), analytical models [31,32,33] and laboratory experiments [22,34,35,36]. Although these studies have improved our knowledge of the dynamics of tidal channel networks, our understanding of their evolution is still incomplete [3] and calls for new insights in order to adopt strategic and sustainable management actions. Improving our understanding of the hydrodynamic and morphodynamic features of tidal channel networks is, in fact, a key step to understand and manage the dynamics of the landscapes they cut through (e.g., [3,23,37]).

Sediment dynamics in tidal channels are strongly governed by the marked tidal asymmetry and distortion of tidal currents that develop as a consequence of non-linear tidal effects as the tidal wave propagates from the inlets of the tidal system towards its inner portions (e.g., [16,38,39,40,41,42]). The non-linear friction effects produce a time-lag in sea-level changes, being the flood periods shorter and stronger than ebbs [39,40,43]. Moreover, the non-linear effects, which are more important for areas with complex bathymetries, play an important role in the morphodynamic behavior, promoting a tidally averaged net flow of sediment (residual transport) [44,45].

As an example of a complex tidal creek system we consider here the case of the Sancti-Petri and Carraca Creeks (SPC and CC, respectively, hereinafter), whose evolution is strongly influenced by the complex tidal hydrodynamics of the Cádiz Bay, heavily altered by human interventions during last decades [46,47]. Several researchers have focused their efforts to analyze hydrodynamics [48,49,50,51,52], sediment dynamics [53] and tidal energy [17] of the outer and inner basins of the Cádiz Bay. Nevertheless, to the authors knowledge no specific studies have been conducted to analyze the tidal and suspended-sediment transport characteristics at SPC and CC that directly connect the inner bay with the open ocean. The singularity of these tidal creeks lies on their convergence in an inflow–outflow channel which extends from the inner basin of the Cádiz Bay directly to the Atlantic Ocean. These two different bodies of water have different hydrodynamic characteristics, with tidal waves propagating through both creeks with different amplitudes and phases until their intersection [54], an unexplored critical point [51] for which the morphodynamic evolution can play a major role on the interaction between the two waves and, hence, on hydro-morphodynamics of the whole system.

In such complex systems, numerical models are used to describe tidal and sediment transport dynamics in order to plan and design effective restoration activities and evaluate their impact. Some examples of such 2D and/or 3D models, which are usually coupled to sediment transport models [55], are MIKE21 [56], SWAN [57], TELEMAC [58], SHYFEM [59] and Delft3D. However, these models present limitations to simulate the combined effect of tidal currents and wind waves on the morphological evolution in very shallow domains with highly irregular bathymetries [46]. Here we used a numerical model (Wind-Wave Tidal Model, see Section 2.3) specifically developed by [60] and [61] for shallow and irregular domains characterized by the presence of salt marshes periodically experiencing wetting and drying transition. This model has been widely applied to study different study sites such as the Venice lagoon [62,63,64], the lagoons of the Virgina Coastal Reserve [65] and the connection between the inner and outer basins of Cádiz Bay [46].

The main objective of this manuscript is to provide new insights on the tidal dynamic processes of complex tidal creek networks with double connection with the open sea, focusing on the role of tidal asymmetries and human interventions in driving sediment transport processes which can lead to irreversible morphological changes. The analysis is made by means of a hydro-morphodynamic numerical model developed by [61,62], which is applied to a complex tidal creek network at the Cádiz Bay as an existing example of this type of coastal systems. A field survey was employed to calibrate and test model results, which are used to discuss the impacts of (1) tidal asymmetries along the creeks and (2) human interventions carried out in the area. The paper is organized as follows: Section 2 describes the study site, the numerical model and its implementation on the study site. Section 3 presents the main results of the characterization of the creek hydrodynamics (tidal range and currents) and their sediment dynamics. Section 4 discusses the impacts of tidal asymmetries and human interventions on the creeks. Finally, the last section summarizes the main conclusions derived from this study.

2. Material and Methods

2.1. The Sancti-Petri (SPC) and Carraca (CC) Creeks

Cádiz Bay is a $7000\mathrm{ha}$, marsh-dominated, coastal plain estuary in southern Spain (${36}^{°}31$${}^{\prime }$ N, $6^{°}{17}^{\prime }$ W, Figure 1) that supports a surprising wealth of wildlife despite the important increase in human activities, such as infrastructure development or channel deepening. These interventions have modified the ecosystem of the bay [66] causing the migration of many species [67]. SPC and CC constitute the most important portions of the channel network that cuts through the coastal marsh system of the Cádiz Bay Natural Park (at the South of Cádiz (yellow)-Figure 1), which encompasses a 105.22 ${\mathrm{km}}^2$ flat landscape of sandy beaches, marshes, salt pans, freshwater lakes and tidal inlets, as well as the two natural areas of Trocadero Island and the creeks of SCC and CC (Figure 1b).

The natural configuration of SPC and CC has been frequently modified by human interventions, such as the construction of different bridges and roads along the creeks, salt extraction at SPC and aquaculture in their surroundings, among others. The distance between the mouths of SPC and CC, and the point where the two creeks meet is about 15 km and 7.5 km, respectively. The entire creek system drains an area of $170\mathrm{ha}$ that comprises intertidal flats around the fringes, subtidal flats (up to 2 m deep) and tidal channels (6–8 m deep), discharging into the Atlantic Ocean across SPC. The tide in the area is mesotidal and semidiurnal, with a mean neap tidal range of ≈ 1.0 m, while the mean spring tidal range is ≈4.0 m. The water column is typically well-mixed. The highest values of suspended sediment concentrations measured at SPC are close to $13\mathrm{mg}/\mathrm{L}$ predominantly composed by muds, since this area is protected from winds and waves [68].

The mouths of the two creeks, as well as the morphology of the area, are quite dynamic as shown in Figure 2. This Figure shows the morphological changes observed since 1956: the coastline has been characterized by rapid changes that increased in magnitude during the last ten years. It shows areas of horizontal progradation, especially at the northern and western area of the creeks, and eastern locations where the retreat of the coastline due to erosion was observed. Although CC has experienced smaller modifications, relevant deposition processes have been observed at the connection between the two creeks (Figure 2). The results presented throughout this work show that this effect is due to the increase in sediment availability in this area, and also that it has important consequences on the hydrodynamic circulation and sediment dynamics of the whole system.

2.2. Field Survey

A field survey carried out between February and March (2012) was used to determine the characteristics of the tide in the study site, with two tidal gauges (T2 and T3 in Figure 1) deployed to measure the water level every $15\min$ within the two creeks and four current meters (ADCPs) deployed within the bay (I1–I4 in Figure 1). These data were used to monitor the hydrodynamic field (Figure 3) and provide data for model calibration and testing. The ADCPs collected data for 40 consecutive days with a sampling time of $15\min$ and an average interval of $120\mathrm{s}$ [46]. Tidal elevation data were obtained from the water level time series [69].

2.3. The Hydro-Morphodynamic Model

The numerical model used in the present study is a hydro-morphodynamic model made up of two modules sharing the same computational grid: (i) a hydrodynamic module coupled with a Wind-Wave Tidal Module (hereinafter WWTM), developed by [61,62]; and (ii) a Sediment Transport and Bed Evolution Module (hereinafter STABEM) developed by [63].

The hydrodynamic module solves the 2-D shallow water equations using a semi-implicit staggered finite element method based on Galerkin’s approach. The governing equations are modified for reproducing flooding and drying processes in very shallow and irregular domains (see [25,60] for a detailed description of the equations and of the numerical scheme). Considering the case of a turbulent flow over a rough wall, the hydrodynamic module computes the bottom shear stress induced by tidal currents using the Strickler formula.

The wind-wave module solves an approximation of the wave action conservation equation [70] that is parameterized using the zero-order moment of the wave action spectrum in the frequency domain [71]. The wind-wave module assumes that the direction of wave propagation instantaneously adjusts to the wind direction and compute: the average wave energy, the significant wave height (derived using linear wave theory) and the peak wave period. The wave period is estimated using an empirical correlation function, which relates the mean peak wave period to the local wind speed and water depth following the approach suggested by [72]. We refer the reader to [62] for a detailed description. The wind-wave module computes the bottom shear stress induced by wind waves as a function of the maximum horizontal orbital velocity at the bottom, related to the significant wave height through the linear theory.

Considering the nonlinear interaction between waves and currents within the boundary layer, the total bottom shear stress resulting from the combined effect of tidal currents and wind waves is computed by WWTM on the basis of the empirical formulation suggested by [73].

The sediment transport and bed evolution module (STABEM) deals with mixtures of sand and mud (comprised by clay and silt [63]) at each computational element, which can be defined as cohesive (sand), non-cohesive (mud) or a mixture of both homogeneously distributed depending on the mud content which varies in space and time. To evaluate the sediment resuspension rate, both the total bottom shear stress and the critical shear stress are treated as random variables. The stochastic approach gives a transport parameter which varies gradually at near-threshold conditions and improves the description of the initiation of the erosion process.

The sediment content is continuously updated by the model in response to the local entrainment and deposition. For that, STABEM solves the advection diffusion equation and the Exner equation for each fraction of sediment and is coupled with the WWTM following an in line approach [74], i.e., bed elevation and composition are updated at each computational time step. In addition to the calibration and testing for the Cádiz Bay [46] WWTM and STABEM have been widely tested [62,63] and used to describe the hydrodynamics and the sediment dynamics in the Venice lagoon [64,75,76] and in the lagoons of the Virginia Coast Reserve, VA, USA [65].

2.4. Model Application to the C ádiz Bay

The simulations presented herein were carried out using a mesh covering the Cádiz Bay, its tidal creek systems and a portion of the Atlantic coast of Spain close to the Gibraltar strait (Figure 1). The mesh was locally refined over the SPC-CC system on the basis of an accurate bathymetric survey and of a precise topography, which were not available for the model setup presented in [46]. The offshore bathymetry data were provided by the Instituto Hidrográfico de la Marina (Spanish Ministry of Defense, Madrid, Spain), while a detailed multi-beam bathymetry of the bay, provided by the Port Authority of Cádiz (Cádiz, Spain), was integrated by the Digital Elevation Model of the Regional Government (Junta de Andalucía, Spain) at a resolution of $10\times 10{\mathrm{m}}^2$. In areas where multiple data sets overlapped, preference was given to the highest resolution multi-beam data [77]. The computational grid consists of 31,400 nodes and more than 62,100 triangular elements. The finest grid size ($10\times 10{\mathrm{m}}^2$) was used for describing the central section of the bay, SCP and CC, together with the tidal flat portions of the bay (Figure 1b). The coarser cell sizes (approximately $80\times 80{\mathrm{m}}^2$), were used to discretize the open ocean and the wide flat areas of the bay (Figure 1b). The time step for the computation, performed for a period of 21 days to include a complete spring-neap tidal cycle, was set equal to $2\mathrm{s}$ to ensure numerical stability.

According to the STABEM model, sediments over the numerical domain were classified in two types: non-cohesive sand and cohesive mud, which is the sum of clay and silt. The actual bed composition is described in the model defining, at each element of the computational grid, the percentage of mud per unit volume which can be 0, 50 or 100%. Figure 4 shows the initial sediment distribution used for the Cádiz Bay: the deepest channel, the tidal flats and the marshes are almost completely composed by mud (Figure 4), whereas the coast and the surroundings of the marshes are composed by fine sand (Figure 4). The mean grain diameter for pure sand and mud were defined as 200 and 20 $\mathsf{\mu }$m, respectively. This distribution was contrasted with previous works [78,79].

Although the WWTM module was calibrated and tested in the connection between the outer and inner basins of the bay by [46], here we focused our efforts on the calibration and testing of the hydrodynamics at SPC and CC. A four-step approach [52] was followed to ensure accurate results. First, the seaward boundary of the computational grid was divided into five portions and tidal levels were imposed at each portion reconstructed considering the amplitudes and phases of the twelve dominant components provided by the Oregon State Tidal Prediction Software (OTPS, [80]), and using the tools presented by [81] to reproduce the tidal levels from harmonic constituents. Second, the system was forced assuming a spatially uniform wind field characterized with data collected from Buoy 2342 (see Figure 1). Third, the model was calibrated considering a temporal subset of the flow velocity data collected at stations T2–T3, which enabled us to identify the more suitable value for the main model parameters, including a spatially variable bottom friction coefficient (Table 1 and Table 2, for a detailed description of the parameters the reader is referred to [62,63]. Finally, the model performance was tested by comparing the numerical results with the remaining (i.e., not used for calibrating) observed tidal elevations at T2–T3.

A quite good agreement between observed and simulated water levels was achieved (correlation coefficient ranges between $R^2=0.99$ and $R^2=0.88$ and skill coefficient [13,82] ranges between $0.99$ and $0.84$, see Table 3). The regression coefficient $R^2$ for tidal currents was observed to be lower, but the agreement remains quite good ($R^2=0.94$). Model results have been compared with extensive field studies (i.e., [15,83,84]) to verify the accuracy of the hydrodynamic results, being within the same range of values (0.8–0.99). In particular, Figure 5 shows the scatter plots of the computed vs observed water levels at the T2 and T3 stations, which are located within the creek system. In the plots we distinguish between data used for calibration and testing as well as between data referring to the ebb and flood periods. Better agreements were obtained during ebbs than during floods (red points in Figure 5). The differences, which are only relevant in the case of the T2 station, can be ascribed to the complex and dynamically active bed levels in the area, which can be hardly captured with the bathymetric data: despite the high-resolution bathymetry and topography data used, the accuracy in these very shallow areas where the data acquisition is complex and the bed levels vary rapidly is significantly lower. This limitation is more important in T2 that is located in the innermost part of the creek system.

3. Results

3.1. Hydrodynamics

Tide propagation and transformation within the SPC and CC network were analyzed by means of a harmonic analysis of both water levels and velocities [81], performed for the period between February 22 and March 14. Hereinafter, coordinates $(x,y)$ represent the cross-sectional and the along channel coordinates, respectively. The residual currents were calculated by averaging the flow over one M2 tidal cycle [85].

Model results provide tidal amplitudes at both stations (T2 and T3) that match quite well the observations, with differences of the O (1 cm) (Figure 5). Within the confidence intervals, the maximum (relative) difference between both sets of variables is of about 2% (Figure 5). Tidal phases agree even better with observations, with maximum (relative) differences of approximately 1%.

3.1.1. Surface Elevations

Tidal asymmetry, produced by non-linear terms associated with friction, advection, and tidal amplitude, plays an important role in coastal environments since it is responsible for the morphodynamic changes within these areas. This asymmetry is due to the production of even harmonics from the semi-diurnal and quarter-diurnal principal lunar tidal harmonics (M2 and M4 tides, respectively) for which Figure 6 provides maps showing tidal amplitudes and phases obtained with the numerical model.

Results show that the tidal amplitude of the M2 component (Figure 6a) decreases from the inlets of SPC and CC ($\sim 0.9\mathrm{m}$) towards the midpoint between (MP hereinafter), where the lowest amplitude of the M2 component is observed (up to 0.5 m in the main channel and 0.1 m over the tidal flats). The change in phase (Figure 6b) from the mouth of SPC to MP point is of about ${50}^{°}$, which represents a tidal delay of about 1.72 h, while the phase lag is about 1.03 h between the mouth of CC and MP. Considering the distance between the mouth of SCP and CC to MP point, 15 km and 7.5 km, respectively, it can be concluded that the tide travels quicker through SPC than CC. These phase lags are mainly due to dissipative friction effects within the both creeks until MP, thus suggesting that the wave takes longer when propagating through the creeks from their mouths, due to the larger length, the shallower bottoms and the presence of the tidal flat. Panels e and f of Figure 6 provide information on the M2 amplitude and phase, respectively, along the path shown in panels (a–d) (white dashed line). It emerges that the M2 amplitude (phase) reaches a minimum (maximum) value in MP. Hence, the tidal wave propagates faster in the main channels than over the tidal flats, in agreement with results of previous studies [42,51,86].

The M4 tidal constituent is locally generated in Cádiz creeks (Figure 6c) and further amplified from both mouths (0.01 m) up to MP (0.05 m). This induces a greater distortion of the tidal wave as it propagates upstream (downstream) starting from SPC (CC). The maximum amplitudes of the M4 component values are located in correspondence of the tidal flats (0.12 m) due to the great influence of friction in these areas. The behavior of the M4 phases (d—Figure 6) along the two creeks is very different due to the different widths of the channel cross-sections, the different changes in depths and the limited local influence. In the case of SPC, the phase increases quickly (from ${220}^{°}$ to ${330}^{°}$ in 4 km) until the narrowing of the creek is observed. However, this value is reduced up to 60% from the narrowing until MP. The phase at CC decreases from its mouth (${100}^{°}$) until MP. This behavior is clearly shown in panels g and h of Figure 6, where the maximum value of M4 is reached at MP (shallower areas).

Hence, the reduction of the tide amplitude and the amplification in the phase when the wave penetrates through both inlets generate a critical and complex point which corresponds to the convergence zone of the tidal wave.

3.1.2. Currents

Figure 7 shows the magnitude and spatial distribution of the flow rates ($q_x$, panels a and b), and $q_y$, panels c and d) along the x and y coordinates (across and along the channel, respectively), during flood and ebb, for spring and neap tides. Maximum values of $q_x$ and $q_y$ are observed for spring tides (first and third column) within the channel, whereas, over the tidal flats the flow rates are negligible. The east-west component of the velocities ($q_x$, panels a1–d1—Figure 7) ranges between 6${\mathrm{m}}^2/\mathrm{s}$ (flood period) and $-6{\mathrm{m}}^2/\mathrm{s}$ (ebb period). The maximum value is found at the mouth of CC and SPC. As with the North-South component of the flow rates per unit width ($q_y$, panels a2–d2—Figure 7), the maximum flood/ebb values are reached at the mouth of SPC (±6${\mathrm{m}}^2/\mathrm{s}$), whereas a decrease up to ±3${\mathrm{m}}^2/\mathrm{s}$ is observed at the mouth of CC. A strong reduction in both components, $q_x$ and $q_y$, is observed at MP, where flow rates almost vanish.

The results of the tidal-current harmonic analysis for the M2 component are shown in Figure 8. The maximum flow (a—Figure 8) are found in both the mouths of the creeks (0.6 m${}^2$/s). This is likely due to the decrease in channel width and depth from the mouth to point MP. There is a steady decrease in the magnitude of the M2 tidal flow from the mouths to the MP (0.1 m${}^2$/s), whereas the flow in the tidal flat is negligible. The maximum of the M2 phase (b—Figure 8) is observed at MP (80${}^{°}$) due to the convergence of the two waves that propagate from the mouths towards the inner creek region, in agreement with the behavior described in Section 3.1.1. Panel c in Figure 8 shows the difference in phase between tidal elevations and current velocities. Water levels and flow velocities are in quadrature (± 90${}^{°}$) within the deeper portions of the channels whereas the difference vanishes over the tidal flats. Furthermore, tides at SPC and CC can be considered to be two progressive waves of the same period traveling in opposite directions, from the mouths of the main channels to wards MP. The first wave propagates from the open sea (from SPC to MP); the second wave propagates from the inner bay towards MP (from CC to MP). At MP, the two waves have the same phase. Thus, a quasi-stationary character of the tidal wave is observed in the creek network, in analogy with previous observations [51].

Figure 9 shows the tidal ellipse parameters of the M2 constituent, namely the semi-major axis (a), the semi-minor axis (b), the inclination of the semi-major axis (c) of the ellipse, and the Greenwich phase (d). The semi-major axis (a—Figure 9) decreases towards MP from $3{\mathrm{m}}^2/\mathrm{s}$ to $0.1{\mathrm{m}}^2/\mathrm{s}$ and its phase (d—Figure 9) decreases from $\pm {340}^{°}$ to $\pm {20}^{°}$. The sharp reduction in M2 currents indicates that bottom friction exerts a large drag on the major astronomical tidal wave entering this shallow creek. The semi-minor axis of M2 (b—Figure 9) exhibits values close to 0. The motion of the tidal flow velocity is almost rectilinear in the main channel of both the creeks, so the semi-minor axis of the flow ellipse is smaller than the semi-major axis. However, important changes in the direction of the rotation of the ellipse are observed (a positive value indicates a cyclonic or anti-clockwise rotation) at both the mouths and in the proximity of the tidal flats due to the water depth changes. Figure 9c represents the inclination of the ellipse; the patchy areas corresponding to areas (tidal flats) where the magnitude of semi-minor axis increases. These areas can be related to the presence of circulation cells and the production of vorticities due to the change in water depths.

Thus, the same complex point as for water levels is found after the analysis of the currents. The flood-ebb dominance in each tidal creek generate overtides, which can be related with the sediment fluxes.

3.2. Sediment Dynamics within the Tidal Creek Network

The tidal flow patterns allow us to correlate the hydrodynamic changes with the sedimentary patterns. The high sediment transport rates are generally associated with large residual tidal flows. Figure 10 shows the residual flow magnitudes, which were estimated by time-averaging over each tidal cycle. The residual flow is typically an order of magnitude smaller than the most energetic tidal flow. The largest values ($\approx 3$ ${\mathrm{m}}^2/\mathrm{s}$) are observed at both mouths where the bathymetry deepens and has abrupt changes. Values of the residual flow close to zero are observed in the surroundings of MP (convergence zone). The residual circulation is more complex along the main channel due to the irregularity of the margins.

The presence of alternating ebb and flood channels is a general feature of creeks, bays and estuarine systems [87,88]. In particular, at SPC and CC strong tidal currents continually resuspend and rework sediment, specifically during spring tides. Ebb and flood tides can follow mutually evasive channels, and currents may be intense enough to cause scouring at the channel base. Both SPC and CC lack relevant amounts of fluvial sediments and are characterized mainly by the presence of allocutions marine sands and autochthonous mud.

Panels b–d in Figure 10 show the concentration of suspended sand (SSC hereinafter). The highest SSC values at MP are reached in correspondence of the minimum depths ($8\mathrm{mg}/\mathrm{L}$). Moving from MP to CC2, the SSC almost vanish. Panels e–g in Figure 10 show the Mud Suspended Concentration (MSC hereinafter). At MP, the highest values are also reached in the shallower areas ($250\mathrm{mg}/\mathrm{L}$). Moving to CC, the MSC is much higher ($120\mathrm{mg}/\mathrm{L}$) than for the transition to SPC, when this is presumably strongly tidally dependent.

The connection between SPC and CC experienced extensive tidal flat erosion and salt marsh deterioration during the last century [89], as emphasized by the analysis of different bathymetric surveys. Here we analyze the sediment dynamics on the basis of model results (Figure 11), displaying the net evolution of bottom elevations at the end of the numerical simulation. The overall net deposition pattern indicates that SPC and CC receive sediments from ocean and bay, respectively. Most of the sediments are deposited within the main channel. The highest deposition is observed in proximity of MP, where there is a decrease in tidal flow velocity and where the two progressive waves meet. Furthermore, at this location the two progressive waves meet and there is a decrease in tidal flow velocity. Sediments eroded from the mouth of the SPC are transported and deposited within the creeks during the ebb phases (inward), where there is convergence of the residual fluxes. Sediment eroded at the mouth of CC are instead deposited inside the creek during the flood phases (seaward).

4. Discusion

This work investigates how double inlets in tidal creek networks affect the dynamic behavior of these systems through the assessment of surface elevations, currents and sediment transport, which are closely related with tidal asymmetry. The analysis is performed by means of a numerical hydro-morphodynamic model which is able to simulate the effect of both tidal and wind waves. The obtained results show that wind waves are not relevant for the dynamics of the tidal creek networks due to the sheltering effect of the constriction between the outer and the inner basins of the Cádiz Bay (Figure 1) and the shoals present at the mouth of SPC (Figure 11). Furthermore, since human interventions can affect the sediment balance of tidal inlet systems [90,91], here it is also discussed how channel deepening has altered the sediment fluxes, which in turn can impact the ecosystem dynamics.

4.1. Tidal Asymmetry

Non-linear terms associated with friction, advection, and amplitudes play an important role in coastal environments and in particular within the creeks [44,92,93]. As the tidal wave propagates within the bay and through the creek network, it becomes asymmetric because the wave trough is more delayed than the wave crest, as a consequence of the smaller water depth relate to the stronger effect of bottom friction. This asymmetry, that strongly depends on the bathymetry, can produce flood periods shorter than ebb periods, and flood velocities larger than ebb velocities [42,86,94]. The areas with the highest variation in tidal amplitude are frequently the ones where geometric variations associated with marsh loss are more pronounced [95].

The tidal asymmetry can be addressed by analyzing the ratio between the amplitudes of the M2 and M4 constituents (a${}_{M4}$/a${}_{M2}$) (a—Figure 12). This ratio is related to the morphological characteristics and configuration of creeks [96]. [40] showed that estuaries with (1) low freshwater input, (2) significant ratio of tidal amplitude to channel depth, and (3) long narrow channels, the ratio amplitude to tidal asymmetry (a${}_{M4}$/a${}_{M2}$) increases as the tide propagates toward the inner estuary. Thus, a larger value of a${}_{M4}$/a${}_{M2}$ involves a greater relevance of tidal asymmetry. For the SPC–CC system, the minimum values for ratio between amplitudes is observed at the mouth of both creeks (≈0.02 at SPC and ≈0.04 at CC), whereas the maximum is obtained at MP (≈0.12) where the effect of the tidal asymmetry is stronger and areas of significant deposition were found (Figure 11). The relationship between hydrodynamic and morphodynamic variables was shown by [17], where the highest values of a${}_{M4}$/a${}_{M2}$ corresponded also with an increase of sedimentation.

The relative phase difference $2{\varphi }_{\mathrm{M}2}-{\varphi }_{\mathrm{M}4}$ is another important parameter to define tidal asymmetry [44] (b—Figure 12). This parameter is smaller than ${90}^{°}$ at CC and smaller than ${200}^{°}$ at SPC, thus suggesting that the tide has ebb-dominant characteristics in CC and flood-dominant characteristics in SPC. The tidal asymmetry is considered positive if the flow direction is from the ocean towards the inner basin, and negative otherwise (Figure 13). These results are consistent with those emerging from the harmonic analysis (see Section 3.1.1) suggesting the occurrence of a relevant tidal asymmetry. Thus, the current distribution is mostly tide-induced. As a consequence, the local dominance of ebb or flood currents triggers the residual sediment flux [17,38,39,42].

Finally, tidal currents can also be used to analyze tidal asymmetry. The magnitude of the asymmetry is defined as:

(1)${\mathrm{Ass}}_{qi}={\displaystyle \frac{max\mid q_i{\mid }_{flood}-max\mid q_i{\mid }_{ebb}}{max\mid q_i{\mid }_{flood}+max\mid q_i{\mid }_{ebb}}}$

with $i=x,y$. This tidal asymmetry ratio was computed along the creeks considering either spring or neap tides (Figure 14). The color-coded representation of the asymmetries suggests that the largest changes of the flow rates per unit width occur during neap tides, whereas asymmetries vanish during spring tides. The ebb dominance prevails at MP and CC for neap tides (blue color in panels b1–b2 in Figure 14). However, the flood dominance prevails at SPC (red color in panels b1–b2 in Figure 14).

4.2. Impact of Channel Deepening

Human interventions, such as infrastructure construction or channel deepening, among others, can disrupt the sediment balance of tidal inlet systems. Here, it is shown that the channel deepening carried out in SPC can potentially enlarge the problem of the future closure of the creeks at MP showed in Section 3.2. For that reason, a scenario with the deepening of the mouth of SPC, which has been a recurring intervention at the area during last decades for navigational reasons, was simulated. The depth in this area was changed from 6 up to 8 m, which are the typical values for this intervention. Figure 15 shows that the original evolution of sediment change significantly after the channel deepening. A problem of great concern is related to the possible closure of the creeks at MP by an accumulation of sediments. This effect was identified both with digital images (Figure 2) and the results of the numerical simulations (Figure 11). According to the results, the dredging of SPC mouth implies morphodynamic variations far away from the dredged area, enhancing the infilling of the SPC–CC system especially at MP. This point is located almost 10 km away along a complex creek network. This infilling may induce a potential negative impact on the ecological and morphological evolution of the Cadiz bay Natural Park (Figure 15).

Hence, bathymetric changes can modify the hydrodynamics of tidal creeks triggering changes in their physical characteristics, such as tidal amplitude or sediment fluxes. These variations can potentially alter ecosystems through variations in salt intrusion [91] or vertical exchange of oxygen [97].

5. Conclusions

Improving current understanding of hydrodynamics and sediment circulation of creek networks located in complex tidal systems is of major importance to face current and future challenges derived from climate change and an increasing human pressure. The aim of this work is to deepen the knowledge of the hydro-morphodynamics of complex creek networks that connect basins with very different characteristics, identifying their morphodynamic trends and the potential impacts of human interventions such as channel deepening. In this case, we selected the SPC–CC system, that connects directly the Atlantic Ocean with the inner basin of the Cádiz Bay, to study the interactions between the tidal waves that penetrate the creeks from these two different water bodies, analyzing the tidal asymmetry and the morphodynamic trend of the system. Towards this goal, we set up a hydro-morphodynamic model specifically developed for areas with very shallow and complex channels [62]. Our results can be extended to other study sites where these types of creek networks are present, such as Ria Formosa, Portugal [98].

The importance of the analyses carried out in the present study lies in the particular configuration of areas such as SPC and CC, which are influenced by the open ocean and an inner bay due to their double connection that makes tidal waves to propagate from both ends with different amplitudes, phases and flow velocity. Both creeks import, export and deposit fine sediment sourced from both landward and seaward.

Our analyses show a decrease in the M2 tide amplitude and an amplification in the phase when the wave penetrates through the creeks. A critical and complex point from a hydro-morphodynamic point of view appears at MP, where a maximum in the phases of the tide constituents occurs. This location corresponds to the convergence zone of the tidal wave.

Our results further show that tidal currents at SCP and CC are distorted from the sinusoidal form of their astronomical forcing, thus suggesting a transfer of energy from M2 to M4. The main consequence of this generation of overtides is the strengthening of flood currents compared to ebb currents in the case of the CC, whereas the opposite behavior is observed in the SPC.

This flood-ebb asymmetry plays an important role in sediment circulation in SPC and CC. Also, the decrease/increase of the tidal current have strong effect in the sediment transport. Thus, a decrease (increase) in tidal currents and a strong (weak) presence of tidal asymmetry cause sedimentation (erosion), which is focused on MP (bank of the creeks).

At both creeks tidal current variations are influenced by the creek geometry. Sediment is also redistributed within the tidal creeks on a regular basis by tidal currents. During the flooding periods large amounts of sediment are transported into the creeks. Two different types of creeks can be identified regarding their behaviors: (1) deep creeks along the main channel; and (2) shallow channels characterized by a landward decreasing depth, where the tidal flats predominate promoting the residual import of fine sediment. The first type is observed at SPC, where flood dominance develops and water fluxes propagate landward starting from the channel mouth. The second case is more common at CC, with water flowing from the creeks to the inner bay. In the connection between SPC and CC (MP), there is a decrease in strength of the tidal flow from the main channel to the shallower area. In addition, the tidal flat can induce a seaward flux of coarse suspended sediment.

Finally, the results also show that although the effect of wind waves on the tidal creek network dynamics is negligible, the channel deepening carried out periodically in SPC can potentially enlarge the problem of the future closure of the creeks at MP, since it tends to enhance the infilling of the SPC–CC system causing a potential negative impact on the ecological and morphological evolution of this complex system.

Author Contributions

Conceptualization, A.D. and L.C.; Data curation, C.Z.; Methodology, C.Z., A.D., L.C.O. and A.L.-R.; Supervision, A.D., L.C., M.D.-M. and M.O.-S.; Writing – original draft, C.Z. and A.L.-R.; Writing—review & editing, A.D., L.C., A.L.-R. and M.O.-S.

Funding

This work was funded by the Cádiz Bay Port Authority, the Department of Innovation, Science and Business of the Andalusian Regional Government (Project P09-TEP-4630), the Spanish Ministry of Economy and Competitiveness through the Projects CTM2009-10520-MAR, CTM2017-89531-R (PIRATES), BIA2015-65598-P (VIVALDI) and PCIN-2017-108, and by the “Programa Iberoamericano de Ciencia y Tecnología para el Desarrollo”, CYTED (project PROTOCOL 917PTE0538). The work of the first author was partially funded by the Andalusian Regional Government, Research Grant RNM-6352.

Acknowledgments

The authors would thank to the staff of the Environmental Fluid Dynamics research group (GDFA, University of Granada) for their support during the field campaign.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

CC	Carracas Creeks
MP	Middle Point
SPC	Sancti Petri Creeks

Figures and Tables

Figure 1. (A) Location of the Cádiz Bay in Spain; (B) Computational grid and color coded bathymetry of the Cádiz Bay with the locations of the stations (labels I1–I4 correspond to current profile meters (marked with white crosses), and labels T1–T3 to tidal gauges (marked with white circles) where data were collected). Dashed color lines correspond to the five boundaries of the mesh. The Cádiz Bay Natural Park includes the square yellow (more detail in panel (C) and the Trocadero Island (blue square); (C) Detailed bathymetry of the creek network.

Figure 2. Morphological changes that have affected the CC and SPC networks from 1956 to 2014. Portions in which landscape erosion and deposition have been observed are in different colors. T2 and T3 refer to tidal gauges.

Figure 3. Water levels and East and North velocities measured during the field survey with the current meters (I1–I4) and water levels measured with the tidal gauges (T1–T3).

Figure 4. Distribution of the sediment along the SPC and CC.

Figure 5. Observed vs modelled tidal water levels, for stations T2 (a) and T3 (b). Both water level data used for calibration (empty dots) and testing (solid circles) are shown, for the two different phases (red symbols refer to the ebb phase, blue symbols correspond to the flood phase).

Figure 6. Results of the tidal harmonic analysis. Panels a and b show a color-coded representation of the modelled amplitude (a) and phase (b) fields for the M2 component, whereas panels (c,d) portray the same fields for the M4 component. The panels (a1–d1) show tidal amplitudes and phases for the flood and ebb phase along the path represented by the with dashed lines in panels (a–d). For the sake of clarity water depths along the white dashed path are also shown in green in the same panels.

Figure 7. (a,b) Along channel variations of the flow rates per unit width qx and qy along the x,y coordinates, respectively, when the maximum flood and ebb conditions during the simulation time are attained during spring and neap tides (qx is positive towards the East, qy is positive towards the North). Spatial distribution of qx during maximum flood for spring (a1) and neap (b1) tides, and spatial distribution of qy during maximum flood for spring (a2) and neap (b2) tides. (c1,d1; c2,d2) same as above for the maximum ebb tide. Orange, blue, red and green dashed lines correspond to flood spring, flood neap, ebb spring and ebb neap, respectively.

Figure 8. Results of the tidal-current harmonic analysis for the M2 component. The magnitude of qM2 (m2/s) and the phase of qM2 (°) are labeled in the panels (a,b), respectively. Panel (c) corresponds to the difference between the phase associated with the velocities and the water levels for M2 (°).

Figure 9. Tidal ellipse parameters of the flow rate. The semi-major axis MaM2 (m2/s), semi-minor axis maM2 (m2/s), inclination φM2 (0° E–90° N) and phase ϕM2 (°) fields are shown in panels (a–d), respectively.

Figure 10. (a) Color coded representation of the residual flow magnitudes (m2/s) at maximum flood period. The sand suspended concentration (mg/L; panels (b–d)) and the mud suspended concentration (mg/L; panels (e–g)) are shown for 3 cross-sections (North, central and south, defined in panel (a). Vertical axes on panels (b–g) refer to the along-channel position of 10 points equally distributed along the areas highlighted in panel (a).

Figure 11. Net evolution of the bottom elevation (mm) computed by the bed evolution model at the end of the simulation.

Figure 12. Hydrodynamic variables along the navigation channel: (a) aM4/aM2 and (b) 2ϕM2+ϕM4. Left panel: transect N-S (dashed line).

Figure 13. Criteria of sign in Cádiz creeks. Considering positive (negative) if the flow direction is from the ocean to the basin (from the basin to the ocean).

Figure 14. Spatial distribution of the asymmetry magnitude Assqx (panels (a1,b1)) and Assqy (panels (a2,b2)) along the creeks for spring (a) and neap (b) tides, respectively.

Figure 15. Net evolution of the bottom elevation (mm) computed by the bed evolution model at the end of the simulation considering the deepening of the navigation channel at SPC. Insets represent the bathymetry before and after the dredging.

Table 1

Distribution of the Strickler bed roughness coefficient K${}_s$ along the creeks.

	Navigation Channel	Shallower Area	Flood Plain
K${}_s$ (m${}^{1/3}$/s${}^{-1}$)	10	25	40

Table 2

Parameters used in the WWTM and STABEM simulations.

Parameter Description	Value
Mean grain diameter for pure sand	200 $\mathsf{\mu }$m
Mean grain diameter for pure mud	20 $\mathsf{\mu }$m
Sediment porosity	0.4
Critical shear stress for pure sand	0.4 Pa
Critical shear stress for pure mud	0.8 Pa
Critical shear stress for deposition	1 Pa
Specific entrainment parameter for pure mud	$1\times {10}^{-5}$ gm/s
Specific entrainment parameter for non-cohesive mixture	$2\times {10}^{-5}$
Initial active layer thickness	0.02 m
Wind drag coefficient	0.003
Settling velocity for sand	0.025 m/s
Settling velocity for mud	0.00034 m/s

Table 3

Root mean square errors ($RMSE$), correlation coefficients (R) and skill coefficients (S), for the calibration period of the elevations ($I1$, $I2$, $T2$ and $T3$) and velocities ($I1$ and $I2$).

	Tidal Level $\left(\mathbf{m}\right)$		East Velocity $(\mathbf{m}/\mathbf{s})$		North Velocity $(\mathbf{m}/\mathbf{s})$
	$\mathit{R}$	$\mathit{S}$	$\mathit{R}$	$\mathit{S}$	$\mathit{R}$	$\mathit{S}$
$I_1$	0.99	0.99	0.94	0.83	0.96	0.91
$I_2$	0.99	0.99	0.95	0.91	0.95	0.91
$T_2$	0.88	0.84	-	-	-	-
$T_3$	0.99	0.97	-	-	-	-