During the past 2 weeks I have:
The stand-alone program provides a convenient way to test new meandering code and observe results without having to write C++ code and compile it into Child, and then wait for an entire run to complete before viewing graphical results afterward. Parameters can be entered into the control window, and results can be viewed as the model is running. The meandering calculations are performed in a separate execution thread so that the GUI does not block. The main program provides for udating and displaying the river coordinates, and for interpolating new points as they are required, either manually or automatically.
Output to, within, and from the meandering routine can be displayed on the terminal:
time = 15, stations = 2212, stnserod = 2212
i x y xs dels flow rerody lerody slope width depth diam lambda
----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
0 -15.163632 1002.873743 0.000000 1.484104 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
1 -15.441774 1001.415935 1.484104 1.483674 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
2 -15.698652 999.954668 2.967778 1.498002 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
3 -15.930872 998.474775 4.465780 1.494417 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
4 -16.149647 996.996459 5.960196 1.649278 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
5 -16.401656 995.366548 7.609474 1.660638 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
6 -16.637211 993.722702 9.270112 1.951864 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
7 -16.867574 991.784479 11.221976 1.969085 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
8 -17.020755 989.821362 13.191061 2.181093 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
9 -17.061157 987.640643 15.372153 2.195039 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
10 -16.955219 985.448163 17.567192 2.375829 25.000000 1.000000 1.000000 0.001000 10.000000 2.500000 0.000000 0.000000
...
mrate=1.000000 mexp=0.750000 ldist=5.000000 sdist=10.000000
i s curvature curvature_e speed dspeed acc mass force
----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
0 0.000000 0.009782 0.009782 1.000000 1.048910 0.009782 10000.000000 97.819940
1 1.484104 0.009782 0.009782 1.000000 1.048910 0.009782 10000.000000 97.819940
2 2.967778 0.012319 0.010628 1.000000 1.053138 0.010628 10000.000000 106.276548
3 4.465780 0.005831 0.009428 1.000000 1.047142 0.009428 10000.000000 94.284858
4 5.960196 -0.004121 0.006719 1.000000 1.033593 0.006719 10000.000000 67.185576
5 7.609474 0.006692 0.006714 1.000000 1.033571 0.006714 10000.000000 67.141454
6 9.270112 0.013302 0.007655 1.000000 1.038277 0.007655 10000.000000 76.553165
7 11.221976 0.020619 0.009203 1.000000 1.046017 0.009203 10000.000000 92.034551
8 13.191061 0.028596 0.011820 1.000000 1.059099 0.011820 10000.000000 118.198202
9 15.372153 0.030526 0.015936 1.000000 1.079678 0.015936 10000.000000 159.356902
10 17.567192 0.036522 0.022710 1.000000 1.113548 0.022710 10000.000000 227.095144
...
i delx dely rightdepth leftdepth
----- -------------- -------------- -------------- --------------
0 -1.016661 0.201633 2.504991 2.495009
1 -1.020546 0.186371 2.504991 2.495009
2 -1.022769 0.171160 2.505422 2.494578
3 -1.020802 0.159939 2.504810 2.495190
4 -1.016932 0.153410 2.503428 2.496572
5 -1.015729 0.146291 2.503426 2.496574
6 -1.020982 0.129021 2.503906 2.496094
7 -1.030591 0.095968 2.504696 2.495304
8 -1.043707 0.045872 2.506031 2.493969
9 -1.060519 -0.019813 2.508130 2.491870
10 -1.078162 -0.101876 2.511586 2.488414
...
One interesting feature of the program is that the interpolation is performed using circular arcs, rather than straight line segments, so as to preserve the curvature of the stream in new points. Thus new points are created with a curvature similar to the surrounding points (rather than zero), and therefore the curvature does not increase erroneously at existing points. As an example, consider a stream consisting of a single circular arc. If, as the stream lengthens, points are interpolated using straight segments, the stream departs from circularity and some points along the perimeter accumulate excess curvature:
If circular arc interpolation is used, circularity (and curvature at each point) is preserved:
Arcs are computed by fitting a circle to every contiguous set of 3 points. Line segments are used only if the points are colinear. This method is also used to compute the curvature at each point for the meandering calculations (i.e. 1 / radius of fitted circle), rather than the angular change in discrete tangential vectors.
Another nice feature of the program is that the state of the stream is automatically saved whenever the simulation is paused. Any number of intermediate states can be accumulated, and it is possible to go back to previous states by clicking a button. The simulation can then be resumed from any past state (possibly using new parameters), or can be stepped forward again without having to recompute the points.
This program has been used to test both the JP meandering model and other, simpler, models. The outcome of these tests indicate that it is both possible and fairly easy to generate the characteristic meandering shape of a river (e.g. 'Kinoshita curves'), including the upstream asymmetry of their lobes, using a broad variety of methods and parameters, not just a complex physical model such as JP. All models are derived from curvature (both 'local' and 'effective' (which may be the same)). Derived models can use either force or speed (or both), or a first or second-order integral (as in JP). Nevertheless, all produce roughly the same general behavior.
There appear to be 3 general requirements for meandering, irrespective of the specific model used:
The last 2 requirements can be satisfied by having the function f() perform a spatial integration (e.g. by averaging the upstream curvature along some interval of s). (There is also a fourth requirement for stability that I will discuss below that may require a second temporal integration of the output.)
The JP implementation I have done follows the equations in Lancaster 1998 (PhD thesis) and in Lancaster and Bras 2002:
with
as in the Topographic Steering model, rather than the equations from Johannesson and Parker 1989:
although with the appropriate change of variables they amount to pretty much the same thing. Following Lancaster, the equation for sigma_s has not been used, with the 'effective' curvature being the same as the local curvature (although there is support in the code for performing a second integral if necessary). Integration is via the trapezoid method. Here is the Python code to implement this model.
The simulation shown above used a simple a * f(C(s))^b model where f() averaged the upstream curvature to yield the effective curvature. Here is an example of a JP simulation:
At present the program does not detect or remove neck or chute cutoffs.
Two problems I am having with the JP model at present are:
i s U C_f A_s K curvature integrand integral u_b1 acc
----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
0 0.000000 1.000000 0.024500 22.430862 16.066791 0.017810392 1.781039e-02 0.000000e+00 0.224002 0.017810
1 0.645995 1.000000 0.024500 22.430862 16.066791 0.017810392 1.803733e-02 1.157873e-02 0.284691 0.017810
2 1.291991 1.000000 0.024500 22.430862 16.066791 0.017259162 1.770180e-02 2.312239e-02 0.336738 0.017259
3 1.937981 1.000000 0.024500 22.430862 16.066791 0.018912779 1.964499e-02 3.518522e-02 0.417675 0.018913
4 2.583972 1.000000 0.024500 22.430862 16.066791 0.021938347 2.307805e-02 4.898457e-02 0.523097 0.021938
5 3.230195 1.000000 0.024500 22.430862 16.066791 0.017271713 1.840058e-02 6.238680e-02 0.528071 0.017272
6 3.876419 1.000000 0.024500 22.430862 16.066791 0.014512966 1.565860e-02 7.339171e-02 0.543604 0.014513
7 4.522633 1.000000 0.024500 22.430862 16.066791 0.011815059 1.291021e-02 8.262250e-02 0.549970 0.011815
8 5.168847 1.000000 0.024500 22.430862 16.066791 0.008401341 9.297083e-03 8.979783e-02 0.536402 0.008401
9 5.817849 1.000000 0.024500 22.430862 16.066791 0.007798342 8.740269e-03 9.565096e-02 0.551095 0.007798
10 6.466850 1.000000 0.024500 22.430862 16.066791 0.007240014 8.218382e-03 1.011541e-01 0.564081 0.007240
...
3273 2246.248259 1.000000 0.024500 22.430862 16.066791 0.017952528 2.369064e+17 1.044139e+20 42.226106 0.017953
3274 2247.099988 1.000000 0.024500 22.430862 16.066791 0.014091143 1.890808e+17 1.045953e+20 41.553978 0.014091
3275 2247.918975 1.000000 0.024500 22.430862 16.066791 0.015566595 2.122591e+17 1.047596e+20 40.977630 0.015567
3276 2248.737962 1.000000 0.024500 22.430862 16.066791 0.015670320 2.171310e+17 1.049355e+20 40.396884 0.015670
3277 2249.554868 1.000000 0.024500 22.430862 16.066791 0.015877484 2.235524e+17 1.051155e+20 39.828822 0.015877
3278 2250.371773 1.000000 0.024500 22.430862 16.066791 0.016157294 2.311638e+17 1.053012e+20 39.271787 0.016157
3279 2251.180293 1.000000 0.024500 22.430862 16.066791 0.016044347 2.332145e+17 1.054889e+20 38.724683 0.016044
3280 2251.988812 1.000000 0.024500 22.430862 16.066791 0.016022536 2.366176e+17 1.056789e+20 38.187020 0.016023
3281 2252.797036 1.000000 0.024500 22.430862 16.066791 0.016029808 2.405048e+17 1.058717e+20 37.658332 0.016030
3282 2253.605260 1.000000 0.024500 22.430862 16.066791 0.016029808 2.443451e+17 1.060676e+20 37.137879 0.016030
I am currently working to solve these problems. Any suggestions or comments based on your own use of this method would be appreciated.
Another problem appears to be either an inherent or a numerical instability in the simulation that causes regions of low curvature, or where the curvature changes sign, to degenerate into oscillations and chaos:
(Click for larger image)
I don't know if this phenomenon is related to the JP method itself (possibly because I am not calculating the integral for sigma_s, and have therefore turned a second-order system into a first-order system), or if it is endemic to the overall simulation. Often when an instability of this kind occurs, it means that the integration method (which computes x(t) from x(t - 1), dx, and dt) is not good enough. I am currently updating x and y just by adding delx and dely (i.e. Euler's method), and which is, I think, also the way that Child does it. If this is the culprit, then simply increasing the number of points or decreasing dt (i.e. the overall migration rate) will not fix it. Instead, I may need to lowpass filter the output (via temporal integration), or, better yet, use a multi-pass integration method (such as Runge-Kutta 4) in order to correctly represent this second-order system. I will be working on this next.