Resolution of the shape of a Beach by Bruun 64 & Dean 91 model

    It is believed that the shore cross profile is a power law Z=Ax^n where x is distance offshore (x<100m), 0<Z<10m is water depth, averaged measurements give n=0.66 (in fact 0.03<x<1.4 see Pilkey 94, we did himself a mistake in the equation). Bruun obtained that from fits of measurements.

    To establish this, Dean supposes that there is a uniform wave energy dissipation constant per unit volume D_*, and that within the breaking zone the flux of energy is steady and dissipated at a rate D_*: \displaystyle D_*= \frac{1}{h}\frac{d}{dx}(E C_g) where E and C_g are the energy density and group velocity. The waves come from the right to the left, we suppose that the mean sea level defines level 0, as E \propto Z^2 and C_g \propto \sqrt{Z}, hence the equation without dimension is : \displaystyle \frac{d}{dx}(Z^{5/2})= Z we solve it with a relaxation as: \displaystyle \frac{\partial}{\partial x}Z + \frac{\partial}{\partial x}(Z^{5/2})= Z the method is similar to which explains the notions of advection, testing the flux, coded with Basilisk (see the end, maybe $Z^2 $ is better).

    Starting by a slope and a flat bottom, we wil compute evolution to the steady power law
    set xlabel "x"
    set ylabel "depth"
    set arrow from 3,-(3.*3./5)**(2./3) to 3,0
    set label "depth -Z(x)" at  3.5,-1
    set key left bottom
    p[0:10][]  -0.25*(x<1? x:1) t'initial', -(3*x/5)**(2./3) t'expected' w l  linec 3 


    Note that here Z is positive, in regular Shallow Water configuration we use z_b which is -Z, as here \eta=0.


    mandatory declarations:

    #include "grid/cartesian1D.h"
    #include "run.h"

    definition of the field Z, the flux F, Boundary conditions

    scalar Z[];
    scalar F[];
    Z[left] = dirichlet(0);
    Z[right] = dirichlet(5);;  
    double m=5./2;  // change for tests, m=1 advection, m=2 Burgers
    double flux(double z)
      return pow(fabs(z),m);
    double celerity(double z)
      return fmin(m*pow(fabs(z),m-1.),10000);

    Main with definition of parameters, note that time step is small

    int main() {
      L0 = 10.;
      X0 = 0;
      N = 256;    
      DT = (L0/N)/16/4; 

    initial elevation: a “slope and a flat bottom”

    event init (t = 0) {
        Z[] = .25*fmin(x,1);
      boundary ({Z});

    begin the time loop, print data, in practice a max time of 5 is enough.

    event printdata (t += .5; t <=5)  
        fprintf (stdout, "%g %g %g  \n", x, Z[], t);
      fprintf (stdout, "\n\n");

    integration step, at each time step

    event integration (i++) {
      double dt = DT;
      double cDelta;

    finding the good next time step

      dt = dtnext(dt);

    the algorithm is based on a split in time, the flux and the source term.

        cDelta = (celerity(Z[])+celerity(Z[-1]))/2.;
        F[] = (flux(Z[])+flux(Z[-1]))/2.  - cDelta *(Z[]-Z[-1])/2;
        boundary ({F});

    split: explicit step update \displaystyle Z_i^{n+1}=Z_i^{n} -{\Delta t} \dfrac{F(Z_{i+1})-F(Z_{i})}{\Delta x}

        Z[] +=  - dt* ( F[1] - F[] )/Delta;
      boundary ({Z});

    split: implicit resolution of \dfrac{\partial Z}{\partial t} = Z : \displaystyle Z_i^{n+1}=Z_i^{n} +{\Delta t} Z_i^{n+1}

        Z[] = Z[]/(1 - 1.0 * dt); //Z[]*(1 + 1.0 * dt);

    Boundary condition

      boundary ({Z});


    Then compile and run:

    qcc  -g -O2 -DTRASH=1 -Wall  bruundean.c -o bruundean ;./bruundean > out


    The analytical solution solves the equilibrium flux/dissipation \displaystyle dZ(x)^{m}/dx = Z so that Z=(\frac{(m-1)}{m} x)^{\frac{1}{m-1}} with m=5/3, so that we obtain Z(x)= (3 x/5)^{2/3} the Bruun’s solution.

    Note that there is obviously a problem of B.C. at the right:

     set xlabel "x"
     set ylabel "-Z"
     set key left bottom
     p[:][-3.5:0]'out' u ($1):(-$2)t'num'w l,-(3*x/5)**(2./3) t'exact' w l linec 3 


    which is -Z(x,t) plotted here for t=0 .5 1 1.5 2 …

    time evolution with splot:
     set xlabel "x"
     set ylabel "t"
     set zlabel "Z"
     sp [:][:][0:5]'out' u 1:3:2 w l,(3*(x)/5)**(2./3) t'exact' w l    



    Solve now with \partial_t Z^2 instead of \partial_t Z, to have an energy like equation in Z^2

    \displaystyle \frac{\partial}{\partial x}Z^2 + \frac{\partial}{\partial x}(Z^{5/2})= Z

    which is as well, m=5/2

    \displaystyle \frac{\partial}{\partial x}Z + \frac{\partial}{\partial x}(\frac{m}{2 (m-1)}Z^{m-1})= \frac{1}{2}


    • Lagrée P-Y “Equations de Saint Venant et application, Ecoulements en milieux naturels” Cours MSF12, M1 UPMC

    • Robert G. Dean “Equilibrium Beach Profiles: Characteristics and Applications” Journal of Coastal Research, Vol. 7, No. 1, 1991

    • O. Pilkey, Jr. “Mathematical Modeling of Beach Behavior Doesn’t Work”