/** # 1D test case for the Shallow Ice Approximation diffusive equation: $$\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: $$\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](https://doi.org/10.1029/JC086iC11p11065)) holds for a flat substratum ($z_b = 0$) and a zero mass balance ($M=0$): $$\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](https://doi.org/10.1029/JC086iC11p11065) 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: \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, $$\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](https://en.wikipedia.org/wiki/Beta_function) $\Beta\!\left(\frac{3n+1}{2n+1},\frac{n}{n+1}\right)$, so we find the condition for $a$: $$V=\frac{2n}{n+1}(aV)\ \Beta\!\left(\frac{3n+1}{2n+1},\frac{n}{n+1}\right)$$ hence $$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