# 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 http://basilisk.fr/sandbox/M1EMN/BASIC/advecte1.c 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.

## Code

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;
run();
}

initial elevation: a “slope and a flat bottom”

event init (t = 0) {
foreach(){
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)
{
foreach()
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.

  foreach()
{
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}

  foreach()
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}

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

Boundary condition

  boundary ({Z});
}

## Run

Then compile and run:

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

## Results

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    

## Exercise

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}

## Bibliography

• 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”