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.
In addition, the initial erosion rate can be automatically set to a given value by varying the migration rate using a binary search. Here is typical output of setting the initial rate to 15 ha/yr/km:
0: mrate = 1.000000, erate = 0.026246 1: mrate = 2.000000, erate = 0.052491 2: mrate = 4.000000, erate = 0.104982 3: mrate = 8.000000, erate = 0.209963 4: mrate = 16.000000, erate = 0.419921 5: mrate = 32.000000, erate = 0.839818 6: mrate = 64.000000, erate = 1.679544 7: mrate = 128.000000, erate = 3.358720 8: mrate = 256.000000, erate = 6.715983 9: mrate = 512.000000, erate = 13.426230 10: mrate = 1024.000000, erate = 26.830305 11: mrate = 768.000000, erate = 20.130939 12: mrate = 640.000000, erate = 16.779264 13: mrate = 576.000000, erate = 15.102919 14: mrate = 544.000000, erate = 14.264617 15: mrate = 560.000000, erate = 14.683779 16: mrate = 568.000000, erate = 14.893351 17: mrate = 572.000000, erate = 14.998136 18: mrate = 574.000000, erate = 15.050527 19: mrate = 573.000000, erate = 15.024332 20: mrate = 572.500000, erate = 15.011234 21: mrate = 572.250000, erate = 15.004685 22: mrate = 572.125000, erate = 15.001410 23: mrate = 572.062500, erate = 14.999773 24: mrate = 572.093750, erate = 15.000592 25: mrate = 572.078125, erate = 15.000182
In this case, the migration rate was changed from 1.0 to a final value of 572.078125 after 25 steps. This is quick to accomplish, as only the meandering and erosion calculations need be performed each step, without boundary or cutoff checks, or reinterpolation.
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:
# MNRRSim command file # load base model load base.mnrr # change settings runSteps = 100 method = Johannesson-Parker showPoints = True showCenter = False displayEvery = 10 checkCutoffsEvery = 10 # 0 ldist = 0.25 sdist = 0.25 initializeErosion 15.0 run saveImage test000.png save test000.mnrr # 1 toBeginning ldist = 0.25 sdist = 0.50 initializeErosion 15.0 run saveImage test001.png save test001.mnrr # 2 toBeginning ldist = 0.25 sdist = 0.75 initializeErosion 15.0 run saveImage test002.png save test002.mnrr ... # 33 toBeginning ldist = 1.50 sdist = 1.00 initializeErosion 15.0 run saveImage test033.png save test033.mnrr # 34 toBeginning ldist = 1.50 sdist = 1.25 initializeErosion 15.0 run saveImage test034.png save test034.mnrr # 35 toBeginning ldist = 1.50 sdist = 1.50 initializeErosion 15.0 run saveImage test035.png save test035.mnrr quit
This generated short (100 step) simulations of all combinations of ldist and sdist from 0.25 to 1.50 at increments of 0.25 river widths (limited here to 500 m). Rather than show all results, I will present a few samples:
ldist = 0.50
sdist = 0.25
ldist = 0.75
sdist = 0.75
ldist = 1.00
sdist = 1.50
ldist = 1.25
sdist = 1.50
These results demonstrate that it is possible to get large features (e.g. 'loops') without using long integration distances (I was previously using ldist = 0.5 and sdist = {2.0 - 8.0}), which are physically reasonable in that we expect the upstream effect of geometry on meandering to be confined to about 1 river width.
Of these 36 results, 11 were selected for full 1000 step runs. Here are some samples of the new results:
ldist = 1.00
sdist = 1.00
ldist = 1.00
sdist = 1.50
ldist = 1.25
sdist = 1.25
ldist = 1.50
sdist = 1.75
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:
wmax = 800.00
flow = 1200.00
wmax = 600.00
flow = 1200.00
wmax = 400.00
flow = 1200.00
wmax = 400.00
flow = 1600.00
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:
depth = 3.50
diam = 0.0010
depth = 3.00
diam = 0.0010
depth = 3.00
diam = 0.0015
depth = 2.50
diam = 0.0015
Both increasing grain size and decreasing depth appear to have the effect of making the river more 'active', although decreasing depth also magnifies the differences in activity between smaller and larger grain sizes.
An additional change I made next was to 'soft clip' the river width to a maximum value, rather than simply truncating all widths greater than the cutoff value. This has the effect of maintaining variation in the river width in its widest parts, rather than replacing the widest parts (a significant fraction of the overall river) by a smooth channel of constant width. I don't know if this will have any significant effect on the simulations, but it seemed like a 'good idea'. Here is a plot of the 'soft clipping' falloff for a maximum width of 800 m (about the widest that the JP method can handle):
Here are some 800 m results based on the previous depth x diameter study:
ldist = 1.00
sdist = 1.00
wmax = 800.00
flow = 1200.00
depth = 3.25
diam = 0.001
ldist = 0.75
sdist = 0.75
ldist = 0.875
sdist = 0.875
flow = 1100.00
depth = 3.125
diam = 0.00125
ldist = 1.5
sdist = 1.5
flow = 1200.00
depth = 3.125
diam = 0.00125
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:
ldist = 1.00
sdist = 1.00
wmax = 600.00
flow = 1200.00
depth = 2.50
diam = 0.0010
10 years
30 years
50 years
70 years
80 years
90 years
100 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.
Nevertheless, it now appears clear that the end is in sight: