Towards a probabilistic assessment of process-based, morphodynamic models

Coastal Engineering 75 (2013) 52–63

Contents lists available at SciVerse ScienceDirect

Coastal Engineering journal homepage: www.elsevier.com/locate/coastaleng

Towards a probabilistic assessment of process-based, morphodynamic models M. van der Wegen a,⁎, B.E. Jaffe b, 1 a b

UNESCO-IHE, PO Box 3015, 2601 DA Delft, The Netherlands USGS Pacific Coastal & Marine Science Center, 400 Natural Bridges Drive, Santa Cruz, CA 95060, United States

a r t i c l e

i n f o

Article history: Received 4 September 2012 Received in revised form 28 December 2012 Accepted 3 January 2013 Available online 21 February 2013 Keywords: Estuarine geomorphology Morphodynamic modeling Probabilistic approach Model uncertainty San Pablo Bay

a b s t r a c t Recent advances in the development of numerical, process-based models have led to remarkable performance in reproducing measured decadal morphodynamic developments. The advantage of this type of models is that they have a detailed output allowing for a close analysis of relevant processes. Drawback is that the output is associated with a high level of presumed uncertainty, because of the large number of processes involved and the high quality level of input data required. This study aims to explore possibilities to assess uncertainty levels associated with process-based morphodynamic modeling. In a probabilistic approach we consider the outcome of an ensemble of runs including variations of model input parameters and forcing schematizations. We propose to evaluate model performance by both a skill criterion (How well does the model reproduce observed patterns?), a confidence criterion (How sensitive are model results to uncertain input?) as well as a combination of these criteria. This methodology provides an objective assessment of the performance of process-based morphodynamic models. In addition, it can determine which input parameters cause largest uncertainty in the model outcome. The San Pablo Bay case study shows that 60% of the modeled volume meets the skill and confidence criteria for the depositional period (1856–1887) and 46% for the erosional period (1951–1983). Approximately 50% of the volume allocation meets the confidence criterion for a 30 year morphodynamic forecast (1983–2013). Model results are sensitive to model input variations only to a limited extend. We attribute this to the high impact of the San Pablo Bay plan form and bathymetry. The forecast shows continuous erosion of the channel and on the northern shoals and a continuous accretion of the channel slopes, albeit more concentrated in the western part of the channel than in preceding decades. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Estuaries are valuable areas. Breeding fish and migratory birds are only examples of the estuarine ecosystem richness. Local fishery, aquaculture and tourism are important economical sectors that profit from the estuarine environment. Numerous ports are situated along estuaries, forming the logistical link between ocean trade and the hinterland. It is of utmost importance to understand estuarine morphodynamic processes so that impact of human interference (e.g. dredging of access channels to ports) and long-term processes (e.g. sea level rise, changing discharge regime) can be estimated and evaluated. Tidally fluctuating water levels, wind, wind waves and the interaction between salt and fresh water drive the hydrodynamics in estuaries to a great extent. The interplay between hydrodynamic forcing,

⁎ Corresponding author. Tel.: +31 152151811; fax: +31 152122921. E-mail addresses: [email protected] (M. van der Wegen), [email protected] (B.E. Jaffe). 1 Tel.: +1 831 427 4748. 0378-3839/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.coastaleng.2013.01.009

associated sediment transport and the estuarine geometry leads to the development of morphological features that range from small scale ripples to sandy channels systems, tidal mud flats and vegetated salt marshes (Perillo, 1995; Hibma et al., 2004). Many of the hydrodynamic and morphodynamic processes involved are highly non-linear, which makes morphodynamic understanding and prediction a challenge. Small differences in roughness (for example, caused by ripple formation) may play a significant role in the development of major channels. Although equilibrium of channel shoal systems seems to exist on short timescales of days to weeks, significant development may occur on a decadal timescale. Adequate prediction of morphological developments requires that models cover a range of interacting spatial and time scales (De Vriend et al., 1993a,b). Werner (2003) further introduces the concept of hierarchical modeling which assumes a hierarchy in processes that depend on each other, but do not interact. Parameterized small-scale processes can thus act as input for modeling larger-scale and longer-term processes. For any modeling effort, this raises the question on the level of process parameterization that is required. Some modeling strategies start from a highly schematized setting focusing on assumed

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

dominant forcing or empirical relationships, whereas other strategies include many detailed, non-linear, physical process descriptions and their interactions on beforehand (De Vriend et al., 1993a,b; Werner, 2003; Murray, 2003). This study focuses on the latter type of models further referred to as process-based models. Examples are Delft3D (Lesser et al., 2004) or ROMS (Warner et al., 2008). Over the last decades, process-based models have evolved into powerful tools that approach the complexity of reality. The models generate appealing, detailed and high resolution output, which allows for a close analysis of processes and forcing mechanisms. The drawback is that high quality model input is required which cannot always be guaranteed by measurements. Also, the amount of processes and input data required raises questions on the models' output uncertainty. The appropriateness of morphodynamic models for estuarine management is essentially related to the models' uncertainty (Vreugdenhil, 2006). Important aspects that contribute to this uncertainty are the parameterization and schematization of physical processes, the forcing and the initial conditions. For example, a morphodynamic model needs data on sediment characteristics, bed composition and sediment concentration boundary conditions. The more processes are included and the higher the level of process detail, the more data are required. Sometimes the process descriptions require such a high level of detailed data input that they cannot be covered by economically or logistically feasible measurement campaigns. Schematization (best-guess estimates) of model input is then necessary (Ganju et al., 2009; Ganju and Schoellhamer, 2010; Van der Wegen et al., 2011a,b). Literature review on morphodynamic modeling shows only limited attention to the uncertainty associated with process-based morphodynamic models. Baart et al. (2011) pointed to confidence intervals around morphological forecasts as the result of the propagation of uncertainty through a chain of models from weather (wind and pressure) to hydrodynamic (waves and tides) and finally morphodynamic predictions. Vreugdenhil (2006) and Pinto et al. (2006) presented examples on the uncertainty associated with different sediment transport formulations. Fortunato et al. (2009) evaluated ensemble averaged morphodynamic predictions in an estuarine environment and conclude that uncertainty increases with time and with larger sediment transports. Ruessink and Kuriyama (2008), Plant and Holland (2011a,b) and Van der Wegen et al. (2011b) showed that ensemble runs in morphodynamic predictions are suitable to indicate uncertainty levels, although only the latter apply this to an estuarine environment. The aim of the current work is to explore possibilities to assess uncertainty levels associated with process-based morphodynamic modeling with a focus on decadal morphodynamic developments in San Pablo Bay, California. 1.1. Approach In this study we apply a process-based model (Delft3D) to San Pablo Bay in California, USA. The reason to select this area is twofold. A unique dataset is available that describes bathymetric evolution in this area over the last 150 years on a 30 year interval. This dataset forms a rare opportunity to compare decadal morphodynamic modeling results with measured data. During the 150 year period, San Pablo Bay was subject to a large sediment pulse as a result of hydraulic mining and a sudden drop in sediment supply due to landward dam construction. The modeling results can thus be assessed under varying forcing conditions. Secondly, the area is subject to many phenomena (wind waves, salt-fresh water interaction, mud and sand transports) that need to be defined and which require process schematization and estimates of model input parameters. We assume that the associated model outcome uncertainty will be significant, so that the study gives a good impression on the value of this type of modeling efforts. Earlier morphodynamic modeling work in the San Francisco Estuary describes the decadal morphodynamic calibration and hindcast of

