sandbox/nlemoine/groundwater/Meuvaines/meuvaines.c

    Return to my homepage

    Hydrodynamic modeling of the Meuvaines / Ver-sur-Mer coastal wetland (Calvados, France)

    Picture : Conservatoire du Littoral

    Picture : Conservatoire du Littoral

    #include "grid/quadtree.h"
#include "run.h"
#include "diffusion.h"
#include "input.h"
#include "output.h"
#include "terrain.h"
#include "nlemoine/vector-geom.h"
#include "view.h"
#include "nlemoine/view-utils.h"
#include "nlemoine/groundwater/Meuvaines/colormaps.h"

    Definitions and global variables

    #define LEVEL 8
#define MINLEVEL 5
#define HE 1.

    We first define some space / time constants

    // Chart datum and mean sea level (RAM_PACK data from SHOM)
#define ZH_ref_Arromanches -4.019
#define NM_Arromanches 4.43
#define ZH_ref_Courseulles -3.99
#define NM_Courseulles 4.59

    1-hour runtime limit on the server does not allow to simulate more than about 18 month with LEVEL = 8. Just change the tfin value (in days) to simulate a longer period, e.g. #define tfin 1462. to simulate from 1963-01-01 until 1966-12-31 and validate against observed streamflow in 1965-1966 (a two-year warm-up period covering 1963-1964 is necessary to remove the effet of the (unknown) initial head distribution.

    #define tfin 518.
#define DAYS_PER_YEAR 365.25
#define DAYS_PER_MONTH 30.4375

// Some reference points (RGF Lambert 93, EPSG:2154)
double XlimW,XlimE,YlimS,YlimN;
coord SegNeumannW = {430330.,6912160.};
coord SegNeumannE = {435380.,6911570.};
coord n_int;
double ZMER;
double datenum0 = 1614.; // initial date expressed as elapsed days since 1958-08-01, hence 1963-01-01

    A small hack to convert elapsed simulation time into year & month (avoids bringing out the big guns of “time.h”) Assumes that t+datenum0 is the number of (fractional) days elapsed since the begining of SAFRAN forcings on 1958-08-01:

    int hack_date(double tt, int * _year, int * _month)
{
  int elapsed_years_since_1958 = (int) floor((tt+datenum0+212)/DAYS_PER_YEAR);
  int elapsed_months_since_1958 = (int) floor((tt+datenum0+212)/DAYS_PER_MONTH);
  *_month = (elapsed_months_since_1958 % 12) + 1; 
  *_year = 1958 + elapsed_years_since_1958;
  return(0);
}

typedef struct {
  coord Site;
  char BSS_ID[12];
  int layer;
} Piezo;

struct TimeSeries Forcings;
double CumP,CumPE,P_dt,PE_dt;
double SMA0,SMA1,SMA2,Excess0,Excess1,Excess2;
double maxStorage0; // marls soil storage
double * Forc_t;
double oneday = 1. [0,1];
double sum_dt_month;

double PARAM[11] = {200., 200., 2.e-4, 1.e-5, 0.08, 3.e-4, 150., 1.5e-5, 1.e-4, 0.025, 1.e-5};
int nparam;

scalar isNeumann[], isDirichlet[];
scalar zb[];
scalar msk1[],msk2[];

struct Polygon PolyDirichlet;
struct Polygon BV_Gronde;
struct Polygon BV_Provence;
double QM_Gronde, QM_Provence;
Piezo * PIEZO;
int npiezo;
double * HPIEZO;

scalar maliere_mur[], marnesPB_mur[], bathonien_mur[];

    All scalars are indexed (1) for the semi-confined aquifer in Bajocian limestone, and (2) for the unconfined aquifer in mid-Bathonian limestone.

    scalar zbottom1[],zbottom2[]; // elevation of the base of each aquifer layer
scalar ztop1[],ztop2[];
scalar h1[],h2[]; // groundwater head
face vector T1[], T2[]; // transmissivity
scalar S1[],S2[]; // storage
scalar r1[],r2[]; // recharge
scalar beta1[],beta2[];
scalar seep1[],seep2[];
scalar isConfined1[],isConfined2[]; 
scalar hasThickness1[],hasThickness2[];
scalar recharge_layer[];

double dt;
mgstats mgd;

