<<

. 13
( 14)



>>

setupgridneighbors();
topo=matrix(1,lattice_size_x,1,lattice_size_y);
toposave=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
depth=matrix(1,lattice_size_x,1,lattice_size_y);
slope=matrix(1,lattice_size_x,1,lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fscanf(fp1,"%f",&topo[i][j]);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{toposave[i][j]=topo[i][j];
calculatechannelslope(i,j);
if (slope[i][j]<0.01) slope[i][j]=0.01;}
for (k=1;k<=niterations;k++)
{printf("iteration = %d/%d\n",k,niterations);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{flow[i][j]=pow(runoff,1.66667)*sqrt(slope[i][j])*deltax/manningsn;
topovec[(j-1)*lattice_size_x+i]=topo[i][j];}
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
A2.3 SUCCESSIVE FLOW ROUTING WITH MFD METHOD 241


if (i==lattice_size_x) j--;
mfdflowroute(i,j);
depth[i][j]=pow(flow[i][j]*manningsn/
(deltax*sqrt(slope[i][j])),0.6);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
if ((toposave[i][j]+depth[i][j])>topo[i][j])
topo[i][j]=toposave[i][j]+depth[i][j]/niterations;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]-toposave[i][j]);
}
Appendix 3




Codes for solving the advection equation

A3.1 Coupled 1D bedrock-alluvial channel evolution

The following code models the coupled 1D evolution of a bedrock-alluvial channel longitudinal
pro¬le described in Section 4.5.


main()
{ FILE *fp1,*fp2;
float transport,c,D,*hold,time,factor,max,totsed,deltah,*h,sum;
int bedrock,xc,i,lattice_size,check;

fp1=fopen("appalachiansspace1","w");
fp2=fopen("appalachianstime1","w");
c=0.1;
D=1.0;
lattice_size=128; /* delta= 1 km, so basin is 128 km in length */
bedrock=16; /* bedrock-alluvial transition 16 km from divide */
h=vector(1,lattice_size);
hold=vector(1,lattice_size);
xc=bedrock;
for (i=1;i<=lattice_size;i++)
if (i>bedrock)
{h[i]=0.0;
hold[i]=0.0;
/* plot the initial condition first */
fprintf(fp1,"%f %f\n",i/(float)(lattice_size),h[i]);}
else {h[i]=1.0;
hold[i]=1.0;
fprintf(fp1,"%f %f\n",i/(float)(lattice_size),h[i]);}
factor=1.0;
time=0.0;
check=10000;
while (time<100000.0)
{if (time>check)
{sum=0;
A3.2 MODELING THE DEVELOPMENT OF TOPOGRAPHIC STEADY STATE IN THE STREAM-POWER MODEL 243


check+=10000;
for (i=1;i<=lattice_size;i++)
{sum+=h[i];
fprintf(fp1,"%f %f\n",i/(float)(lattice_size),h[i]);}
fprintf(fp2,"%f %f\n",time,sum/lattice_size);
fprintf(fp1,"%f %f\n",xc/(float)(lattice_size),h[xc]);}
max=0;
totsed=0;
deltah=c*factor*(1/(float)(lattice_size))*(h[1]-h[2]);
if (fabs(deltah)>max) max=fabs(deltah);
h[1]-=deltah;
totsed=deltah;
for (i=2;i<=lattice_size-1;i++)
{deltah=c*factor*(i/(float)(lattice_size))*(hold[i]-hold[i+1]);
totsed+=deltah;
transport=D*factor*(i/(float)(lattice_size))*(hold[i]-hold[i+1]);
if ((transport>(totsed+0.01))&&(i<=bedrock))
{xc=i;
if (fabs(deltah)>max) max=fabs(deltah);
h[i]-=deltah;}
else
{deltah=D*factor*((i-1)/(float)(lattice_size))*(hold[i-1]-hold[i])-
D*factor*(i/(float)(lattice_size))*(hold[i]-hold[i+1]);
if (fabs(deltah)>max) max=fabs(deltah);
h[i]+=deltah;}}
if (max<0.001)
{for (i=1;i<=lattice_size;i++)
hold[i]=h[i];
time+=factor;}
else
{h[i]=hold[i];
factor/=3;}
if (max<0.0001) factor*=3;}
fclose(fp1);
fclose(fp2);
}



A3.2 Modeling the development of topographic steady state in the
stream-power model

The following code implements the stream power model on a landscape with m = 1/2 and n = 1 for a
uniformly uplifting square mountain block. In order to apply the stream-power model, we must create
some initial landscape. In this example, the initial landscape is created by ¬rst generating a random
white noise with very minor relief and then diffusing that landscape as it is uplifted by 1 m. This
process creates a low-relief landscape that drains towards the edges but also has microtopographic
variability that establishes initial drainage pathways in an unstructured manner (Figure A3.1a). After
the initial landscape is created, the block is uplifted at a rate U = 1 m/kyr and the stream-power
244 CODES FOR SOLVING THE ADVECTION EQUATION



(a) (b)




steady state
initial uplift phase

Fig A3.1 Output of the code in Section A3.2 that
implements the stream-power model in a uniformly uplifting
mountain block. (a) Early in the simulation, drainages grow
headward into a diffusing plateau with initially random
topography. (b) After development of steady state.



model with K w = 1 kyr’1 is imposed using upwind differencing. On hillslopes, a simple avalanching
rule is used: any slope above 30—¦ sheds mass to restore the 30—¦ angle of stability. The 30—¦ stability
criterion is imposed in routine avalanche as a threshold relief between adjacent pixels.

#define sqrt2 1.414213562373
#define oneoversqrt2 0.707106781186
#define fillincrement 0.001

