Interaction of fluctuations on different scales

Examples of two scales of fluctuations superimposed are calculated.
Revision of stochasflow8, where now the steps have equal lengths, and the step height parameters are modified to ensure equal traversal times. The equal step lengths makes it possible to make the two scales completely commensurate. At the same time, some changes in notation and definitions as set out in section 21 are implemented here.

In particular, 
i) We now have steps on a uniform grid with spacing L. This gives the length of a single fluctuation as 3L (rather than 3L/2, with L defined as in stochasflow 7-9).
ii) Instead of the time parameter τ0 used previously we now introduce the scale length Λ = -V1 t0 right from the start and hence have the fluctuation length scale parameter λ = L/Λ.
iii) Θ as used here is a dimensionless version of the step time θ used elsewhere.

Stir[x_] := (((x - 1))^(x - 1)) 
        Exp[1 - x] (Sqrt[2  π (x - 1)])

Single fluctuation formulas

Implementing equations 21.1-21.12

A single 3-step fluctuation is characterised by 2 parameters: its length scale λ and step size of first step Δ. Stochasflow8.nb shows that sizes of the other steps for the equal step length case are determined by the average speed condition as :

Δ1[Δ_] = {Δ, -((2 Δ)/(1 + Δ)), (Δ/(1 + 2 Δ))}

The dimensionless times for the 3 steps are:

Θ1[Δ_] = {1, (2/(1 + Δ)), 
        3}

The formulas below also require the starting and ending times, where there are no actual steps. That will not be necessary for the case of multiple steps below, where use is made of a cyclical definition; but when we get to the combined scales case, it will become necessary not only to add those two times but also to extend definitions for several more different cases that arise. To keep everything consistent, it seems best to introduce an extended times list as a function as follows:

Θ1ex[Δ_, i_] := 
      Which[i == 0, 0, 
        i ≤ 3, Θ1[Δ][[i]], 
        i == 4, 4]

The α-parameters for the i-th step is given by the formula

α1[λ_, Δ_, i_] := 1

α1[λ_, Δ_, 
        i_] := ((λ Θ1ex[Δ, i] + 1)/(λ Θ1ex[Δ, i - 1] + 1))

The step enhancement factor is

a1[λ_, Δ_, i_] := 
      Together[α1[λ, Δ, 
            i + 1] (α1[λ, Δ, 
                    i] + Δ1[Δ][[i]])^2/(α1[λ, Δ, 
                    i + 1] (α1[λ, Δ, 
                          i] + Δ1[Δ][[i]])^2 - 4  α1[λ, Δ, 
                        i]^2 Δ1[Δ][[i]])]

and so for the fluctuation as a whole

A1[λ_, Δ_] := 
      Simplify[a1[λ, Δ, 1]*
          a1[λ, Δ, 2]*
          a1[λ, Δ, 3]]

A1[λ, Δ]

Check consistency with previous formulas in stochasflow8.nb using τ0. Note the relation between them is τ0 = -1/(2λ) :

Simplify[
      A1[λ, Δ] == 
        Aeq[Δ, -1/(2  λ)]]

A1[5, 0.15]

