MNRR JP and variable-width tests

None of the Johannesson-Parker implementations will work using the MNRR input data as-is:

width = variable (min=218.901 max=1480.894 mean=803.771 in 200 pt 19877.41 m subset)
depth = 3.429
flow = 899.33
slope = 0.0001918
diam = 0.001

All produce fatal numerical errors in the calculations. Here is the output for implementation 3 (JP89 integral eqns with partial summation, probably the best):

    i              U            C_f          chi_1            chi         chi_20          delta        C_tilde            A_s              A            phi          sigma        sigma_s     integrand1      integral1     integrand2      integral2           u_1b
----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
    0       0.180396       0.198057       0.173020      -0.160314       2.818374      -3.474884       0.002655      -0.000237       3.147522       0.000000       1.930374       1.930374   1.930374e+00   0.000000e+00   1.930374e+00   0.000000e+00       1.930374
    1       0.177104       0.205488       0.169862      -0.163471       2.878720      -3.184809       0.002655      -0.000244       3.147522      99.929554       1.966254 -3160604853.992569   3.127006e+05   1.562411e+07  -5.026425e+14  -2.511442e+16 -29784391179.655685
    2       0.177628       0.204277       0.170365      -0.162968       2.868916      -3.229242       0.002303      -0.000243       3.147522      99.608831       1.699916 -663616290109375872.000000   2.424762e+05   1.207648e+07  -9.465832e+22  -4.714402e+24 -6196875914419073024.000000
    3       0.181984       0.194615       0.174543      -0.158790       2.790280      -3.624717       0.000035      -0.000233       3.147522      99.938426       0.025294 -563471415042582015060738048.000000   2.136592e+03   1.068488e+05  -4.759665e+31  -2.378367e+33 -5029452071087756093497016320.000000
    4       0.191777       0.175247       0.183935      -0.149398       2.630410      -4.721768       0.001716      -0.000203       3.147522      99.974190       1.173588 -16920977129793365379595106484762968064.000000   3.217321e+04   1.608247e+06  -4.638786e+41  -2.318795e+43 -136052681449621237668582504876839272448.000000
    5       0.208259       0.148605       0.199744      -0.133589       2.405373      -7.606553       0.001073      -0.000121       3.147522      99.886018       0.675837   -3.37807e+51   3.888858e+03   1.942799e+05  -1.943790e+55  -9.707871e+56   -2.30124e+52
    6       0.225136       0.127160       0.215931      -0.117402       2.219601     -13.601221       0.000937       0.000009       3.147522      99.921461       0.545722   -2.60977e+73   9.026013e+02   4.512838e+04  -4.316452e+76  -2.156531e+78    -1.5219e+74
    7       0.238514       0.113296       0.228762      -0.104571       2.097129     -25.382153       0.000749       0.000150       3.147522      99.990238       0.411920  -6.83215e+109   3.050599e+02   1.527879e+04 -5.059751e+112 -2.529629e+114  -3.55238e+110
    8       0.250164       0.102990       0.239935      -0.093398       2.004812     -63.146341       0.000313       0.000303       3.147522      99.982178       0.164000  -1.54118e+192   6.655432e+01   3.347715e+03 -6.254375e+194 -3.126630e+196  -7.28419e+192
Exception in thread stepThread:
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/threading.py", line 486, in __bootstrap_inner
    self.run()
  File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/threading.py", line 446, in run
    self.__target(*self.__args, **self.__kwargs)
  File "/Users/sky/Projects/Lancaster/MNRRSim/MNRRSimData.py", line 629, in stepCalcs
    params)
  File "/Users/sky/Projects/Lancaster/MNRRSim/MeanderJP3.py", line 225, in meander
    edpr = exp(delta * phi / r)
FloatingPointError: overflow encountered in exp

Here is the full output for implementations 1 (Lancaster 98), 3, and 4 (JP89 differential eqns). These all fail because:

1. While u_1b starts at a reasonable value of 1.9, it then immediately becomes unbounded, because
2. The integrals for sigma_s and sigma are unbounded, because
3. A_s and delta are negative (and A_s is way too small in magnitude due to the wide river (it should be comparable to A)), because
4. chi is negative (also note that chi_20 is way too high (it should be about 1.0)), because
5. C_f is way too big (it should be about 1e-3), because
6. Either U (flow) is too low or H (depth) and/or I (slope) are too large.