53

Suisun Bay (Ganju et al., 2008, 2009) and San Pablo Bay (Van der Wegen et al., 2011a,b) as well as a forecast of developments under different climatic change scenarios (Ganju and Schoellhamer, 2010). The current work continues the work of Van der Wegen et al. (2011b) who hindcasted morphodynamic development for the 1856–1887 period in San Pablo Bay. The settings of their model will be used to hindcast morphodynamic development of a period in which net erosion occurred (1951–1983) as well as to forecast 30 year morphodynamic development from the latest (in 1983) measured bathymetry onwards. We calibrate the modeled volumes and patterns by comparison to measured bathymetries and use the optimal model settings for the forecast. Sensitivity analysis reveals uncertainty bounds due to model input parameter variation. Finally, we propose a methodology to quantify model uncertainty in a single indicator. 1.2. Study site Literature provides detailed background on the physical and biological characteristics of the Bay-Delta system (i.e. Kimmerer, 2004). Here we focus on the main aspects concerning morphodynamic development. The Delta is known as the area where rivers draining the Central Valley and Sierras of California, including the Sacramento and San Joaquin Rivers, meet before discharging into the northeastern end of the San Francisco Estuary (Fig. 1). The river discharge is low (~200–500 m3/s) during the relatively dry season (summer/autumn) and much higher (up to 18,000 m3/s) during the wet season (winter/spring), although flow distribution and magnitude vary considerably during a year as well as over different years. The San Joaquin River is responsible for about 10–15% of the discharge and the Sacramento River for about 80%. The remaining discharge originates from minor tributaries. Magnitude and timing of high river flow are subject to human interferences by means of managed reservoir releases in the watershed and water export from the Delta for the purpose of fresh water supply to Southern California. The tide near the Golden Gate is highly irregular with strong diurnal components. Main tidal constituents are M2, O1, and K1. During high river flow no tidal influence is present at Sacramento (155 km upstream from the Golden Gate), whereas minor tidal fluctuations are observed during low river flow. The tidally averaged salt intrusion up the estuary depends primarily on fresh water flow and to a lesser extent on spring-neap tidal variations. The 2 psu isohaline is found most often in Suisun Bay. Dam water release during dry seasons is managed in such a way that that salt does not intrude farther landward. High river flows may cause salinities to decrease to less than 5 psu east of San Pablo Bay. Stratification and gravitational circulation occurs in particular during neap tides, when tidally-driven mixing processes are weak. The location of the largest salinity gradients varies depending on the river flow. San Pablo Bay, a sub-embayment in the northern Estuary, is circular with an area of about 250 km 2 and an average tidal range of about 1.5 m. It is rather shallow (average depth b2 m below mean sea level depth of 4 m close to the main channel) and muddy apart from the main sandy channel that transects the Bay from East to West conveying the river flow seaward. The Bay-Delta system has been subject to a dramatic change in sediment loads and morphodynamics over the last 150 years (Gilbert, 1917; Porterfield, 1980; Cappiella et al., 1999; Jaffe et al., 1998, 2007). Hydraulic mining during the Gold Rush in the mid 19th century caused excessive supply of sediment from the Sacramento River watershed towards the Delta and San Francisco Estuary. At the end of the 19th century hydraulic mining was put to a hold and reservoir dam construction took place in the watershed. As a result sediment loads decreased considerably. Wright and Schoellhamer (2004) reported a decrease of about 50% in the Sacramento River of measured sediment load in the 1957–2004 period. Bathymetric surveys from the 1850s to 1980s on a 30-year interval show that the morphology of San Pablo Bay has changed markedly

54

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

Fig. 1. Location of San Francisco Estuary and San Pablo Bay. Dashed (dotted) lines in the right figure represent ~9.1 m (~1.8 m) contour lines of the 1983 bathymetry.

(Jaffe et al., 1998, 2007). Deposition of more than a quarter billion cubic meters of hydraulic gold mining debris reduced the average depth of San Pablo Bay by 85 cm in the middle and late 1800s. In the late 1900s the intertidal flats narrowed and the major channel in the Bay deepened as more sediment was lost to the sea than entered from rivers. Processes of sediment redistribution caused the main channel to become narrower as well, a trend observed over the last 150 years. Jaffe et al. (2007) report minor dredging activities on the northern shoal to maintain the channel towards Petaluma River before 1922. An area of dredging in the SPB main channel after 1922 is delineated by the straight section of the 30-foot contour lines (Fig. 1). Between 1955 and 1983 on average about 0.2 * 10 6 m 3/year was dredged and released in San Pablo Bay (Ogden Beeman and Associates and Ray Krone and Associates, 1992). These volumes are about 25% of the net erosion volume averaged over the 1951–1983 erosional period. Since exact timing and location of the dredging and dumping activities were not reported and since the dredged volumes were deposited again within San Pablo Bay, we did not include these in the model runs.

2. Model description Lesser et al. (2004) provide details on the Delft3D morphodynamic model. Here we present a brief overview of the model setup described in more detail by Van der Wegen et al. (2011a,b) for the 1856–1887 depositional period.

