sandbox/nlemoine/ice/halfar1D.c

    Return to my homepage

    1D test case for the Shallow Ice Approximation diffusive equation:

    \displaystyle \frac{\partial\eta}{\partial t} + \nabla\cdot(h\mathbf{U}) = M where \eta is the ice surface elevation, z_b the elevation of the substratum, h = \eta - z_b the ice thickness, \mathbf{U} the depth-average velocity and M the local mass balance (difference between accumulation and ablation).

    In the SIA, the flux density \mathbf{q} = h\mathbf{U} boils down to a diffusive flux with a nonlinear diffusivity D which depends on ice thickness and ice surface slope: \displaystyle \mathbf{q} = h\mathbf{U} = -\underbrace{\textstyle\frac{2A (\rho g)^n}{n+2} h^{n+2}\big\|\nabla \eta\big\|^{n-1}}_{D} \nabla\eta n is the exponent of Glen flow law, n\approx 3

    Halfar 1D similarity solution (Halfar, 1981) holds for a flat substratum (z_b = 0) and a zero mass balance (M=0): \displaystyle \eta(x,t)=h(x,t)=\sqrt{aV}f(t)\,g\!\left(\frac{x f(t)}{\sqrt{aV}}\right)

    //#include "grid/multigrid1D.h"
#include "grid/bitree.h"
#include "run.h"
#include "diffusion.h"

#define LEVEL         8
#define MINLEVEL      6
#define ETAE          20 // error on surface elevation
#define DE            1e7 // error on diffusivity (m2/yr)
#define sec_per_year  31557600.
#define tfin          100000.
#define nGlen	      3.
#define AGlen	      5.0e-24

scalar zb[],eta[];
face vector D[];
double Volume,aV;
double t0,Gamma;

double dt, tini;
mgstats mgd;