It's probably more complicated than this, as there is likely a thin n-dimensional subspace of acceptable values whose relationships cannot be easily partitioned by simple ranges. I have noticed this before: the JP method only accepts inputs in a fairly narrow range of values. In general, it likes very flat, shallow, streams with fairly high flows. In order to get something working with the MNRR geometry, I have tried starting from a different set of initial conditions:

width = 100.0
depth = 1.0
flow = 1000.0
slope = 0.0001
diam = 0.001

Here is the JP1 simulation:

Here is the JP3 simulation:

Here is the JP4 simulation:

In the interest of time, I have decided to continue to work only with the JP3 implementation, as I think that stands the best chance of producing something useful. The initial output of the simulation shown above (at middle) was the following:

   i              U            C_f          chi_1            chi         chi_20          delta        C_tilde            A_s              A            phi          sigma        sigma_s     integrand1      integral1     integrand2      integral2           u_1b
----- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- --------------
    0      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.002655       3.523061       6.499205       0.000000       0.132775       0.132775   1.327748e-01   0.000000e+00   1.327748e-01   0.000000e+00       0.132775
    1      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.002655       3.523061       6.499205      99.929554       0.132775       0.133013   1.330351e-01   1.328114e+01   1.332733e-01   1.329304e+01       0.135014
    2      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.002303       3.523061       6.499205      99.608831       0.115130       0.130623   1.153545e-01   1.235794e+01   1.308787e-01   1.314295e+01       0.154782
    3      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.000035       3.523061       6.499205      99.938426       0.001755       0.110377   1.758526e-03   5.840805e+00   1.105936e-01   1.205342e+01       0.269462
    4      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.001716       3.523061       6.499205      99.974190       0.085814       0.095156   8.598263e-02   4.385754e+00   9.534263e-02   1.028334e+01       0.186351
    5      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.001073       3.523061       6.499205      99.886018       0.053665       0.088262   5.377057e-02   6.971289e+00   8.843460e-02   9.169066e+00       0.219634
    6      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.000937       3.523061       6.499205      99.921461       0.046845       0.078613   4.693703e-02   5.026171e+00   7.876683e-02   8.344864e+00       0.227289
    7      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.000749       3.523061       6.499205      99.990238       0.037461       0.069288   3.753424e-02   4.218560e+00   6.942360e-02   7.401094e+00       0.237326
    8      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.000313       3.523061       6.499205      99.982178       0.015643       0.058118   1.567359e-02   2.656244e+00   5.823152e-02   6.374823e+00       0.259524
    9      10.000000       0.000010      24.596748      24.263414       1.000110     297.786749       0.000709       3.523061       6.499205      99.999972       0.035461       0.050287   3.553068e-02   2.558679e+00   5.038547e-02   5.425148e+00       0.239977
...

This is much more reasonable, as C_f, chi_20, A_s, and u_b1 are all within their proper magnitude ranges. By experimenting with the parameters for width, depth, flow, slope, and diameter, while attempting to bring them closer to their values for the MNRR, I have produced some reasonable looking simulations:

width = 100.0, depth = 2.5, flow = 100.0, slope = 0.0001, diam = 0.001

width = 200.0, depth = 2.5, flow = 150.0, slope = 0.0001, diam = 0.001

width = 300.0, depth = 2.5, flow = 300.0, slope = 0.00019, diam = 0.001

width = 300.0, depth = 3.0, flow = 600.0, slope = 0.00019, diam = 0.001

width = 300.0, depth = 3.0, flow = 800.0, slope = 0.00019, diam = 0.001

width = 600.0, depth = 3.0, flow = 800.0, slope = 0.00019, diam = 0.001

width = 800.0, depth = 3.0, flow = 800.0, slope = 0.00019, diam = 0.001

The wider simulations are somewhat less interesting, in part because the cutoff detection routine is removing new bends a bit sooner than it should, as new points in the bends are closer together than a river width. However, these simulations suggest that it might be possible to get close to the input values for the MNRR and to produce realistic output from the JP3 method.

