Large-scale coastal and fluvial models constrain the late Holocene evolution of the Ebro Delta

Abstract

The distinctive plan-view shape of the Ebro Delta coast reveals a rich morphologic history. The degree to which the form and depositional history of the Ebro and other deltas represent autogenic (internal) dynamics or allogenic (external) forcing remains a prominent challenge for paleo-environmental reconstructions. Here we use simple coastal and fluvial morphodynamic models to quantify paleo-environmental changes affecting the Ebro Delta over the late Holocene. Our findings show that these models are able to broadly reproduce the Ebro Delta morphology, with simple fluvial and wave climate histories. Based on numerical model experiments and the preserved and modern shape of the Ebro Delta plain, we estimate that a phase of rapid shoreline progradation began approximately 2100 years BP, requiring approximately a doubling in coarse-grained fluvial sediment supply to the delta. River profile simulations suggest that an instantaneous and sustained increase in coarse-grained sediment supply to the delta requires a combined increase in both flood discharge and sediment supply from the drainage basin. The persistence of rapid delta progradation throughout the last 2100 years suggests an anthropogenic control on sediment supply and flood intensity. Using proxy records of the North Atlantic Oscillation, we do not find evidence that changes in wave climate aided this delta expansion. Our findings highlight how scenario-based investigations of deltaic systems using simple models can assist first-order quantitative paleo-environmental reconstructions, elucidating the effects of past human influence and climate change, and allowing a better understanding of the future of deltaic landforms.

Introduction

The Ebro Delta, Spain, with its distinctive plan-view shape, has experienced significant morphologic changes over the last millennia caused by the growth and reworking of different delta lobes (Fig. 1) (Canicio and Ibáñez, 1999). While autogenic processes might have caused some of these morphological changes, others aspects could be attributable to past climate changes or anthropogenic activities within the drainage basin. Many different scenarios leading to the modern morphology have been proposed, including high-frequency (centennial scale) sea-level fluctuations (Somoza et al., 1998), human-induced sediment load changes in the Ebro River (Guillén and Palanques, 1997a), and climate fluctuations affecting river discharge (Xing et al., 2014).

Reconstructed morphologic development of the Ebro Delta, modified from Canicio and Ibáñez (1999).

[Figure omitted. See PDF]

Many deltas around the world have experienced substantial morphologic changes over the last millennia due to anthropogenic factors such as river damming, land-use change, and climate change (Anthony et al., 2014; Giosan et al., 2012; Maselli and Trincardi, 2013; Syvitski and Saito, 2007). The Ebro Delta lends itself particularly well to quantitative reconstructions because it is morphologically constrained (Nelson, 1990), it displays a distinctive plan-view shape (Fig. 1), and its environment is relatively well-studied (Cearreta et al., 2016; Maldonado, 1975). Here, we use a coastline evolution model and a river profile evolution model to quantitatively constrain the style, timing, and rate of morphologic change and the associated fluvial transport to the Ebro Delta during the late Holocene.

Our goal in this paper is to investigate the general evolution of a wave-influenced delta–river system using “scenario-based” and quantitative model experiments. We do not attempt to capture the precise morphology or geochronology of any one segment of the Ebro Delta, but rather approximate delta paleo-morphodynamics to assess the potential physical mechanisms that could have formed this delta plain. Our scenario-based approach has two objectives: (i) to investigate whether we can reproduce the broad morphology of the Ebro Delta plain with simple models and available data on fluvial and wave climate histories, and, if possible, then (ii) use simple models to quantify first-order sediment fluxes and timescales. As a test of the suitability of the delta and the river models, we compare the model predictions to observed deltaic (Jiménez and Sánchez-Arcilla, 1993) and fluvial change (Vericat and Batalla, 2006) over the last century. Overall, our approach allows us to test existing hypotheses of environmental changes that may have influenced the development of the Ebro Delta.

Background

Ebro River

The Ebro River reached the Mediterranean Sea, after an endorheic phase, sometime between 13 and 5 million years ago (Babault et al., 2006; Garcia-Castellanos et al., 2003). Its modern drainage basin extends over 85 530 km${}^{\mathrm{2}}$, covering a large portion of the Pyrenees, the Cantabrian mountains, and the Iberian Massif (Mikeš, 2010). The average channel width in the lower course of the river is $\sim$ 150 m, with a bankfull flow depth of $\sim$ 5 m (Guillén and Palanques, 1997a). Average (pre-dam) discharge has been estimated at about 500 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ (Batalla et al., 2004). The fluvial sediment flux during the Holocene highstand, based on radiometric dating of Ebro continental shelf deposits, is estimated to be ca. 200 kg s${}^{-\mathrm{1}}$ (6.3 Mt yr${}^{-\mathrm{1}})$ (Nelson, 1990). The suspended load consists mostly of clay and silt (Muñoz and Prat, 1989), while the bedload is predominantly sand and gravel (Vericat and Batalla, 2006).

Ebro Delta

At the Ebro River outlet to the Mediterranean Sea, fluvial sediment deposition over the course of millions of years has expanded the Ebro continental shelf and constructed successive deltas (Babault et al., 2006; Nelson, 1990). During the Holocene, strong waves and limited coarse-grained sediment input have shaped the Ebro coast towards a wave-dominated deltaic morphology with a smooth shoreline and single thread distributary network (Jiménez et al., 1997). The Ebro nearshore zone consists mostly of sand-size sediment (Maldonado, 1975; Somoza et al., 1998) to a depth of $\sim$ 12 m, transitioning into muds farther offshore (Guillén and Palanques, 1997b). Two distinctive features on the Ebro Delta plain are the spits to the north (El Fangar) and south (La Banya) of the current river mouth, considered to be formed by wave reworking of abandoned delta lobes (Fig. 1) (Maldonado, 1975).

Ebro Delta Holocene evolution

Similar to many other deltas around the world, Holocene sea-level rise led to the transgression of the previous Pleistocene delta deposit (Maldonado, 1975). The maximum flooding surface of the Ebro Delta is dated to about 6900 years BP, with its landward extent near the town of Amposta (Lario et al., 1995; Somoza et al., 1998). Several studies have interpreted historical references to suggest that the Ebro was still an estuary $\sim$ 2000 years ago (Guillén and Palanques, 1997a; Maselli and Trincardi, 2013); however, radiocarbon dating of relict, arcuate, beach ridges on the delta plain (Canicio and Ibáñez, 1999) and recently dated sandy beach shells from boreholes (Cearreta et al., 2016) indicate that the delta was already formed by $\sim$ 6000 years BP.

Dated beach ridges show that the Ebro Delta coast was small, cuspate, and wave-dominated at least until 3000 years BP (Canicio and Ibáñez, 1999). The same study suggested that sometime between 1400 and 1000 years BP the Riet Vell lobe had grown rapidly and extended approximately 20 km into the Mediterranean Sea, although no confirming dates are currently available. This increase in progradation rate, at least 2 to 3 times faster than previous and initiating sometime after 3000 but before 1400 years BP, is commonly ascribed to an increase in fluvial sediment supply (Canicio and Ibáñez, 1999).

What could have caused this increase in fluvial sediment supply? Thorndycraft and Benito (2006) examining periods of fluvial flooding in Spain suggest, although with limited data, that extensive fluvial flooding occurred prior to 9000 yr BP and after 3000 years BP, with most records from the period in between $\sim$ 9000 and $\sim$ 3000 years BP pointing to floodplain forestation and low-energy floodplain deposition. Floodplain alluviation of Spanish rivers after 3000 years BP suggests three periods of intense flooding: 2710–2320, 2000–1830, and 910–500 years BP (Benito et al., 2008). The first of these three periods has been associated with large-scale climate variability causing increased flooding. The last period of floodplain aggradation, however, is not in phase with paleoflood records, which Benito et al. (2008) therefore attributed to anthropogenic modifications such as deforestation that increased the Ebro River sediment load. Other deltas around the Mediterranean and the Black Sea, whose hinterlands have comparable observed land-use change histories, also show periods of increased progradation in response to human activities (Anthony et al., 2014; Giosan et al., 2012; Maselli and Trincardi, 2013).

Xing et al. (2014) used the long-term fluvial discharge and sediment supply model HydroTrend (Kettner and Syvitski, 2008) to quantify the effects of anthropogenic and climate change on fluvial suspended sediment supply to the Ebro Delta. HydroTrend uses empirical relations between basin area, land-cover, drainage basin relief, temperature, and precipitation and is calibrated using modern sediment transport records to simulate river sediment load. The model results of Xing et al. (2014) suggest that discharge variation was mostly a result of climatic variability, whereas forest clearing likely contributed to changes in suspended sediment load. Their study estimated a 40 % increase in the fluvial suspended sediment load in response to deforestation. Other studies, such as Nelson (1990) and Guillén and Palanques (1997a), who reconstructed a sediment budget from delta plain and shelf aggradation rates, have estimated a greater fluvial sediment flux increase of 350 %.