float **flow1,**flow2,**flow3,**flow4,**flow5,**flow6,**flow7,**flow8,**flow;
float **topo,**topoold,**topo2,**slope,deltax,*ax,*ay,*bx,*by,*cx,*cy,*ux,*uy;
float *rx,*ry,U,K,D,**slope,timestep,*topovec,thresh;
int *topovecind,lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void hillslopediffusioninit()
{ int i,j,count;
float term1;

ax=vector(1,lattice_size_x);
ay=vector(1,lattice_size_y);
bx=vector(1,lattice_size_x);
by=vector(1,lattice_size_y);
cx=vector(1,lattice_size_x);
cy=vector(1,lattice_size_y);
ux=vector(1,lattice_size_x);
uy=vector(1,lattice_size_y);
A3.2 MODELING THE DEVELOPMENT OF TOPOGRAPHIC STEADY STATE IN THE STREAM-POWER MODEL 245


rx=vector(1,lattice_size_x);
ry=vector(1,lattice_size_y);
D=10000000.0;
count=0;
term1=D/(deltax*deltax);
for (i=1;i<=lattice_size_x;i++)
{ax[i]=-term1;
cx[i]=-term1;
bx[i]=4*term1+1;
if (i==1)
{bx[i]=1;
cx[i]=0;}
if (i==lattice_size_x)
{bx[i]=1;
ax[i]=0;}}
for (j=1;j<=lattice_size_y;j++)
{ay[j]=-term1;
cy[j]=-term1;
by[j]=4*term1+1;
if (j==1)
{by[j]=1;
cy[j]=0;}
if (j==lattice_size_y)
{by[j]=1;
ay[j]=0;}}
while (count<5)
{count++;
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
topo2[i][j]=topo[i][j];
for (i=1;i<=lattice_size_x;i++)
{for (j=1;j<=lattice_size_y;j++)
{ry[j]=term1*(topo[iup[i]][j]+topo[idown[i]][j])+topoold[i][j];
if (j==1)
ry[j]=topoold[i][j];
if (j==lattice_size_y)
ry[j]=topoold[i][j];}
tridag(ay,by,cy,ry,uy,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=uy[j];}
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
topo2[i][j]=topo[i][j];
for (j=1;j<=lattice_size_y;j++)
{for (i=1;i<=lattice_size_x;i++)
{rx[i]=term1*(topo[i][jup[j]]+topo[i][jdown[j]])+topoold[i][j];
if (i==1)
rx[i]=topoold[i][j];
if (i==lattice_size_x)
246 CODES FOR SOLVING THE ADVECTION EQUATION


rx[i]=topoold[i][j];}
tridag(ax,bx,cx,rx,ux,lattice_size_x);
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=ux[i];}}
}

void avalanche(i,j)
int i,j;
{
if (topo[iup[i]][j]-topo[i][j]>thresh)
topo[iup[i]][j]=topo[i][j]+thresh;
if (topo[idown[i]][j]-topo[i][j]>thresh)
topo[idown[i]][j]=topo[i][j]+thresh;
if (topo[i][jup[j]]-topo[i][j]>thresh)
topo[i][jup[j]]=topo[i][j]+thresh;
if (topo[i][jdown[j]]-topo[i][j]>thresh)
topo[i][jdown[j]]=topo[i][j]+thresh;
if (topo[iup[i]][jup[j]]-topo[i][j]>(thresh*sqrt2))
topo[iup[i]][jup[j]]=topo[i][j]+thresh*sqrt2;
if (topo[iup[i]][jdown[j]]-topo[i][j]>(thresh*sqrt2))
topo[iup[i]][jdown[j]]=topo[i][j]+thresh*sqrt2;
if (topo[idown[i]][jup[j]]-topo[i][j]>(thresh*sqrt2))
topo[idown[i]][jup[j]]=topo[i][j]+thresh*sqrt2;
if (topo[idown[i]][jdown[j]]-topo[i][j]>(thresh*sqrt2))
topo[idown[i]][jdown[j]]=topo[i][j]+thresh*sqrt2;
}

main()
{ FILE *fp1;
float deltah,time,max,duration;
int printinterval,idum,i,j,t,step;

fp1=fopen("streampowertopo","w");
lattice_size_x=250;
lattice_size_y=250;
U=1; /* m/kyr */
K=0.05; /* kyr^-1 */
printinterval=100;
deltax=200.0; /* m */
thresh=0.58*deltax; /* 30 deg */
timestep=1; /* kyr */
duration=200;
setupgridneighbors();
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topo2=matrix(1,lattice_size_x,1,lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
slope= matrix(1,lattice_size_x,1,lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
A3.2 MODELING THE DEVELOPMENT OF TOPOGRAPHIC STEADY STATE IN THE STREAM-POWER MODEL 247


flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
{topo[i][j]=0.5*gasdev(&idum);
topoold[i][j]=topo[i][j];
flow[i][j]=1;}
/*construct diffusional landscape for initial flow routing */
for (step=1;step<=10;step++)
{hillslopediffusioninit();
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{topo[i][j]+=0.1;
topoold[i][j]+=0.1;}}
time=0;
while (time<duration)
{/*perform landsliding*/
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=0;
while (t<lattice_size_x*lattice_size_y)
{t++;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
avalanche(i,j);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{flow[i][j]=1;
topovec[(j-1)*lattice_size_x+i]=topo[i][j];}
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
248 CODES FOR SOLVING THE ADVECTION EQUATION


{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{topo[i][j]+=U*timestep;
topoold[i][j]+=U*timestep;}
/*perform upwind erosion*/
max=0;
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{calculatealongchannelslope(i,j);
deltah=timestep*K*sqrt(flow[i][j])*deltax*slope[i][j];
topo[i][j]-=deltah;
if (topo[i][j]<0) topo[i][j]=0;
if (K*sqrt(flow[i][j])>max) max=K*sqrt(flow[i][j]);}
time+=timestep;
if (max>0.3*deltax/timestep)
{time-=timestep;
timestep/=2.0;
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
topo[i][j]=topoold[i][j]-U*timestep;}
else
{if (max<0.03*deltax/timestep) timestep*=1.2;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];}
if (time>printinterval)
{printinterval+=250;
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
{fprintf(fp1,"%f\n",topo[i][j]);}}}
fclose(fp1);
}




A3.3 Knickpoint propagation in the 2D sediment-¬‚ux-driven bedrock
erosion model

The following code solves for the ¬‚uvial bedrock incision of an uplifted square block according to the
sediment-¬‚ux-driven bedrock erosion model. In the model, tectonic uplift at a rate of U = 1 m/kyr
occurs for 1 Myr. After that initial period, isostatic rebound is calculated using the 2D Fourier-¬ltering
A3.3 KNICKPOINT PROPAGATION IN THE 2D SEDIMENT-FLUX-DRIVEN BEDROCK EROSION MODEL 249




Fig A3.2 Grayscale maps of example output of the code in
Section A3.3. This code implements the sediment-¬‚ux driven
model in a mountain block that is uniformly uplifted for 1 Myr
and then undergoes erosional decay with isostatic rebound.
(a) Topography, and instantaneous (b) erosion rate, (c)
sediment ¬‚ux, and (d) ¬‚exural rebound at the end of a 40 Myr
simulation.



method (see Chapter 5). The code accepts an input DEM named sedfluxdriventopoinitial (available
from the CUP website). The initial topography is divided by a factor of ¬ve to create a low-relief
surface with a drainage network geometry already formed. Figure A3.2 illustrates example output
from the model.
250 CODES FOR SOLVING THE ADVECTION EQUATION


#define sqrt2 1.414213562373
#define oneoversqrt2 0.707106781186

float delrho,alpha,*load,**slope,**erode,S,averosionrate,avreboundrate,avelevation;
float oneoverdeltax,oneoverdeltax2,**U,K,D,X,duration,timestep;
float *topovec,thresh,**deltah,**channel,**area,**sedflux,**sedfluxold,**flow1;
float **flow2,**flow3,**flow4,**flow5,**flow6,**flow7,**flow8,**flow,**topo,**topoold;
float **topo2,deltax,*ax,*ay,*bx,*by,*cx,*cy,*ux,*uy,*rx,*ry;
int *nn,*topovecind,lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void computeflexure()
{ int i,j,index;
float fact;

for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{load[2*(i-1)*lattice_size_y+2*j-1]=erode[i][j];
load[2*(i-1)*lattice_size_y+2*j]=0.0;}
fourn(load,nn,2,1);
load[1]*=1/delrho;
load[2]*=1/delrho;
for (j=1;j<=lattice_size_y/2;j++)
{fact=1/(delrho+pow(alpha*PI*j/lattice_size_y,4.0));
load[2*j+1]*=fact;
load[2*j+2]*=fact;
load[2*lattice_size_y-2*j+1]*=fact;
load[2*lattice_size_y-2*j+2]*=fact;}
for (i=1;i<=lattice_size_x/2;i++)
{fact=1/(delrho+pow(alpha*PI*i/lattice_size_x,4.0));
load[2*i*lattice_size_y+1]*=fact;
load[2*i*lattice_size_y+2]*=fact;
load[2*lattice_size_x*lattice_size_y-2*i*lattice_size_y+1]*=fact;
load[2*lattice_size_x*lattice_size_y-2*i*lattice_size_y+2]*=fact;}
for (i=1;i<=lattice_size_x/2;i++)
for (j=1;j<=lattice_size_y/2;j++)
{fact=1/(delrho+pow(alpha*sqrt((PI*PI*i*i)/(lattice_size_x*lattice_size_x)
+(PI*PI*j*j)/(lattice_size_y*lattice_size_y)),4.0));
index=2*i*lattice_size_y+2*j+1;
load[index]*=fact;
load[index+1]*=fact;
index=2*i*lattice_size_y+2*lattice_size_y-2*j+1;
load[index]*=fact;
load[index+1]*=fact;
index=2*lattice_size_x*lattice_size_y-2*(i-1)*lattice_size_y-2*(lattice_size_y-j)+1;
load[index]*=fact;
load[index+1]*=fact;
index=2*lattice_size_x*lattice_size_y-2*lattice_size_y-
2*(i-1)*lattice_size_y+2*(lattice_size_y-j)+1;
load[index]*=fact;
A3.3 KNICKPOINT PROPAGATION IN THE 2D SEDIMENT-FLUX-DRIVEN BEDROCK EROSION MODEL 251


load[index+1]*=fact;}
fourn(load,nn,2,-1);
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
{U[i][j]=(load[2*(iup[i]-1)*lattice_size_y+2*j-1]+
load[2*(idown[i]-1)*lattice_size_y+2*j-1]+
load[2*(iup[i]-1)*lattice_size_y+2*jdown[j]-1]+load[2*(i-
1)*lattice_size_y+2*jup[j]-1])/(4*4*lattice_size_x*lattice_size_y);}
}

void setupmatrices()
{ int i,j;

ax=vector(1,lattice_size_x);
ay=vector(1,lattice_size_y);
bx=vector(1,lattice_size_x);
by=vector(1,lattice_size_y);
cx=vector(1,lattice_size_x);
cy=vector(1,lattice_size_y);
ux=vector(1,lattice_size_x);
uy=vector(1,lattice_size_y);
rx=vector(1,lattice_size_x);
ry=vector(1,lattice_size_y);
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topo2=matrix(1,lattice_size_x,1,lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
sedflux=matrix(1,lattice_size_x,1,lattice_size_y);
sedfluxold=matrix(1,lattice_size_x,1,lattice_size_y);
area=matrix(1,lattice_size_x,1,lattice_size_y);
slope=matrix(1,lattice_size_x,1,lattice_size_y);
deltah=matrix(1,lattice_size_x,1,lattice_size_y);
load=vector(1,2*lattice_size_x*lattice_size_y);
erode=matrix(1,lattice_size_x,1,lattice_size_y);
U=matrix(1,lattice_size_x,1,lattice_size_y);
channel=matrix(1,lattice_size_x,1,lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
}

main()
252 CODES FOR SOLVING THE ADVECTION EQUATION


{ FILE *fp0,*fp1,*fp2,*fp3,*fp4,*fp5;
float change,maxelevation,capacity,time,max,hillslopeerosionrate;
int printinterval,i,j,t;

fp0=fopen("sedfluxdriventopoinitial","r");
fp1=fopen("sedfluxdriventopo","w");
fp2=fopen("sedfluxdrivenerosion","w");
fp3=fopen("sedfluxdrivensedflux","w");
fp4=fopen("sedfluxdrivenrebound","w");
lattice_size_x=256;
lattice_size_y=256;
deltax=500; /* m */
delrho=0.28; /* (rho_m-rho_c)/rho_m */
oneoverdeltax=1.0/deltax;
oneoverdeltax2=1.0/(deltax*deltax);
timestep=1.0; /* kyr */
alpha=100000; /* m */
hillslopeerosionrate=0.01; /* m/kyr */
K=0.0005; /* m^1/2/kyr */
D=1.0; /* m^2/kyr */
X=0.005; /* m^-1 */
alpha/=deltax;
setupmatrices();
setupgridneighbors();
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp0,"%f",&topo[i][j]);topo[i][j]/=5;
if ((i==1)||(j==1)||(i==lattice_size_x)||(j==lattice_size_y)) topo[i][j]=0;
topoold[i][j]=topo[i][j];
flow[i][j]=deltax*deltax;
area[i][j]=0;
U[i][j]=0;
erode[i][j]=0;
channel[i][j]=0;
sedflux[i][j]=0;
sedfluxold[i][j]=0;}
time=0;
nn=ivector(1,2);
nn[1]=lattice_size_x;
nn[2]=lattice_size_y;
duration=40000; /* 40 Myr */
printinterval=40000;
while (time<duration)
{for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];
/*compute uplift rate*/
if (time>1000) computeflexure(); /* compute isostatic uplift */
else
A3.3 KNICKPOINT PROPAGATION IN THE 2D SEDIMENT-FLUX-DRIVEN BEDROCK EROSION MODEL 253


{/* or prescribe a tectonic uplift rate */
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
U[i][j]=1.0; /* m/kyr */}
avelevation=0;avreboundrate=0;maxelevation=0;
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
{avelevation+=topo[i][j];
if (topo[i][j]>maxelevation) maxelevation=topo[i][j];
avreboundrate+=U[i][j];
topoold[i][j]+=U[i][j]*timestep;
topo[i][j]+=U[i][j]*timestep;
deltah[i][j]=0;
channel[i][j]=0;}
avelevation/=(lattice_size_x-2)*(lattice_size_y-2);
avreboundrate/=(lattice_size_x-2)*(lattice_size_y-2);
printf("%f %f %f %f %f\n",time,avelevation,averosionrate,avreboundrate,maxelevation);
/* hillslope erosion can occur by prescribing a uniform rate, as done here
or using the ADI technique to solve the diffusion equation on hillslopes */
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
if (channel[i][j]==0) topo[i][j]-=hillslopeerosionrate*timestep;
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
erode[i][j]=(topoold[i][j]-topo[i][j])/timestep;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
/*route water from highest gridpoint to lowest*/
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
flow[i][j]=deltax*deltax;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
area[i][j]=flow[i][j];
/*perform upwind differencing */
max=0;
254 CODES FOR SOLVING THE ADVECTION EQUATION


for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{calculatealongchannelslope(i,j);
capacity=slope[i][j]*sqrt(area[i][j]);
if (capacity>X)
{change=timestep*K*sqrt(fabs(sedflux[i][j]))*deltax*slope[i][j];
deltah[i][j]+=change;
erode[i][j]+=change/timestep;
if (deltah[i][j]<0) deltah[i][j]=0;
channel[i][j]=1;}
topo[i][j]-=deltah[i][j];
if (topo[i][j]<0) topo[i][j]=0;
if (K*sqrt(fabs(sedflux[i][j]))*deltax>max)
max=K*sqrt(fabs(sedflux[i][j]))*deltax;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
averosionrate=0;
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{flow[i][j]=erode[i][j];
averosionrate+=erode[i][j];}
averosionrate/=(lattice_size_x-2)*(lattice_size_y-2);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{sedfluxold[i][j]=sedflux[i][j];
sedflux[i][j]=flow[i][j];}
time+=timestep;
if (max>5.0*deltax/timestep)
{time-=timestep;
timestep/=2.0;
for (i=2;i<=lattice_size_x-1;i++)
for (j=2;j<=lattice_size_y-1;j++)
{sedflux[i][j]=sedfluxold[i][j];
topo[i][j]=topoold[i][j]-U[i][j]*timestep;}}
else
{if (max<0.5*deltax/timestep) timestep*=1.2;}
A3.3 KNICKPOINT PROPAGATION IN THE 2D SEDIMENT-FLUX-DRIVEN BEDROCK EROSION MODEL 255


if (time>=printinterval)
{printinterval+=40000;
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
{fprintf(fp1,"%f\n",topo[i][j]);
fprintf(fp2,"%f\n",erode[i][j]);
fprintf(fp3,"%f\n",sedflux[i][j]);
fprintf(fp4,"%f\n",U[i][j]);}}}
fclose(fp0);
fclose(fp1);
fclose(fp2);
fclose(fp3);
fclose(fp4);
fclose(fp5);
}
Appendix 4




Codes for solving the ¬‚exure equation

A4.1 Fourier ¬ltering in 1D

The following code uses the Numerical Recipes routine realft to Fourier-transform the loading
function input from ¬le load1dandes (two-column format with elevation as the ¬rst column and
distance along the transect as the second column; ¬le available from the CUP website), ¬lter it using
the kernel of the ¬‚exure equation, and inverse-Fourier-transform it back to real space. The pro¬le is
read into the ¬rst half of a vector of length 4096; the second half of the vector is zeroed (to avoid
periodic boundary effects). In this example, the loading function is constructed from a topographic
cross section through the Andes Mountains with a resolution of 0.52 km/pixel. In the code, only
elevation values greater that 1000 m a.s.l. contribute to the load; elevation values lower than 1000 m
are made equal to zero before the Fourier-¬ltering procedure begins.
main()
{ float deltax,alpha,*w,dum,delrho, L,k;
int lattice_size_x,i;
FILE *fp1,*fp2;

fp1=fopen("load1dandes","r");
fp2=fopen("load1ddeflect50","w");
lattice_size_x=4096;
deltax=0.52; /* km */
L=deltax*2*lattice_size_x;
alpha=50.0; /* (D/(rho_c*g))^0.25 */
delrho=0.27; /* (rho_m-rho_c)/rho_m */
w=vector(1,2*lattice_size_x);
for (i=1;i<=2*lattice_size_x;i++)
{if (i<=lattice_size_x/2) fscanf(fp1,"%f %f",&w[i],&dum);
else w[i]=0;
if (fabs(w[i])<1000) w[i]=0.0;}
realft(w,2*lattice_size_x,1);
w[1]=w[1]/delrho;
k=PI/deltax;
w[2]=w[2]/(delrho+pow(k*alpha,4.0));
for (i=3;i<=2*lattice_size_x;i++)
{k=((i-1)/2)*PI/L;
A4.2 INTEGRAL SOLUTION IN 2D 257


w[i]=w[i]/(delrho+pow(k*alpha,4.0));}
realft(w,2*lattice_size_x,-1);
for (i=1;i<=2*lattice_size_x;i++)
fprintf(fp2,"%f %f\n",i*deltax,w[i]/(lattice_size_x));
fclose(fp1);
fclose(fp2);
}


A4.2 Integral solution in 2D

The integral solution method uses point source solutions to the ¬‚exure equation. The point source
solution has the form of the Kelvin function kei(r ). The functions given below provide approximate
expressions for the kei function in terms of series expansions for the ber and bei functions.

#define PI 3.141592653589793

float ber(r)
float r;
{ float x,z;

x = pow(pow(r/8.0,2),2);
z = 1.0 + x * (-64.0 + x * (113.77777774 +
x * (-32.36345652 + x * ( 2.64191397 +
x * ( -0.08349609 + x * ( 0.00122552 +
x * -0.00000901))))));
return z;
}

float bei(r)
float r;
{ float x,y,z;
y = pow(r/8.0,2);
x = pow(y,2);
z=y* (16.0 + x * (-113.77777774 +
x * (72.81777742 + x * ( -10.56765779 +
x * ( 0.52185615 + x * ( -0.01103667 +
x * 0.00011346))))));
return z;
}

float kei(r)
float r;
{ float x,y,z;

if (r == 0) r = 1e-50;
y = pow(r/8.0,2);
x = pow(y,2);
z = (-log(0.50 * r) * bei(r)) - (0.25*PI*ber(r)) +
258 CODES FOR SOLVING THE FLEXURE EQUATION


y * ( 6.76454936 + x * (-142.91827687 +
x * (124.23569650 + x * ( -21.30060904 +
x * ( 1.17509064 + x * ( -0.02695875 +
x * 0.00029532))))));
return z;
}


A4.3 Fourier ¬ltering in 2D

The following code uses the Numerical Recipes routine fourn to Fourier-transform the loading func-
tion input from ¬le load2dandes (in two-column format with elevation as the ¬rst column and
distance as the second column; available from the CUP website), ¬lter it using the kernel of the
¬‚exure equation, and inverse-Fourier-transform it back to real space. In this example, the load-
ing function is constructed from a DEM of the central Andes Mountains with a resolution of
1.1 km/pixel. In the code, only elevation values greater that 1000 m a.s.l. contribute to the load;
elevation values lower than 1000 m are made equal to zero before the Fourier-¬ltering procedure
begins.

main()
{ float dum,delta,alpha,fact,*w,delrho;
int lattice_size_x,lattice_size_y,*nn,i,j;
FILE *fp1,*fp2;

fp1=fopen("load2dandes","r");
fp2=fopen("load2dandesdeflect50","w");
lattice_size_x=2048;
lattice_size_y=4096;
delta=1.1; /* km */
alpha=50.0; /* (D/(rho_c*g))^0.25 */
delrho=0.27; /* (rho_m-rho_c)/rho_m */
nn=ivector(1,2);
nn[1]=lattice_size_x;
nn[2]=lattice_size_y;
w=vector(1,2*lattice_size_x*lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f",&dum);
if (dum<1000) w[2*(i-1)*lattice_size_y+2*j-1]=0;
else w[2*(i-1)*lattice_size_y+2*j-1]=dum;
w[2*(i-1)*lattice_size_y+2*j]=0.0;}
fourn(w,nn,2,1);
w[1]*=1/delrho;
w[2]*=1/delrho;
for (j=1;j<=lattice_size_y/2;j++)
{fact=1/(delrho+pow(alpha*j*PI/lattice_size_y,4.0));
w[2*j+1]*=fact;
w[2*j+2]*=fact;
w[2*lattice_size_y-2*j+1]*=fact;
A4.4 ADI TECHNIQUE APPLIED TO GLACIAL ISOSTATIC RESPONSE MODELING 259


w[2*lattice_size_y-2*j+2]*=fact;}
for (i=1;i<=lattice_size_x/2;i++)
{fact=1/(delrho+pow(alpha*i*PI/lattice_size_x,4.0));
w[2*i*lattice_size_y+1]*=fact;
w[2*i*lattice_size_y+2]*=fact;
w[2*lattice_size_x*lattice_size_y-2*i*lattice_size_y+1]*=fact;
w[2*lattice_size_x*lattice_size_y-2*i*lattice_size_y+2]*=fact;}
for (i=1;i<=lattice_size_x/2;i++)
for (j=1;j<=lattice_size_y/2;j++)
{fact=1/(delrho+pow(alpha*sqrt((PI*PI*i*i)/(lattice_size_x*lattice_size_x)+
(PI*PI*j*j)/(lattice_size_y*lattice_size_y)),4.0));
w[2*i*lattice_size_y+2*j+1]*=fact;
w[2*i*lattice_size_y+2*j+2]*=fact;
w[2*i*lattice_size_y+2*lattice_size_y-2*j+1]*=fact;
w[2*i*lattice_size_y+2*lattice_size_y-2*j+2]*=fact;
w[2*lattice_size_x*lattice_size_y-2*(i-1)*lattice_size_y-
2*(lattice_size_y-j)+1]*=fact;
w[2*lattice_size_x*lattice_size_y-2*(i-1)*lattice_size_y-
2*(lattice_size_y-j)+2]*=fact;
w[2*lattice_size_x*lattice_size_y-2*lattice_size_y-
2*(i-1)*lattice_size_y+2*(lattice_size_y-j)+1]*=fact;
w[2*lattice_size_x*lattice_size_y-2*lattice_size_y-
2*(i-1)*lattice_size_y+2*(lattice_size_y-j)+2]*=fact;}
fourn(w,nn,2,-1);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",w[2*(i-1)*lattice_size_y+2*j-1]
/(lattice_size_x*lattice_size_y));
fclose(fp1);
fclose(fp2);
}


A4.4 ADI technique applied to glacial isostatic response modeling

The following code uses the ADI technique to solve the ¬‚exure equation with spatially variable
rigidity. The code accepts two ¬les: westtelcomb3km (a 3 km/pixel grid of elastic thickness values) and
westelamask3km (a mask of equilibrium lines recti¬ed to the same coordinate system as the elastic
thickness grid). These ¬les can be downloaded from the CUP website. The ADI method is repeated
5000 times in this example (the method converges slowly). In practice, the number of iterations
should be determined by the required level of accuracy.
#define PI 3.141592653589793
#define adiiterations 5000

int idum,*iup,*idown,*jup,*jdown,lattice_size_y,lattice_size_x;
float **D,delrho,deltax,deltay,**load,**deflect,**deflect2,te;
260 CODES FOR SOLVING THE FLEXURE EQUATION


void solveimpl()
{
int i,j,count;
float d,**ax,**ay,*xx,*xy,*rx,*ry,**alx,**aly,tot;
unsigned long *indx,*indy;

indx=lvector(1,lattice_size_x);
indy=lvector(1,lattice_size_y);
ax=matrix(1,lattice_size_x,1,5);
ay=matrix(1,lattice_size_y,1,5);
xx=vector(1,lattice_size_x);
xy=vector(1,lattice_size_y);
rx=vector(1,lattice_size_x);
ry=vector(1,lattice_size_y);
alx=matrix(1,lattice_size_x,1,2);
aly=matrix(1,lattice_size_y,1,2);
for (count=1;count<=adiiterations;count++)
{for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
deflect2[i][j]=deflect[i][j];
tot=0;
for (j=1;j<=lattice_size_y;j++) {
for (i=1;i<=lattice_size_x;i++)
if ((i>2)&&(i<lattice_size_x-1))
{rx[i]=-load[i][j];
ax[i][1]=D[i][j];
ax[i][2]=-4*D[i][j];
ax[i][3]=12*D[i][j]+delrho;
rx[i]+=-D[i][j]*(deflect2[i][jup[jup[j]]]+deflect2[i][jdown[jdown[j]]]-
4*(deflect2[i][jup[j]]+deflect2[i][jdown[j]]));
ax[i][4]=-4*D[i][j];
ax[i][5]=D[i][j];}
else
{rx[i]=0;
ax[i][1]=0;
ax[i][2]=0;
ax[i][3]=1;
ax[i][4]=0;
ax[i][5]=0;}
bandec(ax,lattice_size_x,2,2,alx,indx,&d);
banbks(ax,lattice_size_x,2,2,alx,indx,rx);
for (i=1;i<=lattice_size_x;i++)
deflect[i][j]=rx[i];}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
deflect2[i][j]=deflect[i][j];
for (i=1;i<=lattice_size_x;i++) {
for (j=1;j<=lattice_size_y;j++)
if ((j>2)&&(j<lattice_size_y-1))
A4.4 ADI TECHNIQUE APPLIED TO GLACIAL ISOSTATIC RESPONSE MODELING 261


{ry[j]=-load[i][j];
ay[j][1]=D[i][j];
ay[j][2]=-4*D[i][j];
ay[j][3]=12*D[i][j]+delrho;
ry[j]+=-D[i][j]*(deflect2[iup[iup[i]]][j]+deflect2[idown[idown[i]]][j]-
4*(deflect2[iup[i]][j]+deflect2[idown[i]][j]));
ay[j][4]=-4*D[i][j];
ay[j][5]=D[i][j];}
else
{ry[j]=0;
ay[j][1]=0;
ay[j][2]=0;
ay[j][3]=1;
ay[j][4]=0;
ay[j][5]=0;}
bandec(ay,lattice_size_y,2,2,aly,indy,&d);
banbks(ay,lattice_size_y,2,2,aly,indy,ry);
for (j=1;j<=lattice_size_y;j++)
deflect[i][j]=ry[j];}
}
}

main()
{ FILE *fp1,*fp2,*fp3;
int i,j;

fp1=fopen("westtelcomb3km","r");
fp2=fopen("westsnowlinemask3km","r");
fp3=fopen("westdeflectcombnew","w");
lattice_size_x=626;
lattice_size_y=642;
deltax=3; /* km */
deltay=3; /* km */
delrho=6000; /* kg/m^3 */
deflect=matrix(1,lattice_size_x,1,lattice_size_y);
deflect2=matrix(1,lattice_size_x,1,lattice_size_y);
load=matrix(1,lattice_size_x,1,lattice_size_y);
D=matrix(1,lattice_size_x,1,lattice_size_y);
setupgridneighborsperiodic();
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f",&te);
if ((i<5)||(i>lattice_size_x-5)||(j<5)||(j>lattice_size_y-5)) te=0;
D[i][j]=te*te*te*70000000/(deltax*deltax*deltax*deltax)/
(12*(1-0.25*0.25));
fscanf(fp2,"%f",&load[i][j]);
load[i][j]*=0.82; /* rho_c/rho_m */
if ((i<5)||(i>lattice_size_x-5)||(j<5)||(j>lattice_size_y-5)) load[i][j]=0;
deflect[i][j]=load[i][j];}
262 CODES FOR SOLVING THE FLEXURE EQUATION


solveimpl();
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp3,"%f\n",deflect[i][j]);
fclose(fp1);
fclose(fp2);
fclose(fp3);
}
Appendix 5




Codes for modeling non-Newtonian ¬‚ows

A5.1 2D radially symmetric lava ¬‚ow simulation

The following code implements the 2D radially symmetric gravity ¬‚ow model with temperature-
dependent viscosity described in Section 6.3.


main()
{ FILE *fp1;
int lattice_size_x,*iup,*idown,i;
float deltar,nu,ro,rn,*height,*heightold,*temp,*tempold,*flux;
float simulationlength,outputinterval,timestep,elapsedtime,maxdelt,maxdelh,del;

fp1=fopen("viscousflownu100","w");
lattice_size_x=260;
deltar=0.02;
nu=100.0;
ro=1.0;
rn=1.1;
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
height=vector(1,lattice_size_x);
heightold=vector(1,lattice_size_x);
temp=vector(1,lattice_size_x);
tempold=vector(1,lattice_size_x);
flux=vector(1,lattice_size_x);
simulationlength=10.0;
outputinterval=0.1;
for (i=1;i<=lattice_size_x;i++)
{temp[i]=exp(-pow(i*deltar/rn,20));
tempold[i]=temp[i];
if (i*deltar<=ro) height[i]=pow(4.0/(1+nu)*log(rn/ro),0.25);
else if (i*deltar<=rn) height[i]=
pow(4.0/(1+nu)*log(rn/(i*deltar)),0.25); else height[i]=0.001;
heightold[i]=height[i];
flux[i]=0;
264 CODES FOR MODELING NON-NEWTONIAN FLOWS


iup[i]=i+1;idown[i]=i-1;}
iup[lattice_size_x]=lattice_size_x;idown[1]=1;
timestep=0.00001;
elapsedtime=0.0;
while (elapsedtime<simulationlength)
{maxdelt=0;maxdelh=0;
flux[1]=1;
for (i=1;i<=lattice_size_x;i++)
if (i*deltar<ro) flux[i]=1;
else flux[i]=i*(1+nu*tempold[i])*pow(0.5*(heightold[idown[i]]+heightold[i]),3)*
(heightold[idown[i]]-heightold[i]);
for (i=1;i<=lattice_size_x-1;i++)
{del=timestep/(i*deltar*deltar)*(flux[i]-flux[iup[i]]);
height[i]+=del;
if (fabs(del)>maxdelh) maxdelh=del;
del=0;
if (heightold[i]>0)
{del=timestep*(flux[i]*(tempold[idown[i]]-tempold[i])/
(heightold[i]*i*deltar*deltar)-
tempold[i]/(heightold[i]*heightold[i]));
if (i*deltar<ro) del=0;
temp[i]+=del;} else del=0;
if (fabs(del)>maxdelt) maxdelt=fabs(del);
if (temp[i]<0) temp[i]=0;
if (temp[i]>1) temp[i]=1;}
for (i=1;i<=lattice_size_x;i++)
if (i*deltar<ro) height[i]=height[(int)(ro/deltar)];
elapsedtime+=timestep;
if ((maxdelt>0.1)||(maxdelh>0.1))
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<=lattice_size_x;i++)
{height[i]=heightold[i];
temp[i]=tempold[i];}}
else
if ((maxdelt<0.01)&&(maxdelh<0.01)) timestep*=1.2;
for (i=1;i<=lattice_size_x;i++)
{heightold[i]=height[i];
tempold[i]=temp[i];}}
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%d %f %f\n",i,height[i],temp[i]);
fclose(fp1);
}



