The Numerical Modeling of River Deltas Eric W.H. Hutton, James P.M. Syvitski, and Yusuke Kubo Environmental Computation and Imaging Group, INSTAAR University of Colorado at Boulder Boulder CO, 80309-0450, USA Phone: (303) 492-6233 Fax: (303) 492-6388 Eric.Hutton@colorado.edu An important way to investigate the interwoven complexity of earth surface processes is through the judicious application of surface dynamic models. Communitydeveloped models (HydroTrend and SedFlux, for example) are upgraded continuously in their representation of physical processes. Their general use and accuracy, however, often depend on the availability of appropriate numerical boundary conditions, such as the time dependency of spatial climatic factors (e.g. sediment supply, rainfall, and temperature), and geomorphic basin properties (e.g. digital elevation models (DEMs), geology, and distribution of ice-fields and lakes). Uncertainty We see three main sources of uncertainty in the numerical modeling of river-delta systems, and for geologic systems in general. These are uncertainty in the identification and parameterization of the dominant sediment transport processes, errors in validation data, and the boundary conditions that constrain the model simulation. When developing a numerical model that simulates the formation of stratigraphy within a river delta, one must identify the major processes that act to move sediment. The relative importance of these processes will vary from region to region and so one must continually evaluate the applicability of a model when applying it to new areas. A model that distributes sediment solely using a surface plume may work well for one area, but will work miserably when applied to a region with high wave energy that is able to rework the sediment completely. Models that have identified the major processes can also become non-portable when these processes are being modeled incorrectly. That is, the models give the right results for the wrong reasons (the term equifinality describes this concept). In this case, the model may work well in a specific location and agree with field data there, but because the process is being modeled incorrectly, when moved to a new location the model may no longer be useful. Typically, models that exhibit equifinality are those that are based on empirical relationships. Uncertainty in the numerical modeling of geologic systems also arises due to uncertainty in the field measurements against which the modeler verifies his or her model. This type of error is more easily characterized than the previously mentioned source as oftentimes field measurements come with an estimate of error. We can only be as certain in the model as our certainty of the verification data. The final source of error that we have identified are the initial conditions used for a model simulation, and will be the topic of this paper. Uncertainty in the starting conditions (initial bathymetry), and driving forces (sea level, or sediment supply, for example) will propagate through a model simulation. The characterization of this error is not well understood and is currently a hot research topic that deserves attention. The first-order forces that drive delta building are sediment supply and accommodation. Most delta-building models recognize these two factors to be boundary conditions that the model user will specify. That is, these models will not generate their own sea-level curve or sediment supply (although it may modify it). This is a fine assumption since the building of a delta will not affect, in any significant way, climate that will in turn affect sea level or sediment supply. However, delta-building models should be able to modify these driving forces. For instance, compaction of sediments and isostatic subsidence due to sediment loading will alter accommodation. In addition, deposition and erosion of the evolving delta plain will be able to modify the sediment supply to the sea. Obtaining boundary conditions Numerical models oftentimes require detailed boundary conditions that are location specific, or do not correspond to direct field measurements. These types of parameters may be difficult to establish. However, because of the recent release of various gridded global data, estimates of the larger, first order boundary conditions are typically readily available. Sediment Supply The amount of sediment that is given to a delta is critical to predicting accurately both its internal and external structure. Currently many river gauges throughout the world provide detailed time series of sediment flux. Unfortunately, many more remain ungauged and even those with gauges typically do not provide data for more than a few decades. The paleo-modeler is left with even fewer data. Where detailed discharge data are unavailable, the modeler is now able to look toward globally gridded data sets. Recently, NOAA/NCDC and GHCN/University of Delaware have released global discharge data from GRDC and University of New Hampshire (1960-1994). This dataset provides globally gridded data of discharge with a .5º spatial resolution. In addition, the relations developed by Syvitski et al. (2003) allow the combination of global grids of basin area (or discharge), basin relief, and temperature to produce new grids of sediment flux. Accommodation While the total sediment supply establishes the volume of the delta, the amount of space that is available for the sediment to accumulate (accommodation) shapes the delta. Typically, one thinks of accommodation as the space available for the sediment to accumulate underwater. However, many delta systems trap large amounts of sediment on the delta before ever reaching the sea (Meade, 1996) and so this too can be though of as accommodation. Accommodation is a balance between sea level, bathymetry, and tectonics. Digital Elevation Models (DEMs) are widely available for land and range in resolution from tens of meters to kilometers. For global bottom bathymetry, one turns to the 2º by 2º ETOPO2 data set (Smith and Sandwell, 1997). Although this resolution will be too coarse for the modeling of many river deltas, it is still valuable as it is able to provide accommodation estimates. As with bathymetry, global sea-level curves (such as that obtained from the GRIP ice-core record) exist that may provide estimates of changing accommodation. However, because sea level curves will vary from region to region, site-specific sea-level curves should be used whenever possible. Example: The Po River Delta In our first example, we demonstrate the effect of accommodation on the delta system of the Po River since the last glacial maximum and compare it to the interpretation of the region by Cattaneo and Trincardi (1999). For this simulation, daily sediment fluxes were taken from Kettner and Syvitkski (2004), sea level from Asioli et al. (1999), and present-day bathymetry was used to approximate that of the Adriatic at the last glacial maximum. Figure 1 compares the sedflux2D model simulation with an interpretation derived from seismic data from the region. The general character of the final deposit is similar in both simulations. A large delta complex is built during the low stand of the last glacial maximum followed by a series of smaller deltas during sea level rise, and concluding with a second delta that is the Po River delta of today. Although the two larger deltas correspond well in size, shape, and location to those of the interpreted profile, the intermediate deltas seem to be shifted laterally along the profile. Figure 2 explores the possible impact of the major boundary conditions of sediment supply and accommodation on the deposit. First, to examine the effect that a changing rate of sediment input has on the deposit we have repeated the simulation using a constant flux of sediment to the basin. In this case (Figure 2, upper left), we note that similar deltas are built and so conclude that their general character are a result not of sediment supply, but of accommodation. The upper-right and lower-left figures of Figure 2 show the same simulation conducted again with constant sediment input but now the bottom bathymetry is set to a constant slope (Figure 2, upper right), and the sea level curve simplified to a constant rise (Figure 2, lower left). In the case of constant sea-floor slope, we note that a series of deltas build in similar water depths as those of the original run. These deltas then, are a result of variations in sea level change. Similarly, the bottom bathymetry of our profile is responsible for the delta system of Figure 1. Our sea level curve seems to produce deltas are at the correct water depth, but our bathymetry causes them to shift position laterally. Example: The Eel River Delta Our second example demonstrates the effect of the timing of sediment delivery to the Eel River delta in Northern California. The drainage basin of the Eel River is relatively small in the sense that convective storms are able to rain over the entire basin and so produce large floods. Unlike larger rivers, the Eel River moves most of its sediment over a small percentage of time. In fact, 90% of the sediment that leaves the Eel River is transported over 5.7% of the time. Here we have conducted two sedflux2D simulations to examine the effect that the timing of sediment transport events have on the stratigraphy and final structure of a delta. Both simulations use similar boundary conditions; the one difference being that the first simulation has used annually averaged river discharge data. Although the same amount of sediment is delivered to each of the systems, the way that the sediment is delivered has changed. Figure 3 shows a cross-section of grain size for the two model simulations. We note that the annually averaged data has caused the delta to prograde further into the sea than the daily time step did. When the daily discharge data are used, high-energy events introduce most of the sediment to the system. Because of an increase in river velocity, these high-energy flooding events are able to transport their sediment load further into the basin. As a result, the delta becomes thinner and longer. To illustrate this point further, Figure 4 shows synthetic cores of grain size from each of the simulations at 15m water depth. There are two main differences to notice here. The first is that the daily simulation has produced sediment layering that is significantly more variable. The second is that the overall core has coarsened. Again, when daily river data introduce sediment into the basin, the high-energy flood events distribute nearly all of the sediment and so move coarser grained sediments to deeper water. Conclusion Globally gridded data sets exist that can provide source terms for delta building models. Many times though, these data sets are of course resolution and may not provide the accuracy that the modeler desires. However, they are useful tools that provide estimates of total sediment flux to larger delta systems. We note that this sediment flux controls to first-order the size of the delta but that the timing of the flooding events can play a significant role in modifying a delta’s shape. Similar globally gridded data sets exist for seafloor bathymetry which, when coupled with a sea-level curve, give the boundary conditions that control accommodation. On a global scale, these data sets also are generally too coarse for the modeling of most delta systems. For these regions, they provide first-order estimates of accommodation. Finally, many factors contribute to the uncertainty of the numerical modeling of river deltas, not the least of which is uncertainty in sediment supply and accommodation. Attention should be paid to the estimation of this uncertainty in geologic numerical models. This is a hot research topic and a consensus of how uncertainty in these models should be measured and presented has not been reached. Certainly, because of the complexity of these models, it does not make sense to provide a single “standard deviation” for each model simulation. However, model results can be presented with a “best guess” simulation as well as “lower” and “upper” bound scenarios. References Asioli, A, Trincardi, F, Lowe, F, Oldfield, F (1999) Short-term climate changes during the Last Glacial-Holocene transition: comparison between Mediterranean records and the GRIP event stratigraphy. Journal of Quaternary Science, 14(4), 373-381 Cattaneo, A. and Trincardi, F., 1999. The late-Quaternary transgressive record in the Adriatic epicontinental sea: basin widening and facies partitioning. In: Isolated Shallow Marine Sand Bodies: Sequence Stratigraphic Analysis and Sedimentologic Interpretation (Eds. K. Bergman and J.W. Snedden). SEPM Spec. Publ. 64, 127-146. Kettner, A.J., and Syvitski, J.P.M., submitted 2004. Predicting Discharge and Sediment Flux of the Po River, Italy since the Late Glacial Maximum. IAS, Special Issue. Meade, R.H., 1996. River-sediment inputs to major deltas. In: Sea-Level Rise and Coastal Subsidence. J. Milliman and B. Haq (eds.). Kluwer, London. pp. 63-85. Morehead, M. D., Syvitski, J. P. M., and Hutton, E. W. H., 2001. The link between abrupt climate change and basin stratigraphy: a numerical approach. Global and Planetary Change, v. 28, pp. 107-127. Smith, W. H. F., and Sandwell, D. T., 1997. Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings. Science, v. 277, pp. 1956-1962. Syvitski, J.P.M., Morehead, M., Nicholson, M., 1998a. HYDRO-TREND: a climatedriven hydrologic-transport model for predicting discharge and sediment load to lakes or oceans. Com-put. Geosci. 24, 51–68. Syvitski, J.P., Morehead, M.D., 1999. Estimating river-sediment discharge to the ocean: Application to the Eel margin, northern California. Mar. Geol. 154, 13–28. Syvitski, J.P.M., Hutton, E.W.H. 2001. 2D SEDFLUX 1.0C: An advanced processresponse numerical model for the fill of marine sedimentary basins. Computers and Geoscience 27(6): 731-754. Syvitski, J.P.M., Peckham, S.D., Hilberman, R.D. and T. Mulder, 2003. Predicting the terrestrial flux of sediment to the global ocean: A planetary perspective. Marine Geology 162, 5-24. Syvitski, J.P.M., Kettner, A.J., Vörösmarty, C.J. and Green, P.A., 2004. Climate and Climate Variability as Input to Surface Dynamic Models. AGU Fall meeting, San Francisco, USA. Figure 1 -- Comparison of interpreted paleo-deltas of the Adriatic (modified from Cattaneo and Trincardi, 1999) a sedflux2D simulation of the Po River, Italy since the last glacial maximum. Cross−section of Grain Size 0 0 Cross−section of Grain Size 50 50 Depth (m) 150 Depth (m) 100 100 150 200 200 250 0 100 200 300 Y−Distance (km) 400 500 600 250 0 100 200 300 Y−Distance (km) 400 500 600 5 5.5 6 6.5 7 7.5 Grain Size( φ) 8 8.5 9 5 5.5 6 6.5 7 7.5 Grain Size( φ) 8 8.5 9 Cross−section of Grain Size 0 50 Depth (m) 100 150 200 250 0 100 200 300 Y−Distance (km) 400 500 600 5 5.5 6 6.5 7 7.5 Grain Size( φ) 8 8.5 9 Figure 2 – three experiments to examine the effect of sediment supply (upper left), sea level (upper right), and initial bathymetry (lower left) on the formation of deltas in the Po River, Italy since the last glacial maximum. Cross−section of Grain Size 0 20 40 0 20 40 Cross−section of Grain Size Depth (m) 60 80 100 120 0 2 4 6 8 Y−Distance (km) 10 12 Depth (m) 60 80 100 120 0 2 4 6 8 Y−Distance (km) 10 12 6.2 6.4 6.6 6.8 7 7.2 7.4 Grain Size( φ) 7.6 7.8 8 5.5 6 6.5 Grain Size( φ) 7 7.5 Figure 3 – Grain-size cross-sections for two simulations of the Eel River, California. Comparison of deltas built using (left) annual averaged sediment input, and (right) daily values for sediment input. The annually averaged data causes an increase in the progradation rate. 15 4 km 20 15 3.3 km 20 Water Depth (m) 30 Water Depth (m) 5.5 6 6.5 Grain Size (phi) 7 7.5 25 25 30 35 35 40 40 45 5 45 5 5.5 6 6.5 Grain Size (phi) 7 7.5 Figure 4 – Synthetic grain-size cores taken at 15m water depth for two simulations of the Eel River, California. Comparison of deltas built using (left) annual averaged sediment input, and (right) daily values for sediment input. The annually averaged data show a fining in grain size, and decrease in variability.