The increase in the Ebro Delta plain progradation first formed the Riet Vell lobe (Fig. 1) (Canicio and Ibáñez, 1999). From relict channel deposits on the delta plain (Maldonado, 1975), published maps, and historical evidence, Canicio and Ibáñez (1999) and Somoza and Rodriguez-Santalla (2014) suggest that the progradation of the Riet Vell lobe stopped after 1000 years BP but prior to 600 years BP, when the avulsion of the main channel started the new Sol de Riu lobe to the north (Fig. 1), which also prograded rapidly. Subsequently, the Riet Vell lobe was reworked into the La Banya spit to the south. After a second river avulsion about 300 yr ago to the Mitjorn-Buda lobe in between the previous active channels, the Sol de Riu lobe was also abandoned and reworked into the northern El Fangar spit (Fig. 1).

Recent changes

As of today, over 187 dams have been built in the Ebro that have highly regulated its discharge and currently impound 57 % of the mean annual runoff (Batalla et al., 2004). The average fluvial water discharge based on hydrographic records before dam construction was approximately 500 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, while post-dam discharge has averaged about 340 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ (Batalla et al., 2004).

Prior to the construction of the major dams in the Ebro, peak discharge was about 50 % higher than modern flows (Batalla et al., 2004). As a consequence, while bedload-transporting river flows ($\sim$ 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}})$ were previously exceeded 15 % of the time, dams reduced the exceedance frequency of these floods to just 4 % of the year and thereby lowered the bedload sediment flux at the delta outlet (Vericat and Batalla, 2006). Additionally, reservoirs behind dams trap about 90 % of the upstream suspended sediment load and 100 % of the upstream bedload (Vericat and Batalla, 2006).

Accurate measures of pre-dam fluvial sediment flux are challenging, but from early 20th century sediment concentration and discharge measurements, Guillén and Palanques (1992) obtained an annual average suspended load estimate of $\sim$ 600 kg s${}^{-\mathrm{1}}$ (20 Mt yr${}^{-\mathrm{1}})$. Syvitski and Saito (2007) used Bagnold's (1966) equation to estimate a pre-dam bedload flux of 71 kg s${}^{-\mathrm{1}}$ (2.2 Mt yr${}^{-\mathrm{1}})$.

Post-dam measurements taken 50 km upstream of the delta (25 km downstream of the Flix Dam, the last major dam in the main river channel) estimate a modern total sediment load of about 28 kg s${}^{-\mathrm{1}}$ (0.9 Mt yr${}^{-\mathrm{1}})$, of which 40 % is transported as bedload (Vericat and Batalla, 2006). Using predictive sediment transport formulae from van Rijn (1984) combined with discharge measurements, Jiménez (1990) estimated the modern sand (bedload) transport at the mouth of the delta at 1.6 kg s${}^{-\mathrm{1}}$ (0.05 Mt yr${}^{-\mathrm{1}})$. A comparison of estimates of pre-dam to post-dam bedload transport to the delta suggests a reduction of about $\sim$ 95 %. Evidence of this sediment deficit include scours of the lower course of the channel bed and the formation of armored layers. Immediately downstream of the Flix Dam, the channel bed surface consists of coarse gravel ($D_{\mathrm{50}}=$ 38 mm) while the subsurface consists of mixed sand and gravel ($D_{\mathrm{50}}=$ 17 mm) (Vericat et al., 2006).

The 20th century fluvial sediment flux reduction has also led to morphologic changes of the delta at the coast. While for much of the last millennia the Ebro Delta's mouth likely exhibited, at least periodically, close to a river-dominated morphology, the sediment supply reduction has led to a more wave-dominated form of the modern Mitjorn-Buda lobe (Jiménez et al., 1997).

River damming may not be the only cause for large-scale coastal changes in the future. A bathtub-style estimate projected that subsidence and sea-level rise may submerge about 40 % of the delta plain by 2100 (Ibáñez et al., 1997). However, the projected effects of sea-level rise on coastal change up to 2050 are negligible compared to ongoing change resulting from alongshore sediment transport gradients (Sánchez-Arcilla et al., 2008). These gradients have caused retreat rates of 50 m yr${}^{-\mathrm{1}}$ near the river mouth, and have resulted in spit accretion at rates of approximately 10 m yr${}^{-\mathrm{1}}$ between 1957 and 1992 (Jiménez and Sánchez-Arcilla, 1993).

Modeling wave-influenced deltas

Many numerical models have been developed over the last decades to quantitatively reproduce, predict, and understand the dynamics of deltaic systems. Complex “simulation models” such as Delft3D typically are used to reproduce a particular well-constrained natural environment (e.g., van der Wegen et al., 2011) or to parameterize poorly understood physical processes (e.g., Nienhuis et al., 2016a). Simple “exploratory models” of “reduced complexity”, on the other hand, are designed to capture the essential feedbacks leading to an observed phenomenon (Murray, 2003). Because, in the long term, the millennial- to centennial-scale development of the Ebro Delta is poorly constrained, we apply exploratory models of wave-influenced delta dynamics to capture the essential physical mechanisms affecting the evolving morphology of the Ebro Delta using scenario-based approaches.

The plan-view shape of the Ebro Delta, like other wave-dominated deltas, is governed by wave-driven alongshore sediment transport (Bakker and Edelman, 1964; Bhattacharya and Giosan, 2003; Jiménez and Sánchez-Arcilla, 1993). Modeling of wave-dominated delta shape is therefore usually performed with coastline models (e.g., Ashton and Giosan, 2011; Bakker and Edelman, 1964). By assuming the cross-shore profile maintains a constant shape, gradients in alongshore transport can be linearly related to accretion or erosion of a single contour line, typically the coastline. Such one-contour-line models calculate alongshore sediment transport based on surf-zone-averaged equations such as the CERC formula (Komar, 1971), which relate the relative wave angle and height to a sediment transport flux. The cuspate coastline shape typical of wave-influenced deltas arises when adding a point source of (fluvial) sediment to an otherwise straight sandy coast (Grijm, 1960).

By comparing fluvial and wave-driven sediment fluxes, Nienhuis et al. (2015) quantified when deltas would be expected to attain a wave-dominated versus a river-dominated shape. If the fluvial sediment supply is larger than the maximum potential alongshore sediment transport away from both delta flanks, waves cannot transport fluvial sediment delivered at the river mouth alongshore and a delta would be expected to be river-dominated. Their study defined a river dominance ratio $R$ as the fluvial sediment flux ($Q_{\mathrm{r}})$ divided by the maximum alongshore sediment transport flux away from the river in both directions ($Q_{\mathrm{s},\max })$. For $R$ > 1, the delta is river-dominated. If $R$ < 1, there is an equilibrium plan view delta flank orientation such that the fluvial sediment flux ($Q_{\mathrm{r}})$ equals the wave-driven sediment flux ($Q_{\mathrm{s}})$ away from the river mouth along both flanks. The amount of shoreline deflection at the river mouth of a wave-dominated delta is therefore an indicator of its wave dominance, with flatter coasts being more wave-dominated.

Ashton and Giosan (2011) showed that for very obliquely approaching waves, wave-dominated deltas can become asymmetrical and develop downdrift migrating sand waves and spits on the downdrift flank. These shoreline instabilities can form on growing deltas in which case they are oriented roughly parallel to the delta flank. Nienhuis et al. (2013) later showed that prominent recurved spits can develop from the reworking of delta lobes. These recurved spits develop after a reduction in fluvial sediment supply to a delta lobe (e.g. due to avulsion or dam construction) only if one or both flanks of the delta grew past the maximum in alongshore sediment transport. These recurved spits are generally not oriented parallel to the delta coastline. Rather, the orientation of free spits is controlled by the wave climate and the rate of delta lobe retreat (Ashton et al., 2016; Nienhuis et al., 2013).

Modeling fluvial sediment supply

The sand-sized sediment feeding the Ebro Delta is supplied as bedload and suspended load through the Ebro River, interacting with the alluvial river bed (Jiménez et al., 1990). In alluvial rivers, channel–bed interaction sets up an equilibrium between the along-stream slope, river discharge, and sediment supply (Lane, 1955). One of the first attempts to numerically model fluvial sediment transport was by Hirano (1971), who combined the depth-averaged, one-dimensional Saint-Venant equations for fluid flow with a simple formulation for sediment transport. Their model resulted in a typical concave up longitudinal river profile for a scenario of gradually increasing water discharge downstream (Hirano, 1971; Snow and Slingerland, 1987).

If one assumes that changes in bed elevation are less pronounced compared to changes in flow characteristics, the flow can be approximated as steady (de Vries, 1965). If it is further assumed that the flow is locally uniform (spatial changes are small compared to the flow), then the steady flow becomes quasi-normal and an alongstream momentum balance relates bed shear stress to water depth and bed slope. Combined with an Exner equation for sediment conservation and a Chezy or Manning coefficient for form drag, the Saint-Venant equations for fluvial flow can then be reduced to a simple analytical expression for longitudinal river profile shape and equilibrium sediment transport rates (Parker, 1978).

Methods

Coastline Evolution Model

We study the morphologic evolution of the Ebro Delta plain using the Coastline Evolution Model (CEM), an exploratory, process-based one-contour-line model (for a full description see Ashton and Murray, 2006). In this model, the plan-view coastal zone is discretized into 50 m square cells that are either filled (land), empty (water), or partially filled (coastline), the latter allowing for a smooth, continuous shoreline. Incoming deep-water waves are refracted and shoaled across parallel contours from the toe of the shoreface up to the breaking wave depth. Breaking wave characteristics are then used to compute alongshore sediment transport ($Q_{\mathrm{s}}$, kg s${}^{-\mathrm{1}})$ with the CERC formula (Komar, 1971). Alongshore sediment transport is adjusted with a factor 0.6 based on calibration studies of Jiménez and Sánchez-Arcilla (1993).