A5.2 Sandpile method for ice-sheet and glacier reconstruction

The following code implements the sandpile method of ice-sheet or glacier reconstruction using the
Greenland ice sheet as an example. The code accepts the input ¬les greenlandsmallbed (a 10 km/pixel
A5.2 SANDPILE METHOD FOR ICE-SHEET AND GLACIER RECONSTRUCTION 265


DEM of the bed topography) and greenlandsmallmask (a 10 km/pixel mask grid de¬ning locations of
ice coverage) from the local directory and outputs the ¬le greenlandsmallsurf. The routine checkmin
is just a simple version of fillinpitsandflats that only applies to the masked area (i.e. the area
with ice cover).

#define fillincrement 0.1

float taugam,**topo,block,*heightvec,**mask,**height,delta;
int *heightvecind,*iup,*jup,*idown,*jdown,lattice_size_x,lattice_size_y,change;

void addiflessthanthreshold(i,j)
int i,j;
{ float slope,slopex,slopey;

height[i][j]+=block;
slopex=height[i][j]-height[iup[i]][j];
if (height[i][j]-height[idown[i]][j]>slopex) slopex=height[i][j]-height[idown[i]][j];
slopey=height[i][j]-height[i][jup[j]];
if (height[i][j]-height[i][jdown[j]]>slopey) slopey=height[i][j]-height[i][jdown[j]];
slope=sqrt(slopex*slopex+slopey*slopey)/delta;
/* this version implements slope-dependent threshold
tau(S)=15S^0.55 */
if (pow(slope,0.45)*(height[i][j]-topo[i][j])>taugam*15)
height[i][j]-=block;
else change=1;
}