// Definition of hydrodynamic attributes (uniform for the moment, but could be scalar [])
// They will be attached to the scalars zbottomN[]

attribute {  
  double maxStorage;
  double Ksat; // Saturated hydraulic conductivity
  double Tf; // fixed transmissivity component
  double Ss; // Spectific storage (confined case)
  double omega_d; // porosité efficace (cas libre)                           
} 
 
FILE * resfile, * flowfile;

    Main

    #include "nlemoine/groundwater/Meuvaines/aquifer_ML.h"

int main ()
{
  if(pid()==0){
    
    system("wget https://dropsu.sorbonne-universite.fr/s/Zg24YkSA3w2e7wk/download && mv download data.zip");
    system("unzip data.zip");
    
    resfile = fopen("results.csv","wt");
    flowfile = fopen("streamflow.dat","wt");
  }

  run(); 

  fclose(resfile);
  fclose(flowfile);
  return(0);
}

    Initialization

    event init (i = 0)
{
  // Mean sea level

  ZMER = 0.5*((NM_Arromanches+ZH_ref_Arromanches)+
              (NM_Courseulles+ZH_ref_Courseulles));
    
  // Max resolution: (2^12) * (2^12) cells of 5 m 

  L0 = 5. * (double)(1 << 12);

  // Initialize grid

  XlimW = 429162.5; // East of Port-en-Bessin
  XlimE = 449102.5; // Courseulles
  YlimS = 6911432.5;
  YlimN = YlimS+L0;

  size (L0);
  X0 = (XlimW+XlimE)/2. - L0/2;
  Y0 = YlimS;
  origin (X0, Y0);
  N = 1 << LEVEL;
  init_grid (1 << LEVEL);

  // Inward normal to the domain at the south border,
  // defined by a segment approximately following the N13 road south of Bayeux
  // between Aure river (SegNeumannW) and Seulles river (SegNeumannE)

  n_int = (coord){-(SegNeumannE.y-SegNeumannW.y),+(SegNeumannE.x-SegNeumannW.x)};
  double norm = sqrt(sq(n_int.x) + sq(n_int.y));
  n_int.x /= norm;
  n_int.y /= norm;

  // Read polygon specifying Dirichlet conditions (Aure and Seulles rivers)

  (void) ReadPolygon ( "PolyDirichlet_v2.dat", & PolyDirichlet);
  printf("%d vertices read.\n",PolyDirichlet.nv);

  // Read contours of gauged catchments

  (void) ReadPolygon ( "Gronde_SGR_simplif.dat", & BV_Gronde);
  printf("%d vertices read.\n",BV_Gronde.nv);

  (void) ReadPolygon ( "Provence_SGR_simplif.dat", & BV_Provence);
  printf("%d vertices read.\n",BV_Provence.nv);

  // Read locations of observation boreholes
    
  (void) read_piezo_sites ( & PIEZO, "BSS_Basilisk.dat", & npiezo);

  HPIEZO = (double *) malloc(npiezo*sizeof(double));

  if(pid()==0)
  {
    fprintf(resfile,"# PARAM :");  
    for(int kp=0;kp<nparam;kp++)
      fprintf(resfile," %g",PARAM[kp]);
    fprintf(resfile,"\n");
          
    fprintf(resfile,"t,P,PE,Q_Gronde,Q_Provence");
    for(int kp=0;kp<npiezo;kp++)
      fprintf(resfile,",%s",PIEZO[kp].BSS_ID);
    fprintf(resfile,"\n");
  }
    
  // Load bathy-topographic data

  terrain (zb, "bessin_topo", NULL);
  boundary({zb});

  // Load interface elevations

  terrain (maliere_mur, "maliere_mur", NULL);
  terrain (marnesPB_mur, "marnesPB_mur", NULL);
  terrain (bathonien_mur, "bathonien_mur", NULL);

  (void) update_geom_stack();

    We set the hydrodynamic properties of each aquifer layer according to the content of the PARAM vector. Specific storage (for the confined case) is a volume of water released by a unit volume of porous medium in response of a unit head drop, hence in \textrm{L}^3\cdot\textrm{L}^{-3}\cdot\textrm{L}^{-1} i.e. \textrm{L}^{-1}.

        // Soil storage over marls (0)

    maxStorage0 = PARAM[0]; // in mm
    
    // Bajocian semi-confined aquifer (1)
    
    zbottom1.maxStorage = PARAM[1];
    zbottom1.Ksat = PARAM[2]*86400.; // convert m/s => m/day
    zbottom1.Tf = PARAM[3]*86400.;  // convert m2/s => m2/day   
    zbottom1.omega_d = PARAM[4];
    zbottom1.Ss = PARAM[5]; // specific storage in m^-1
    
    // mid-Bathonian unconfined aquifer (2)
    
    zbottom2.maxStorage = PARAM[6];
    zbottom2.Ksat = PARAM[7]*86400.; // convert m/s => m/day
    zbottom2.Tf = PARAM[8]*86400.; // convert m/s => m/day
    zbottom2.omega_d = PARAM[9];
    zbottom2.Ss = PARAM[10]; // specific storage in m^-1
    
  // Initialized isDirichlet[] and isNeumann[] scalars

  (void) update_BC_cells();

  // Initialize groundwater heads

  scalar zbmin = zb.dmin;
  scalar msk[];
  scalar msk_base_maliere[], msk_base_marnesPB[], msk_base_bathonien[];
  scalar msk_outcrop_maliere[], msk_outcrop_marnesPB[], msk_outcrop_bathonien[];

  foreach()
  {
    msk[] = (isDirichlet[]>0.) | (isNeumann[]>0.) ? -1. : 1. ;
    double drainage_level = zbmin[]<=zb[] ? fmax(zbmin[],ZMER) : fmax(zb[],ZMER);
      
    h1[] = fmax(zbottom1[],drainage_level);
    seep1[] = 0.;
    
    h2[] = fmax(zbottom2[],drainage_level);
    seep2[] = 0.;

    msk_base_maliere[]    = (maliere_mur[]  <=zb[])      && (msk[]>0.) ? 1. : -1.;
    msk_base_marnesPB[]   = (marnesPB_mur[] <=(zb[]-3.)) && (msk[]>0.) ? 1. : -1.;
    msk_base_bathonien[]  = (bathonien_mur[]<=zb[])      && (msk[]>0.) ? 1. : -1.;

    msk_outcrop_maliere[]    = (maliere_mur[]  <=zb[])      && (marnesPB_mur[] >=(zb[]-3.)) && (msk[]>0.) ? 1. : -1.;
    msk_outcrop_marnesPB[]   = (marnesPB_mur[] <=(zb[]-3.)) && (bathonien_mur[]>=zb[])      && (msk[]>0.) ? 1. : -1.;
    msk_outcrop_bathonien[]  = msk_base_bathonien[];
  }

  // Illustrate interfaces
  float qview_x[4],qview_z[4],qview[4];
  (void) gl_axis_to_quat ((float[]){1,0,0}, 0.42*PI, qview_x); // angle w.r.t vertical
  (void) gl_axis_to_quat ((float[]){0,0,1}, 0.72*PI, qview_z); // rotation in horizontal plane
  gl_add_quats(qview_x, qview_z, qview);
  
  char str[120];

  view (fov = 13.5, quat = {qview[0],qview[1],qview[2],qview[3]},
        sx = 1., sy = 1., sz = 40.,
        width = 1200, height = 768);

  translate(x = -X0-L0/2.+0.12*L0,y = -Y0-L0/2.+0.12*L0,z=+40.){
    NCCLASS = 25;
    masked_squares("maliere_mur", linear = true, z = "maliere_mur", mask = msk_base_maliere, min = -170., max = 80., map = discrete_blues);
    masked_squares("marnesPB_mur", linear = true, z = "marnesPB_mur", mask = msk_base_marnesPB, min = -140., max = 110., map = discrete_grays);
    masked_squares("bathonien_mur", linear = true, z = "bathonien_mur", mask = msk_base_bathonien, min = -120., max = 130., map = discrete_reds);
    surf_cells(zb,mask=msk_outcrop_maliere,lc={0.,0.,0.7});
    surf_cells(zb,mask=msk_outcrop_marnesPB,lc={0.4,0.4,0.4});
    surf_cells(zb,mask=msk_outcrop_bathonien,lc={0.7,0.,0.});

    sprintf (str, "                             Blue: Base (solid) and outcrops (wireframe) of Bajocian limestone");
    draw_string (str, 1, size = 100, lc = {0,0,1}, lw = 1);  // from top-left

    sprintf (str, "\n                             Gray: Base (solid) and outcrops (wireframe) of lower-Bathonian marls");
    draw_string (str, 1, size = 100, lc = {0,0,0}, lw = 1);  // from top-left

    sprintf (str, "\n\n                             Red:  Base (solid) and outcrops (wireframe) of mid-Bathonian limestone");
    draw_string (str, 1, size = 100, lc = {1,0,0}, lw = 1);  // from top-left

    sprintf (str, "Z-stretch x 40 ");
    draw_string (str, 3, size = 64, lc = {0,0,0}, lw = 1); // bottom-right
  }
  save("interfaces.png");
  clear();

    3D view of the geological structure

    The coastal wetland is at the edge of the outcrops of both Bajocian limestone and mid-Bathonian limestone. Geological layers have a tilt oriented towards the North-East direction.

    3D geological structure

    3D geological structure

      // Initialize T & S maps

  (void) update_TS(zbottom1,ztop1,h1,T1,S1,isConfined1,hasThickness1,seep1);
  (void) update_TS(zbottom2,ztop2,h2,T2,S2,isConfined2,hasThickness2,seep2);
  boundary(all);

  // Load cumulative P & PE time series
  const char * separators = " ,;|";
  (void) load_timedata(& Forcings,"safran.csv",1,separators,datenum0); 
    
  // Initialize cumulative P & PE
  Forc_t = (double *) malloc(2*sizeof(double));
  (void) interpolate_timedata(& Forcings,0.,&Forc_t);
  CumP = Forc_t[0];
  CumPE = Forc_t[1]; 
   
  SMA0 = 0.;
  SMA1 = 0.;
  SMA2 = 0.;
  sum_dt_month = 0;
  QM_Gronde = 0.;
  QM_Provence = 0.;
}

    Final event (assemble movies)

    event stop (t = tfin)
{  
  system ("for f in seepage-*.png; do convert $f ppm:- && rm -f $f; done | "
	  "ppm2mp4 seepage.mp4");
  system ("for f in gw_levels-*.png; do convert $f ppm:- && rm -f $f; done | "
	  "ppm2mp4 gw_levels.mp4");
  fprintf(stderr,"Done.\n");
  return 1;
}

    Time integration

    event integration (i++)
{
  scalar S1back[],S2back[];
  scalar r1back[],r2back[];
  
  dt = dtnext(2.);  // fixed time step = 2 days
    
  (void) update_geom_stack();  
  (void) update_TS(zbottom1,ztop1,h1,T1,S1,isConfined1,hasThickness1,seep1);
  (void) update_TS(zbottom2,ztop2,h2,T2,S2,isConfined2,hasThickness2,seep2);
  (void) apply_head_conditions();

  (void) interpolate_timedata(& Forcings,(t+dt)/oneday,&Forc_t);
  P_dt = Forc_t[0]-CumP;
  PE_dt = Forc_t[1]-CumPE;  
  CumP = Forc_t[0];
  CumPE = Forc_t[1];

    The soil moisture accounting (SMA) procedure is taken from the GR4J model. It is computed in the three zones defined, namely over impervious marls (0), over outcropping Bajocien limestone (1), and over outcropping mid-Bathonian limestone (2):

      (void) SMA_GR(& SMA0, maxStorage0,P_dt,PE_dt,(dt/oneday),false,& Excess0); 
  double excess_rate0 = (1.e-3)*Excess0/dt;

  (void) SMA_GR(& SMA1, zbottom1.maxStorage,P_dt,PE_dt,(dt/oneday),false,& Excess1); 
  double excess_rate1 = (1.e-3)*Excess1/dt;

  (void) SMA_GR(& SMA2, zbottom2.maxStorage,P_dt,PE_dt,(dt/oneday),false,& Excess2); 
  double excess_rate2 = (1.e-3)*Excess2/dt;
      
  foreach()
  {
    r1[] = recharge_layer[]==1. ? excess_rate1 : 0.; 
    r2[] = recharge_layer[]==2. ? excess_rate2 : 0.; 
    beta1[] = 0.;
    beta2[] = 0.;
    S1back[] = S1[];
    r1back[] = r1[];
    S2back[] = S2[];
    r2back[] = r2[];
  }

    Then we solve the diffusion equation over the timestep. Note that the r, \beta and \theta fields will be modified by the solver.

      mgd = diffusion(h1,dt,T1,r1,beta1,S1);
  mgd = diffusion(h2,dt,T2,r2,beta2,S2);
  
  foreach()
  {
      S1[] = S1back[];
      r1[] = r1back[];
      S2[] = S2back[];
      r2[] = r2back[];
  }
    
  (void) update_seepage();
  
  // Sum seepage fluxes in catchments to estimate discharge
    
  double Q_Gronde = 0.;
  double Q_Provence = 0.;
    
  foreach(reduction(+:Q_Gronde) reduction(+:Q_Provence))
  {
     coord P = (coord){x,y};  
     int inGronde = isInPolygon(P, & BV_Gronde );
     int inProvence = isInPolygon(P, & BV_Provence );
 
     if(inGronde>0.)
     {
        if(seep1[]>0.) Q_Gronde += seep1[]*sq(Delta)*oneday/86400.;
        if(seep2[]>0.) Q_Gronde += seep2[]*sq(Delta)*oneday/86400.;
        if(recharge_layer[]==0.) Q_Gronde += excess_rate0*sq(Delta)/86400.;   
     }       

     if(inProvence>0.)
     {
        if(seep1[]>0.) Q_Provence += seep1[]*sq(Delta)*oneday/86400.;
        if(seep2[]>0.) Q_Provence += seep2[]*sq(Delta)*oneday/86400.;
        if(recharge_layer[]==0.) Q_Provence += excess_rate0*sq(Delta)/86400.;   
     }
  }

  QM_Gronde += Q_Gronde * dt/oneday;
  QM_Provence += Q_Provence * dt/oneday;
  sum_dt_month += dt/oneday;

  // Groundwater head at boreholes

  for(int kp=0;kp<npiezo;kp++)
  {
    ‘Point point’ is obsolete, use ‘foreach_point/region’ instead 
         Point point = locate ( PIEZO[kp].Site.x, PIEZO[kp].Site.y );
     if (point.level < 0 ) // borehole is not in the process subdomain
       HPIEZO[kp] = -9999.;
     else
     {
        if( PIEZO[kp].layer == 1 ) HPIEZO[kp] =  interpolate_linear(point,h1,PIEZO[kp].Site.x, PIEZO[kp].Site.y);   
        if( PIEZO[kp].layer == 2 ) HPIEZO[kp] =  interpolate_linear(point,h2,PIEZO[kp].Site.x, PIEZO[kp].Site.y);   
      }
  }
    
@if _MPI
    MPI_Allreduce (MPI_IN_PLACE, HPIEZO, npiezo, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
@endif    

  if (pid() == 0) {
    fprintf(resfile,"%f,%.2f,%.2f,%.1f,%.1f",t+oneday*datenum0,P_dt*oneday/dt,PE_dt*oneday/dt,1000.*Q_Gronde,1000.*Q_Provence);  
  
    for(int kp=0;kp<npiezo;kp++)
      fprintf(resfile,",%.2f",HPIEZO[kp]);
    fprintf(resfile,"\n");
    fflush(resfile);
  }

  boundary(all);
}

    Grid adaptivity settings

    int adapt() {

#if TREE
  scalar h1w[],h2w[],maxseep[],indic_MPB[],indic_cote[];
  
  foreach()
  {
    h1w[] = (isDirichlet[]>0.) | (isNeumann[]>0.) ? ZMER : h1[];
    h2w[] = (isDirichlet[]>0.) | (isNeumann[]>0.) ? ZMER : h2[];
    maxseep[] = (isDirichlet[]>0.) | (isNeumann[]>0.) ? 0. : fmax(seep1[],seep2[]);
    indic_MPB[] = (zb[]>=marnesPB_mur[]) && (zb[]<=bathonien_mur[]) ? 1. : 0.;
    indic_cote[] = zb[]>=ZMER ? 1. : 0.;
  }
  
  boundary ({h1w,h2w,maxseep});
  double SE = 5.e-4;
  double IE = 0.1;
    
  astats s = adapt_wavelet ({h1w,h2w,maxseep,indic_MPB,indic_cote}, (double[]){HE,HE,SE,IE,IE},
			    LEVEL, MINLEVEL);

  // Updating isNeumann[] and isDirichlet[] to newly adapted grid
  (void) update_BC_cells();

  boundary(all);
    
//  fprintf (ferr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
  return s.nf;
#else // Cartesian
  return 0;
#endif
}

event do_adapt(t+=8.)
{
  (void) adapt();
}

    Movie output

    event movie (t+=2.)
{
  scalar seepmsk1[],seepmsk2[],msk[],l[];
  scalar h1draw[], h2draw[];
  scalar zbpos[],colorcode[],seasurface[],msksea[];
  
  foreach()
  {
//    msk[] = (isDirichlet[]>0.) || (isNeumann[]>0.) || (zb[]<-4.) ? -1. : 1.;
    msk[] = (isDirichlet[]>0.) || (isNeumann[]>0.) ? -1. : 1.;
    msk1[] = (msk[] < 0.) || (hasThickness1[] < 1.) ? -1. : 1. ;
    msk2[] = (msk[] < 0.) || (hasThickness2[] < 1.) ? -1. : 1. ;
    seepmsk1[] = 0.5*(1.+msk1[]) * seep1[] ;
    seepmsk2[] = 0.5*(1.+msk2[]) * seep2[] ;
    l[] = level;
    h1draw[] = h1[];
    h2draw[] = msk2[] > 0 ? fmax(h1[]+0.3,h2[]) : h2[];
    seasurface[] = ZMER;    
    msksea[] = (zb[] >= (ZMER+2.)) ||  (msk[]<0.) ? -1. : 1.; 

    colorcode[] = 0.3;
    if(marnesPB_mur[]<(zb[]-3.))colorcode[] = 0.5;
    if(bathonien_mur[]<zb[])colorcode[] = 0.7;
    if(seep1[]>0. && zb[]>= (ZMER-5.)) colorcode[] = 0.1;
    if(seep2[]>0. && zb[]>= (ZMER-5.)) colorcode[] = 0.9;
  }

  boundary({zbpos});
    
  char str[100];
  int year, month;
  (void) hack_date(t/oneday,& year, & month);

  float qview_x[4],qview_z[4],qview[4];
  (void) gl_axis_to_quat ((float[]){1,0,0}, 0.36*PI, qview_x); // angle w.r.t vertical
  (void) gl_axis_to_quat ((float[]){0,0,1}, 0.7*PI, qview_z); // rotation in horizontal plane
  gl_add_quats(qview_x, qview_z, qview);

  view (fov = 11., quat = {qview[0],qview[1],qview[2],qview[3]},
        sx = 1., sy = 1., sz = 25.,
        width = 1200, height = 768);

  translate(x = -X0-L0/2.+0.05*L0,y = -Y0-L0/2.+0.12*L0){
    masked_squares_checkerboard("colorcode", linear = true, z = "zb", mask = msk, min = 0., max = 1., map = topocolormap);
    surf_cells (zb,mask = msk);
    masked_squares("seasurface", linear = true, z = "seasurface", mask = msksea, min = 0., max = 1., map = uniform_gray);
    surf_cells (seasurface,mask = msksea);

    draw_polygon_on_surface ( & BV_Gronde, zb, lc = {0,1,0}, lw = 3.5, z_offset = 2.5);
    draw_polygon_on_surface ( & BV_Provence, zb, lc = {0,1,0}, lw = 3.5, z_offset = 2.5);

    sprintf (str, "%4.4d-%2.2d",year,month);
    draw_string (str, 0, size = 32, lc = {0,0,0}, lw = 4); // bottom-left

    sprintf (str, "                                             Blue:  Bajocian limestone w/ semi-confined aquifer");
    draw_string (str, 1, size = 98, lc = {0,0,1}, lw = 1);  // from top-left

    sprintf (str, "\n                                             White: Lower-Bathonian marls (impervious layer)");
    draw_string (str, 1, size = 98, lc = {0,0,0}, lw = 1);  // from top-left

    sprintf (str, "\n\n                                             Red:   Mid-Bathonian limestone w/ unconfined aquifer");
    draw_string (str, 1, size = 98, lc = {1,0,0}, lw = 1);  // from top-left

    sprintf (str, "\n\n\n\n\n   Satured color = seepage from underlying aquifer");
    draw_string (str, 2, size = 120, lc = {0,0,0}, lw = 1);  // from top-left

    sprintf (str, " Green outlines:");
    draw_string (str, 1, size = 98, lc = {0,1,0}, lw = 2); // top-left
    sprintf (str, "\n Gronde & Provence catchments");
    draw_string (str, 1, size = 98, lc = {0,1,0}, lw = 2); // top-left
   
    sprintf (str, "Z-stretch x 25 ");
    draw_string (str, 3, size = 64, lc = {0,0,0}, lw = 1); // bottom-right
  }
  char fname[200];
  sprintf(fname,"seepage-%4.4d.png",(int)t);  
  save (fname);
  clear();


  view (fov = 11., quat = {qview[0],qview[1],qview[2],qview[3]},
        sx = 1., sy = 1., sz = 25.,
        width = 1200, height = 768);

  translate(x = -X0-L0/2.+0.05*L0,y = -Y0-L0/2.+0.12*L0){
    NCCLASS = 7;
    masked_squares("h1", linear = true, z = "h1draw", mask = msk1, min = 0., max = 70., map = discrete_blues);
    surf_cells (h1draw,mask = msk1,lc={0.,0.,0.2});
    masked_squares("h2", linear = true, z = "h2draw", mask = msk2, min = 0., max = 70., map = discrete_reds);
    surf_cells (h2draw,mask = msk2,lc={0.2,0.,0.});

    sprintf (str, "%4.4d-%2.2d",year,month);
    draw_string (str, 0, size = 32, lc = {0,0,0}, lw = 4); // bottom-left

    sprintf (str, "                                           Blue: Head in Bajocian semi-confined aquifer");
    draw_string (str, 1, size = 90, lc = {0,0,1}, lw = 1);  // from top-left

    sprintf (str, "\n                                           Red:  Head in mid-Bathonian unconfined aquifer");
    draw_string (str, 1, size = 90, lc = {1,0,0}, lw = 1);  // from top-left

    sprintf (str, "Z-stretch x 25 ");
    draw_string (str, 3, size = 64, lc = {0,0,0}, lw = 1); 
  }
  sprintf(fname,"gw_levels-%4.4d.png",(int)t);  
  save (fname);
  clear();

}

     Function to compute Gronde & Provence average monthly flow ‘on the fly’

    This would be painfull with gnuplot. The function also writes observed data in face of each simulated value for years 1965 & 1966.

    event monthly_discharge(t=DAYS_PER_MONTH;t+=DAYS_PER_MONTH)
{
  // Monthly flow observations for 1965 & 1966, in L/s (see )
  double Gronde_obs[24] = {580.,689.,113.,91.,44.,14.,2.9,0.5,1.7,1.6,25.,408.,455.,178.,105.,85.,89.,32.,11.,3.,0.7,90.,201.,224.};
  double Provence_obs[24] = {64.,68.,24.,32.,34.,31.,20.,21.,29.,23.,27.,79.,83.,73.,57.,53.,53.,37.,33.,31.,20.,38.,68.,86};
  double S_Gronde = 20.30; // in km2
  double S_Provence = 9.48; // in km2

  int year,month;
  QM_Gronde *= 1000./sum_dt_month; // average in L/s
  QM_Provence *= 1000./sum_dt_month; // average in L/s

  hack_date(t/oneday-DAYS_PER_MONTH/2.,& year, & month);
  if( pid()==0 && year>=1965 && year<=1966)
  {
    int ix = 12*(year-1965)+month-1;
    fprintf(flowfile,"%d %d %g %g %g %g %g %g %g %g\n",year,month,
      QM_Gronde,Gronde_obs[ix],
      QM_Provence,Provence_obs[ix],
      DAYS_PER_MONTH * 0.0864 * QM_Gronde/S_Gronde,DAYS_PER_MONTH * 0.0864 * Gronde_obs[ix]/S_Gronde,
      DAYS_PER_MONTH * 0.0864 * QM_Provence/S_Provence,DAYS_PER_MONTH * 0.0864 * Provence_obs[ix]/S_Provence);
    fflush(flowfile);
  }

  sum_dt_month = 0.;
  QM_Gronde = 0.;
  QM_Provence = 0.;
}

    Animation of the solution (simulation timestep: 2 days)

    Saturated areas

    Groundwater head in both aquifers. Filled contours have 10-meter interval. The head in both aquifers is initially set to the topographic level, so that a warm-up period is needed in order for it to drop to a deeper, more realistic level