<< . .

. 46
( : 51)



. . >>

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

<< . .

. 46
( : 51)



. . >>