double Beta(double xx, double yy){
 return exp(lgamma(xx)+lgamma(yy)-lgamma(xx+yy));
}

    Diffusivity function

    int update_diffusivity(scalar zb, scalar eta, face vector D)
{
  foreach_face()
  {
    double hm = (eta[]+eta[-1]-zb[]-zb[-1])/2.; // ice thickness interpolated at face center
    hm = max(hm,0.);
    double sn = (eta[]-eta[-1])/Delta; 				// slope component along face normal
    D.x[] = sec_per_year*Gamma*pow(hm,nGlen+2.)*pow(fabs(sn),nGlen-1.);
  }

 return(0);
}

    Halfar 1D similarity solution (Halfar, 1981)

    In Halfar’s paper the solution is scaled using total volume V multiplied by a constant a. Noting W(t) the half-width at time t, the explanation is as follows: \displaystyle \begin{aligned} V(t) & = 2\int_0^{W(t)=\frac{\sqrt{aV}}{f(t)}}\sqrt{aV}f(t)\,g\!\left(\frac{x}{W(t)}\right)dx\\ & = 2(aV)\int_0^1 g(\xi)d\xi \end{aligned} Moreover, \displaystyle \int_0^1 g(\xi)d\xi = \int_0^1 \left(1-\xi^{\frac{n+1}{n}}\right)^{\frac{n}{2n+1}}d\xi = \frac{n}{n+1}\int_0^1 \left(1-\xi'\right)^{\frac{n}{2n+1}}\xi'^{\frac{n}{n+1}-1}d\xi' The latter integral is the Beta function \Beta\!\left(\frac{3n+1}{2n+1},\frac{n}{n+1}\right), so we find the condition for a: \displaystyle V=\frac{2n}{n+1}(aV)\ \Beta\!\left(\frac{3n+1}{2n+1},\frac{n}{n+1}\right) hence \displaystyle a=\frac{n+1}{2n\ \Beta\!\left(\frac{3n+1}{2n+1},\frac{n}{n+1}\right)}

    double Halfar1D(double xx, double tt)
{
    double ft = pow(t0/tt,1./(3.*nGlen+2.)); // eq. (12)
    double argg = fabs(xx*ft/sqrt(aV));	// eq. (10)
    double res = argg < 1. ? sqrt(aV)*ft*pow( (1.0 - pow(argg,1.+1./nGlen)) , nGlen/(2.*nGlen+1.) ) : 0. ; // eq. (11)
    return res  ;
}

    Main

    int main (int argc, char * argv[])
{
  N = 1 << LEVEL;
  L0 = 2.0e6;
  size (L0);
  origin (-L0/2.,-L0/2.);

  Volume = 2.0e9;  // Total volume (m2)

  double bet = Beta((3.*nGlen+1.)/(2.*nGlen+1.),nGlen/(nGlen+1.));
  aV = Volume*(nGlen+1.)/(2.*nGlen*bet);
  Gamma  = 2. * AGlen * pow(9.81*910.,nGlen)/(nGlen+2.) ;
  t0 = pow((2.*nGlen+1.)/(nGlen+1.),nGlen)*pow(aV,-nGlen/2.)/(3.*nGlen+2.); // eq. (13)
  t0 = t0/Gamma/sec_per_year;

  tini = 1.;

  run();
  return(0);
}

event init (i=0)
{
  foreach()
  {
    zb[] = 0.;
    eta[] = Halfar1D(x,tini);
  }
  boundary({zb,eta});
}

event profile (t=0. ; t+=100.)
{
  FILE * fp;
  char name[80];
  double xp,eta_num,eta_ana;

  sprintf(name,"transect-%d.dat",(int)t);
  fp = fopen(name,"w");

  scalar Dc[];
  foreach()
     Dc[] = (D.x[]+D.x[1])/2.;

  for(int np = 0;np<N;np++)
  {
    xp = -L0/2.+(L0/N)*((float)np+0.5);
    eta_num = interpolate (eta, xp);
    eta_ana = Halfar1D(xp,t+tini);
    double Dp = interpolate (Dc, xp);
    Point point = locate(xp);
    fprintf (fp, "%g %g %g %g %d\n",xp/1000.,eta_num,eta_ana,Dp,point.level);
  }
  
  fclose(fp);
}

event stop (t = tfin);

    Time integration

    The SIA is a basically a diffusion equation, so we use the Poisson solver for integration. Since this solver is time-implicit, there isn’t any stability problem; however, since the diffusivity depends on the solution, setting a too large time step will cause spurious features to develop in the solution if the diffusivity field is not updated often enough. So we set the time step to a ‘’reasonable’’ (up to 10) multiple of the explicit time step \displaystyle \Delta t_\textrm{explicit} = \frac{1}{2}\frac{\Delta x^2}{D_\textrm{max}}

    event integration (i++)
{
  (void) update_diffusivity(zb,eta,D);

  stats s = statsf (D.x);
  double dtExplicit = 0.5*L0*L0/N/N/s.max;
  dt = dtnext(4.0*dtExplicit);

  mgd = diffusion(eta,dt,D);

  // preserve positivity of ice tickness
  foreach()
   eta[] = max(zb[],eta[]);
}

// Adaptivity

int adapt() {
#if TREE
  scalar Dc[];
  foreach()
     Dc[] = (D.x[]+D.x[1])/2.;

//  astats s = adapt_wavelet ({eta,Dc}, (double[]){ETAE,DE},
  astats s = adapt_wavelet ({eta}, (double[]){ETAE},
			    LEVEL, MINLEVEL);
//  fprintf (ferr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
  return s.nf;
#else // Cartesian
  return 0;
#endif
}

// And finally we set the event that will apply our *adapt()* function at every timestep.

event do_adapt (i++) adapt();
       reset
   set xlabel "x (km)"
   set ylabel "Ice thickness (m)"
   ydim = 800 
   xdim = 480
   plot 'transect-0.dat' u 1:3 w lines lc rgb "black" lw 3 title 't = 0', \
        'transect-100.dat' u 1:3 w lines lc rgb "red" title 't = 0.1 ka (analytical)', \
        'transect-100.dat' u 1:2 w point pointtype 6 ps 1 lc rgb "red" title 't = 0.1 ka (numerical)', \
        'transect-1000.dat' u 1:3 w lines lc rgb "yellow" title 't = 1 ka (analytical)', \
        'transect-1000.dat' u 1:2 w point pointtype 6 ps 1 lc rgb "yellow" title 't = 1 ka (numerical)', \
        'transect-10000.dat' u 1:3 w lines lc rgb "green" title 't = 10 ka (analytical)', \
        'transect-10000.dat' u 1:2 w point pointtype 6 ps 1 lc rgb "green" title 't = 10 ka (numerical)', \
        'transect-100000.dat' u 1:3 w lines lc rgb "blue" title 't = 100 ka (analytical)', \
        'transect-100000.dat' u 1:2 w point pointtype 6 ps 1 lc rgb "blue" title 't = 100 ka (numerical)'
    (script)

    (script)