sandbox/lopez/guadalquivir.c

Estuary of the river Guadalquivir

The river Guadalquivir is an old river with a lot of history. Despite its low depth the river has been always navigable. In fact, it was used by the vikings to atack Sevilla in the IX century. Today, the river is still navigable up to Sevilla where a port is in use. Ships with large drafts must ride the tide to sail up to the port of Sevilla.

The Estuary of the Guadalquivir

The river discharges in the Atlantic sea. The Doñana National Park is located on the right riverbank. The village of Sanlucar de Barrameda is on the left bank and is famous for its wines.

Simulation

The simulation of the estuary is performed using an explicit scheme for the Saint-Venant shallow water equations with adaptation. That means to use the most general model hydro.c assuming a single layer and that the vertical momentum equation reduces to the hydrostatic relationship. These assumptions are the default one. So nothing has to be set in addition. Details of the layered model can be found in Popinet, 2020.

#include "grid/multigrid.h"
#if !SV
#include "layered/hydro.h"
#else
#include "saint-venant.h"
#endif
#include "terrain.h"

riemann.h contains kurganov() which is required in discharge.h.

#if !SV
#include "riemann.h"
#endif
#include "discharge.h"

#define MINLEVEL 4
int maxlevel = 8;
scalar river[];

The sea goes up and down sinusoidally because of a semidiurnal moon tide M2 of amplitude equal to 0.7 m. The period of this component is of 12.42 h (or 745.2 min). We assume that the origin of times is in phase with the tide. We have decided to measure the lengths in meters and the times in minutes. There is no problem with it meanwhile we were consistent with the system of units elected.

#define OM  (2*pi/745.2) // frequency (The M2 period is 12.42h)
#define PHASE 0. // Time delay (in case the time origin is shifted)
#define AMP 0.7 // Amplitude

u.n[bottom] = - radiation(AMP*sin(OM*t + PHASE));
u.n[top]    = + radiation(AMP*sin(OM*t + PHASE));
u.n[left]   = - radiation(AMP*sin(OM*t + PHASE));

In the main() we define the left bottom corner of the domain (with origin()) and the size of the computational domain of 21000 x 21000 m^2. Note we are using a cartesian UTM proyection (ED50 H 30N).

int main()
{
  size (21000);
  origin (189000, 4071000);
  G = 9.81*sq(60); //The time is measured in minutes 
  init_grid (1 << maxlevel);
  run();
}

In the init event we initialize the bathymetry zb using the database guadalquivir. The bathymetry must go beyond the simulation domain. Un particular must be availabe even for the coarsest ghost cells. If not a fake ground is set in the boundary (run the present test with minlevel equal to 2 and see). The folder where the database locates must be well specified. Also, conserve_elevation() serves to keep the depth of the fluid instead of the fluid volume when coarsening/refining the grid. If not, gravity waves are induced in every change of grid.

event init (i = 0)
{
  terrain (zb, "/home/basilisk/terrain/guadalquivir", NULL);

  if (restore (file = "dump"))
    conserve_elevation();
  else {
    conserve_elevation();
    
    foreach() {
      h[] = max(0., - zb[]);
      river[] = (y > 4087935  && y < 4088357);
    }
    boundary ({h});
  }
  output_ppm (zb, file = "zb.png");
}

The bathymetry is initially described with a grid size equal to \Delta = L0/2^N \simeq 82 m where N=8 is the level of refinement used.

Bathymetry

The course of the river is very well described because of the high number of points in the database at that location. Not the same happens with some dry parts. The bathymetry also shows that the low curse of the Guadalquivir goes through a marshland.

An inflow enters through the right. The inflow depends on time. Over a constant flow rate Qo = 150 m^3/s a flood having the form of a gaussian bell of height Qp = 2000 m^3/s is set. The peak of flood enters the domain with a delay of 500 minutes after the first peak of the tide. If the riverbed is not well distinguished from the riverbank, the inflow could spread along the right boundary as it trys to compute a water elevation \eta compatible with the inflow. To avoid this unwanted mechanism we force the inflow to enter through the region for which river is equal to 1. That is, through the narrow gap going from y_{low} = 4087935 up to y_{up} = 4088357.

#define gaussian(t,to,b) (exp(-sq((t-to)/(2.*b))))
#define Qo (150*60)
#define Qp (2000.*60)

