<< . .

. 45
( : 51)



. . >>

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;

<< . .

. 45
( : 51)



. . >>