1/4/09 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! * The river has a tendency to 'stick' to the valley walls, and to increase the length of contiguous contact over time. Still to do are: * Accumulate erosion area and coverage statistics (e.g. 2d joint distribution) over time and plot river location with respect to grid, map, and satellite images. * * Formulate strategy and methodology for coordinating automated runs and handling results from multiple computers (well, at least 2). Full simulations yield from 1-5 Mb of data (depending on number of cutoffs), so that a full 'population' of 1000 runs might require several Gb of storage and management to easily view the final or intermediate results of any of the simulations for side-by-side comparison, and for generating illustrations for the final report and website. The search for the one true simulation parameter combination continues, but now it can be perfomed automatically using prepared command files (i.e. 'batch files'). I've added the ability for the simulation program to read and parse commands to load and initialize files, change parameters, perform runs, save the results as both binary data and images, reset the simulation, change parameters, and continue with more runs, etc... This makes it possible, for example, to perform 2d grid searchs of parameter spaces to see which yield the 'best' looking results. (There is even a 'goto' command that allows a batch run to be stopped and resumed on a particular simulation at a later time.) Final river geometries must still be reviewed manually after all runs have completed. Subsets of results which appear promising can then be used to refine the parameter search. Ideally this iterative process will yield a reasonable starting point for the large set of runs to be started next month. A programmed aggregate test was performed to examine the relationship between geometry and feature size with respect to the two integration parameters ldist (lag distance) and sdist (sum distance), which are used to calculate the effective curvature at each point. A small Python program was used to generate a command file of 36 simulations which tested 6 by 6 value combinations of each of these parameters. The resulting command file looks like this: Of these 36 results, 11 were selected for full 1000 step runs. Here are some samples of the new results: Next, 9 new runs using combinations of maximum width and flow having values of {800, 600, 400} m x {800, 1200, 1600} m^3/s were generated using ldist = 1.00 and sdist = 1.00. Here are some samples: Next, using ldist = 1.00, sdist = 1.00, wmax = 600.00, flow = 1200.00 as a base, 9 new runs using combinations of depth and grain diameter having values of {3.5, 3.0, 2.5} m x {0.0005, 0.001, 0.0015} m were generated. Although not much is published about the bed composition of the MNRR, the equations we are using for A (also known as K) are appropriate for a 'sandy' bed rather than a 'gravel' bed, so we would expect the grain sizes to be on the order of 1 mm. Here are some samples of these results: Nevertheless, it now appears clear that the end is in sight: 1. The search for a set of initial simulation parameters which yield meandering behavior similar to the actual MNRR is converging, and 2. That by slightly varying each of the 6 individual parameters {ldist, sdist, wmax, flow, depth, diam} in combination, a distribution of simulations about a common base can be generated. For example, by choosing 3 nearby values of each of the parameters, 729 slightly different simulations can be performed. During the past 2 weeks I have run additional sets of N x M parameter combinations in order to narrow down the search for a starting point and set of parameter values to use in the large set of runs. I think I've come up with such a set, all of which produce 'reasonable' or 'possible' looking configurations of the MNRR, although this evaluation is based primarily on my own subjective interpretation of the resulting geometry and dynamics. Using runs performed with paired combinations of 5 x 5 values (25 combinations) of {ldist and sdist}, {width and depth}, and {flow and diameter} (i.e. 3 independent 2-dimensional subspaces of the 6-dimensional parameter space), I've decided on the following discrete values: This will yield 1x3x4x4x4x4x3 = 2304 total runs, requiring an estimated 33.6 days to complete, if all goes well. Assuming some unanticipated problems, this should realistically allow me to finish the runs by early- to mid-November, allowing sufficient time for the analysis and write-up. While the runs are in progress, I will be doing other things (see below). Here are some representative results. Except as noted, all use the base values of: 3.4.9 Common area coverage and definition of migration corridor ---------------- Finally, the program can now accumulate coverage arrays for the entire boundary area, currently at a resolution of 0.5 x 0.5 km. These arrays show both the maximum extent of the river in 2d, and also the total amount of time spent by the river (although not necessarily contiguous time) at any point in the valley. Here is the final 100 yr configuration and a series of coverage extents for intermediate times during a particular simulation. Cutoffs are omitted from the intermediate displays, as well as values near the inlet, which would otherwise tend to dominate the plot. Blue is shortest residence time (0.1 year), red longest (up to 97.5 years near inlet), with a mean value of 5.12 years: It is hoped that the accumulation of many (e.g. 1000) such arrays will enable the calculation of time and location statistics (mean and standard deviation) at each 3d point which can be used to generate probability curves (contours) which define the overall 'migration corridor' as a function of time. However, bear in mind that all of the simulations shown above have been initialized at 10x the current erosion rate so as to better visualize the differences in geometry. Simulations run at a nominal rate will probably not be as wide-ranging. Contour plots I've created some contour plots of the final coverage and erosion, calculated areas covered by different fractions of all simulations 3.4.10 Target areas in domain -------------------------------------------------- 'Target' areas within the river valley can be defined which correspond to towns, roads, airstrips, parks, and other geographic landmarks: The 100 year coverage array can be overlayed on these targets to see which of them are encroached upon by the river, and to what degree: The database can also find examples of individual rivers which reach and cover several designated targets at some time during the simulation, but which miss others. In the following plots, two towns (in green and red at upper left of valley) are eventually covered by the same river in 4 different examples: This information, along with the 2d coverage maps shown above, can be used to assess the relative danger from the river to specific locations within the valley at different times during the next 100 years. For example, here are contour plots showing the percentage (95%, 75%, 50%, 25%, 5%, 1%) of all simulated rivers reaching points in the valley. Towns lying on or inside these lines have a 1/100, 5/100, etc..., chance of being overrun by the river during 100 years: Using the 0.5 x 0.5 km coverage grid, various regions ('targets') can be defined corresponding to geographic areas on the maps and satellite images. For example, in the following plots 13 such areas are defined:
3.4.10.1 Location -------------------------------------------------------------- This information, along with the 2d coverage maps shown above, can be used to assess the relative danger from the river to specific locations within the valley at different times during the next 100 years. For example, here are contour plots showing the percentage (95%, 75%, 50%, 25%, 5%, 1%) of all simulated rivers reaching points in the valley. Towns lying on or inside these lines have a 1/100, 5/100, etc..., chance of being overrun by the river during 100 years: Using the 0.5 x 0.5 km coverage grid, various regions ('targets') can be defined corresponding to geographic areas on the maps and satellite images. For example, in the following plots 13 such areas are defined:
3.4.10.2 Percent coverage ------------------------------------------------------ 3.4.10.3 Time of invasion ------------------------------------------------------ 3.4.11 Comparison of 2 meandering models --------------------------------------- 85% of the river simulations were performed using the Johannesson-Parker meandering model of 1989, which is based on the solution of 2 integral equations. Most of the simulated rivers shown in the previous pages have been produced using this method. However, as a comparison, 540 simulations were perfomed using a much simpler (and less 'physical') meandering model which is based on the difference in circumference of two circles of different radii (e.g. around river bends). Shown below are 4 examples from this set of runs --the longest (5061), the shortest (4351), covering the most area (5278), and a typical 'average' river (5373): 4.0 Results -------------------------------------------------------------------- 4.1 Summary -------------------------------------------------------------------- A total of 3612 river simulations have been performed, using slightly different input parameters (e.g. width, depth, flow, etc...) each time, so as to generate a range of river shapes and behaviors. During each simulation, the location of the river was tracked and accumulated in several 2d arrays. After all simulations were complete, the 2d coverage arrays for all rivers could be summed up and plotted on either the topographic maps or satellite images of the valley. These plots show the percentage of all simulations which passed through various points in the valley, where those points which were occupied at least once by all rivers are shown in red. This area is designated as the probable '100 year migration corridor' for the river: Initial run results During the past 11 days I have completed 1064 simulations of the river, using the combination of parameters described previously. Runs are performed on 2 Apple Macintosh computers having 2.0 and 2.4 GHz Intel dual-core CPUs running OS X 10.5.5. When I am not doing anything else with the computers, 2 simultaneous simulations (as Unix processes) can be run, taking less combined time than 2 sequential simulations. Because the runs are progressing well, I have decided to perform an additional 768 simulations using the same combinations of the 5 variables {width, depth, diam, flow, sdist}, but with an initial erosion rate of 1.8 ha/yr/km, so that 4 values of initial erosion will be used: {1.2, 1.8, 2.4, 3.6}. This is partly because many of the simulations decrease their erosion rate over time (so that some end with erate << 1.2), but also because after viewing the initial results I feel that an intermediate value between 1.2 and 2.4 is needed, whereas there is less final difference between runs started at 2.4 and 3.6 ha/yr/km. This will yield a new total of 3072 rather than 2304 distinct simulations. At the current rate of 96/day, this should take a total of 32 days, allowing all runs to be completed by the end of October, barring unforseen problems. I have also written a simple database to collect, display, and manipulate the results, in both numerical and graphical form. This database has a programming and query syntax similar to some other common programs, so I have elected to call it 'RiverSquirrel', after 'my squirrel' (MySQL) or 'post grey squirrel' (PostgreSQL). A database creation program peruses all simulation files and accumulates several numerical results, creating a large table of values. The database query program can perform the following operations on the table: * Print entries and fields in any order. * Sort entries about any field. * Select subsets of entries using numerical relationships and combinations of 'and' and 'or' clauses. * Plot the river for any parameter combination. * Accumulate and plot common coverage and average erosion as 2d arrays for subsets of entries. Currently the query interface is a programming one (from Python), but it could be modified to either read text command files, or to take interactive input from the terminal. At present it is sufficient to perform rudimentary searches and display subsets of results as text and images. It is expected that this interface and its capabilities will evolve as necessary over the next few months as the data are analyzed. Intermediate run results As of this morning, I have completed 2615 out of 3072 simulations (85%). I have also added new features to the database and plots, and have performed some new analyses of the results. Highlights of these are presented below. If all continues to go well, all simulations should be completed well before the end of October. So far there are no indications that any runs will need to be performed again, but I am keeping some time available at the end of the month in case it is needed. The project schedule calls for all additional analyses to be performed during November (as well as any other experiments such as variable-width or multiple-thread simulations if time permits), and for December to be used only for writing the final report and completing the website. Completed run results As of the end of October, all river simulations --a total of 3612-- have been completed. This is considerably in excess of the 1000 runs I had hoped to perform several months ago. 3072 simulations were performed using the Johannesson-Parker meandering model, and another 540 were then performed using the simpler Circumferential Speed method (in which the excess speed fraction is equal to sigma in the JP model) for use as a comparison to the JP results. 4.2 Examples of simulations ---------------------------------------------------- Here are several fields from the table showing the 6 independent, discrete, input parameters (sorted by width->depth->diam->flow->sdist->erate) and 5 dependent, continuous, output values:
Using the database, one can plot, say, the longest and shortest final river:
Or the ones with the most cutoffs (also the longest river) and with the longest average cutoffs:
Or the 'wildest' and 'tamest' rivers (by final erosion rate):
4.2.1 Extrema ------------------------------------------------------------------ By looking at the combinations of input parameters in the previous printout, one can see that the 'wilder' rivers tend to be narrower, shallower, have more discharge, larger grain size, depend on more upstream geometry, and have higher initial erosion rates (pretty much as expected). As before, the database can show the longest and shortest rivers:
In addition, rivers can be sorted by total length = current length plus lengths of all cutoffs:
And by erosion rate:
As before, the database can be used to select specific rivers --here the longest and shortest, longest and shortest including cutoffs, highest and lowest erosion rates, and greatest and least area coverage (many of which are the same river):
4.2.2 Nominal ------------------------------------------------------------------ Here are some 'more nominal' results (for length < 195 Km and erate < 1.95):
Here are some more nominal rivers which are within 0.05 standard deviations of the mean length and erosion rate:
Two near-average runs are:
4.3 Histograms and 2/3d regressions -------------------------------------------- Relationships between all combinations of the output values can also be plotted (in the following, the single-point outlier in many plots is run 143 shown above): Although these plots show only about 35% of all runs (and only 2 initial widths so far), they do inspire a few observations/hypotheses: 1. There appears to be a downward trend in river length, coverage and eroded area, and final erosion rate with an increase in river width. 2. There appears to be a positive trend in both coverage and eroded area with respect to length, superimposed on a more random variation. Limiting these plots to specific subsets of the initial parameters might show where this transition occurs. However, although the relationship between river length and erosion rate is also positive, it appears that almost all lengths can be generated from fairly low values of erosion. This may be due to the increased number of cutoffs which are also associated with increased erosion rate, and which limit river length. 3. As expected, there are strong positive correlations between coverage and eroded area, and between both and the final erosion rate. Further sorting and selecting of combinations of input variables (rather than all combinations lumped together as shown above) might show differences in these plots and allow one to see where specific trends form and change (i.e. for fixed subsets of the other variables, do variations in a single or pair of input variables lead to specific variations in output values?). Here are plots of pairs of output variables showing best-fit lines and r values (Pearson product-moment correlation coefficients of Y x X): Plots can now be colored to show the different locations of points with respect to a third variable (in this case the average river width). Note the differences in erosion rate, length, coverage area, and eroded area with respect to narrow vs. wide rivers: If length is replaced by total length, the distributions and fits shown above become much better. Thus, the scattered 'cloud' in the previous plots is due to cutoffs reducing the length of the river. Here the third axis is flow: Here are some histograms of a few variables across all runs. Note the differences between initial and final erosion rates, and between river length and total length (including cutoffs): 4.4 Common coverage & erosion -------------------------------------------------- For this project, the most important result might be the cumulative common coverage of the river (i.e. showing the fraction of all simulations which pass through 2d points in the valley), and the average cumulative total erosion over 100 years (also at 2d points in the valley): Notes: * The common coverage shows the fraction of all simulated rivers that occupy points in the valley at least once over 100 years. This goes from red (points where all 1064 rivers pass through), to green (where 50% of all rivers go), to dark blue (where only one river goes), to black (where no river goes). These values could also be interpreted as probabilities of occupancy by at least one hypothetical river over the course of 100 years. The limit (contour) of some specific probability value might then be called the 'migration corridor'. * Total erosion is plotted as hectares per (0.5 km)^2. Note that 0.25 km^2 is 25 ha, so that points showing more than 25 ha total erosion were, on average, eroded by the river more than once. Also note that this value is the total continuously integrated (each 0.25 year) erosion at a point, not the difference between final and initial land areas (as are reported in EJ06). Thus, the river could stay mostly in one place and 'wiggle' slightly, and it would continue to erode surface area, possibly beyond the capacity of the landscape to continue to provide sediment. In this analysis, it is assumed that surface runoff, or some other mechanism, will continue to provide sediment necessary to satisfy any cumulative erosion. 4.4.1 As 2d image and contours ------------------------------------------------- Here are plots of this coverage area at 0, 20, 40, 60, 80, and 100 years: Here is a plot of the final average erosion, in hectares per half kilometer squared: Finally, here is the common coverage plot (to within 500 m x 500 m) for all 2615 runs: By interpolating and smoothing this array by 4x, I can plot the coverage to within a 125 m x 125 m grid. This does not really contain any additional information, but is more visually appealing: Here is the final 100 year coverage for the JP runs. This array can be overlayed on either the topographic map or satellite image mosaics, and can be drawn with any opacity between 0 and 1: Contour plots I've created some contour plots of the final coverage and erosion, calculated areas covered by different fractions of all simulations, and tried to clean up the graphics a bit to make all elements easier to see. I don't know how this looks on your display (I have an Apple LCD), so let me know if you have problems seeing the results. Here are just the contour limits for coverage by 95%, 75%, 50%, 25%, 5%, and 1% of all rivers: Here are contours overlayed on the maps and satellite images. The minimum non-zero value of the blue area outside of the 1% line is 0.000326 (1 out of 3072 rivers, although not necessarily the same one everywhere). Note that several towns are outside of the 1/100 danger line, while Mission Hill is at 5/100, etc...: 4.4.2 Coverage wrt time -------------------------------------------------------- Here are the coverages at 0, 20, 40, 60, 80, and 100 years: 4.4.3 Statistics --------------------------------------------------------------- Here are approximate areas between the contour lines at both the final 4x interpolation (125 x 125 m) and at the original accumulated coverage (500 x 500 m). The original area of the full-width digitized river (avwidth = 747.84 m, length = 95353.75 m) is around 71.31 km^2, although all the width-limited simulations start with less than this value:
4.4.4 Average erosion ---------------------------------------------------------- Here is the final average erosion (in ha per 0.25 km^2): Here are contours for the total eroded area (in ha per 0.25 km^2): 4.4.5 Coverage with respect to parameter subsets ------------------------------- Since the coverage for all runs contains fairly active river examples, it might be useful to accumulate and plot coverage for specific subsets of these runs. For example, all runs having less than or equal to the mean length and erosion rate (n=971): All runs which do not cross the valley to Vermillion (n=989): All runs which do not have any cutoffs (n=212): By providing coverage examples for different subsets of rivers, it might be possible to better define the '100 year migration corridor' under different conditions and assumptions. However, while different plots show differences in coverage values somewhat less than 1.0, they do not significantly change the shape of the 'core' set of points reached by all rivers (shown in red). Coverage for input parameter subsets This update shows differences in 100 year common coverage for subsets of runs having fixed values of the 6 different input parameters, arranged in increasing sequences of these values. Thus, it does not represent new information about the simulations, but is a filtering of the full output so as to look for trends in the data, in this case indicated by the 95% to 1% contour lines of the coverage. Shown are subsets of all runs for both the Johannesson-Parker (u1b) and the simpler Circumferential Speed (sigma) simulations for comparison. Comments are interspersed with the plots, discussion at bottom. Numerical values indicating cumulative and differential coverage area between contour lines are also shown, along with statistics about when subsets reach various targets (e.g. towns) in the valley. The first plot shows differences in coverage (JP method) with respect to maximum river width. Note that the total coverage decreases with increasing river width, although the 95% area increases slightly with width. The 50% cumulative area is virtually the same for all widths: The next plots shows coverage with respect to river depth. In this case, the total coverage and the 95% to 50% areas are greatest for the shallowest rivers: The next plot shows coverage with respect to bed particle diameter. In this case, all areas increase with increasing particle size: The next plot shows coverage with respect to flow. As expected, the total coverage increases with increasing flow, but this is not necessarily the case for the fractional areas. The 95% to 50% areas appear to be greatest for a moderate flow, and then decrease with larger flows: The next plot shows coverage with respect to upstream summation distance (in local river widths). This clearly indicates that the coverage increases, and is somewhat sensitive to, the amount of upstream river geometry that affects the lateral migration at any point: Finally, as expected, the major contributing parameter to all coverage areas appears to be the initial (set by program to specific value) erosion rate: During the simulation, the erosion rate is not adjusted, and may end up either higher or lower than its initial value for specific rivers. However, since this seems to be the primary variable responsible for variations in coverage, it might be worthwhile to repeat the previous analyses for a single initial erosion rate (e.g. 1.8 ha/yr/km) to see if the observed trends remain. Here is a sample of the numerical output of the coverage areas and target acquisition for the cases shown above. Click the links below for the full files:
4.5 Target invasion ------------------------------------------------------------ The database can track the course of every individual river over time to see when, and how much area of, each target it will reach and cover: Rivers can be sorted according to whether they reach certain points in the valley. In this case, specific points (0.5 km x 0.5 km areas) near the Vermillion shore: 4.5.1 Percent coverage ---------------------------------------------------------
Information about which rivers occupy grid cells contained in any of the targets at any time during the simulation can be accumulated by the database, and an 'average incursion' or 'probability of invasive danger' over the entire 100 year period can be calculated for each area:
Here is a plot showing the final coverage overlayed with the target areas: Note that 4 out of 13 targets are completely missed by the rivers. One of these is Meckling, right in the middle of the valley, although the coverage map shows that at least one river grazes this target to within 0.5 km. Searching for this river yields number 191 from above, which is also the longest and most active instance (and possibly the least realistic): So, while Meckling escapes being flooded, it doesn't do so by much (and may offer some potential scenic riverside property investment opportunities :-) ). 4.5.2 Time of invasion ---------------------------------------------------------
The database can also accumulate information about when targets are invaded, both by specific rivers and by all rivers on average:
A -1 signifies that the target is not reached in 100 years. Note that, as defined, Ponca state park is partially within the initial river, so that the incursion time is 0. However, the differential erosion of park grid area could be tracked over time, if desired. 4.5.3 With respect to parameter subsets ---------------------------------------- The runs can then be divided into sets which (for example) either reach or do not reach the Vermillion shore. Here are the shortest (total length) river which crosses the valley, and the longest river which does not:
Here is a sample of the numerical output of the coverage areas and target acquisition for the cases shown above. Click the links below for the full files:
4.5.4 Example rivers reaching specific targets --------------------------------- The runs can then be divided into sets which (for example) either reach or do not reach the Vermillion shore. Here are the shortest (total length) river which crosses the valley, and the longest river which does not:
The database can also find (for example) all rivers that reach Mission Hill or Gayville, sort them by arrival time, and plot several examples that reach both targets:
4.6 Comparison of JP and CS results -------------------------------------------- These runs, when combined, produce a coverage plot which is similar to that of the Johannesson-Parker simulations: Here is the final coverage and erosion for the Circumferential Speed runs: Here are some extreme selections from the CS database (longest and shortest, longest and shortest including cutoffs, highest and lowest erosion rates, and greatest and least area coverage):
Note in these cases that the longest, most active, and greatest area rivers are now distinct, and that the CS rivers do not suffer from the oscillation instability which is common to the wider, less active, JP rivers. Two near-average runs are:
Here is the coverage (CS method) with respect to maximum river width. Unlike the JP results, the CS coverage is not as sensitive to width, and does not always respond in a monotonic fashion: Here is the coverage with respect to river depth. Unlike the JP results, in this case the coverage is relatively insensitive to depth, as depth is not a primary variable determining migration in this method: Here is the coverage with respect to flow. Again, coverage is relatively insensitive to flow, as flow is not a primary migration variable in this method: The CS method is based on water speed, which varies as the ratio of flow / depth, so that a similar range of values can be generated by holding one parameter constant while varying the other. The CS method does not use bed particle diameter as a parameter. The next input parameter is upstream summation distance. In this case, the CS coverage response is the same as for JP (in fact, the two methods share the same upstream summation code): The CS method also has a similar response to the JP method with respect to initial erosion rate, although this is an actual result of the two meandering methods, rather than because of similar code. (Note that the CS method was not run using an initial erosion rate of 3.6, as 1.2 - 2.4 ha/yr/km were deemed sufficient): Here is a sample of the numerical output of the coverage areas and target acquisition for the cases shown above. Click the links below for the full files:
Full output
Full output Some comparisons between the two meandering methods are noteworthy: 1. In general, the JP total coverage area is greater than for CS, with the following exceptions: (1) flows of 800 m^3/s, and (2) all initial erosion rates (1.2 - 2.4 ha/km/yr). 2. However, in all cases, the CS 95% area was larger than the corresponding JP 95% area. This suggests that most of the CS rivers tend to stay closer to the original river, even though total coverage area of all CS rivers is greater than for JP rivers at the same initial erosion rate. A possible explanation for this is suggested by plots showing the initial and final distributions of erosion rates for the JP runs (top) and the CS runs (bottom): Although both sets of runs start with discrete erosion rates, the JP method seems to redistribute the erosion according to an exponential or Poisson process, yielding a final histogram which has a characteristic unimodal peak, and where some rivers have increased their erosion rates while others have decreased. The CS runs, however, all increase their rates, although it is unclear what redistribution mechanism is at work. 4.6.1 Instabilities ------------------------------------------------------------ You may have noticed that some river plots show instabilities (a short sequence of oscillations), usually in the beginning of the river near Yankton. I don't know the exact reason for these, but they appear to be a product of the JP method possibly resonating with the local width and/or discretization and smoothing distances. Instabilities occur only for a small subset of parameter combinations, and usually only in the beginning of the river (I have noted previously that the JP method needs to 'get up to speed' with increasing s distance). Rather than try to debug this situation (if possible), I have elected to leave the instabilites as they are for 2 reasons: 1. They don't seem to affect the downstream shape of the river, and 2. Since this is 'science', I feel a duty to report on all aspects of the work, including those which are surprising in an annoying way. The JP method is very 'finicky', and I have already put in considerable effort to make it work at all, including ameliorating the other problems I have noted in previous blogs (e.g. 041308 at bottom). Here is a histogram of the maximum lengths of instabilities (sequential oscillations) found in all simulations, and a plot of instability length vs. erosion rate. Note that instabilities occur predominantly at the lowest erosion rates: Here are the rivers with the longest instabilities:
Here are the fraction of all runs with instabilities up to a given length. Note that there are no rivers with zero instabilities. Two rivers having instabilities of length 4 (which are virtually invisible) are plotted. Since these account for more than 90% of all runs, I'm not going to worry about it too much:
5.0 Discussion ----------------------------------------------------------------- 5.1 Applicability of approach used in study ------------------------------------ Everything considered, I am looking for simulations that resemble the river as it was 100 years ago, before human intervention and management began to change and constrain its shape: A total of 3612 river simulations have been performed, using slightly different input parameters (e.g. width, depth, flow, etc...) each time, so as to generate a range of river shapes and behaviors. During each simulation, the location of the river was tracked and accumulated in several 2d arrays. After all simulations were complete, the 2d coverage arrays for all rivers could be summed up and plotted on either the topographic maps or satellite images of the valley. These plots show the percentage of all simulations which passed through various points in the valley, where those points which were occupied at least once by all rivers are shown in red. This area is designated as the probable '100 year migration corridor' for the river: Next month I will use the geometry generated this month to initialize some meandering simulations, including the JP89 method from April. I will still need access to depth, flow, and bed composition data. Robb Jacobson has provided several references to hydraulic data, and additional historical figures appear in the 2006 USGS report. Initially, depth can be set to a constant, although Robb has also mentioned that he has access to bathymetric data and models. Ideally, some cross-stream depth profiles would be useful, as it is apparent from the maps and satellite photos that the MNRR is definitely not an idealized trapezoidal channel. Furthermore, in many places the calculated width is likely to be considerably in excess of the actual flow width, due to the presence of bars and other shallow shoals, and the thalweg may be far from the calculated centerline. It may be possible to partly compensate for this based on the bar and vegetated bar percentage estimates in the 2006 report. In any case, these are all variables with which to experiment next month. 1. The JP method does not really produce results that look like the MNRR. It produces 'river-like' loops and all, but they don't have the same 'flavor' as the MNRR bends. Other meandering methods (for example, the variable-width method I tried briefly several months ago) produce potentially more similar results: Nevertheless, by virtue of the different parameter combinations, the JP method produces a wide variety of geometries which are probaby sufficient for the 'proof of concept' of this project. 2. The MNRR hasn't changed all that much in 100 years. For examples, see figures 29 and 31 of EJ06. Similarly, the simulations performed at the nominal erosion rate of 1.2 ha/yr/km are pretty boring, and have few if any new loops or cutoffs. Only those simulations with exaggerated migration rates are of visual interest. Therefore, making use of such exaggerated simulations in this project might be called a 'dramatization' of the evolution of the river over the next 100 years. Or, what might happen if the circumstances conspired to create greater than average migration which was not properly managed by human intervention. By running many different simulations having different combinations of width, depth, flow, erosion rate, etc... it also becomes possible to delineate at what point the evolution of the river becomes dangerous to established human habitation and other property. 5.2 Appropriateness of meandering methods chosen ------------------------------- The 'simple' models I am using cannot generate the specific detail which is unique to this river (including multiple channels, 'braiding', chutes, islands, and sand bars), but by performing a large number (3000+) of simulations using slightly different input parameters, I can generate a diversity of river shapes and dynamics that, all together, approximate what I think the river might do during the next 100 years: Finally, the single-thread meandering model used in this study is clearly not adequate to accurately represent a river such as the Missouri, primarily because: 1. It reduces the geometry of the river to a single idealized channel defined by a centerline and width and having a trapezoidal cross-section, and 2. It does not track the transport of sediment within the water from where it is eroded from the banks or bottom to where it is eventually deposited downstream. * Continuing with both the JP and other 'physical' methods. I'm still not completely happy with the JP method, and I have had doubts about it since I first read the paper in January. On one hand, I'm happy that a simple method such as centripetal speed works well, but I also acknowledge the need for a more realsitic method that makes use of some physics. Perhaps there is a way to combine the two. My intuition is that the best practical method will be a simple physical model which is also second order (in either integration or differentiation), but my personal bias (perhaps given my education) is toward a mass/acceleration/force model such as Topographic Steering (which isn't such a bad model) rather than one based on bank speed. Now that I have a good and an easy way to test different models, I would like to try some others (suggestions welcomed). However, I suspect that in the end, for this project at least, it will be less the case that we attempt to find the 'one true' meandering model which is best above all others, and more the case that we try to find or modify one sufficiently to fit the actual data for the MNRR, no matter what its basis might be. In conclusion, I would predict that the JP89 method will not provide a very realistic simulation for this river. It is clear from the maps and photos that the MNRR violates both the requirements and the basic principle of the JP method (i.e. that the flow is faster on the outside radius of curves, and that material is eroded there and deposited on the inside radius). The MNRR is generally very wide, has very unsymmetric banks, a very non-linear bed depth profile, and often a shallower bed depth at the outside radius of curves, where material has accumulated rather than eroded. This may follow partly from the hypothesis that the primary flow is confined to a narrow off-center region which is not well delineated by the actual banks, especially in wide segments of the river. Ideally, I would like to develop a meandering model which can take these issues into consideration, by both (a) changing the geometries (erosion/accumulation) of each bank independently of one another, rather than as a single idealized centerline/width/flow, and (b) modeling the cross-stream depth and flow profiles as other than linear gradients so as to better represent both shallow, slow-moving flows near the outer banks and faster, deeper flows away from the bars. I have also tried to create some simulations with deliberately bigger 'loops'. As I said a few weeks ago, getting big loops using my meandering routines is nothing special; mostly it involves picking the proper distances used for curvature averaging. One downside to bigger loops seems to be that the river may then 'hug' the valley wall for protracted periods and distances, as there is little or no downstream variation in curvature to make it do anything else (see images below). Getting big loops with Johannesson-Parker, however, is still an unkown, and a question of picking the proper combination of other physical parameters. I'll be looking into that next week (including using much larger flows). Nevertheless, I could probably use my curvature averaging code with JP, and bigger loops would 'magically' appear for any set of physical parameters you wanted to use (provided they didn't crash the JP calculations). In the simulations below, I have also 'clipped' the river width to a maximum value of 400 m and below, as I think that is more representative of the actual meandering width, especially in segments where I have disregarded islands, shallows, and multiple channels (e.g. at 'Goat Island'). You and I seem to agree that the water surface and bank width are a deceptive kind of 'event horizon' which hide the fast-flowing, narrow, submerged conduits which are actually responsible for sediment exchange and river migration. A consequence of this hypothesis is that there may be a planform beneath the visible planform, and that the surface width of the river may be just the 'hull' or 'envelope' of the active submerged planform. I would still like to investigate this possibility later this year, if time permits. Simple multi-thread simulations As I have stated in previous updates, we believe that a single-thread simulation --even one with a variable width-- is not a good representation of a river like the Missouri, for the following reasons: * The MNRR is a wide river, and varies considerably in width, from a minimum of 186.79 m to a maximum of 1647.97 m, with a mean of 747.84 m, as digitized from 1957-1996 USGS topographic maps. * The 59-mile segment of the MNRR from Gavins Point Dam to Kensler's bend contains numerous islands, bifurcations, and chutes, none of which can be accomodated by a single-thread model (and which partly account for the large width). * In addition, the 'wild' river (prior to human intervention --see 1890 maps) is considerably braided in areas. * Depending on water level, the MNRR contains many submerged or surface bars, and shallow areas near the banks and away from the actively flowing channel(s). * We hypothesize that the width suggested by the banks is misleadingly large, and that the actual active, meandering, channel is much narrower, possibly on the order of from 400-800 m or narrower. * It is clear that the MNRR bed does not have an idealized trapezoidal cross-stream profile, as is expected by the Johannesson-Parker and other meandering models. * We further hypothesize that there may be several narrow, eroding, and sediment carrying 'cores' submerged within the observable banks and water surface of the MNRR, and that these multiple threads are responsible for the actual geometry and dynamics of the river. 5.3 Material properties of simulation domain ----------------------------------- 5.4 Indications for future research -------------------------------------------- Finally, some estimate is needed of any of the following equivalent information: * The bed particle size and density distribution, or * The various Shields stress values (possibly from surveys of similar rivers), or * Empirical measurements of several bed slopes wrt centerline curvature, or * Some other way to estimate the 'scour factor' used in the equations for u_1b. 5.4.1 Separability of simulation framework, meandering methods, and database --- 5.4.2 Variation of input, domain, and boundary parameters ---------------------- 5.4.3 Alternate meandering methods --------------------------------------------- * Modifying the meandering program to represent the river as a macroscopic entity and system rather than a purely 1d geometric object, incuding support for variable width and depth(s), potential depth or speed profiles, and the possibility of allowing islands in the stream or a small number of short alternate paths (e.g. 'chutes'). * Integrating the meandering with the surrounding topography. 5.4.3.1 Multiple threads ------------------------------------------------------- Finally, the single-thread meandering model used in this study is clearly not adequate to accurately represent a river such as the Missouri, primarily because: 1. It reduces the geometry of the river to a single idealized channel defined by a centerline and width and having a trapezoidal cross-section, and 2. It does not track the transport of sediment within the water from where it is eroded from the banks or bottom to where it is eventually deposited downstream. As a logical next step to the current study, I have begun work on a simple 'multiple-thread' simulation which: * Makes use of several independently meandering threads of flow, all of which are submerged within the observable banks and surface of the river, * Represents the river bed by a 3d surface showing the depth of the water at each grid point, and * Estimates sediment erosion and re-deposition by a simple exponential mechanism. Here are examples of 3 and 5 thread simulations, with color and contours showing the depth of the water at each point (red = shallow/ground, blue = deepest): Although this is an initial, simplistic attempt, it is sufficient to demonstrate how the combined, synergetic effect of several discrete threads of flow might produce the geometry and dynamics which are closer to what are seen in the real river, including: * Multiple channels, bifurcation and merging of channels, * 'Braiding' (a complex and changing interwoven surface pattern), * Islands and chutes, * More naturally sculpted 3d basins, * Deposition of sediment to form shallows and bars. Further work on multiple-thread simulations will be continued in the next project. Simple multi-thread simulations As I have stated in previous updates, we believe that a single-thread simulation --even one with a variable width-- is not a good representation of a river like the Missouri, for the following reasons: * The MNRR is a wide river, and varies considerably in width, from a minimum of 186.79 m to a maximum of 1647.97 m, with a mean of 747.84 m, as digitized from 1957-1996 USGS topographic maps. * The 59-mile segment of the MNRR from Gavins Point Dam to Kensler's bend contains numerous islands, bifurcations, and chutes, none of which can be accomodated by a single-thread model (and which partly account for the large width). * In addition, the 'wild' river (prior to human intervention --see 1890 maps) is considerably braided in areas. * Depending on water level, the MNRR contains many submerged or surface bars, and shallow areas near the banks and away from the actively flowing channel(s). * We hypothesize that the width suggested by the banks is misleadingly large, and that the actual active, meandering, channel is much narrower, possibly on the order of from 400-800 m or narrower. * It is clear that the MNRR bed does not have an idealized trapezoidal cross-stream profile, as is expected by the Johannesson-Parker and other meandering models. * We further hypothesize that there may be several narrow, eroding, and sediment carrying 'cores' submerged within the observable banks and water surface of the MNRR, and that these multiple threads are responsible for the actual geometry and dynamics of the river. To begin to explore the ideas listed above, I have created a very simple multiple thread simulation program derived from the MNRR simulation used previously. It combines the vector-based meandering approach we have been using with a cell or grid-based background which represents the depth of the river at every point (or, in this case, in every 50x50 m square): Threads can meander as before (although in the simulations which follow, only simple curvature is used to produce thread migration --see discussion below). At present, threads interact only through the common depth of the water (again, see discussion below). The migration rate of each point on a thread is based on the depth of the water there. When a thread is surrounded by water of the same depth, it can migrate at the maximum rate. When the depth to one side (based on the depth of the next 'cell' in the direction of meandering) is less than that of the thread (possibly dry land), the migration rate is attenuated by a non-linear function having the following shape: I must admit that I just made this function up; it didn't come from any well-defined theory. However, its shape is based on the following ideas: * The function should not be linear. * Migration should be easiest when the thread is surrounded by water at the same depth, and hardest (but not zero) when there is dry land to one side. * Migration should be easier when the depth to one or both sides is similar to the thread depth, and progressively harder as the depth decreases. Depths near to zero should retard migration nearly as much as dry land. In any case, it is easy to substitute another function mapping depth to migration rate, to see what the effect of different functions would be. I am also not modeling sediment transport. Migrating threads 'scour' the river bed to their depth in each cell they occupy. Depth of cells not occupied by a thread 'decay' at a constant exponential rate, typically 0.99 (i.e. depth(t + dt) = depth(t) * 0.99). This serves to simulate sedimentation without having to perform all the computational bookkeeping otherwise required. Here is an example of a single thread meandering unconstrained. The color of the background indicates depth (blue = 0, red = 2.5 m, although the thread is also drawn in blue but is always at maximum depth). In this case, the thread goes everywhere at the maximum migration rate, and the depth remains constant over the extent of its motion. Loop cutoffs are removed as they happen, but are not accumulated since they occur in water. Both endpoints of the thread are fixed. Each small square of the grid is 100 m. Depths are saved for every 50 x 50 m sub-square. The thread is 25 m wide. Actually, the distances are unimportant, since the particular meandering model used does not make use of scaled physical quantities (e.g. flow volume), but only geometry. After each motion step, the thread is re-interpolated at 50 m resolution using circular arcs. In general, the length of the thread continues to increase during the simulation, except when cutoffs are removed: Here is an example of an unconstrained thread meandering against a background whose depth decays at a rate of 0.99: Here are two different examples of a single thread whose migration rate is constrained at each point along its length by the function shown above, depending on the depth of the water in the cells toward which the thread wishes to meander cross-stream. Thus, it is easier (faster) for the thread to meander toward deeper water, and harder (slower) for the thread to meander toward shallower water, or dry land. The depth of each cell occupied by the thread is set to 2.5 m, while the depths of all unoccupied cells decrease to depth(x, y, t) * 0.99 every dt, simulating the effect of re-sedimenting the river bed. Currently no effort is made to keep track of or balance the amount of sediment removed by the thread, or replaced by the decay effect: Here is the example at left above shown as a color DEM, monochrome DEM with contour lines every 0.25 m, and as a 3d wire frame. Each DEM has been interpolated by 4x from the original 50 m grid. Note that the colors are reversed from the simulation program, so that red is ground level and blue is the deepest water. This is a more visualy appealing way to display the resulting river: A single thead with depth-restricted motion and re-sedimentation will in general produce a gently meandering channel like this. Note that there are 2 kinds of migration going on: (1) The migration of the submerged thread, which in this case is given as a constant times the local curvature of the thread, where the value of the constant changes according to the depth of the water, and (2) the overall migration of the eroded chanel, which is in general both wider and considerably less sinuous than that of the thread. The meandering of the channel in (2) is not the result of any specifically programmed meandering model, but is an indirect, aggregate, emergent process resulting from the underlying action of the thread. Here are two steps from the same simulation in which 3 threads of equal depth were initialized 'mostly' straight with slight random perturbations, 500 m apart. The final step configuration is also shown with depth contours and in 3d: The threads migrate independently of one another, and interact only through the depth background. Nevertheless, there is a tendency for threads to synergistically merge and separate due to their combined effect on the depth of the bed, and to create a braided pattern which changes over time. Note again that the overall effect of several threads is to create a wide system of braids or multiple channels and islands which itself gently meanders across the landscape. Simulations can also be performed with 5 (first series of plots) and 7 (second series of plots) independent threads having the same depth: The threads don't need to be the same depth, and we expect that multiple threads in a river might be distributed in {x, y, z} space, so that they can pass one another without touching (in the current simulations, thread pair intersections are simply ignored). In the following 7 (left and right images) and 5 (center) thread simulations, threads are created with the following depths: {0.625 m, 1.25 m, 1.875 m, 2.5 m, 1.875 m, 1.25 m, 0.625 m} and {0.833 m, 1.666 m, 2.5 m, 1.666 m, 0.833 m}. Threads maintain a constant depth, and migrate at the maximum rate if the current depth is equal to or greater than their own. This adds an additional level of synergy, as a shallow thread can 'piggyback' over a deeper thread and move at the maximum rate while the deeper thread is attenuated by scouring the bottom. The threads have also been created 250 m apart (10x their width) rather than at 500 m. Simulations using threads at different depths have, in general (a) a single discernable thalweg, and (b) a more naturally sculpted basin: Discussion The simulations shown above are the simplest I could think of that make use of some of the work I have already done (and that I could reasonably accomplish in only 2 weeks). Obvious improvements to these might include the following: * The threads need to meander by some physical mechanism other than simple curvature. Actually, there should probably be two such mechanisms: one when the thread is submerged and another when it is eroding a bank. The latter might be one of the speed-based methods such as JP89, etc... However, the mechanism for meandering within the submerged flow should probably be based on force rather than speed, in that a 'tube of water flowing within water' which creates expanding loops must: 1. Perform work as its loops increase in size (as the changing radius entails a centripetal force which acts through a distance (the changing radius)). The energy for this expansion could come from 3 sources: 1. The speed (kinetic energy) of the flow could slow, which, for an incompressible flow such as water, would require an increase in thread radius (cross-sectional area), moving against a depth pressure which would require an additional loss of speed, and cause the thread to 'die out' quite quicky. Since we are positing long-distance-stable threads, this is not a tenable mechanism, and we are 'axiomatically' requiring the thread speed to remain constant throughout its length. 2. The thread could become colder (thermal energy). Actually, I suspect that this does happen, and I would like to see some relevant observations of flow speed, temperature, and depth. 3. Primarily, the energy for thread motion comes from the gravitational (potential energy) difference between the river inlet and outlet. All of the simulations shown above were performed on a completely flat level plane, and assumed that the specified rate of water flow could be maintained by whatever means necessary. 2. Provide a change in linear momentum as a thread of water changes direction, which must be accompanied by an equivalent change in momentum of the surrounding water (i.e. it affects the flow field downstream). 3. Provide a change in angular momentum of a curved, expanding thread loop which must be accompanied by a torque at right angles to the momentum: i.e. the thread must also be helical along any length which is curved and accelerating (in direction), and changing in radius. * The threads should interact with one another. Threads could cutoff one another, or they could 'crossover' like chromosomes at their intersections (in which segments of threads are exchanged). I know of no hydrodynamic reason they might attract or repel one another like electrical currents flowing through wires, but perhaps there is a way for them to close distance or merge along short segments. Threads must also be allowed to bifurcate and join. * The threads should very their depth over time and distance, based on some mechanism which is different from that of {xy} meandering. Currently threads maintain a specific depth, without reason. Actual threads should be able to change their depth depending on bed geometry and other conditions (e.g. temperature and pressure). However, I suspect that the mechanism for varying thread depth in z is different from that of thread migration in x and y. Nevertheless, a suitable physical process should be sought. * Thread velocity should be used to extrapolate a 'flow field' over the entire depth domain, which should in turn be used to calculate sediment pickup, transport, and downflow accumulation. For example, previous and unrelated work (by Sky Coyote) has been used to calculate the flow of the atmosphere on Venus from Earth-based telescope images: A similar procedure could be used to calculate the vector flow (direction and magnitude) at every {x, y} or {x, y, z} point in the basin, both near to and away from the individual threads. With this field, additional calcultions concerning localized sediment buildup and removal, vorticity, and departure from laminar flow could be performed, or a solution to the Navier-Stokes equations could be performed. 5.4.3.2 Sediment transport ----------------------------------------------------- * Sediment eroded by threads (or provided at inlets) should be tracked and used to 'decay' the downstream depth by exponential accumulation which is potentially different at each point. This might be based on the speed of flow, as determined by a manner discussed next. Simply, flow speed above an average or threshold would pick up sediment from the bed or banks, while flow below the average speed or another threshold would deposit a fraction of whatever sediment the flow is carrying at that point. 5.4.3.3 Shape parameters ------------------------------------------------------- 1. Digitizing the 1890 maps of the MNRR. 2. Trying to create a quantitative 'metric' for comparing the geometries of the river in the 1890s and the 1990s, and for ranking and sorting the resulting simulations based on similarity to these two configurations. 3. If possible, deriving a 'vector space' of river configurations, in which the 1890 and 1990 points (and the 2304 new variations) form a line or cluster in 6d which can be extrapolated to predict other potential configurations at different times in the past or future of the river. * Digitize 1890 maps. * Find 3d space of 1890-1990 river shape parameters. * Least-squares and other similarity metrics between rivers. * Look for numerical shape parameters describing river geometry on different scales. * Attempt to compare (and sort) different river geometries according to some set of shape parameters or orthogonal basis vectors. 6.0 Conclusion ----------------------------------------------------------------- This update represents the last of the 'real work' on the MNRR project. During December I will work only on the final report and website. Updates of work on the final report will be posted at the parent site as they are available. Nevertheless, I think this is a good conclusion to the project, as it clearly indicates what should be done next. This has been both an interesting and a worthwhile project to work on, and I hope that someone makes use of this work 7.0 References -----------------------------------------------------------------