The model grid covers an area from the seaward end of San Pablo Bay to the east side of Suisun Bay where it connects to the Delta area. It is a curvi-linear grid with grid sizes varying from 100 m to 600 m with 15 sigma layers over the vertical with higher resolution towards the bed and the open water surface to account for expected high velocity gradients (Fig. 2). With a sounding resolution of 250 m in 1856 increasing to 60 m in 1983 (Jaffe et al., 2007), most of the grid cells in San Pablo Bay have multiple bathymetric data. The cell bed level was determined by grid cell averaging or triangulation in case of limited data. Hydraulic boundary conditions were derived by running a larger Delft3D model (see supplementary material) covering the entire San Francisco Estuary, extending about 50 km into the ocean and including the (highly schematized) Delta area up to the area where hardly tidal influence is present. The large model allows defining landward boundary conditions by river flow and seaward boundary conditions by observed tidal constituents. We carried two large model runs with different landward boundary conditions representing a wet season with high river flow and a dry season with low river flow. The large model runtime would have been too long to carry out the sensitivity runs as applied in the current study. The hydraulic boundary conditions in the current study describe 1 month of high river flow (5000 m 3/s) followed by 11 months of relatively low flow (350 m 3/s). The high level of schematization reduced river inflow to three practical ‘tuning’ parameters (i.e. duration of high discharge compared to low discharge and the magnitudes of the high and low river discharges). Salt concentration is set constant

Fig. 2. Model grid.

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

at zero at the landward boundary and at 31 psu at the seaward boundary. Prevailing wind conditions are schematized by a diurnal sinusoidal signal varying from 0 at midnight to 7 m/s at noon uniformly distributed over the domain based on the wind climatology described by Hayes et al. (1984) and additional analysis of wind data from http:// wwwcimis.water.ca.gov/cimis/welcome.jsp. For the wet season and 6 months of dry season the wind comes from the west and for the remaining 5 months of dry season the wind comes from the south east mimicking the seasonal variations in wind field. Every hour the SWAN model uses wind and hydrodynamic data to generate a wave field and returns resulting adapted hydrodynamic parameters to the flow module. The model includes 3 sand fractions and 5 mud fractions as well as a model for bed composition administration. The transport of cohesive mud is modeled by the Partheniades–Krone formulations (Ariathurai, 1974; Krone, 1962, 1993). For the transport of non-cohesive sediment, Van Rijn's (1993) approach is followed. The non-cohesive sediment transport magnitude and direction are corrected for bed slope effects (Bagnold, 1966; Ikeda, 1982). Following the methodology by Van der Wegen et al. (2011a), we obtained an initial distribution of sediments across the model domain by running the model for 1 month without bed level changes, but allowing for a bed composition update. Sediment concentration boundary conditions were derived following the method proposed by Ganju et al. (2008). They related sediment loads to the estimated historical discharges bþ1

Q s ¼ aQ r

ð1Þ

In which Qs Qr a, b

sediment load, kg/s river discharge, m 3/s calibration coefficients with units depending on each others' values.

Ganju et al. (2008) suggest a constant value of b = 0.13 and a time varying value of a = 0.02 (before hydraulic mining) to 0.13 (when hydraulic mining stopped in 1884), with the parameter a slowly decreasing to approximately 0.03 for more recent decades. To enhance morphodynamic development a Morphological Factor (MF) multiplies every time step the bed level changes calculated from the divergence of the sediment transport field. This approach requires that the upscaled bed level changes (per time step) are small compared to the water depth, so that they have negligible influence on the hydrodynamics. Details on this methodology can be found in Roelvink (2006) and Van der Wegen and Roelvink (2006). For example, to reproduce 30 years of morphological development, the model calculates first 1

55

month of high river discharge with a MF of 30 followed by 4 months of low river discharge with a MF of 82.5 (smaller transports during the dry season allow for a slightly higher MF). A typical run takes about 1.5 days on a 3.0 GHz, 3.46 Gb RAM PC. 3. Model results 3.1. Description To hindcast the morphodynamic developments over the 1951– 1983 erosional period, we initially applied slightly modified model parameter settings as suggested by Van der Wegen et al. (2011b) for the 1856–1887 depositional period. To reflect the impact of the dam construction works in the watershed and the consolidation of mud in San Pablo Bay, we adjusted the river discharge during the wet season (from 5000 to 4000 m 3/s), the depth-averaged, land boundary mud sediment concentrations following suggestions by Ganju et al. (2008) (from 300 mg/l to 30 mg/l) and the erosion coefficient M (from 1 * 10 −4 to 5 * 10−5 kg/m2/s) to account for mud consolidation. We performed an ensemble of runs in which we systematically varied a number of model input parameters. While other parameters were kept constant at the standard settings, the sensitivity analysis shows the impact of varying a single model parameter value. Variations were determined within reasonable and possible range. These parameters were landward boundary sediment concentration (c) in mg/l, critical mud erosion shear stress (τcr) in N/m 2, erosion coefficient (M) in kg/m 2s, the mud fall velocity (w) in m/s, horizontal eddy diffusivity (D) in m 2/s, roughness parameter Manning (s/m 1/3) or Chézy (m 1/2/s), river discharge (Qr) in m 3/s, and wind speed in (m/s). There are several other model parameters that could have been included in this sensitivity analysis, but this work focuses on a methodology to deal with uncertainty rather than reporting on the impact of varying a ‘complete’ set of model parameters. Model results are compared in terms of measured volumes in Fig. 3 and in terms of the Brier Skill Score (BSS, see Appendix A and Sutherland et al. (2004) for clarification) as a proxy for the quality of erosion and sedimentation patterns in Fig. 4. In the latter figure we applied a value of 0.2 m 3 for δ which is a conservative estimate to account for the measurement error. Jaffe et al., 2007 suggest that bias errors of the bathymetric soundings (from inaccuracies in determining the relation of MLLW datums for different surveys and from grid representation differing from the sounding values) and random error (associated with sounding inaccuracy) may locally add up to (0.04 + 0.07 = ) 0.11 m. A BSS of 1 reflects perfect modeling results, whereas negative values suggest that the model results are worse than referring to the initial bathymetry (it would have been better not to model at all). The BSS is determined by shifts in phase (same patterns but on a different location), deviating amplitude (more or less sedimentation

Fig. 3. Modeled deposition and erosion volumes (in m3/year) of (a) 1856–1887 depositional period, (b) 1951–1983 erosional period with model parameters of depositional period, (c) 1951–1983 period with adapted model parameter settings, (d) 1983–2015 period. Thick black lines indicate measured values.

56

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

Fig. 4. Brier Skill Scores (BSS, see Appendix A for definition) of (a) 1856–1887 depositional period, (b) 1951–1983 erosional period with model parameters of depositional period, (c) 1951–1983 period with adapted model parameter settings. The dotted line includes the effect of measurement uncertainty (see Eq. (A.2) with δ = 0.2 m3).

or erosion) or different total volume (average bed level change). In an extreme case, a perfectly reproduced pattern may thus lead to a negative BSS if the average bed level change is not modeled correctly or when the patterns' location is slightly different. A drawback of the BSS is its sensitivity to small modeled deviations from measurements in areas where the measured volumes are also small. As a result, the BSS covering the full domain may be poor, despite a good BSS in areas of large volumetric changes. Although interpretation of the BSS needs careful attention, it offers an objective indicator that appears especially useful when comparing the performance of different model runs and not focusing too much in its actual value. For the erosional period modeled volumes systematically underestimate measured developments (Fig. 3b), whereas Fig. 4b shows a negative BSS for all runs of the erosional period. To improve the results we carried out further, extensive sensitivity analysis. Combining the parameter settings of the best performing runs for the BSS (for example by doubling the fall velocity and increasing the critical shear stress) would not automatically lead to an improved BSS. Also, improved performance in BSS could lead to a worse predicted volume. Selecting the best runs, we adjusted the different parameters one by one and performed an ensemble run including again variations of the other model parameters. Performing a few of these cycles leads to improved performance with parameter settings shown in Table 1. We used these settings for the 1951–1983 hindcast period and the 1983–2013 forecast period. The sand fractions (only present in the channel) are defined coarser and the lateral bed slope factor (accounting for downhill sediment movement on a slope perpendicular to the flow direction) is lower. We attribute this to the observation that channel morphodynamics have a higher impact on the overall results during the erosional period. Probably, the new parameter settings would have improved the deposition period performance as well, albeit to a limited extent. Since shoal erosion by wind waves accounts for a large part of the morphodynamic development during the erosional period, model results are also more sensitive to the wind schematization. Fig. 3c shows that the modeled erosion volumes with improved parameter settings fit better to the measured values than Fig. 3b and Fig. 4c now shows largely positive BSS's for the 1951–1983 period. Fig. 5 shows an overview of measured and modeled erosion and sedimentation patterns for the 1856–1887 period, the 1951–1983 period and the 1983–2013 period with the standard model parameter settings. Major observed and measured deposition occurred along the channel slopes leading to a narrowing of the channel during the depositional, erosional and future periods. During the erosional period observed erosion on the North Western shoals and in the edges of the main channel is reproduced quite well. Channel slope accretion appeared mainly at the eastern parts in the depositional period, but