I have also tried a variable-width meandering method, although at present this is implmented only with the centripetal speed method (which is equivalent to the unmodified sigma calculation from JP). The basis of this method is straightforward, and is derived from the class of 'bank speed' methods we have been using. The basic idea is the following:

• In the various speed methods, the outer bank speed is calculated as a proportional excess of the centerline speed.
• The erosion at the outer bank is based on this speed excess.
• The sediment removed from the outer bank is cycled via secondary flow to the inner bank (where there is a speed deficit), where it is deposited and accumulates.
• The ongoing transfer of sediment from the outer bank to the inner bank maintains a constant width to the river as it migrates.

However:

• This might be true if there was no primary flow.
• Once the outer bank sediment is in the water, the primary flow carries it downstream.
• The secondary helical flow may continue to move the sediment toward the inside bank, but it will eventually reach a point on that bank which is downsteam of the point where it left the outer bank.
• Therefore, the outer and inner banks of a bend move independently of one another. The outer bank migrates outward via erosion, as described by the speed methods. However the inner bank moves only as a response of the accumulation of sediment already in the water from upstream. Thus, the outer and inner banks of a bend may move at different rates (but always in the same direction, in this model).

The variable-width routine extends the centripetal speed method by:

1. The effective curvature is calculated as before, as the average of the (signed) curvature around a 'lag distance' (which can be positive or negative), summed over a 'sum distance'.
2. A speed excess (sigma or u_1b) and lateral motion (along the left perpendicular vector) is calculated for the outer banks along the river.
3. The inner bank motion is then calculated as the average of the absolute value of the outer bank motion, again averaged over the same 'lag' and 'sum' distances upstream (or downstream if lag +/- sum > s).
4. The centerline is then moved at the average rate of the outer and inner banks at each point, and the width changes by their difference.

As a simplification, both the curvature and sediment distances are the same. Actually, sediment entering the water at a point s should be available to all downstream points > s --not just those at the lag distance--, and should be removed and accumulated by an exponential decay process. However, the current implementation is sufficient to test the concept. There are also other possible modifications that I will discuss below.

Here are some examples of variable-width tests. Note that the initial condition may involve a constant width along the river (which will change), or a variable width. In all cases, the minimum width of any evolved segment was limited to at least 25 m:

width = 100.0, depth = 3.429, flow = 899.33, ldist = 150.0, sdist = 300.0

width = 100.0, depth = 3.429, flow = 899.33, ldist = 110.0, sdist = 550.0

width = 100.0, depth = 3.429, flow = 909.33, ldist = 110.0, sdist = 550.0

width = actual river, depth = 3.429, flow = 899.33, ldist = 110.0, sdist = 550.0

In the following progression of simulations, two additional modifications were added:

1. The lag and sum distances were based on the current river width at the calculated point (i.e. as proportions: 1.0 * width, 1.5 * width, etc...).
2. The outer bank proportional speed (e.g. 1.15, etc...) was raised to the power 1.5 before calculating lateral motion. Since this factor was always ≥ 1, its magnitude was therefore increased non-linearly wrt curvature. (Any exponent can be used as an input parameter.)

width = 400.0, depth = 3.0, flow = 800.0, ldist = 220.0, sdist = 440.0, exp = 1.0

width = actual river, depth = 3.429, flow = 899.33, ldist = 0.5w, sdist = 1.0w, exp = 1.0

width = actual river, depth = 3.429, flow = 899.33, ldist = 0.49w, sdist = 1.0w, exp = 1.5

width = actual river, depth = 3.429, flow = 899.33, ldist = 0.50w, sdist = 1.0w, exp = 1.5

I am particularly interested in the latter simulations because:

• They work well at larger widths.
• They create both upstream and downstream loops.
• All features slowly migrate downstream, similar to those in the MNRR.

Next week I will continue working on these methods in the attempt to:

1. Get them working using the actual MNRR input parameters.
2. Produce migration behavior that is in general similar to the MNRR, as seen in the maps and photos.
3. Produce loops and other features having similar amplitudes, wavelengths, and downstream distributions as are seen in the MNRR.

I will also spend at least one day working on a braided meandering test.