double etar;
event inflow (i++) {
  double qlaw = Qo + Qp*gaussian (t, 500, 100);
  etar = eta_b (qlaw,right, river, 1);

  h[right] = (river[] ? max (etar - zb[], 0.) : 0.);
  eta[right] = (river[] ? max (etar - zb[], 0.) : 0.) + zb[];
}

Friction slow down the fluid according to,

\displaystyle \frac{d \mathbf{u}}{d t} = - \frac{\lambda}{4 r_h} |\mathbf{u}| \mathbf{u} = - \frac{c_f}{r_h} |\mathbf{u}| \mathbf{u}

where \lambda is the Darcy friction factor. c_f = \lambda/4 is the Fanning friction factor. The hydraulic radius is in this case r_h = h. We set c_f = 10^{-3}. The above equation is linearized and solved implicitly,

\displaystyle \frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} = - \frac{c_f}{h} |\mathbf{u}|^n \mathbf{u}^{n+1}

#define fa 1e-3 // The Fanning friction factor

event friction (i++) {
  foreach() {
    double a = h[] < dry ? HUGE : 1. + fa*dt*norm(u)/(h[]);
    foreach_dimension()
      u.x[] /= a;
  }
  boundary ((scalar *){u});
}

Adaption

We adapt the grid after the water elevation.

#if TREE
event adapt (i++) {
  scalar etap[];
  foreach() 
    etap[] = h[] > dry ? eta[] : 0.;
  boundary ({etap});
  adapt_wavelet ({etap}, (double[]){3e-4}, maxlevel, MINLEVEL);
}
#endif

Outputs

We save each 100 time step a dump file in case we need to restart the simulation.

event save_dump( i += 100) {
  dump("dumpfile");
}

Also, we save in ASCII ESRI format the inundation area at instant t = 600.

event inundation(t = 600) {
  scalar hp[];
  foreach()
    hp[] = h[] > dry ? h[] : nodata;
  
  FILE * fp = fopen ("inundation.asc", "w");
  output_grd(hp, fp);
  fclose(fp);
}

Also, we keep a time record (with a time step of 1 min) of the amount of water in the domain as well as we monitorize the relevant variables near the “Salina of Monte Algaida”.

event logfile(t++) {
  double vol = 0;
  foreach(reduction(+:vol))
    vol += h[]*sq(Delta);
  Point point = locate(203000, 4089170); //Near "salina del monte algaida"
  fprintf(stderr, "%g %g %g %g %g\n",
	  t, vol, h[], eta[], sign(u.x[])*sqrt(sq(u.x[])+ sq(u.y[])));
  fflush(stderr);
}

om=2.*pi/745.2
gaussian(x)=exp(-((x-500)/(2.*100))**2)
set multiplot layout 3, 1
set arrow from 0.,0. to 1200., 0. nohead front lc rgb "black" lw 2  dashtype "-"
plot 'log' u 1:($5/60) w l lw 2 t 'velocity (m/s)', '../guadalquivir-sv/log' u 1:($5/60) w l lw 2 t 'velocity (SV)'
plot 'log' u 1:4 w l lw 2 t 'eta', '../guadalquivir-sv/log' u 1:4 w l lw 2 t 'eta (SV)'
set xrange [0 : 1200]
plot sin(om*x) lw 2 t 'tide', gaussian(x) lw 2 t 'flood'
set xlabel "t (min)"
unset multiplot

top: velocity at Salina del monte Algaida, middle: water elevation low: Time dependece of the tides and the flood (script)

Finally, we record a couple of movies illustrating the time evolution of the elevation \eta and the norm of the velocity |u|.

#include "view.h"
event snapshot(t++; t < 1200) {
  char str[40];
  sprintf (str, "t = %g ", t);
  
  scalar etap[], vmax[];
  foreach() {
    etap[] = h[] > dry ? eta[] : nodata;
    vmax[] = h[] > dry ? norm(u) : nodata;
  }
  boundary({etap, vmax});

To center the graph in the canvas we set the translation tx = -X0/L0 - 0.5 and ty = -Y0/L0 - 0.5 where X0 and Y0 son the coordinates of the left bottom corner.

  view (fov = 22, tx = -X0/L0 - 0.5, ty = -Y0/L0 - 0.5);
  draw_string (str, 2, size = 25, lw = 1); 
  squares("etap", linear = false);
  save ("eta.mp4");
  clear();
  squares("vmax", linear = false);
  vectors("u", scale = 2);
  draw_string (str, 2, size = 25, lw = 1);
  save ("vel.mp4");
}