moved to the western part of San Pablo Bay during the erosional period and the future period. For a further discussion on the processes governing the modeled patterns reference is made to Van der Wegen et al. (2011b, submitted). Improving the model results appeared to be a time demanding task since combining best scoring parameters does not necessarily lead to an increase in model skill. We do not claim that we found the most realistic or the best possible parameter settings. For example, the wind amplitude magnitude seems quite high and may be compensated in runs with a relatively low value for the erosion coefficient M. We found that results for the 1951–1983 period were particularly sensitive to the wind direction and magnitude. However, detailed wind data are not readily derived from nearby measurement stations because of the large impact of geological features around San Pablo Bay redirecting the wind. Furthermore, model results may be improved by more detailed forcing conditions, such as including extreme events (storms, river discharge). Confirming the work of Van der Wegen et al. (2011b), varying the model parameters within reasonable limits leads to limited volume and BSS variations. The sensitivity analysis distinguishes dominant

Table 1 Overview of model parameter settings. Period Sand characteristics Sand fractions diameter (μm) Lateral bed slope factor (αbn) Mud characteristics Inflow concentrations (mg/l) Critical shear stress (N/m2) Erosion coefficient (kg/m2/s) Sediment fall velocity (mm/s) Hydrodynamic forcing Horizontal eddy diffusivity (m2/s) Roughness n (s/m1/3) Manning C (m1/2/s) Chézy River discharge (m3/s) Wind Amplitude (m/s) Direction Hydrodynamic run period (months) Morphological factor (−)

Wet season Dry season Wet season Dry season

1856–1887

1951–1983 and 1983–2013

150, 300, 500 100

150, 700, 1600 25

300 0.3–0.8 1 *10−4 0.16–0.4

30 0.36–0.84 5 * 10−5 1

1 0.2

0.5

55 5000 7

2100 9

West and South East ~1 4 30 82.5

South West and South East 0.5 ~2 80 160

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

57

Fig. 5. Erosion (−) and sedimentation (+) patterns in meter in which warm colors indicate deposition (a,b) measured and (c,d,e) modeled for (a,c) the 1856–1887 depositional period, (b,d) the 1951–1983 erosional period and (e) the 1983–2015 period. Contour lines represent the 1856 (a,c), 1951 (b,d) and 1983 (e) measured bathymetry at a 5 meter interval.

from secondary forcing mechanisms. For example, the suspended sediment concentration and river discharge at the landward boundary have significant effect for the 1856–1887 period, whereas they do not change model results for the 1951–1983 period and the 1983–2013 period. The latter periods have become much more sensitive to the erosion factor M and the definition of the wind field. Visual comparison of the erosion and sedimentation patterns for the variations in (a,b) the erosion coefficient (c,d) wind amplitude (Fig. 6) shows that patterns remain strikingly similar, although erosion and deposition volumes may vary considerably. Apparently, the model covers an important governing factor in the morphodynamic development, namely the interaction between hydrodynamics, the bathymetry and the plan form of San Pablo Bay. 3.2. Probabilistic interpretation of model results In this section we further elaborate a methodology to address the model outcome uncertainty. For a number of runs we derive the mean (ensemble averaged) model outcome (μ) as well as its standard deviation (σ) per model grid cell. An uncertainty criterion (for example μ/σ ≥ 3) reflects when variations are small compared to the mean result. Combining this “model uncertainty score” with a BSS threshold (for example BSS ≥ 0.5) gives an indication of the model performance in terms of model uncertainty and its skill. A “model performance indicator” (MPI) reflects the percentage of modeled volume that is likely to be deposited with an acceptable skill. The modeled bed level change in a grid cell averaged over a number of runs [i= 1 ..n] (μmod) is defined by

μ mod

n 1X ¼ Δzmod;i n i¼1

with a standard deviation (σmod) defined by 2

σ mod ¼

ð3Þ