Following the one-contour-line approach, the divergence of alongshore sediment transport is related to shoreline accretion or erosion up to the shoreface depth using the shoreline Exner equation, $${\displaystyle \frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{1}}{D_{\mathrm{sf}}}}{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-p)\cdot {\mathit{\rho }}_{\mathrm{s}}}}{\displaystyle \frac{\mathrm{d}{Q}_{\mathrm{s}}}{\mathrm{d}x}},$$ where d$\mathit{\eta }/\mathrm{d}t$ is erosion or progradation of the shoreline (m s${}^{-\mathrm{1}})$, $D_{\mathrm{sf}}$ is the shoreface depth (m), and d$Q_{\mathrm{s}}/\mathrm{d}x$ is the alongshore gradient in alongshore sediment transport (kg s${}^{-\mathrm{1}}$ m${}^{-\mathrm{1}})$. The density of sediment is ${\mathit{\rho }}_{\mathrm{s}}$ (kg m${}^{-\mathrm{3}})$ and $p$ is the dry mass void fraction.

As Eq. (1) shows, the shoreface depth $D_{\mathrm{sf}}$ represents an important scaling parameter for coastline change rates. In a study of short-term (decadal) coastal change of the Ebro Delta, Jiménez and Sánchez-Arcilla (1993) suggest a 7 m shoreface depth based on cross-shore profile variability. Research suggests, however, that across longer timescales, the shoreface depth increases as a result of lower frequency (storm) events (Hands, 1983). As potential indicator, Guillén and Palanques (1997b) found that the sand–mud transition of the Ebro Delta is located at approximately 12 m water depth based on bed-surface samples. In our centennial-timescale modeling of the Ebro Delta we take advantage of a recent quantitative analysis of shoreface evolution (Ortiz and Ashton, 2016), which suggests morphological response rates may set the effective shoreface closure depth. For 1 m wave heights, a 100-year depth of closure is approximately 40 % deeper than a 10-year timescale depth of closure. In our model, we therefore choose a shoreface depth of 10 m.

The characteristic shoreface slope (0.01) and shelf slope (0.002) are set based on the geometry offshore of the Ebro (Guillén and Jiménez, 1995; Jiménez and Sánchez-Arcilla, 1993). The inclusion of a shelf slope in CEM makes the delta plain prograde slower as it builds out into deeper water further from the coast (Ashton and Murray, 2006). Shoreline retreat maintains a minimum shoreface depth of 10 m. This approach results in a more realistic mass balance, yet does not fully capture potential long-term shoreface dynamics; the latter would be difficult without appropriate centennial-scale measurements of shoreface dynamics.

An advantage of the CEM is its ability to produce arbitrarily sinuous shoreline shapes such as spits. When shoreline erosion causes a neck of a spit to reach a critical width, overwash occurs and sediment is transported from the shoreface to the backbarrier to maintain a minimum width (Jiménez and Sánchez-Arcilla, 2004). Overwash allows spits and barriers to retreat without disconnecting from the rest of the coastline (Ashton and Murray, 2006). Following observations of Jiménez and Sánchez-Arcilla (2004) of the La Banya spit, we set the critical barrier width to 250 m. The overwash depth is determined geometrically assuming a shoreface slope (0.01) and an overwash volume. Even though this is obviously a simplification that could result in overwash depths that are unrealistically deep, it avoids the need for a complicated assessment of backbarrier elevations coastwide.

CEM does not have the ability to incorporate base-level changes in its shoreline change estimates. Surface elevation tables on delta topset deposits indicate a relatively high relative sea-level rise rate of about $\sim$ 3 mm yr${}^{-\mathrm{1}}$ (Ibáñez et al., 1997). Although relative sea-level rise rates in the sandy delta foreset deposits were likely significantly lower, we cannot rule out their potential effects on Ebro Delta change.

Application to the Ebro Delta

We have adapted the CEM to model growth and reworking of the different Ebro Delta lobes. Rather than growing perpendicularly to the initial coastline, we force individual channels to grow along channel paths that we choose based upon the paleo- and modern channels of the Ebro Delta (Fig. 2; Maldonado, 1975). The first lobe builds out at 5${}^{\circ }$ from shore normal and represents the growth of the Riet Vell. The second (Sol de Riu) lobe grows $-$45${}^{\circ }$ from shore normal, and the modern Mitjorn-Buda lobe is oriented at $-$20${}^{\circ }$. For all these lobes, the river channel is highly simplified and only represented as the location alongshore where the littoral-grade portion of the fluvial sediment is deposited. By modeling the mass balance this way, we assume that fine-grained fluvial sediment is winnowed by waves and eventually deposited largely offshore beyond the shoreface (Guillén and Palanques, 1997b).

Schematic of modeling scenario, highlighting the succession and orientation of Ebro Delta lobes, shown on top of the modern morphology (NASA Landsat image) and the inferred paleo-channels (dotted lines, from Maldonado, 1975). In the model, the straight reference coastline is assumed to be non-erodible. Names refer to the spits and the lagoons on the Ebro Delta. Numbers refer to the (1) Riet Vell, (2) Sol de Riu, and the (3) Mitjorn-Buda lobes.

[Figure omitted. See PDF]

As a second modification to the original model, we disable alongshore sediment transport out of a cell that is part of the initial coastline. This modification accounts for the fact that the Ebro Delta juts out of the rocky coastline of Mediterranean Spain, and is not connected to an updrift littoral sediment source (Fig. 2).

Even though the Ebro Delta channel orientations are likely in part determined by wave climate, fluvial sediment supply, and alongshore sediment bypassing of the river mouth (Nienhuis et al., 2016b), we choose to impose channel orientations directly to constrain model variability. Similarly, channel avulsion has been suggested to be controlled by backwater length and channel filling timescales (Chatanantavet et al., 2012). To limit model sensitivity we do not allow autogenic river avulsions in our model, instead we model avulsions at their geologically inferred locations (Maldonado, 1975). Using the channel orientations of the existing delta, we then run 42 simulations with varying avulsion times to determine which scenarios of avulsion timing can best match the current delta plain morphology.

It is important to note that we are not explicitly simulating the history of the Ebro Delta plain; rather we use simple models to constrain fluvial sediment fluxes and delta growth in a broadly representative wave-dominated environment. Because both fluvial and coastal models are exploratory and there is no feedback, we do not couple the two models directly. We run scenarios of different fluvial sediment supply rates to investigate how sediment delivery rates affect delta morphology, including the growth of spits. We also run scenarios of different channel avulsion timings and compare the resulting modeled delta shape to the modern Ebro Delta plain shape to constrain Ebro Delta chronology. See Table 1 for an overview of the model parameters.

Overview of the model parameters chosen for the fluvial and the coastal models.

Name	Value	Units	Note
CEM (Ashton and Murray, 2006)
$K_{\mathrm{1}}$	0.035	m${}^{\mathrm{3}/\mathrm{5}}$ s${}^{-\mathrm{6}/\mathrm{5}}$	deep-water alongshore sediment transport coefficient
${\mathit{\rho }}_{\mathrm{s}}$	2650	kg m${}^{-\mathrm{3}}$	sediment density
$p$	0.4		dry mass void fraction
$H_{\mathrm{s}}$	0.8	m	significant wave height, from the Cap Tortosa wave buoy (Bolanos et al., 2009)
$T_{\mathrm{p}}$	4.1	s	peak wave period, from the Cap Tortosa wave buoy (Bolanos et al., 2009)
$D_{\mathrm{sf}}$	10	m	shoreface depth
Shoreface slope	0.01
Shelf slope	0.002
Critical barrier width	250	m	width below which storm overwash occurs
Cell width	50	m
Time step	1	day	coastal model time step
River Profile Model (Parker, 2004)
$D$	10	mm	$D_{\mathrm{50}}$ of transported bed material (Vericat et al., 2006)
$B$	150	m	channel width
$S$	5.8 $\times$ 10${}^{-\mathrm{4}}$		modern profile gradient
$I$	0.3		flow intermittency
$k_c$	100	mm	bed roughness, 3$\times$ bed material $D_{\mathrm{50}}$ (Vericat et al., 2006)
Profile length	450	km	along channel distance of the upstream boundary condition
Time step	0.2	yr	fluvial model time step
${\mathit{\alpha }}_r$	8.1		dimensionally homogenous Manning–Strickler coefficient (Parker, 2004)
${\mathit{\alpha }}_{\mathrm{t}}$	8		sediment transport coefficient (Parker, 2004)
$n_{\mathrm{t}}$	1.5		sediment transport exponent
${\mathit{\tau }}_{\mathrm{c}}$*	0.045		critical Shields stress for sediment motion
$\mathit{\phi }$	1		skin friction fraction

Wave climate

