# 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[]);
}

#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)'