in which Δzmod,i is the realization (modeled bed level change) of a particular run. We apply a Confidence Index (CI) to assess the level of uncertainty of the model outcome in a grid cell slightly modifying the concept of Fortunato et al. (2009): CI ¼ 1−

σ 2mod μ 2modabs

ð4Þ

In which (μmodabs) reflects the mean of the absolute modeled bed level changes to correct for low mean values as a result of both positive and negative realizations: μ modabs ¼

n   1X Δz  mod;i n i¼1

ð5Þ

A negative CI value implies high uncertainty (the standard deviation is larger than the mean value). Positive values mean that uncertainty is relatively low, with 1 being the perfect score (no uncertainty). A Confidence Performance Index (CPI) can be defined as the ratio of the displaced volume that meets a certain confidence criterion and the total modeled displaced volume: j¼m X

CPI ¼ ð2Þ

n  2 1X μ mod −Δzmod;i n i¼1

j

j

j

μ modabs ⋅a ⋅CP

j¼1 j¼m X j¼1

( where

j j μ modabs ⋅a

CPj ¼ 1 for CI≥ CIth CPj ¼ 0 for CIb CIth

ð6Þ

58

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

Fig. 6. Modeled erosion (−) and sedimentation (+) patterns (in meter) for 1951–1983 erosional period in which warm colors indicate deposition (a) 50% decreased erosion coefficient M (b) doubled erosion coefficient M (c) 25% wind amplitude decrease (d) 25% wind amplitude increase. Contour lines represent the 1951 measured bathymetry at a 5 meter interval.

Where (j) refers to the grid cell number, (m) refers to the total number of grid cells in the domain, (a) refers to the area of a grid cell, and CIth is a confidence threshold value. Following the concept of the BSS (Sutherland et al., 2004) we define a Skill Index (SI) in a grid cell as follows: SI ¼ 1−

ðμ mod −Δzmeas Þ2 Δz2meas

ð7Þ

j¼m X

in which the suffix “meas” refers to a measured quantity. The SI differs from the BSS in the sense that the SI considers a local (grid cell) value, whereas the BSS considers a domain averaged score. Eqs. (2)–(4) refer to variations of modeled values only, whereas Eq. (7) refers to a comparison of modeled values and measurements. A Skill Performance Index (SPI) can be defined as the ratio of the displaced volume that meets a certain skill criterion and the total modeled displaced volume: j¼m X

SPI ¼

j

j

j

μ modabs ⋅a ⋅SP

j¼1 j¼m X

( where

j

μ modabs ⋅a

j

volumes. Similar to Eq. (A.2) we could include a correction for the measurement error in Eq. (7) leading to slightly larger values for the SI and SPI (see also the larger BSS in Fig. 4 when the measurement error is taken into account). By combining both CPI and SPI a “Model Performance Index” (MPI) is defined by the displaced volume that meets both the CI and SI criteria divided by the total modeled volume that is displaced.

SPj ¼ 1 for SI≥ SIth j SP ¼ 0 for SIb SIth

ð8Þ

j¼1

In which SIth is a skill threshold value. The BSS and SI are very sensitive to relatively small modeled deviations from measurements in areas of small measured differences. Interpreting model results in terms of BSS and SI implies that every location has an equal weight. The advantage of the SPI is that it includes a correction for volumes. As a result, areas with weak SI performance and small measured volumes have less weight in the SPI in favor of areas with large measured

MPI ¼

j

j

j

μ modabs ⋅a ⋅MP

j¼1 j¼m X

( and

j j μ modabs ⋅a

j

MP ¼ 1 for CI≥ CIth ∧SI≥ SIth j MP ¼ 0 for CIb CIth ∨SIb SIth

ð9Þ

j¼1

A higher MPI value indicates a higher model performance because of a smaller model uncertainty and/or a better fit to measurements given criteria for model uncertainty and skill. It is noted that a lower uncertainty in model performance accomplished by a new set of runs does not automatically lead to a higher MPI given the same performance criterion, since it does not necessarily imply an improvement of the BSS (and vice versa). A second remark concerns the “perfect” model with a BSS of 1. This model would have a MPI smaller than 1 due to the uncertainty associated with the model parameter variations. Indeed, there may be one run (a mean performance) that has the perfect mix of input parameters to reach the ultimate BSS. However, the MPI presumes uncertain model input and probabilistic output, so that a BSS of 1 should be regarded as a matter of coincidence. In the absence of measured data (e.g. for future predictions), the BSS and the MPI are not suitable. In such a case only the uncertainty part of the model performance is relevant.

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

59

Fig. 7. (a) 1856–1887 measured erosion (−) and sedimentation (+) pattern in meters (warm colors indicate deposition); (b) 1856–1887 mean (ensemble averaged) modeled erosion and sedimentation pattern in meters (warm colors indicate deposition); (c) standard deviation of ensemble averaged erosion and sedimentation pattern in meter, truncated above 1; (d) CI, truncated below 0, see Eq. (3); (e) SI, truncated below 0, see Eq. (6); (f) Area (in white) that meets both confidence criterion (CI ≥0 see (d)) and the skill criterion (SI ≥ 0.35 see (e)). Contour lines represent the 1856 measured bathymetry at a 5 meter interval.

Fig. 8. (a) 1951–1983measured erosion (−) and sedimentation (+) pattern in meters (warm colors indicate deposition); (b) 1951–1983mean (ensemble averaged) modeled erosion and sedimentation pattern in meters (warm colors indicate deposition); (c) standard deviation of ensemble averaged erosion and sedimentation pattern in meter, truncated above 1; (d) CI, truncated below 0, see Eq. (3); (e) SI, truncated below 0, see Eq. (6); (f) Area (in white) that meets both confidence criterion (CI≥0 see (d)) and the skill criterion (SI≥0.35 see (e)). Contour lines represent the 1951 measured bathymetry at a 5 meter interval.

60

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

Based on the above mentioned indicators we can now analyze the performance of our model in terms of skill and uncertainty. Figs. 7 and 8 illustrate the results of this analysis for the 1856–1887 and the 1951–1983 periods respectively by presenting (a) the measured erosion and deposition pattern, (b) the modeled mean patterns, (c) their standard deviations, (d) the CI values and (e) the SI values. Finally, (f) show the area (in white) that meets a confidence criterion (CI ≥ 0) and the skill criterion (SI ≥ 0.35). Since it is a forecast Fig. 9 only includes (b) the modeled mean patterns, (c) their standard deviations and (d) the CI values. A higher standard deviation seems to be related to a larger mean deposition (the eastern part of the channel for the depositional period and the north-western part of the channel for the erosional and future periods). Most of the area has a positive CI, whereas positive SI covers less area. In particular, the model answers skill and uncertainty criteria for the depositional period on the North Eastern shoals and along the South Eastern channel slopes and for the erosional period on the North Western shoals and along the North Western channel slope. It is difficult to draw general conclusions from Figs. 7 and 8, since the performance is quite scattered over the domain. Good performance may even be due to locally very small erosion or sedimentation volumes. Therefore, apart from the performance over the area, model results are evaluated in terms of displaced volumes as well. Fig. 10(a) and (b) presents SPI and CPI for different SI and CI threshold values, whereas (c) and (d) show MPI values for the depositional and erosional period respectively. For example the MPI for the depositional period is obtained by the ratio of the domain

integrated product of Fig. 7(f) and (b) and domain integration of (b). By comparing Fig. 10(a) and (b) it becomes clear that the limiting factor for high MPI scores is the SPI. An SI threshold level of 0 only results in 62% (52%) of the measured volume being covered by the modeled volume for the depositional (erosional) period, whereas low CI threshold values lead to near 100% confidence. This implies that the uncertainty related to the model input has a smaller impact on the model performance than the skill of the model. The depositional period leads to higher SPI and CPI (more skill and confidence). This suggests that the applied model process description and forcing schematization performs less during the erosional period. For a given confidence threshold level, the erosional period and 1983 forecast have a lower CPI, implying that more recent bathymetries (with erosional conditions) lead to more uncertainty in model outcome. The MPI scores given in Fig. 10(c) and (d) indicate that, by applying SI and CI threshold values of 0.35 and 0 respectively, about 60% (46%) of the modeled volume meets the skill and confidence criteria and for the depositional (erosional) period. Given the limited change in observed forcing conditions over the past 50 years, we consider that the forecast based on the 1983 bathymetry will perform similarly as both hindcasts, i.e. that the model will predict approximately 50% of the volume allocation correctly and that uncertainty due to model process and forcing schematizations will be limited. This implies continuous erosion in the channel near Point San Pablo and Carquinez Strait and on the northern shoals and a continuous accretion of the channel slopes, albeit that accretion becomes more concentrated in the western part of the channel (compare Fig. 3(c) and (d) and Fig. 5(d) and (e)).

Fig. 9. (a) 1983 Measured bathymetry in meter; (b) 1983–2015 mean (ensemble averaged) modeled erosion (−) and sedimentation (+) in meters (warm colors indicate deposition); (c) standard deviation of ensemble averaged erosion and sedimentation pattern in meter, truncated above 1; (d) CI, truncated below 0, see Eq. (3); Contour lines represent the 1983 measured bathymetry at a 5 meter interval.

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

61

Fig. 10. (a) Skill threshold value versus SPI. (b) Confidence threshold value versus CPI. Red line reflects the 1856–1887 depositional period, blue line reflects the 1951–1983 erosional period and the black line reflects the 1983–2015 future period. (c) MPI for depositional period versus SIth and CIth. (d) MPI for erosional period versus SIth and CIth.

4. Discussion We present a methodology to assess the quality of morphodynamic model predictions. For the case study in the current research model skill levels appear to have a higher impact on the model performance than the uncertainty of model input. This, however, depends also on subjective threshold criteria applied to assess skill and uncertainty. Another important aspect is the impact of the selection of model parameters that are part of the ensemble. For example, wind direction and magnitude, extreme events of wind and river flow or the impact of the highly irregular tide are not considered extensively. Including these may significantly improve the skill of the model. The question is whether a larger ensemble set will lead to different conclusions. Model input variations within reasonable limits show limited qualitative variations in the reproduced erosion and sedimentation patterns. Apparently the model configuration covers governing morphodynamic processes already. We argue that the interaction between the hydrodynamic model on the one hand and the plan form of San Pablo Bay and its bathymetry on the other hand determines the morphodynamic process to a large extent. For example, the rocky outcrops at the landward and seaward end fix the location of the

main channel to a high degree. The resulting (tide residual) spatial shear stress gradients may be indicative already where sediment settles and where it erodes. Our findings do not necessarily hold for other, freer systems that are less confined by the plan form. Model performance may be improved by a more detailed grid covering the interaction between hydrodynamic processes and shape of San Pablo Bay in a better way. A better performance of the erosional periods may be obtained by a better description of the bed composition and bed erosion processes. The CPI methodology may also be used to assess the impact of different model parameters on the model outcome uncertainty. We determined the mean, standard deviation and CPI of three realizations of a specific input parameter (i.e. the standard, minimum and maximum values). Although this number of realizations is small, it is useful enough to explain the methodology. Fig. 11 shows that wind and river discharge (Qr) variations account for the largest model uncertainty during the depositional period. Variations in wind and the erosion factor (M) cause largest uncertainty for the erosional and future periods. Decreasing the uncertainty of these model parameters (for example by more detailed measurements or a better schematization) will have largest impact on increasing the confidence in the model outcome.

Fig. 11. Confidence threshold values versus CPI for different input parameter variations for (a) depositional period, (b) erosional period and (c) future period.

62

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

Table 2 CPI0.8 values per period and model input parameter.

c τ M w D n Qr Wind

Depositional

Erosional

Future

0.980 0.983 0.981 0.955 0.986 0.962 0.918 0.951

1.000 0.989 0.498 0.945 0.993 0.949 1.000 0.551

0.999 0.973 0.499 0.903 0.981 0.920 0.998 0.520

the 1983 San Pablo Bay bathymetry will predict approximately 50% of the volume allocation with sufficient skill and that uncertainty due to model process and forcing schematizations is limited. This implies continuous erosion in the channel and on the northern shoals and a continuous accretion of the channel slopes, albeit more concentrated in the western part of the channel than in preceding decades. A better determination of the wind climate and the erosion parameter will have most impact on improving the confidence in the model results. Acknowledgments

A useful standard to compare the sensitivity of model parameters could be CPI0.8 indicating the CPI value of a certain parameter with CIth =0.8. Table 2 gives an overview of these values. A suggested classification of the impact of the model results by the specific parameter uncertainty is 1 > CPI0.8 >0.9 — very limited, 0.9 >CPI0.8 > 0.7 — weak, 0.7 > CPI0.8 > 0.5 — reasonable, 0.5 > CPI0.8 > 0.3 — high and 0.3 > CPI0.8 > 0 — very high. The CPI value spread for the different model parameters also becomes larger in the erosional and future periods compared to the depositional period. This confirms the larger variations in BSS observed in Fig. 4 for the erosional and future periods. We attribute this effect to changing morphodynamic conditions. The depositional area is governed by a major sediment supply that governs the morphodynamic development. The morphodynamics during the erosional and future periods are more subtle and subject to multiple processes and parameters that were, apparently, of less importance during the depositional period. For example, San Pablo Bay has become net erosional so that morphodynamics more and more depend on unknown bed sediment characteristics. Furthermore, the erosion and sedimentation volumes become smaller so that the effect of model limitations and bias becomes relatively large. However, not all high BSS variations are reflected by the CPI curves as well. Compare for example the high BSS variations and the high CPI of the critical shear stress (τcr). This is probably due to the large impact of τcr-variations in areas that have only limited volumetric changes. 5. Conclusions Literature review shows that considerable progress has been made in process-based morphodynamic modeling last decades. Modeling studies of different case studies show that measured (decadal) developments have been modeled with significant skill. However, process-based models depend heavily on process description, model input and forcing schematization leading to an inherent uncertainty of the models' outcome. Only limited attention has been paid to address and assess this uncertainty. We present a methodology to assess model skill compared to measured developments (Skill Performance Index, SPI) and to evaluate model performance in terms of uncertainty (Confidence Performance Index, CPI) at the same time. The Model Performance Index (MPI) combines both indices and quantifies the modeled volume percentage that meets a skill criterion as well as an uncertainty criterion. The described case study of San Pablo Bay, California, shows a MPI of 60 to 46% for a decadal depositional and erosional period respectively. MPI variations related to model parameter uncertainty are limited compared to the impact of the model skill. This suggests that further performance improvement should focus on improvement of skill and that input uncertainty has a limited effect. Determining the CPI as the result of variations of a single model input parameter reveals which parameters have highest impact on the model uncertainty. Such an analysis may act as a guideline which model parameters should be determined more adequately to increase confidence in the model results. Given the limited change in observed forcing conditions over the past 50 years, we have confidence that a 30 year forecast based on

The research is part of the US Geological Survey CASCaDE climate change project (CASCaDE contribution 32). The authors acknowledge the US Geological Survey Priority Ecosystem Studies and CALFED for making this research financially possible. Friendly reviews by Dano Roelvink (UNESCO-IHE, Deltares and TU Delft), Ali Dastgheib (UNESCOIHE), Ilgar Safak (USGS) and two anonymous reviewers considerably improved the quality of the work. Appendix A In order to assess the skill of morphodynamic models Sutherland et al. (2004) suggest the use of the Brier Skill Score (BSS). For the current study the BSS is defined as follows: D

2

ðΔvolmod −Δvolmeas Þ   BSS ¼ 1− Δvol2meas

E ðA:1Þ

in which Δvol mod meas

volumetric change compared to the initial bed, (m 3) modeled quantity, measured quantity

and the 〈〉 denote an arithmetic mean (a spatial average in this case). A BSS of 1 is a perfect model, whereas lower values indicate poorer model performance. Sutherland et al. (2004) note that the BSS is unbounded at the lower limit and that the BSS can be extremely sensitive to small changes when the denominator is low. Following work by Murphy and Epstein (1989), Sutherland et al. (2004) further elaborate on a methodology to decompose the BSS in terms of an amplitude error, a phase error, or a deviation from the average value. In order to account for the effect of measurement errors Van Rijn et al. (2003) suggest the following extended BSS:

BSSvR

D E ðjΔvolmod −Δvolmeas j−δÞ2   ¼ 1− Δvolmeas 2

ðA:2Þ

in which δ (m 3) is the volumetric measurement error and in which |Δvolmod − Δvolmeas| − δ is set to zero if |Δvolmod − Δvolmeas| b δ. Van Rijn et al. (2003) proposed a classification of BSS and BSSvR values as presented in Table 3.

Table 3 BSS classification.

Excellent Good Reasonable Poor Bad

BSS

BSSvR

0.5–1.0 0.2–0.5 0.1–0.2 0.0–0.1 b0.0

0.8–1.0 0.6–0.8 0.3–0.6 0.0–0.3 b0.0

M. van der Wegen, B.E. Jaffe / Coastal Engineering 75 (2013) 52–63

References Ariathurai, C.R., 1974. A finite element model for sediment transport in estuaries. PhD thesis, University of California, Davis. Baart, F., Van Gelder, P.H.A.J.M., Van Koningsveld, M., 2011. Confidence in real-time forecasting of morphological storm impacts. Journal of Coastal Research SI 64, 1835–1839 (Proc. 11th Int. Coastal Symp.). Bagnold, R.A., 1966. An approach to the sediment transport problem from general physics. U.S. Geological Survey Professional Paper 422-I. Cappiella, K., Malzone, C., Smith, R., Jaffe, B.E., 1999. Sedimentation and bathymetric change in Suisun Bay 1876–1990. U.S. Geological Survey Open File Report 9–563. De Vriend, H.J., Zyserman, J., Nicholson, J., Roelvink, J.A., Pechon, P., Southgate, H.N., 1993a. Medium-term 2DH coastal area modeling. Coastal Engineering 21 (1–3), 193–224. De Vriend, H.J., Capobianco, M., Chesher, T., De Swart, H.E., Latteux, B., Stive, M.J.F., 1993b. Approaches to long term modeling of coastal morphology: a review. Coastal Engineering 21 (1–3), 225–269. Fortunato, A.B., Bertin, X., Oliveira, A., 2009. Space and time variability of uncertainty in morphodynamic simulations. Coastal Engineering 56–8, 886–894. http://dx.doi.org/ 10.1016/j.coastaleng.2009.04.006. Ganju, N.K., Schoellhamer, D.H., 2010. Decadal‐timescale estuarine geomorphic change under future scenarios of climate and sediment supply. Estuaries and Coasts 33, 15–29. http://dx.doi.org/10.1007/s12237-009-9244-y. Ganju, N.K., Knowles, N., Schoellhamer, D.H., 2008. Temporal downscaling of decadal sediment load estimates to a daily interval for use in hindcast simulations. Journal of Hydrology 349, 512–523. Ganju, N.K., Schoellhamer, D.H., Jaffe, B.E., 2009. Hindcasting of decadal-timescale estuarine bathymetric change with a tidal-timescale model. Journal of Geophysical Research 114, F04019. http://dx.doi.org/10.1029/2008JF001191. Gilbert, G.K., 1917. Hydraulic mining debris in the Sierra Nevada. U.S. Geological Survey Professional Paper 105 (148 pp.). Hayes, T.P., Kinney, J.J., Wheeler, N.J.M., 1984. California Surface wind climatology. California Air Resources Board, Aerometric Data Division. (73 pp.). Hibma, A., Stive, M.J.F., Wang, Z.B., 2004. Estuarine morphodynamics. Coastal Engineering 51 (8), 765–778. http://dx.doi.org/10.1016/j.coastaleng.2004.07.008. Ikeda, S., 1982. Lateral bed load transport on side slopes. Journal of Hydraulics Division, ASCE 108 (11), 1369–1373. Jaffe, B.E., Smith, R.E., Zink Toressan, L., 1998. Sedimentation and bathymetric change in San Pablo Bay 1856–1983. U.S. Geological Survey Open File Report 98, 0795. Jaffe, B.E., Smith, R.E., Foxgrover, A.C., 2007. Anthropogenic influence on sedimentation and intertidal mudflat change in San Pablo Bay, California: 1856–1983. Estuarine, Coastal and Shelf Science 73, 175–187. http://dx.doi.org/10.1016/j.ecss.2007.02.017. Kimmerer, W.J., 2004. Open water processes of the San Francisco Estuary: From physical forcing to biological responses [online]. San Francisco Estuary and Watershed Science 2 (1) (article 1. [Available at http://repositories.cdlib.org/jmie/sfews/vol2/iss1/art1]). Krone, R.B., 1962. Flume studies of the transport of sediment in estuarial shoaling processes. Final Report Hydraulic Engineering Laboratory and Sanitary Engineering Research Laboratory. University of California, Berkeley. Krone, R.B., 1993. Sedimentation revisited. Nearshore and Estuarine Cohesive Sediment Transport. In: Mehta, A.J. (Ed.), AGU, Coastal and Estuarine Studies, pp. 108–125. Lesser, G.R., Roelvink, J.A., Van Kester, J.A.T.M., Stelling, G.S., 2004. Development and validation of a three-dimensional morphological model. Coastal Engineering 51 (8), 883–915. http://dx.doi.org/10.1016/j.coastaleng.2004.07.014.

63

Murphy, A.H., Epstein, E.S., 1989. Skill scores and correlation coefficients in model verification. Monthly Weather Review 117, 572–582. http://dx.doi.org/10.1175/ 1520-0493(1989) 117. Murray, A.B., 2003. Contrasting the goals, strategies and predictions associated with simplified numerical models and detailed simulations. In: Wilcock, P.R., Iverson, R.M. (Eds.), Prediction in Geomorphology. Geophys. Monogr. Ser., 135. AGU, Washington, D.C, pp. 151–165. Ogden Beeman and Associates and Ray Krone and Associates, 1992. Sediment Budget Study for San Francisco Bay. Final report prepared for U.S.Army Corps of Engineers, San Francisco District, under contract DACW07-89-D-0029, USA. Perillo, G.M.E., 1995. Geomorphology and Sedimentology of Estuaries. Definitions and Geomorphologic Classifications of Estuaries. Development in Sedimentology 53. Pinto, L., Fortunato, A.B., Freire, P., 2006. Sensitivity analysis of non-cohesive sediment transport formulae. Continental Shelf Research 26 (15), 1826–1839. Plant, N.G., Holland, K.T., 2011a. Prediction and assimilation of surf-zone processes using a Bayesian network: Part I: Forward models. Coastal Engineering 58 (1), 119–130. Plant, N.G., Holland, K.T., 2011b. Prediction and assimilation of surf-zone processes using a Bayesian network: Part II: Inverse models. Coastal Engineering 58 (3), 256–266. Porterfield, G., 1980. Sediment transport of streams tributary to San Francisco, San Pablo, and Suisun Bays, California, 1909–1966. U.S. Geological Survey WaterResources Investigations Report, 80-64 (91 pp.). Roelvink, J.A., 2006. Coastal morphodynamic evolution techniques. Coastal Engineering 53, 177–187. Ruessink, B.G., Kuriyama, Y., 2008. Numerical predictability experiments of cross-shore sandbar migration. Geophysical Research Letters 35, L01603. http://dx.doi.org/ 10.1029/2007GL032530. Sutherland, J., Peet, A.H., Soulsby, R.L., 2004. Evaluating the performance of morphological models. Coastal Engineering 51, 917–939. Van Rijn, L.C., 1993. Principles of sediment transport in rivers, estuaries and coastal seas. AQUA Publications, the Netherlands, p. 535. Van Rijn, L.C., Walstra, D.J.R., Grasmeijer, B., Sutherland, J., Pan, S., Sierra, J.P., 2003. The predictability of cross-shore bed evolution of sandy beaches at the time scale of storms and seasons using process-based profile models. Coastal Engineering 47, 295–327. Van der Wegen, M., Jaffe, B.E., Roelvink, J.A., 2011a. Process-based, morphodynamic hindcast of decadal deposition patterns in San Pablo Bay, California, 1856–1887. Journal of Geophysical Research 116, F02008. http://dx.doi.org/10.1029/2009JF001614. Van der Wegen, M., Dastgheib, A., Jaffe, B.E., Roelvink, J.A., 2011b. Bed composition generation for morphodynamic modeling: case study of San Pablo Bay in California, U.S.A. Ocean Dynamics. http://dx.doi.org/10.1007/s10236-010-0314-2. Vreugdenhil, C.B., 2006. Appropriate models and uncertainties. Coastal Engineering 53 (2-3), 303–310. Warner, J.C., Sherwood, C.R., Signell, R.P., Harris, C.K., Arango, H.G., 2008. Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model. Computers & Geosciences 34, 1284–1306. Werner, B.T., 2003. Modeling landforms as self-organized, hierarchical dynamical systems. In: Wilcock, P.R., Iverson, R.M. (Eds.), Prediction in Geomorphology. Geophys. Monogr. Ser., 135. AGU, Washington D.C., pp. 133–150. Wright, S.A., Schoellhamer, D.H., 2004. Trends in the sediment yield of the Sacramento River, California, 1957–2001. San Francisco Estuary and Watershed Science 2 (Article 2).