Wave height and the directional distribution of incoming waves exert a first-order control on wave-influenced delta evolution (Ashton and Giosan, 2011). We compared five different wave climatology sources from nearby the Ebro Delta and investigated their effect on modeled alongshore sediment transport. Wave climates are extracted from two directional wave buoys and three hind-casted wave models (Fig. 3 and Table 2). All sources are located in sufficiently deep water for the waves to be treated as deep-water waves (depth > 1/4 $\mathit{\pi }{T}_{\mathrm{p}}^{\mathrm{2}})$, and all sources show peaks of wave intensity from the east and from the south that affect Ebro Delta alongshore transport. The different wave sources differ particularly in the relative strength of the waves approaching from the south. This could be because the southerly (summer) waves are generated more locally (Jiménez et al., 1997) and therefore their magnitude may be sensitive to buoy location or hindcast methodology.

(a) Comparison of the five different wave roses and their location on a map from NOAA (2015). See Table 2 for an overview of the sources. (b) Angular distribution of alongshore sediment transport potential for the five different sources.

[Figure omitted. See PDF]

Overview of five different sources of wave climate data close to the Ebro Delta. See Fig. 3 for an overview of locations and the angular distribution of alongshore sediment transport potential. Wave height is the effective, yearly averaged wave height weighted by its ability to move sediment alongshore, i.e. ($\mathrm{\Sigma }$ $H_{\mathrm{s}}^{\mathrm{2.4}}{)}^{\mathrm{1}/\mathrm{2.4}}$. The $R^{\mathrm{2}}$ value is the coefficient of determination of the alongshore sediment transport calculated from the wave data versus the measurements of Jiménez and Sánchez-Arcilla (1993).

Name	Type	Lat	Long	Water	Wave	Wave	$Q_{\mathrm{s},\max }$	$R^{\mathrm{2}}$	Data period	Reference
		${}^{\circ }$N	${}^{\circ }$E	depth (m)	height (m)	period (s)	(kg s${}^{-\mathrm{1}}$)		(yr)
Cap Tortosa	buoy	40.7	1.0	60	0.8	4.1	47.9	0.89	1990–2011	Bolanos et al. (2009)
Tarragona	buoy	41.0	1.2	24	1.0	5.5	72.4	0.86	2004–2011	Puertos del Estado (2015)
MedAtlas	model	40.0	1.0	222	0.7	4.0	48.3	0.76	1992–2002	Gaillard et al. (2004)
Hipocas	model	40.8	1.0	68	1.1	4.9	71.1	0.87	1958–2001	Sotillo et al. (2005)
Wavewatch III®	model	40.8	0.8	63	0.7	4.9	31.1	0.86	1979–2009	Chawla et al. (2013)

Testing the alongshore sediment transport model assumptions

Jiménez and Sánchez-Arcilla (1993) used aerial photographs from 1957 to 1989 and beach profile measurements between 1988 and 1992 to calculate Ebro coastline change. Their study found sustained multi-decadal rates of erosion of up to 50 m yr${}^{-\mathrm{1}}$ close to the river mouth, and progradation of about 10 m yr${}^{-\mathrm{1}}$ along the spits (Jiménez and Sánchez-Arcilla, 1993). These measured recent shoreline changes allow us to test the one-line shoreline assumptions underlying the delta evolution model. We use the back-refracted CERC formula (Ashton and Murray, 2006) to calculate alongshore sediment transport ($Q_{\mathrm{s}}$, kg s${}^{-\mathrm{1}})$ from deep-water wave characteristics, $$Q_{\mathrm{s}}=K_{\mathrm{1}}\cdot {\mathit{\rho }}_{\mathrm{s}}\cdot \left(\mathrm{1},-,p\right)\cdot H_{\mathrm{s}}^{\mathrm{12}/\mathrm{5}}{T}^{\mathrm{1}/\mathrm{5}}{\cos }^{\mathrm{6}/\mathrm{5}}\left({\mathit{\phi }}_{\mathrm{0}},-,{\mathit{\phi }}_{\mathrm{s}}\right)\sin \left({\mathit{\phi }}_{\mathrm{0}},-,{\mathit{\phi }}_{\mathrm{s}}\right),$$ where $H_{\mathrm{s}}$ is the offshore deep-water significant wave height (m), $T$ is the wave period (s), ${\mathit{\phi }}_{\mathrm{0}}$ is the deep-water wave approach angle (which equals $\mathit{\gamma }$-$\mathit{\theta }$ in a regional setting, Fig. 3a), and ${\mathit{\phi }}_{\mathrm{s}}$ is the local shoreline orientation (Ashton and Murray, 2006). We use a $K_{\mathrm{1}}$ of 0.035 m${}^{\mathrm{3}/\mathrm{5}}$ s${}^{-\mathrm{6}/\mathrm{5}}$ compared to the typical coefficient of 0.06 m${}^{\mathrm{3}/\mathrm{5}}$ s${}^{-\mathrm{6}/\mathrm{5}}$ (Komar, 1998; Nienhuis et al., 2015) based on Ebro Delta calibration studies of Jiménez and Sánchez-Arcilla (1993).

For five different wave sources (Table 2 and Fig. 3), we computed net alongshore sediment transport along the modern Ebro Delta shoreline, extracted from the NOAA shoreline database (NOAA, 2015). We correct for shadowing of certain wave approach angles by other portions of the delta coastline.

The calculated littoral sediment transport trends along the Ebro Delta coastline are similar between the five wave climates (Fig. 4), showing sediment transport is greatest along both spits and close to the modern river mouth. The computed sediment transport magnitude however between the wave climate sources differs by almost a factor of 3. All wave climates except for the MedAtlas have similar correlation coefficients when compared to sediment transport patterns estimated based on observed beach change (Fig. 4b, black markers) (Jiménez and Sánchez-Arcilla, 1993). We choose to use the Cap Tortosa buoy data (described in Bolanos et al., 2009) in the delta evolution model because its 21-year record is sufficiently long, it is located close to the mouth of the modern Ebro River, and its wave height and wave period are bound by the other four wave sources.

(a) The Ebro Delta coastline, colored by the simulated alongshore sediment transport flux from the Cap Tortosa data. (b) Alongshore sediment transport along the Ebro Delta coastline from all five wave climate sources (and assuming no sediment was supplied by the Ebro River). Alongshore transport is positive to the right when looking offshore. Black markers indicate alongshore sediment transport estimates from Jiménez and Sánchez-Arcilla (1993). (c) The Cap Tortosa buoy data recast into shoreline change rates using the one-contour-line approach (Eq. 2) compared to the measured shoreline change rates from Jiménez and Sánchez-Arcilla (1993).

[Figure omitted. See PDF]

From the computed alongshore sediment transport gradients from the Cap Tortosa data, we predict shoreline accretion and erosion using the one-contour-line approach (Eq. 1). For this comparison we use a decadal timescale shoreface depth of 7 m (Jiménez and Sánchez-Arcilla, 1993). In general, the rate of shoreline change is well predicted ($R^{\mathrm{2}}=$ 0.84) by the one-contour-line model and the wave climate from the Cap Tortosa buoy (Fig. 4c).

Aside from testing our model, we can draw two observations from the measurements of Jiménez and Sánchez-Arcilla (1993) about the ongoing coastal changes of the Ebro Delta. First, around the river mouth there is rapid coastal retreat to the south, and deposition further to the north. The field measurements align with the one-contour-line simulations close to the river mouth – these simulation results do not include a fluvial sediment contribution, and are therefore consistent with other studies suggesting negligible modern fluvial sediment supply to the coast (Jiménez and Sánchez-Arcilla, 1993).

Secondly, the sediment transport patterns along the spits can be cast in the framework proposed by Ashton et al. (2016). Along the barrier sections of the Ebro Delta spits, the computed alongshore sediment transport gradients are nearly zero, whereas measured shoreline retreat is approximately 10 m yr${}^{-\mathrm{1}}$ (Fig. 4c). This suggests that in these regions overwash is driving coastline retreat without gradients in alongshore sediment transport. The barrier section (the “neck”) is fed by a sediment source upcoast and is generally erosional up to a fulcrum point, where alongshore sediment transport is maximized and erosion transitions into deposition (Ashton et al., 2016). The measured and simulated shoreline change indicate that the northern and the southern spit are indeed depositional and are prograding at approximately 10 m yr${}^{-\mathrm{1}}$ (Fig. 4c).

River Profile Model

We investigate the response timescales of the river basin to climate and land-use changes using an exploratory 1-D river profile model (Parker, 2004). In this model, sediment is not merely a passive tracer, but interacts with the bed elevation to reach a longitudinal profile in morphodynamic equilibrium (Carling and Cao, 2002). The interaction between flow and topography creates a dynamic model – rivers are not treated as static pipes – which allows us to use the computed longitudinal profiles together with the observed modern longitudinal profile to investigate potential past and present sediment transport conditions. Additionally, by focusing on the interaction of the flow with the channel bed, we can model the bed material load – the sediment that makes up most of the delta shoreface (Maldonado, 1975) – while we ignore the finer grained washload that is mostly deposited farther offshore. In the absence of subsidence or sea-level changes, and if the bed material load and the flow discharge are in equilibrium, the bed slope does not change and the capacity is equal to the supply.

The channel bed in the model is freely erodible and our approach is therefore strictly applicable to alluvial, transport-limited systems (Parker, 2004). A similar 1-D river profile model was recently applied to study timescales of sediment supply decreases in the Mississippi River (Nittrouer and Viparelli, 2014). Their study suggested a long ($O$ 100 years) delay between dam construction $\sim$ 1000 km upstream and sand load changes near the coast.