void checkmin(i,j)
int i,j;
{ float min;

min=10000;
if (mask[i][j]>0.1) {
if (height[iup[i]][j]<min) min=height[iup[i]][j];
if (height[idown[i]][j]<min) min=height[idown[i]][j];
if (height[i][jup[j]]<min) min=height[i][jup[j]];
if (height[i][jdown[j]]<min) min=height[i][jdown[j]];
if (height[i][j]<min) {height[i][j]=min+0.1;
checkmin(iup[i],j);
checkmin(idown[i],j);
checkmin(i,jup[j]);
checkmin(i,jdown[j]);} }
}

main()
{ FILE *fp1,*fp2,*fp3;
int i,j,t;

fp1=fopen("greenlandsmallbed","r");
fp2=fopen("greenlandsmallmask","r");
266 CODES FOR MODELING NON-NEWTONIAN FLOWS


fp3=fopen("greenlandsmallsurf","w");
lattice_size_x=150;
lattice_size_y=280;
delta=10000.0; /* km */
taugam=100000.0/(920.0*9.8); /* tau = 10^5 Pa */
height=matrix(1,lattice_size_x,1,lattice_size_y);
mask=matrix(1,lattice_size_x,1,lattice_size_y);
topo=matrix(1,lattice_size_x,1,lattice_size_y);
setupgridneighbors();
heightvec=vector(1,lattice_size_x*lattice_size_y);
heightvecind=ivector(1,lattice_size_x*lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f",&topo[i][j]);
fscanf(fp2,"%f",&mask[i][j]);
height[i][j]=topo[i][j];}
block=100; /* start at 100 m block size */
while (block>0.9) /* stop when block size < 1 m */
{change=1;
while (change>0)
{change=0;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
heightvec[(j-1)*lattice_size_x+i]=height[i][j];
indexx(lattice_size_x*lattice_size_y,heightvec,heightvecind);
t=0;
while (t<lattice_size_x*lattice_size_y)
{t++;
i=(heightvecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(heightvecind[t])/lattice_size_x+1;
if (mask[i][j]>0.1)
addiflessthanthreshold(i,j);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
checkmin(i,j);}
block=block/3.33;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{if (height[i][j]>0)
fprintf(fp3,"%f\n",height[i][j]);
else fprintf(fp3,"0\n");}
fclose(fp1);
fclose(fp2);
fclose(fp3);
}
Appendix 6




Codes for modeling instabilities


A6.1 Werner™s eolian dune model

The following code implements Werner™s eolian dune model described in Section 7.4 with parameters
l = 5, ps = 0.6, pns = 0.4, and a shadow-zone angle equal to half of the angle of repose. The code was
written by the author based on the description in Werner (1995). The initial condition is a ¬‚at bed
with three sand slabs at each grid point.


int **height,*iup,*jup,*idown,*jdown,lattice_size_x,lattice_size_y;
float **mask,thresh;

void avalanchedown(i,j)
int i,j;
{ float hshadow;
int i2,violate,min,mini,minj;

violate=0;min=height[i][j];mini=i;minj=j;
if (height[i][j]-height[iup[i]][j]>thresh)
{violate=1;
if (height[iup[i]][j]<min)
{min=height[iup[i]][j];mini=iup[i];minj=j;}}
if (height[i][j]-height[idown[i]][j]>thresh)
{violate=1;
if (height[idown[i]][j]<min)
{min=height[idown[i]][j];mini=idown[i];minj=j;}}
if (height[i][j]-height[i][jup[j]]>thresh)
{violate=1;
if (height[i][jup[j]]<min)
{min=height[i][jup[j]];mini=i;minj=jup[j];}}
if (height[i][j]-height[i][jdown[j]]>thresh)
{violate=1;
if (height[i][jdown[j]]<min)
{min=height[i][jdown[j]];mini=i;minj=jdown[j];}}
if (violate==1)
{height[i][j]--;
268 CODES FOR MODELING INSTABILITIES


height[mini][minj]++;
avalanchedown(mini,minj);}
else
{hshadow=height[i][j]-0.5;
i2=iup[i];
while (height[i2][j]<hshadow)
{if (mask[i2][j]<hshadow) mask[i2][j]=hshadow;
i2=iup[i2];hshadow-=0.5;}}
}

void avalancheup(i,j)
int i,j;
{ float hshadow;
int i2,violate,min,mini,minj;

violate=0;min=height[i][j];mini=i;minj=j;
if (height[iup[i]][j]-height[i][j]>thresh)
{violate=1;
if (height[iup[i]][j]>min)
{min=height[iup[i]][j];mini=iup[i];minj=j;}}
if (height[idown[i]][j]-height[i][j]>thresh)
{violate=1;
if (height[idown[i]][j]>min)
{min=height[idown[i]][j];mini=idown[i];minj=j;}}
if (height[i][jup[j]]-height[i][j]>thresh)
{violate=1;
if (height[i][jup[j]]>min)
{min=height[i][jup[j]];mini=i;minj=jup[j];}}
if (height[i][jdown[j]]-height[i][j]>thresh)
{violate=1;
if (height[i][jdown[j]]>min)
{min=height[i][jdown[j]];mini=i;minj=jdown[j];}}
if (violate==1)
{height[i][j]++;
height[mini][minj]--;
avalancheup(mini,minj);}
else
{hshadow=height[i][j]-0.5;
i2=iup[i];
while ((hshadow<mask[i2][j])&&(hshadow>0))
{if (height[i2][j]>=hshadow) mask[i2][j]=0; else mask[i2][j]=hshadow;
i2=iup[i2];hshadow-=0.5;}}
}

main()
{ FILE *fp1,*fp2;
int idum,t,l,ijump,jjump,i,j,duration;
float psand,pbed,p;
A6.1 WERNER™S EOLIAN DUNE MODEL 269


fp1=fopen("wernermodeltopo","w");
fp2=fopen("wernermodelshadowmask","w");
idum=-56;
psand=0.6;
pbed=0.4;
thresh=1;
l=5;
lattice_size_x=300;
lattice_size_y=300;
duration=100;
height=imatrix(1,lattice_size_x,1,lattice_size_y);
mask=matrix(1,lattice_size_x,1,lattice_size_y);
setupgridneighborsperiodic();
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
{height[i][j]=3;
mask[i][j]=0;}
for (t=1;t<=duration*lattice_size_x*lattice_size_y;t++)
{if (t%(lattice_size_x*lattice_size_y)==0) printf("%d\n",t);
ijump=(int)(ran3(&idum)*lattice_size_x)+1;
while (ijump>lattice_size_x) ijump=(int)(ran3(&idum)*lattice_size_x)+1;
jjump=(int)(ran3(&idum)*lattice_size_y)+1;
while (jjump>lattice_size_y) jjump=(int)(ran3(&idum)*lattice_size_y)+1;
while ((height[ijump][jjump]==0)||(mask[ijump][jjump]>0.1))
{ijump=(int)(ran3(&idum)*lattice_size_x)+1;
while (ijump>lattice_size_x) ijump=(int)(ran3(&idum)*lattice_size_x)+1;
jjump=(int)(ran3(&idum)*lattice_size_y)+1;
while (jjump>lattice_size_y)
jjump=(int)(ran3(&idum)*lattice_size_y)+1;}
height[ijump][jjump]--;
avalancheup(ijump,jjump);
ijump=(ijump+l)%lattice_size_x+1;
if (mask[ijump][jjump]>0.1) p=1;
else if (height[ijump][jjump]>0) p=psand; else p=pbed;
while (ran3(&idum)>p)
{ijump=(ijump+l)%lattice_size_x+1;
if (height[ijump][jjump]>0) p=psand; else p=pbed;}
height[ijump][jjump]++;
avalanchedown(ijump,jjump);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fprintf(fp1,"%d\n",height[i][j]);
fprintf(fp2,"%f\n",mask[i][j]);}
fclose(fp1);
fclose(fp2);
}
270 CODES FOR MODELING INSTABILITIES



A6.2 Oscillations in arid alluvial channels

The following code solves the equations for arid channel oscillations described in Section 7.5.

#define PI 3.141592653589793

main()
{ FILE *fp1;
int idum,*iup,*idown,count,printed,outputinterval,length_y,lattice_size_x,i;
float courantwave,courantdiff,maxdelw,del,*h2,*w2,h0,*flux,*deltaw,*deltah;
float w0,simulationlength,deltax,elapsedtime,timestep,channelslope;
float *height,*heightold,*width,*widthold,maxheightchange;

fp1=fopen("channeloscillations","w");
lattice_size_x=100;
deltax=100; /* m */
idum=-76;
channelslope=0.01; /* m/m */
w0=100; /* m */
h0=1; /* m */
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
height=vector(1,lattice_size_x);
heightold=vector(1,lattice_size_x);
deltah=vector(1,lattice_size_x);
width=vector(1,lattice_size_x);
widthold=vector(1,lattice_size_x);
deltaw=vector(1,lattice_size_x);
flux=vector(1,lattice_size_x);
h2=vector(1,lattice_size_x);
w2=vector(1,lattice_size_x);
simulationlength=5000000.0; /* yr */
outputinterval=100000.0;
for (i=1;i<=lattice_size_x;i++)
{width[i]=w0;
widthold[i]=width[i];
height[i]=-0.05*sin((i-1)*2*PI*3/(lattice_size_x-1))+0.03*gasdev(&idum);
heightold[i]=height[i];
flux[i]=0;h2[i]=0;w2[i]=0;
deltaw[i]=0;deltah[i]=0;
iup[i]=i+1;idown[i]=i-1;}
iup[lattice_size_x]=1;idown[1]=lattice_size_x;
timestep=1;
elapsedtime=0.0;
printed=0;
while (elapsedtime<simulationlength)
{printed=elapsedtime/outputinterval;
maxdelw=0;
A6.3 1D MODEL OF SPIRAL TROUGHS ON MARS 271


for (i=1;i<=lattice_size_x;i++)
flux[i]=((heightold[idown[i]]-heightold[i])/deltax+channelslope)/widthold[i];
for (i=1;i<=lattice_size_x;i++)
{del=0.5*timestep/deltax*(flux[idown[i]]-flux[i]);
h2[i]=0.5*(heightold[idown[i]]+heightold[i])+del;
w2[i]=widthold[i]-0.5*timestep*0.5*
(flux[idown[i]]+flux[i]-2*channelslope/w0)/(h0-h2[i]);}
for (i=1;i<=lattice_size_x;i++)
flux[i]=((h2[i]-h2[iup[i]])/deltax+channelslope)/w2[i];
for (i=1;i<=lattice_size_x;i++)
{del=timestep/deltax*(flux[i]-flux[iup[i]]);
height[i]+=del;
width[i]-=0.5*timestep*(flux[iup[i]]+flux[i]-
2*channelslope/w0)/(h0-heightold[i]);
if (fabs(widthold[idown[i]]-widthold[iup[i]])/widthold[i]>maxdelw)
maxdelw=fabs(widthold[idown[i]]-widthold[iup[i]])/widthold[i];
if (width[i]<10) width[i]=10;
if (width[i]>300) width[i]=300;}
elapsedtime+=timestep;
if ((maxdelw>0.5*deltax/timestep)||(timestep>deltax*deltax))
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<=lattice_size_x;i++)
{height[i]=heightold[i];
width[i]=widthold[i];}}
else
if ((maxdelw<0.5*deltax*deltax/timestep/10)&&(timestep<deltax*deltax/10))
timestep*=1.2;
for (i=1;i<=lattice_size_x;i++)
{heightold[i]=height[i];
widthold[i]=width[i];}
if ((int)(elapsedtime/outputinterval)>printed)
{for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%d %f %f\n",i,height[i]/h0,width[i]/w0);
printed=(int)(elapsedtime/outputinterval);}}
fclose(fp1);
}


A6.3 1D model of spiral troughs on Mars

The following code integrates the 1D spiral trough model described in Section 7.7.

#define N 2
#define epssm 0.0001

float *xp,**yp,dxsav,a,eps;
int kmax,kount,*iup,*jup,*idown,*jdown;
272 CODES FOR MODELING INSTABILITIES


void derivs(x,y,dydx)
float dydx[],x,y[];
{
dydx[1]=y[2];
dydx[2]=(y[2]*(1-y[2])*(y[2]-a)-y[1])/eps;
}
main()
{ FILE *fp1;
int nbad,nok,t,lattice_size_x,i;
float *sol,h1=0.1,hmin=0.0,x1,x2,*lap,deltax,fac;
float *u,*v,timestep,kappa;



fp1=fopen("spiralmodel1d","w");
a=0.2;
eps=0.05;
deltax=0.2;
timestep=deltax*deltax/5;
kappa=timestep/(deltax*deltax);
lattice_size_x=512;
u=vector(1,lattice_size_x);
v=vector(1,lattice_size_x);
lap=vector(1,lattice_size_x);
sol=vector(1,2);
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
for (i=1;i<=lattice_size_x;i++)
{iup[i]=i+1;
idown[i]=i-1;}
iup[lattice_size_x]=1;
idown[1]=lattice_size_x;
for (i=1;i<=lattice_size_x;i++)
if ((i>lattice_size_x/2-5)&&(i<lattice_size_x/2-2))
{u[i]=0.25;v[i]=0.0;}
else
{u[i]=0;v[i]=0;}
/* plot initial condition first */
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%f %f %f\n",(i-lattice_size_x/2)/10.0,u[i],-v[i]);
fac=timestep*kappa/deltax/deltax;
x1=0;x2=timestep;
for (t=1;t<=1000;t++)
{for (i=1;i<=lattice_size_x;i++)

<<

. 13
( 14)



>>