<< Пред. стр. стр. 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)

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 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++)