The 1-D river profile model assumes normal flow conditions, such that a width-averaged momentum balance connects bed slope and flow depth to bed shear stress. Flow depth in the channel is determined using a Manning–Strickler formulation for the flow resistance (Parker, 2004). Because of the gravel bed of the Ebro River (Vericat and Batalla, 2006) we use the Meyer-Peter and Muller (1948) equation to calculate fluvial bed-material load (kg s${}^{-\mathrm{1}})$, $$Q_{\mathrm{r}}=I{\mathit{\rho }}_{\mathrm{s}}B\sqrt{\mathrm{RgD}}D{\mathit{\alpha }}_{\mathrm{t}}{\left\{{\left({\displaystyle \frac{Q^{\mathrm{2}}{k}_c^{\mathrm{1}/\mathrm{3}}}{{\mathit{\alpha }}_r^{\mathrm{2}}{\mathrm{gB}}^{\mathrm{2}}}}\right)}^{\mathrm{3}/\mathrm{10}},{\displaystyle \frac{S^{\mathrm{7}/\mathrm{10}}}{RD}},-,{\mathit{\tau }}_{\mathrm{c}}^{*}\right\}}^{\mathrm{2.5}},$$ where $R$ is the submerged specific gravity of the sediment (1.65); $g$ is gravity (m s${}^{-\mathrm{2}})$; $D$ is the median grain size (m); ${\mathit{\alpha }}_{\mathrm{t}}$ and ${\mathit{\alpha }}_r$ are flow and sediment transport coefficients; $Q$ is a representative flood discharge (m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}})$; $k_c$ is the bed roughness (m); $S$ is the channel bed slope; $I$ is the flood discharge intermittency; ${\mathit{\rho }}_{\mathrm{s}}$ is the sediment density (kg m${}^{-\mathrm{3}})$; $B$ is the channel width (m); and ${\mathit{\tau }}_{\mathrm{c}}^{*}$ is a critical Shields stress for sediment motion (Parker, 2004). See Table 1 for an overview of the model parameters.

From Eq. (3) we can observe that the intermittency $I$, the flood discharge $Q$, and the grain size $D$ are sensitive parameters for the fluvial sediment load estimates. The flood intermittency factor $I$ characterizes the fraction of time (generally a year) the river is in flood. Frequently, this factor is scaled with a particular flood discharge to match an observed annual fluvial sediment flux $Q_{\mathrm{r}}$ (Wright and Parker, 2005). However, the pre-dam fluvial sediment flux of the Ebro River is poorly constrained, so here instead we estimate the flood intermittency $I$ directly from flow records from Batalla et al. (2004) to generate an independent estimate of the fluvial sediment flux. To estimate flood intermittency, we first fit a function to the flow-exceedance frequency statistics of Batalla et al. (2004), $$Q\left(e\right)=\mathrm{550}\cdot e^{-\mathrm{0.25}},$$ where $Q$ is Ebro River discharge (m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}})$ and $e$ is the pre-dam exceedance frequency (i.e., $e=$ 0.1 indicates a discharge that is exceeded 36.5 days each year). From Eq. (4), a representative flood intermittency for an annual bed-material load $Q_{\mathrm{r}}$ can be estimated by taking into account all the floods from an extreme flood ($e=$ 0) to a critical exceedance frequency for bed-material-load motion (700 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, or $e_{\mathrm{crit}}\approx$ 0.25) (Vericat and Batalla, 2006). Because sediment transport is nonlinearly related to flow we do not integrate Eq. (4) directly, but rather we scale discharge to sediment flux with an exponent of 1.5 ($Q_{\mathrm{r}}$ $\sim$ $Q^{\mathrm{1.5}}$, Eq. 3). Formalized, the flood intermittency of a particular flood magnitude can be described as $$I\left(e\right)={\displaystyle \frac{\mathrm{1}}{Q{\left(e\right)}^{\mathrm{1.5}}}}\underset{\mathrm{0}}{\overset{e_{\mathrm{crit}}}{\int }}Q(\mathit{\epsilon }{)}^{\mathrm{1.5}}\mathrm{d}\mathit{\epsilon }.$$ In the river profile model we choose a flood discharge of 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ which occurred relatively often with a pre-dam exceedance frequency $e$ of 0.15 (Vericat and Batalla, 2006). For a 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ flood, Eq. (5) evaluates to an intermittency $I$ of approximately 0.3. In other words, the instantaneous bed-material load of a 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ flow roughly corresponds to an annual fluvial bed-material load if we use a 30 % intermittency factor, which we therefore use in the model.

The third sensitive parameter affecting the fluvial profile model is the grain size. The Ebro River is a gravel-bed river (most mobile $D_{\mathrm{50}}$ is $\sim$ 10 mm) (Vericat and Batalla, 2006), so aggradation and erosion due to divergences in the bed-material load should be modeled using gravel size sediments. However, the median sediment size of the Ebro shoreface is sand ($\sim$ 0.2 mm) (Jiménez and Sánchez-Arcilla, 1993). In the coupled fluvial-delta system we should therefore consider the sand load as the representative bed material load volume at the river mouth. To retain the simplicity of a unimodal fluvial profile model we choose a 10 mm median bed-material load grain size to compute the timescales of profile incision and aggradation. Given the relatively constant slope of the Ebro River ($S=$ 5.8 $\times$ 10${}^{-\mathrm{4}}$, Fig. 5), we assume that the bed material load at the Ebro Delta should be roughly similar to the bed-material load further upstream despite the change in the median grain size.

(a) The Ebro River basin showing the main river channel in light blue and larger tributaries in darker blue, colored by elevation. (b) The elevation profile of the Ebro River, with the equilibrium profile model prediction in red dashed line. The black dashed line shows the cumulative fraction of the Ebro pre-dam discharge from Batalla (2004).

[Figure omitted. See PDF]

Applying the model based on the pre-dam fluvial and discharge conditions, the median bed-material load grain size, and the observed slope ($D_{\mathrm{50}}=$ 10 mm, $Q_{\mathrm{flood}}=$ 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, $I=$ 30 %, $S=$ 5.8 $\times$ 10${}^{-\mathrm{4}})$ we find an annual average bed-material load transport rate $Q_{\mathrm{r}}$ of 70 kg s${}^{-\mathrm{1}}$ (2.2 Mt yr${}^{-\mathrm{1}})$. This estimate however is sensitive to the bed roughness ($k_c)$ which we estimate at 100 mm, $\sim$ 3 times the bed material $D_{\mathrm{50}}$ (Vericat et al., 2006).

We model the Ebro drainage basin as a single channel representing an average of its tributaries. This 1-D river profile model requires the choice of an upstream boundary, representing the average location of the fluvial discharge and sediment supply into the drainage basin. The choice of an upstream boundary is important because it acts as a first-order control on fluvial sediment transport timescales from the drainage basin to the delta. To find an appropriate upstream boundary, we calculated the pre-dam morphologic (2-year) flood discharge along the Ebro River relative to the discharge at the delta from existing hydrologic records (Batalla et al., 2004). We set the upstream boundary condition at 450 km upstream of the delta, where the pre-dam morphologic (2-year) flood discharge is 50 % of its final discharge at the delta and a clear discontinuity in the longitudinal profile occurs (Fig. 5). Note that the observed channel slope remains constant upstream of the confluence with the Cinca River even though the flood discharge decreases significantly, a sign of fluvial or sedimentological heterogeneity within the drainage basin. However, a spatially explicit model of the Ebro basin taking into account these heterogeneities would be a significant departure from our exploratory model approach.

The apex of the delta should be considered the downstream boundary of the fluvial profile model, because the normal flow assumption is invalid in the backwater zone near the river mouth, where the channel aggrades and progrades and the flow is nonuniform (Hotchkiss and Parker, 1991). However, as Chatanantavet et al. (2012) recently demonstrated, annual flooding cycles in the backwater zone often create a condition where aggradation during low flow is nearly balanced by erosion during high flow. This (near) balance suggests that in terms of bedload volumes and neglecting subsidence and sea-level rise, delta progradation is significantly larger than channel aggradation and, therefore, that the absence of a backwater zone in our normal flow model only results in a limited overestimation of the fluvial sediment supply to the river mouth when considering centennial timescales. In our simplified river profile model, we therefore assume that all bedload sediment transported to the apex of the Ebro Delta is deposited near the river mouth as delta foreset. Additionally, this prevents the need to couple the coastal and fluvial models directly. Rather, we treat both models as exploratory and inform timescales and sediment fluxes of coastal and fluvial change based on outcomes from each model.

Testing the fluvial profile model

To test the applicability of the river profile model to the Ebro drainage basin, we compare model estimates to recent measured bed elevation and sediment transport changes 25 km downstream of the lowermost Flix Dam for 55 years after its construction in 1948 (Fig. 6) (Vericat and Batalla, 2006). Between 2002 and 2004, Vericat and Batalla (2006) observed an average bedload transport rate of 12 kg s${}^{-\mathrm{1}}$ (0.4 Mt yr${}^{-\mathrm{1}}$), down from pre-dam estimates of around 70 kg s${}^{-\mathrm{1}}$ (2.2 Mt yr${}^{-\mathrm{1}}$). They also observed downstream scour at a rate of about 0.03 m yr${}^{-\mathrm{1}}$ in Mora d'Ebre, although with much variability. To model river profile response to dam construction, we applied a 100 % reduction in sediment supply immediately downstream of the Flix Dam. Concomitantly, following analysis of Vericat and Batalla (2006), we impose a fourfold decrease in the occurrence of bedload transporting floods of 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ (from a 15 to a 4 % exceedance probability, or a 30 to an 8 % intermittency factor).