Plot[{A1[5, Δ]}, {Δ, 0, 0.3}, 
      PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

That is the same curve as obtained previously with τ0. The value λ = 5 is not physically very relevant, but used here and elsewhere in the notebook for comparison with the τ0 = -0.1 value used in the previous notebooks. 
More relevant values of λ and Δ :

Plot[{A1[0.001, Δ] - 1}, {Δ, -0.2, 
        0.2}, PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

The sharp rise for negative Δ (i.e., reverse fluctuation) may be unrealistic as it is connected with a divergence in the rational form for A1 at Δ near -0.3 . But the main thing is that A1 remains >1 for both "orientations " of the fluctuation. Further down it is shown that the divergence is connected with the fact that if Δ = -1/3, the speed V3 that is needed in the last part of the fluctuation becomes infinite. See equation 21.2 for V3. 
 Basically, the thing is that if you slow down first, then if this slowdown is too large it becomes impossible to catch up with the constant flow over a fixed length interval. For a fluctuation that speeds up first, this problem does not exist; one can always slow down enough to let the constant flow catch up. This is the reason why results are not symmetric with the sign of Δ. The asymmetry is more severe in the present case where there is a fixed step length than in the previous one where one compensated for a large slowdown by only allowing this to happen over a short length. Anyway, as long as one stays in a reasonable Δ range, say |Δ| ≤ 0.1 , numeric problems are avoided although some asymmetry remains.

Also look at how A1 changes if the length scale changes:

Plot[{A1[λ, 0.05] - 1}, {λ, 0.001, 5}, 
      PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

So for very short fluctuations, A1 -> 1; then it increases with the fluctuation length and reaches a maximum at about λ = 2 and decreases very slowly beyond that. The location of the maximum shifts to longer lengths if the step height increases.

Plot3D[
      A1[λ, Δ], {λ, 0, 
        1}, {Δ, -0.2, 0.2}, 
      ViewPoint -> {-1.226,  2.999,  0.975}, 
      PlotRange -> {1, Automatic}, 
      ColorOutput -> GrayLevel]

Options[Plot3D]

Multiple identical fluctuations

Implementing equations 21.13-21.17. Because of the cyclic use of the times, direct substitution can be done to give all α's as one of 3 formulas depending on i (eq 21.15):

αm[λ_, Δ_, m_, i_] := 
      Which[i == 
          1, ((3  λ (m - 1) + λ + 1)/(3  λ (m - 1) + 1)), 
        i == 
          2, ((3  λ (m - 1) + ((2  λ)/(1 + Δ)) + 1)/(3  λ (m - 1) + λ + 1)), 
        i == 
          3, ((3  λ (m - 1) + 3  λ + 1)/(3  λ (m - 1) + ((2  λ)/(1 + Δ)) + 1))]

am[λ_, Δ_, m_, i_] := 
      Together[Which[
          i == 
            1, αm[λ, Δ, m, 
              2] (αm[λ, Δ, m, 
                      1] + Δ1[Δ][[1]])^2/(αm[λ, Δ, m, 
                      2] (αm[λ, Δ, m, 
                            1] + Δ1[Δ][[1]])^2 - 4  αm[λ, Δ, m, 
                          1]^2 Δ1[Δ][[1]]), i == 
            2, αm[λ, Δ, m, 
              3] (αm[λ, Δ, m, 
                      2] + Δ1[Δ][[2]])^2/(αm[λ, Δ, m, 
                      3] (αm[λ, Δ, m, 
                            2] + Δ1[Δ][[2]])^2 - 4  αm[λ, Δ, m, 
                          2]^2 Δ1[Δ][[2]]), i == 
            3, αm[λ, Δ, m + 1, 
              1] (αm[λ, Δ, m, 
                      3] + Δ1[Δ][[3]])^2/(αm[λ, Δ, m + 1, 
                      1] (αm[λ, Δ, m, 
                            3] + Δ1[Δ][[3]])^2 - 4  αm[λ, Δ, m, 
                          3]^2 Δ1[Δ][[3]])]]

Am[λ_, Δ_, m_] := 
      am[λ, Δ, m, 1]*
        am[λ, Δ, m, 2]*
        am[λ, Δ, m, 3]

Am[λ, Δ, 1]

Consistency check:

Simplify[
      A1[λ, Δ] == 
        Am[λ, Δ, 1]]

Table[Am[5, 0.15, n], {n, 1, 9}]

Values look quite OK - compare with unequal length table called efacs based on Asm1 in previous notebook.

Combining two fluctuation scales

The small scale fluctuation is still characterised by the grid length λ and now its step height is called Δs. The big scale fluctuation has a length (3N+1)λ where N is a new parameter, and a height Δb. We are working with one repeat unit here, the big fluctuation, within which there is a variable number of small fluctuations. Implement the collection of equations 21.18 as multiple case functions, similar to that used for the single fluctuation above.

Δc[N_, Δs_, Δb_, i_] := 
      Which[i == 0, 0, 
        i ≤ 
          3  N, Δ1[Δs][[Mod[i, 3, 1]]], 
        i == 3  N + 1, Δb, 
        i ≤ 
          6  N + 1, Δ1[Δs][[Mod[
              i - 3  N - 1, 3, 1]]], 
        i == 
          6  N + 2, ((-2 Δb)/(1 + Δb)), 
        i ≤ 
          9  N + 2, Δ1[Δs][[Mod[
              i - 6  N - 2, 3, 1]]], 
        i == 9  N + 3, Δ(b/(1 + 2  Δb)), 
        i == 9  N + 4, 0]

Table[Δc[2, Δs, Δb, i], {i, 
        0, 22}]

Right, now the same thing for the step times.

Θc[N_, Δs_, Δb_, i_] := 
      Which[i == 0, 0, i ≤ 3  N, 
        3  Quotient[i - 1, 
              3] + Θ1[Δs][[Mod[i, 3, 
                1]]], i == 3  N + 1, 3  N + 1, 
        i ≤ 6  N + 1, 
        3  N + 1 + ((1 - Δb)/(1 + Δb)) (3  Quotient[i - 3  N - 2, 
                    3] + Θ1[Δs][[Mod[
                      i - 3  N - 1, 3, 1]]]), 
        i == 6  N + 2, ((2 (3  N + 1))/(1 + Δb)), 
        i ≤ 
          9  N + 2, ((2 (3  N + 1))/(1 + Δb)) + ((1 + 
                  3  Δb)/(1 + Δb)) (3 
                   Quotient[i - 6  N - 3, 
                    3] + Θ1[Δs][[Mod[
                      i - 6  N - 2, 3, 1]]]), i == 9  N + 3, 
        9  N + 3, i == 9  N + 4, 9  N + 4]

Table[Θc[2, Δs, Δb, i], {i, 
        0, 22}]

ListPlot[Table[{i, Θc[2, 0.1, 0.1, i]}, {i, 0, 22}], 
      Prolog -> PointSize[0.02], 
      PlotRange -> {{0, Automatic}, {0, Automatic}}]

Next adapt the definitions of α's and step factors a for the combined fluctuation to the definitions above

αc[N_, λ_, Δs_, Δb_, 
        i_] := ((λ Θc[N, Δs, Δb, i] + 1)/(λ Θc[N, Δs, Δb, i - 1] + 1))

ac[N_, λ_, Δs_, Δb_, i_] := 
      Together[αc[N, λ, Δs, Δb, 
            i + 1] (αc[
                    N, λ, Δs, Δb, 
                    i] + Δc[
                    N, Δs, Δb, 
                    i])^2/(αc[
                    N, λ, Δs, Δb, 
                    i + 1] (αc[
                          N, λ, Δs, Δb, 
                          i] + Δc[
                          N, Δs, Δb, i])^2 - 
                4  αc[
                        N, λ, Δs, Δb, 
                        i]^2 Δc[
                    N, Δs, Δb, i])]

And now finally we introduce an index n for each small fluctuation and find the overall enhancement factor for it from eq 21.19:

As[N_, λ_, Δs_, Δb_, n_] := 
      Which[n ≤ 
          N, (ac[N, λ, Δs, Δb, 
              3 (n - 1) + 1]*
            ac[N, λ, Δs, Δb, 
              3 (n - 1) + 2]*
            ac[N, λ, Δs, Δb, 
              3 (n - 1) + 3]), 
        n ≤ 
          2  N, (ac[N, λ, Δs, Δb, 
              3 (n - 1) + 2]*
            ac[N, λ, Δs, Δb, 
              3 (n - 1) + 3]*
            ac[N, λ, Δs, Δb, 
              3 (n - 1) + 4]), 
        n ≤ 
          3  N, (ac[N, λ, Δs, Δb, 
              3 (n - 1) + 3]*
            ac[N, λ, Δs, Δb, 
              3 (n - 1) + 4]*
            ac[N, λ, Δs, Δb, 
              3 (n - 1) + 5])]

ListPlot[Table[{n, As[2, 5, 0.1, 0.1, n]}, {n, 1, 6}], 
      Prolog -> PointSize[0.02], 
      PlotRange -> {{0, Automatic}, {1, Automatic}}]

The factor associated with the steps making the big fluctuation is

Ab1[
        N_, λ_, Δs_, Δb_] := (ac[
            N, λ, Δs, Δb, 3  N + 1]*
          ac[N, λ, Δs, Δb, 6  N + 2]*
          ac[N, λ, Δs, Δb, 
            9  N + 3])

The combined enhancement is

Acom[N_, λ_, Δs_, Δb_] := 
      Ab1[N, λ, Δs, Δb]*(∏_(n = 1)^(3  N))As[N, λ, Δs, Δb, 
            n]

Table[As[2, 5, 0.1, 0.1, n], {n, 1, 6}]

Ab1[1, 5, 0.1, 0.1]

Reference cases

If there are only small length fluctuations (case B of page 172) i.e. 3N+1 small fluctuations of length 3 each gives

Asmall[
        N_, λ_, Δs_] := (∏_(n = 1)^(3  N + 1))Am[λ, Δs, n]

Table[Am[5, 0.1, n], {n, 1, 7}]

On the other hand, if there is only one big fluctuation of length (3N+1)*3 we adapt the single fluctuation formulas taking into account that we still want to evaluate the factor at the final time (9N+4) rather than a time associated with a position (3N+1) ahead of the last step as would have been done previously:

Θ1big[N_, Δ_, i_] := 
      Which[i == 0, 0, 
        i ≤ 
          3, Θ1[Δ][[i]]*(3  N + 1), 
        i == 4, 9  N + 4]

α1big[N_, λ_, Δ_, 
        i_] := ((λ Θ1big[N, Δ, i] + 1)/(λ Θ1big[N, Δ, i - 1] + 1))

a1big[N_, λ_, Δ_, i_] := 
      Together[α1big[N, λ, Δ, 
            i + 1] (α1big[N, λ, Δ, 
                    i] + Δ1[Δ][[i]])^2/(α1big[N, λ, Δ, 
                    i + 1] (α1big[N, λ, Δ, 
                          i] + Δ1[Δ][[i]])^2 - 4  α1big[N, λ, Δ, 
                        i]^2 Δ1[Δ][[i]])]

A1big[N_, λ_, Δ_] := 
      a1big[N, λ, Δ, 1]*
        a1big[N, λ, Δ, 2]*
        a1big[N, λ, Δ, 3]

Comparisons

See some values:

A1big[2, 0.1, 0.05]

Asmall[2, 0.1, 0.05]

Acom[2, 0.1, 0.05, 0.05]

Table[A1big[N, 0.01, 0.05], {N, 1, 6}]

We need to investigate dependencies on 4 different arguments. So try to select reasonable values for some then let the rest vary. First, we start with λ = 0.001 (values as small as (10^(-5)) were suggested by Peclet number arguments, but this a physical scale of about 0.1 m, which looks reasonable). For Δs and Δb take 0.05 as representative. Let us see what effect the scale ratio parameter N has:

ListPlot[
      Table[Acom[N, 0.001, 0.05, 0.05] - Asmall[N, 0.001, 0.05], {N, 1, 10}], 
      Prolog -> PointSize[0.02], 
      PlotRange -> {{0, Automatic}, {0, Automatic}}]

ListPlot[
      Table[Acom[N, 0.001, 0.05, 0.05] - A1big[N, 0.001, 0.05], {N, 1, 10}], 
      Prolog -> PointSize[0.02], AxesOrigin -> {0, 0}]

Bearing in mind that A1big > Asmall (see below) these suggest that Acom is somewhere in between, and the difference from either increases with the value of N. Take N = 5, i.e. the big fluctuation is 16 times longer than the small one, as a typical case.

Next, first fix value for the small fluctuation step size Δs and let the big one chnge either way by a factor of 2. Because all values are close to 1, just plot the fractional part::

Plot[
      Evaluate[{Acom[5, λ, Δs, r*Δs] - 1, 
            Asmall[5, λ, Δs] - 1, 
            A1big[5, λ, r*Δs] - 
              1} /. {λ -> 0.001, Δs -> 
              0.05}], {r, 0.5, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Here red is the combined curve, green is small fluctuations only, and turquioise one large fluctuation only.

We see that the single big fluctuation gives a larger factor for equal step sizes (r = 1) as expected from the previously established increase as function of λ of a single fluctuation. Also, the value rises approx linearly as the size of the big step increases, again plausible enough. 

The surprising thing is that the combined steps closely follows the value for the uniform scale small fluctuations; i.e. the presence of the long scale fluctuation is almost negligible unless the step lengths are close (see N=1 below) and the step size difference is considerable:

Plot[
      Evaluate[{Acom[1, λ, Δs, r*Δs] - 1, 
            Asmall[1, λ, Δs] - 1, 
            A1big[1, λ, r*Δs] - 
              1} /. {λ -> 0.001, Δs -> 
              0.05}], {r, 0.5, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

This behaviour is not sensitive to the value of λ unless it becomes of order 1 or larger:

Plot[
      Evaluate[{Acom[5, λ, Δs, r*Δs] - 1, 
            Asmall[5, λ, Δs] - 1, 
            A1big[5, λ, r*Δs] - 
              1} /. {λ -> 0.1, Δs -> 
              0.05}], {r, 0.5, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Plot[
      Evaluate[{Acom[5, λ, Δs, r*Δs] - 1, 
            Asmall[5, λ, Δs] - 1, 
            A1big[5, λ, r*Δs] - 
              1} /. {λ -> 1, Δs -> 
              0.05}], {r, 0.5, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Try the reverse, where the long fluctuation is kept fixed, and the small one varies :

Plot[
      Evaluate[{Acom[5, λ, r*Δb, Δb] - 1, 
            Asmall[5, λ, r*Δb] - 1, 
            A1big[5, λ, Δb] - 1} /. {λ -> 
              0.001, Δb -> 0.05}], {r, 0.5, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

The conclusion is really the same and we conclude that when there are fluctuations on different scales present, the long ones lose their effectiveness on dispersion and it should be sufficient to restrict consideration only to the shortest length fluctuations present. Note that this holds in particular when the step sizes of the two scales of fluctuation are the same. If the short fluctuations have larger step sizes than the long ones, their effect is slightly reduced by the presence of the long ones .

Having settled this, we now return to the single scale case and revise the algebraic calculation done previously for the unequal step lengths for the case of fluctuations on a fixed grid.

Algebraic analysis of repeated fluctuations on grid

What we have so far is

Am[λ, Δ, m]

LeafCount[%]

Remove[LogLinearPlot]

<< Graphics`Graphics`

Names["Graphics`Graphics`*"]

First, we reproduce the behaviour found previously in stochasflow7.nb for unequal step lengths, for the grid model and new parameter choice:

LogLinearPlot[{Am[0.001, 0.1, m] - 
          1, (10^10) (m - 1/(2*0.001))}, {m, 1, 10000}, 
      PlotRange -> {0, 0.0002}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

So, once more, there is for low m values a slow decrease (linear, but relative to log m, so more like a slow exponential) of the enhancement factor, but at around m = 1/(2λ) a sharper decline takes place. For this to be visible, of course λ has to be quite small. Anyway, the behaviour is pretty much as found previously, so follow similar steps to find a simplified approximation to the Am formula.

Asermi = 
      Simplify[Series[
          Am[λ, Δ, m] - 1, {m, ∞, 
            2}]]

Simplify[
      Series[Am[λ, Δ, m] - 1, {λ, 0, 
          2}]]

As in the previous case, the two series have common factors in their first order terms which depend only on Δ. Extract that factor:

QQ[Δ_] = SeriesCoefficient[Asermi, 1]

The reason for dividing by two is for consistence with the Q expression used in stochasflow7.nb, the connction will be further explored below.

As before, we also make a double expansion about the "midpoint" of the m-range:
(PS: from later expression, it appears that m = 1/(3λ) is a better midpoint; but that does not give the simple expression as this choice does)

Simplify[
      Series[Am[λ, Δ, m] - 1, {m, 1/(2  λ), 
          2}, {λ, 0, 2}]]

Lo and behold, once more the leading term contains Q, now with a coefficient 6/5. So as before we pull a rational expression from thin air that has a similar relationship between its various expansions:

Simplify[
      Series[((A m (λ^2) + B λ)/(C m λ + (((m λ - 1))^2))), {m, ∞, 2}]]

Simplify[
      Series[((A m (λ^2) + B λ)/(C m λ + (((m λ - 1))^2))), {λ, 0, 2}]]

Simplify[
      Series[((A m (λ^2) + B λ)/(C m λ + (((m λ - 1))^2))), {m, 1/(2  λ), 2}, {λ, 0, 
          2}]]

To get agreement of the 3 expansions we need to take A = Q, B = 3 Q and solve for C:

Solve[((2 (A + 2 B) )/(1 + 2 C)) == ((6 Q)/5) /. {A -> Q, B -> 3  Q}, C]

So now the trial expression simplifies to

Simplify[((A m (λ^2) + B λ)/(C m λ + (((m λ - 1))^2))) /. {A -> Q, B -> 3  Q, 
          C -> 16/3}]

LeafCount[((3 Q λ)/(1 + 3 m λ))]

And as before, a very simple approximation is obtained.  See how well it works:

Apx[λ_, Δ_, m_] := 
      QQ[Δ] ((3 λ)/(1 + 3 m λ))

Apx[0.001, 0.1, 10]

First look at the (λ,Δ) dependence of the approximation:

Plot3D[
      1 + Apx[λ, Δ, 1], {λ, 0, 
        1}, {Δ, -0.2, 0.2}, 
      ViewPoint -> {-1.226,  2.999,  0.975}, 
      PlotRange -> {1, Automatic}, 
      ColorOutput -> GrayLevel]

Looks qquite close to the exact values:

Plot3D[
      A1[λ, Δ], {λ, 0, 
        1}, {Δ, -0.2, 0.2}, 
      ViewPoint -> {-1.226,  2.999,  0.975}, 
      PlotRange -> {1, Automatic}, 
      ColorOutput -> GrayLevel]

To see it more clearly, this is the difference plotted on a similar scale:

Plot3D[
      A1[λ, Δ] - 1 - 
        Apx[λ, Δ, 1], {λ, 0, 
        1}, {Δ, -0.2, 0.2}, 
      ViewPoint -> {-1.226,  2.999,  0.975}, 
      PlotRange -> {0, 1.5}, ColorOutput -> GrayLevel]

So for small λ, the agreement is pretty good; main deviation is for Δ -> -0.2. For positive Δ, agreement within say 10% even at a large λ->1. The deviation is consistently an underestimate of the enhancement, i.e. on the conservative side. 
But main interest is for time evolution of plume, which comes from m-dependence of the enhancement factor. So lets look at that:

LogLinearPlot[{1 + Apx[0.001, 0.1, m], 
        Am[0.001, 0.1, m], (10^10) (m - 1/(2*0.001))}, {m, 1, 
        10000}, PlotRange -> {1, 1.0002}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Red is the approximation, green the exact expression. And once more, the fit is almost perfect! Try large λ values:

LogLinearPlot[{1 + Apx[0.01, 0.1, m], 
        Am[0.01, 0.1, m], (10^10) (m - 1/(2*0.01))}, {m, 1, 
        10000}, PlotRange -> {1, 1.002}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

LogLinearPlot[{1 + Apx[0.1, 0.1, m], 
        Am[0.1, 0.1, m], (10^10) (m - 1/(2*0.1))}, {m, 1, 1000}, 
      PlotRange -> {1, 1.02}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

LogLinearPlot[{1 + Apx[1, 0.1, m], 
        Am[1, 0.1, m], (10^10) (m - 1/(3*1))}, {m, 1, 100}, 
      PlotRange -> {1, 1.05}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Beyond m = 10 the approximation remains excellent throughout. For low m values, it deteriorates as λ increases, but remains qualitatively OK even up to probably unrealistic fluctuation lengths of the order of Λ. In fact, the consistent underestimate of low order terms simply means that a constant correction factor (i.e., independent of M) will be needed in front of the product for large λ's, the behaviour as function of m derived from large m terms remain correct.

A striking thing when comparing the results of this grid model with the unequal step length model in stochasflow7.nb, is that i) Dependence on Δ is separated off into the overall factor Q, and ii) dependence on m and λ is not just simple but in fact identical in both cases when allowing for the parameter change τ0 -> -1/2λ . 

 In fact, all of the difference between the two cases is contained in different expressions for the Q function. Let us see how these compare; call the old one Quneq:

Quneq[Δ_] := ((4 (Δ^2) (3 - 2 Δ + 24 (Δ^2) + 2 (Δ^3) + 5 (Δ^4)))/(3 (((-1 + Δ))^2) (((1 + Δ))^4)))

Plot[{QQ[Δ], 
        Quneq[Δ]}, {Δ, -0.3, 0.3}, 
      PlotRange -> {0, 0.5}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

The grid case Q function (red) is larger but for positive Δ,  quite similart the earlier one.

As remarked when looking at A1, there is asymmetry and this is larger for the grid case (red) than for the unequal length case (green). The reason for this is demonstrated by seeing how Δ values relate to the velocity ratios V2/V1 (red) and V3/V1 (green) in the two cases. First for unequal step lengths:

Plot[{((1 + Δ)/(1 - Δ)), ((1 - Δ)/(1 + Δ))}, {Δ, -0.3, 
        0.3}, PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Now the fixed grid case:

Plot[{((1 + Δ)/(1 - Δ)), ((1 + Δ)/(1 + 3  Δ))}, {Δ, -0.3, 
        0.3}, PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

All in all, the most important observation is that while the details of Q differ depending on how the fluctuation parameters were chosen, it remains true that Q is a small number, around 0.1 for any plausible step size.

The fact that the m- and λ dependence in the two cases are identical, means that in principle the further work in stochasflow7.nb remains valid and only a different numeric Q might be used in the final formula. However, we would also like to update the notation from τ0 to λ so key steps are repeated here in the new notation.

The step size parameter needed to get a Q of 0.1:

QQ[0.14]

The compined enhancment for a total of M steps:

AM[λ_, Q_, 
        M_] = (∏_(m = 1)^M)(1 + ((3 Q λ)/(1 + 3 m λ)))

The Stirling approximation applies for λ << 1 as before and we jump directly to the approximate form by manually replacing 2 τ0 -> -1/λ in the Abigm formula of stochasflow7.nb:

AMapx[λ_, Q_, 
        M_] := (((1 + ((3  Q λ)/(3  M λ + 1))))^M) (((1 + ((3 M λ)/(3 Q λ + 1))))^Q) (((1 - ((9 M Q (λ^2))/((3 M λ + 1) (3 Q λ + 1)))))^((1/2) + (1/(3 λ))))

Check this by plotting, deliberately choosing a rather large λ otherwise you won't see the  difference:

Plot[{AM[0.1, 0.1, M], AMapx[0.1, 0.1, M]}, {M, 1, 10000}, 
      PlotRange -> {1, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Over a longer range of M values:

LogLinearPlot[{AM[0.1, 0.1, M], AMapx[0.1, 0.1, M]}, {M, 1, (10^6)}, 
      PlotRange -> {1, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

LogLinearPlot[{AM[0.001, 0.1, M], AMapx[0.001, 0.1, M]}, {M, 1, (10^6)}, 
      PlotRange -> {1, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

As before, be careful of the shape distortion by the log scale, the value really does go to 1 at M = 0 as shown by the linear one.

This time, do a log-log plot for a realistic λ as is done for the dispersivity later on:

LogLogPlot[{AM[0.001, 0.1, M], AMapx[0.001, 0.1, M]}, {M, 1, (10^4)}, 
      PlotRange -> {1, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Here is a plot in terms of scaled traversal length, as eventually (stochasflow10.nb) it was decided that the derivative thing for the dispersivity is not so good and in fact the enhancement factor is itself just a scaled version of dispersivity:

LogLogPlot[{AMapx[0.001, 
          0.3, (Γ/(3 (0.001)))]}, {Γ, 0.01, 
        1000}, PlotRange -> {{0.01, 1000}, {1, 10}}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

Algebraic expressions in terms of scaled traversal length

Simplify[
      AMapx[λ, Q, (Γ/(3  λ))]]

Simplify[
      AM[λ, 
          Q, (Γ/(3  λ))] /. {Pochhammer[a_, n_] -> 
            Gamma[a + n]/Gamma[a]}]

Dispersivity as function of traversal length

Go through similar steps as in stochasflow7.nb, but for the grid case we have V1 t = 3 M L (differs by factor 2 from L definition in stochasflow7) so eliminating M in favour of t we have

AMt[λ_, Q_, t_] = 
      AMapx[λ, Q, ((- t/t0)/(3 λ))]

Next, define the variance factor and find its derivative. Once the derivative has been done, it is best to eliminate time in favour of the traversal length defined as TL = V1 t . Even better, replace TL by its dimensionless scaled version defined as Γ = TL / Λ; using Λ = -V1 t0. (The notation here is different from stochasflow7 here - Γ is now the scaled traversal length, there it was the actual length) This yields

AMt[λ, Q, 0]

varf[t_] := (t - t0) AMt[λ, Q, t]

Simplify[(∂_t) 
          varf[t] //. {t -> -Γ t0}]

Simplify[AMt[λ, Q, -Γ t0]]

As before, keep the factors in AMt and collect the rest together in a prefactor, applying also the log approximation justified by the small sizes of Q, and of λ compared to Γ :

Pre[λ_, Q_, Γ_] = 
      Simplify[((((((1 + Γ + 3 Q λ)/(1 + Γ + 3 Q λ + 3 Q Γ λ))))^(-1)) (3 λ (2 + 2 Γ + 
                          3 Q λ) + 
                    2 (1 + Γ) (1 + Γ + 
                          3 Q λ) Log[
                        1 + ((3 Q λ)/(1 + Γ))]))/(6 (1 + Γ) λ (1 + 
                    3 Q λ)) //. {Log[1 + x_] -> 
              x}]

To see how well the approximation does, one has to go to unrealistically large λ value = 10 to see any deviation, and then only for small Γ:

LogLinearPlot[
      Evaluate[{Pre[10, 
            0.1, Γ], ((((((1 + Γ + 3 Q λ)/(1 + Γ + 3 Q λ + 3 Q Γ λ))))^(-1)) (3 λ (2 + 2 Γ + 
                            3 Q λ) + 
                      2 (1 + Γ) (1 + Γ + 
                            3 Q λ) Log[
                          1 + ((3 Q λ)/(1 + Γ))]))/(6 (1 + Γ) λ (1 + 
                      3 Q λ)) /. {λ -> 10, 
              Q -> 0.1}, 1}], {Γ, 2, 50000}, 
      PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

For any plausible λ value, the prefactor really just reduces to a constant:

LogLinearPlot[{Pre[0.001, 0.1, Γ], 1}, {Γ, 
        2, 50000}, PlotRange -> {0, Automatic}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

The value of this constant is obtained by taking the λ -> 0 limit:

Simplify[Pre[0, Q, Γ]]

Simplify[Pre[λ, Q, 0]]

So then, the full expression for the dispersion enhancement factor Den is a product of 3 terms:

Den[λ_, Q_, Γ_] := 
      Den0[λ, Q, Γ]*Den1[λ, Q, Γ]*
        Den2[λ, Q, Γ]

Den0[λ_, 
          Q_, Γ_] = (1 + 
              Q) (((((1 + Γ + 3 Q λ)/(1 + Γ + 3 Q λ + 3 Q Γ λ))))^((1/2) + (1/(3 λ))));

Den1[λ_, 
          Q_, Γ_] = (((1 + ((3 Q λ)/(1 + Γ))))^((Γ/(3 λ))));

Den2[λ_, 
          Q_, Γ_] =  (((1 + (Γ/(1 + 3 Q λ))))^Q);

Den[λ, Q, 0]

LogLinearPlot[{Den0[0.0001, 0.1, Γ], 
        Den0[0.1, 0.1, Γ], 
        Den0[1, 0.1, Γ]}, {Γ, 0.001, 100}, 
      PlotRange -> {0.8, 1.6}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

LogLinearPlot[{Den1[0.0001, 0.1, Γ], 
        Den1[0.1, 0.1, Γ], Den1[1, 0.1, Γ], 
        Exp[0.1]}, {Γ, 0.001, 100}, 
      PlotRange -> {0.8, 1.6}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

LogLinearPlot[{Den2[0.0001, 0.1, Γ], 
        Den2[0.1, 0.1, Γ], 
        Den2[1, 0.1, Γ]}, {Γ, 0.001, 100}, 
      PlotRange -> {0.8, 1.6}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

LogLinearPlot[{Den[0.0001, 0.1, Γ], 
        Den[0.1, 0.1, Γ], 
        Den[1, 0.1, Γ]}, {Γ, 0.001, 100}, 
      PlotRange -> {0.8, 1.6}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

While the detailed values are slightly different, the qualitative behaviour is very similar to that found in stochasflow7.nb for the case of unequal steps. 

Replot using log scales on both axes, and choose the range to fit Λ = 100 and then cover a similar range of lengths as the experimental values from Gelhar.

LogLogPlot[{Den[0.001, 0.1, Γ], 
        Den[0.1, 0.1, Γ], 
        Den[1, 0.1, Γ]}, {Γ, 0.01, 1000}, 
      PlotRange -> {1, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

The curve above shows two distinct regions, separated by a transition around Γ = 1. Below this, the enhancement rises "slowly" i..e., because of the log scale, according to a low power law; but above, it changes to a higher power law. This characterisation somewhat distorts the "real behaviour; look at a normal linear graph of the same thing:

Plot[{Den[0.001, 0.1, Γ]}, {Γ, 0.01, 400}, 
      PlotRange -> {1, 2}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

So really the rise is fast in the beginning (as for a fractional power) and then settles down to a slower rise (albeit that for a higher, but still fractional, power law).

Comparison with experimental results

Comparing the graphs above with that of Gelhar (Fetter, fig 2.18) , the rise according to a fractional power law at large traversal lengths agrees well. Also, the separation into two regions agrees. The most obvious inconsistency is that the experimental results finds in effect a higher power law in the low Γ region, while the theoretical one above found a lower power law. As a direct result of this, the overall enhancement achieved over the same length range is only a factor 2, while the experimental ones found an enhancement of about (10^4). This looks like a large discrepancy, but bear in mind that already the factor 2 is itself a huge enhancement compared to the factor contributed by a single fluctuation; the largest of which, the one for the first fluctuation, is for the parameters chosen, only a factor

AM[0.001, 0.1, 1]

The parameter chosen for the theoretical curves so far were mainly based on plausibility arguments, e.g. that the flow velocity changes by not more than say 20% up or down in a fluctuation. So it seems worthwhile to reverse the argument, and see if parameter values could be chosen to represent the observations, and if such values seem reasonable.

To avoid having to use the large number of experimental points in fig 2.18 all the time, we use instead the line α = 0.1 TL which is shown on that figure as a basis of comparison. First, we need to convert the α value to an enhancement factor by dividing by the zero length ("unperturbed") longitudinal dispersivity of (10^(-2)) m on the figure and the length needs to be converted to the scaled parameter Γ by using the appropriate value of Λ. The equation of the line then becomes e = 0.1/0.001 Λ Γ = (10 Λ) Γ and the value Λ = 100 m used so far does not seem unreasonable as the figure seems to show a transition between low range and high range behaviour somewhere between 10 and 100 m. As a rough estimate of the high range power law exponent, Q = 0.3 seems compatible with the experiment. The value of λ is not  sensitive as long as it is << 1; the value 0.001 would mean fluctuation grid scale of 0.1 m which seems OK:

LogLogPlot[{Den[0.001, 0.3, Γ], 
        1000  Γ}, {Γ, 0.01, 1000}, 
      PlotRange -> {{0.01, 1000}, {1, (10^6)}}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

OK, so the high range slope seems fine but the problem is the low range rise that is too small. Suppose we increase the coefficient of Q in the factor Den1 that dominates the low range behaviour. For example

Den1[λ_, 
          Q_, Γ_] = (((1 + ((75 Q λ)/(1 + Γ))))^((Γ/(3 λ))));

LogLogPlot[{Den[0.001, 0.3, Γ], 
        1000  Γ}, {Γ, 0.1, 1000}, 
      PlotRange -> {{0.1, 1000}, {1, (10^6)}}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

The tendency is fine but the transition takes place at Γ = 10 which for the chosen Λ = 100 m means at 1000 m, whereas the experimental one is at about 50 - 100 m . If we look at the plots of Den1's behaviour, it tends to rise up to a value of about Γ = 10 (as a result of the (1+Γ) factor it contains; it saturates when Γ >>1 ). So to fit the data it seems that we must adjust Λ down to about 10 m to get the transition to occurr at the right place:

Λ = 10 m

LogLogPlot[{Den[0.001, 0.3, Γ], 
        100  Γ}, {Γ, 0.1, 10000}, 
      PlotRange -> {{0.1, 10000}, {1, (10^6)}}, 
      PlotStyle -> {Hue[0.0], Hue[0.3], Hue[0.5], Hue[0.8]}]

That appears to fit the obseved behaviour quite well (within the large spread of experimental values), and the values arrived at for Q, Λ and λ all seem reasonable. 
The main change required was to increase the coefficient of Q in the Den1 expression by a factor of 25. However, this was done in an ad hoc way and it is for example unclear if such a change can be obtained in a consistent way from some modifcation of the underlying response of a single fluctuation.So the next step is to see if a plausible explanation of that can be found.