(a) A close-up of the Ebro drainage basin close to the delta (image from Landsat). (b) Modeled response of the Ebro River downstream of the lowermost modern dam, the Flix Dam.

[Figure omitted. See PDF]

Even though the model does not capture processes such as bed armoring and downstream fining, results show reasonable agreement with the field measurements, estimating about 1.5 m of bed degradation at Mora d'Ebre 55 years after dam construction, an incision rate of approximately 0.03 m yr${}^{-\mathrm{1}}$, and a local sediment bed-material load of 12 kg s${}^{-\mathrm{1}}$ (0.4 Mt yr${}^{-\mathrm{1}})$. At Mora d'Ebre, the measurement location of Vericat and Batalla (2006), Eq. (3) predicts that the change in flooding frequency decreased the coarse-grained sediment flux from 70 kg s${}^{-\mathrm{1}}$ (2.2 Mt yr${}^{-\mathrm{1}})$ to 17 kg s${}^{-\mathrm{1}}$ (0.5 Mt yr${}^{-\mathrm{1}})$. The sediment capture in the reservoirs and the subsequent channel bed slope adjustment decreased the coarse-grained sediment flux further from 17 to 12 kg s${}^{-\mathrm{1}}$. Furthermore, model results suggest that the bed response to dam construction has not yet reached the Ebro Delta. At the delta, the model predicts a bed-material load of 16.5 kg s${}^{-\mathrm{1}}$ 55 years after dam construction, compared to a predicted 17 kg s${}^{-\mathrm{1}}$ immediately after dam construction due to the change in flooding frequency. Of the total reduction in bed-material load to the delta, the model therefore predicts that about 99 % is due to changes in the flooding frequency, whereas only 1 % is due to a capturing of the sediment in the reservoirs and a resulting change in the channel bed slope. The model prediction for the Ebro Delta is higher than the estimate of 1.6 kg s${}^{-\mathrm{1}}$ (0.05 Mt yr${}^{-\mathrm{1}})$ of Jiménez et al. (1990) based on the formulae from van Rijn (1984)

Results

Delta response to increased fluvial sediment supply

We investigated whether changes in fluvial sediment supply could explain the rapid growth of the Riet Vell lobe, which potentially occurred sometime between 3000 and 1100 years BP (Canicio and Ibáñez, 1999). Cast in terms of the fluvial dominance ratio $R$ (Nienhuis et al., 2015), the transition from a slowly growing cuspate delta to a rapidly growing pointy (not cuspate) delta should be expected to occur when $R$ > 1 (or $Q_{\mathrm{r}}$ $\sim$ 50 kg s${}^{-\mathrm{1}}$ for the Cap Tortosa wave data; see Table 2). At a pre-dam estimate of 70 kg s${}^{-\mathrm{1}}$ (2.2 Mt yr${}^{-\mathrm{1}}$) (Syvitski and Saito, 2007), this means that during the period of rapid growth, a single thread channel of the Ebro should have been river-dominated or close to a transition to river dominance, with a fluvial dominance ratio $R$ of 1.4.

We also investigated the effect of fluvial sediment supply on plan-view Ebro Delta morphology with the CEM. After 750 model years, for bedload sediment fluxes up to about 35 kg s${}^{-\mathrm{1}}$ (1 Mt yr${}^{-\mathrm{1}})$, the modeled delta plain exhibits a smooth cuspate morphology (Fig. 7a) while prograding at about 6 m yr${}^{-\mathrm{1}}$ (5 km in 800 years, Fig. 7c). A delta supplied double this sand load (70 kg s${}^{-\mathrm{1}}$; 2 Mt yr${}^{-\mathrm{1}})$, however, progrades 5 times more rapidly ($\sim$ 30 m yr${}^{-\mathrm{1}})$, developing shoreline instabilities along the updrift and downdrift flanks.

A modeled delta lobe (a) after 750 years of growth and (b) after 150 years of reworking (950 years of total model time). (c) Contour diagram of the progradation distance versus time as a function of the fluvial sediment flux $Q_{\mathrm{r}}$, or the river dominance ratio $R$.

[Figure omitted. See PDF]

From the same set of model experiments, we can also study the effect of fluvial sediment supply on post-avulsion abandonment and wave reworking. For low pre-abandonment fluvial sediment supply (< 40 kg s${}^{-\mathrm{1}}$; 1.2 Mt yr-1), because the delta remains wave-dominated during growth ($R$ < 1) and with cuspate and continuous pre-abandonment morphology, no spit forms after abandonment (Fig. 7b) (Nienhuis et al., 2013). For high fluvial sediment supply during growth ($Q_{\mathrm{r}}$ > 50 kg s${}^{-\mathrm{1}}$, $R$ > 1), because the delta coast grows with a pointy shape, a spit forms after abandonment (Fig. 7b).

We therefore estimate that the early cuspate morphology (around 3000 years BP; Canicio and Ibáñez, 1999; Cearreta et al., 2016) was formed with a fluvial sediment supply of at most 35 kg s${}^{-\mathrm{1}}$. The latter, more rapidly growing Riet Vell lobe that was reworked into a spit, was formed with a significantly larger fluvial sediment supply, such that $R$ > 1 (likely more than 50 kg s${}^{-\mathrm{1}})$. Extending the progradation trajectory of the Riet Vell lobe (Fig. 7c) and keeping in mind that the modern bathymetry suggests a maximum Riet Vell lobe extent of $\sim$ 20 km (Canicio and Ibáñez, 1999), we estimate a Riet Vell fluvial sediment supply of $\sim$ 70 kg s${}^{-\mathrm{1}}$. Note that these flux estimates are sensitive to model parameters such as the effective shoreface depth, the littoral CERC formula constant, and the wave height, which were estimated based on modern Ebro Delta change as described in Sect. 3.1.

Timescales of change on the delta plain

Aside from linking fluvial sediment fluxes to Ebro Delta morphology, we can also use the Coastline Evolution Model to assess the timescales of Ebro Delta morphologic change. To investigate the timescales, we have simulated the growth and reworking of all three lobes. In 42 different simulations we use the estimated fluvial sediment supply of 70 kg s${}^{-\mathrm{1}}$ and we vary the growth times of the different lobes. For example, in one simulation we grew the Riet Vell lobe for 800 years, then the Sol de Riu lobe for 400 years, and finally the Mitjorn lobe for 300 years (Fig. 8). In another simulation, we used growth times of respectively 1200, 600, and 300 years.

To assess which one of all the 42 simulations best represents the actual history of the Ebro Delta, we measured the radial lengths of the modeled lobes through time. Then, we measured the radial lengths of the lobes on the modern Ebro Delta from the avulsion apex. Both the paleo-channels of the Riet Vell and the Sol de Riu lobe currently extend approximately 10 km from the avulsion apex. The modern active lobe, the Mitjorn, extends about 15 km from the avulsion apex (Fig. 2). We select the best-matched model simulation as the one where the three lobes reach the currently observed lengths of the modern Ebro Delta at the same time. This “reverse engineering” approach yields an estimate of how long each lobe was active and therefore also of the start of Ebro Delta plain's rapid growth. These estimates are made independently of published field studies, using the modern delta plain morphology.

In the case of the one simulation where the Riet Vell lobe grows for 800 years, the Sol de Riu for 400 years, and the Mitjorn for 300 years (Fig. 8), we find that for these growth times the radial extents of both the Riet Vell and the Sol de Riu are less than 10 km when the Mitjorn reaches 15 km (the current observed channels' lengths) because both the simulated Riet Vell and the Sol de Riu shores erode past the modern shore following avulsion.

The best-matched model scenario of the consecutive growth of the three delta lobes has growth times of 1200, 600, and 300 years, respectively (solid lines in Fig. 8), when the modeled Riet Vell and Sol de Riu have eroded back to the modern observed lengths of 10 km, and the Mitjorn has prograded to 15 km. We estimate therefore, based on this best-matched model scenario, that the period of rapid growth of the Ebro Delta plain lasted approximately 1200 $+$ 600 $+$ 300 $=$ 2100 years, placing the time at which rapid delta growth started approximately 2100 years BP (Fig. 8). These growth times would suggest that the second avulsion occurred 300 yr BP, and the first avulsion occurred 900 yr BP. Note, however, that these avulsion times estimates are sensitive to the fluvial sediment supply to the delta (here kept at 70 kg s${}^{-\mathrm{1}})$ and its variability through time. We do not model time-varying fluvial sediment supply in any of our 42 varying avulsion time simulations to limit the number of model variables.

Simulated radial extent of the three different Ebro Delta lobes for a sediment supply of 70 kg s${}^{-\mathrm{1}}$ and forced avulsions after 1200 and 1800 model years (solid lines) and after 800 and 1200 model years (dashed lines). Note that the radial extent can increase without the lobe being active because of littoral sediment transported from adjacent lobes. Three inset deltas show the solid line model run after 700, 1350, and 1900 years. The gray second horizontal axis indicates the real time inferred from the solid line model run and the modern Ebro Delta morphology, where, at the year 2015, lobes 1 and 2 are approximately 10 km long, and the active lobe is 15 km long, measured from the apex (Fig. 2).

[Figure omitted. See PDF]

The best-matched model estimates for the start of rapid delta growth, made purely based on physical constraints set by alongshore sediment transport and fluvial sediment supply, roughly coincides with a simple volumetric estimate based on our assumed shoreface depth ($\sim$ 10 m) and coarse-grained sediment supply ($\sim$ 70 kg s${}^{-\mathrm{1}}$, 2.2 Mt yr${}^{-\mathrm{1}})$, and the modern delta plain area beyond the dated beach ridges of Canicio and Ibáñez (1999) ($\sim$ 280 km${}^{\mathrm{2}})$, $$T_{\mathrm{s}}={\displaystyle \frac{A\cdot D_{\mathrm{sf}}}{{\mathit{\rho }}_{\mathrm{s}}\cdot (\mathrm{1}-p)\cdot Q_{\mathrm{r}}}}\approx \mathrm{2030}\mathrm{years},$$ where $T_{\mathrm{s}}$ is the time since the start of rapid growth and $A$ is the delta plain area (m${}^{\mathrm{2}})$. Our best-matched model also agrees with observations suggesting increased flood plain deposition in the drainage basin ($\sim$ 2000–1800 yr BP; Thorndycraft and Benito, 2006).

The model-estimated avulsion times compare closely with the existing, albeit limited, historical evidence (Canicio and Ibáñez, 1999; Somoza and Rodriguez-Santalla, 2014), at least for the avulsion of the Sol de Riu at $\sim$ 300 years BP. Other model simulation results, such as the maximum extent of the modeled Riet Vell Lobe ($\sim$ 20 km, Fig. 8), approximate earlier estimations of its extent made from bathymetry (Canicio and Ibáñez, 1999). Even though the history of the Ebro Delta was likely more complex than our model simulations, the qualitative agreement between the model scenario and the growth, reworking, and spit formation observed on the Ebro Delta suggests the possibility that the gross morphology of the delta plain can develop without significant sea-level or fluvial sediment supply fluctuations.

Model simulations show the development of spits during both lobe growth and lobe abandonment. However, spits growing during growth and reworking have markedly different orientations (Fig. 8). Ashton et al. (2016) suggest that spit orientation is strongly affected by the updrift shoreline change rate. We speculate that, based upon their more river parallel orientations, the lagoons in the southern region of the modern Ebro Delta plain (e.g. the Encanyissada, Clot, and Tancada lagoons, Fig. 2) formed as they were enclosed by spits created while the delta was growing. On the other hand, the active southern La Banya spit has a different orientation because it was formed as the updrift shoreline retreated during reworking of the Riet Vell lobe.

Wave climate change as a potential cause of delta growth

Investigating the effect of changes in sediment supply on the Ebro Delta, we assumed the wave climate was constant. However, previous studies (Goy et al., 2003; Sabatier et al., 2012) focusing on the western Mediterranean over the last millennia suggest evidence exists of significant changes in wave climate. Goy et al. (2003), studying the cuspate coast of the Gulf of Almeria in southern Spain, correlated beach ridge progradation to periods of negative North Atlantic Oscillation (NAO) because of stronger winds from the southwest that would increase littoral drift to the coast.

To investigate the potential effect of a change in the NAO index on the fluvial dominance ratio $R$, we correlated the monthly NAO index (Jones et al., 1997) with the Hipocas record (Sotillo et al., 2005), the longest wave climate hind-cast record available, spanning 44 years (Table 2). Over this 44-year timespan, there were higher waves from the south during periods of negative NAO (Fig. 9a). For more positive NAO values, average wave height is lower, particularly from the south. Calculating the monthly $Q_{\mathrm{s},\max }$, and comparing it to the NAO index, we find a weak trend from 60 kg s${}^{-\mathrm{1}}$ for strongly negative NAO ($-$4) to 35 kg s${}^{-\mathrm{1}}$ for periods of strongly positive NAO ($+$4) (Fig. 9b).

(a) The angular distribution of alongshore sediment transport potential from Hipocas hind-cast data (Sotillo et al., 2005), separated into periods of negative and positive monthly North Atlantic Oscillation (NAO) index (Jones et al., 1997). Insets show wave roses weighted by alongshore sediment transport potential for positive and negative NAO. (b) The monthly NAO versus the monthly maximum potential alongshore sediment transport $Q_{\mathrm{s},\max }$. Dotted line shows the least-squares best fit line. Inset shows the NAO index distribution for 3000–2000, 2000–600, and 600–0 years BP (Jones et al., 1997; Olsen et al., 2012). (c) Computed distribution of $Q_{\mathrm{s},\max }$ for different time periods.

[Figure omitted. See PDF]

Climate reconstructions suggest that the NAO index since the mid-Holocene can be divided into three distinct periods. Prior to 2000 years BP the NAO index was mostly negative; afterwards, up to about 600 years BP, it changed to become mostly positive. Over the past 600 years, the NAO index has been fluctuating with short but strongly negative periods (Jones et al., 1997; Olsen et al., 2012).

To obtain approximations of $Q_{\mathrm{s},\max }$ of the last 3000 years, we determined distributions of NAO indices from Olsen et al. (2012) and Jones et al. (1997) for each of the three periods (Fig. 9b). We find that extreme NAO indices are rare and that the distributions of NOA indices, even though distinct, also overlap considerably. Therefore, although $Q_{\mathrm{s},\max }$ can vary with changes in the NAO, particularly on a year-to-year basis (Fig. 9b), proxy-record constructions based on the NAO do not suggest significant sustained differences across the previous two millennia (Fig. 9c), also in comparison to suggested increases in the fluvial sediment supply to the Ebro Delta. These suggested fluvial sediment supply increases range from 40 to 350 % and have been sustained up to the 20th century (Guillén and Palanques, 1997a; Xing et al., 2014). Large-scale wave climate changes may have occurred over the past 2000 years; such changes, however, do not jump out of our analysis of NAO cycles, suggesting that one does not necessarily have to appeal to wave climate changes to explain the evolution of the Ebro Delta.

Timescales of environmental change in the fluvial catchment

Results from the CEM in concert with previous records indicative of hydrologic change (Thorndycraft and Benito, 2006), place the start of Ebro Delta plain's rapid growth at approximately 2100 years BP. Additionally, CEM model experiments indicate that roughly a sustained doubling in the sediment flux (from 35 kg s${}^{-\mathrm{1}}$ (1 Mt yr${}^{-\mathrm{1}})$ to 70 kg s${}^{-\mathrm{1}}$ (2 Mt yr${}^{-\mathrm{1}}\left)\right)$ over this period of time could create the observed morphologic changes in growing delta morphology. We have run four different scenarios in the river profile model to estimate the types and timing of drainage basin changes that could explain this increased fluvial sediment supply to the delta from 35 to 70 kg s${}^{-\mathrm{1}}$ starting 2100 years BP and lasting up to the 20th century. The four scenarios are (1) an increase in fluvial sediment supply, (2) an increase in fluvial flood discharge, (3) an increase in fluvial flood discharge and fluvial sediment supply, and (4) an increase in fluvial flood discharge and a 500-year lag in an increase in fluvial sediment supply (Table 3).

Overview of the four river profile model experiments and their final equilibrium slope and bed level change. $Q$ is the fluvial flood discharge, $Q_{\mathrm{r}}$ is the upstream fluvial sediment supply, $i$ is the initial antecedent fluvial environment, and $f$ is the final fluvial environment.

Description	$Q_i$	$Q_f$	$Q_{\mathrm{r},i}$	$Q_{\mathrm{r},f}$	Slope ($i$)	Slope ($f$)	Upstream bed
	(m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$)	(m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$)	(kg s${}^{-\mathrm{1}}$)	(kg s${}^{-\mathrm{1}}$)	($\times$ 10${}^{-\mathrm{4}}$)	($\times$ 10${}^{-\mathrm{4}}$)	level change (m)
Sediment $\times$2	900	900	35	70	4.0	5.8	80
Discharge $\times$1.5	600	900	35	35	5.8	4.0	$-$80
Discharge $\times$1.5 and sediment $\times$2	600	900	35	70	5.8	5.8	0
Discharge $\times$1.5 and sediment $\times$2 w/ delay	600	900	35	70	5.8	5.8	0

In scenario one we change the fluvial sediment supply 450 km upstream from the Ebro Delta from 35 to 70 kg s${}^{-\mathrm{1}}$, with the flood discharge and its intermittency remaining constant. Such a scenario could arise from land clearing that increased sediment supply without altering the discharge. The model experiment shows that the channel bed slowly aggrades to the new sediment supply and that the change in supply signal takes approximately 5000 years to significantly affect the Ebro Delta (Fig. 10). This increase is associated with upstream aggradation of about 80 m. While there are numerous field studies that show large alluvial deposits throughout the Ebro drainage basin that date between 6000 and 2000 years BP (e.g. Benito et al., 2008; Constante-Orrios et al., 2009; Constante et al., 2010; Constante and Peña-Monné, 2009; González-Sampériz and Sopena Vicién, 2002; Gutiérrez-Elorza and Peña-Monné, 1998; Soriano, 1989), the majority of these deposits are on the order of $\sim$ 10 m thick. The unrealistic magnitude of the predicted aggradation is in part caused by the assumption that floodplain width remains constant, although the likely formation of a wider floodplain would not greatly affect the sediment supply to the delta. More importantly, the lack of any observed 80 m thick Holocene deposit makes it unlikely that exclusively a fluvial bedload sediment supply increase occurred in the Ebro drainage basin. Even though subsequent erosion of some deposits is likely, a sustained increase in sediment supply should have been accompanied by a sustained high slope and preserved upstream alluviation (Fig. 10b).

(a) Fluvial sediment flux at the apex of the delta and (b–e) longitudinal river profile evolution from four experiments of the river profile model with an increase in sediment supply (red), flood discharge (blue), sediment supply and flood discharge (orange), and flood discharge with a lagged sediment supply (green). (f–i) Time evolution of the channel bed and water surface elevation through time, 200 km upstream of the delta at the approximate location of the floodplain records from Benito et al. (2008). Note the different scales between (f–g) and (h–i). Expected occurrence of floodplain deposits (period of increasing water surface elevation) shown by the gray bars.

[Figure omitted. See PDF]

In contrast to an increase in fluvial sediment supply, any change in hydrology (flood magnitude and/or flood duration) affects sediment supply to the delta instantaneously. A 50 % increase in the flood magnitude results in a doubling of the fluvial sediment flux delivered to the delta, but would simultaneously cause the channel to start incising upstream (Fig. 10a). Over time, this discharge-driven incision gradually lowers the fluvial sediment flux at the river mouth, returning to the previous value after approximately 5000 years (Fig. 10c). A concave-down river profile would be diagnostic of an ongoing upstream adjustment to a large increase in discharge over the past several thousand years. However, as a concave-down river profile is not observed (Fig. 5b), we find it unlikely that an increase in flood discharge and/or duration is the sole cause of increased Ebro Delta growth.

In a third scenario, we investigated a simultaneous doubling of upstream sediment supply and a 50 % increase in the flood discharge. A combined change in sediment supply and discharge instantly doubles the sediment supply at the delta (Fig. 10a). Over time, incision due to discharge increases is compensated for by the aggradation caused by increased fluvial sediment supply (Fig. 10d).

Lastly, the fourth scenario we tested is also a doubling of the upstream sediment supply and a 50 % increase in the flood discharge, but now including a 500-year lag on the sediment flux. Such a scenario could be result of deforestation, where an instantaneous hydrologic signal is followed by a delayed secondary channel slope signal reaching the main stem of the Ebro River. We find that this fourth scenario has a nearly similar effect on deltaic sediment supply as the simultaneous discharge and sediment supply change (scenario 3). A delay in the increase in fluvial sediment flux has a small and temporary but measurable ($\sim$ 5 m) effect on the fluvial longitudinal profile (Fig. 10e).

Because floodplain aggradation is dependent on the elevation of the channel and water surface with respect to the surrounding floodplain (Heller and Paola, 1996; Schumm and Lichty, 1963), each of the tested scenarios would leave a distinct record in the floodplain deposits. Looking at modeled vertical profile changes 200 km upstream of the Ebro Delta, approximately the location of some of the floodplain records from Benito et al. (2008), the fourth scenario of increased floods (leading to channel incision) and a delayed increase in sediment flux (leading to channel aggradation) shows a double peaked response in floodplain aggradation. Our fluvial profile model suggests that an increase in flood discharge would reflect an initial period of floodplain aggradation that would decrease gradually as the channel starts to incise (Fig. 10i). The second period of floodplain aggradation would result from the increase in fluvial sediment supply. Radiocarbon dating of floodplain aggradation across the entire Iberian Peninsula similarly shows two periods of increased aggradation in the last 2000 years, one between 2000 and 1830 years BP, and one between 910 and 500 years BP (Benito et al., 2008).

In general, the river profile model experiments suggest an increase in either sediment or discharge alone is not responsible for the rapid and sustained growth of the Ebro Delta plain. Instead, a combination of increased flood discharge and increased fluvial sediment supply generates a response that best agrees with our understanding and previous findings of changes on the Ebro Delta plain. The observed channel bed slope appears to be in a long-term equilibrium, with no evidence of thick Holocene floodplain deposits. These model results here show that changing flooding and sediment discharge can cancel each other out, resulting in a sustained signal that can be felt instantaneously at the river delta. Both climate change and human impacts on landscapes such as deforestation can increase both the fluvial flood discharge and the fluvial upstream sediment flux (Syvitski and Milliman, 2007; Xing et al., 2014). Our fluvial profile model is therefore not able to quantify the individual response of either climate or land-use changes. However, the application of this fluvial profile model does illustrate that care should be taken when assuming that any change in the basin can result in an instantaneous and sustained change in sediment delivery to the delta (see also Nittrouer and Viparelli, 2014).

Discussion and conclusions

In this study we used two reduced-complexity models to temporally and physically constrain the late Holocene evolution of the Ebro Delta plain. Where possible, we assumed the simplest possible scenario of environmental change, focusing on the first-order effects on the river and its delta. The Coastline Evolution Model is able to broadly reproduce the size and shape of the Ebro Delta plain using only simplified fluvial and wave climate histories. However, both the delta model and the fluvial profile model are sensitive to a number of less well-constrained parameters, such as the shoreface depth and the fluvial grain size. Therefore, the general agreement of our model outcomes with earlier studies of the Ebro Delta change (e.g., Canicio and Ibáñez, 1999) should not be implied to indicate the absence of complicating factors; rather, it suggests that one does not necessarily have to appeal to complicating factors to explain the large-scale morphology of the Ebro Delta.

Using best-estimate model parameters, we find that an increase in coarse fluvial sediment supply to the delta approximately 2100 years BP is the most likely driver of growth of the modern Ebro Delta plain, whereby the delta prograded approximately 2–3 times faster than before (Cearreta et al., 2016). Additionally, model experiments with the delta evolution model show that Ebro Delta avulsions, where reworking of the abandoned lobes resulted in development of the modern La Banya and El Fangar spits, likely occurred around 900 years BP and 300 years BP, respectively, consistent with previous studies (Canicio and Ibáñez, 1999).

Aside from physically constraining Ebro Delta change, our models also highlight the physical mechanisms responsible for the generation of observed morphology. Simulations point to the formation of spits during delta growth, potentially responsible for delineating the Clot, Encanyissada, and Tancada lagoons, with orientations distinct from large recurved La Banya and El Fangar spits that formed from reworking of abandoned lobes. The suggested changes to the Ebro Delta leading to the formation of the observed spits are possible under a constant sea-level and sediment supply, caused by river avulsions.

Using constraints from the delta evolution model together with a river profile model, we find that a combination of increased fluvial flood discharge and fluvial sediment supply that started approximately 2100 years BP is the most likely cause of a rapid and sustained period of deltaic growth over the last 2100 years. The rapid growth of the Ebro Delta plain is likely not solely caused by an increase in fluvial flood discharge because that would greatly increase fluvial incision. Instead, a combined change in discharge and sediment supply can be felt instantaneously at the river delta while persisting for millennia without a significant channel profile change. A combined change in discharge and sediment supply can also, depending on their respective timing, generate two periods of floodplain aggradation (Fig. 10i).

In this study we have highlighted a few factors that particularly influence the sensitivity of our results. Fluvial sediment supply, wave climate characteristics, and the littoral sediment transport constant all have a first-order effect on gross delta plain shape as reflected by the fluvial dominance ratio $R$. Shoreface characteristics such as the depth of closure and the basin depth determine how the delta plain responds to sediment flux changes. Timescales of the river profile model are particularly sensitive to the median channel bed grain size and the upstream boundary location: the average distance between the delta and environmental change in the drainage basin. In all of the simulations presented here, we have chosen average, representative model parameters frequently mentioned in the literature, with model results showing the broad first-order agreement with other studies of Ebro Holocene evolution.

As the Ebro Delta moves into the 21st century, the effects of sea-level rise and river damming will increasingly manifest themselves in the delta morphology (Sánchez-Arcilla et al., 2008). Even though it is tempting to run our delta model into the 21st century, we emphasize that future delta shoreline predictions should include the effects of sea-level rise (e.g., Sánchez-Arcilla et al., 2008). However, by quantifying potential effects of historical land-use and climate change on historical delta evolution, simple models such as the one discussed here might be able to assess long-term future deltaic change and help guide management decisions (Giosan et al., 2014).

The data for this paper are available upon request from the authors.

The authors declare that they have no conflict of interest.

Acknowledgements

This study was supported by NSF grant EAR-0952146. This article greatly benefitted from constructive comments from reviewers Jorge Guillen, Enrica Viparelli, Edward Anthony, and Jose Jimenez, and from the associate editor Simon Mudd. We are grateful to Carles Ibáñez and Lluis Jornet for discussions on the Ebro Delta. Edited by: Simon Mudd Reviewed by: Enrica Viparelli, Jorge Guillen, Jose A. Jiménez, and Edward Anthony