FUN WITH EROSION Earth's surface topography is largely controlled by erosion processes, and the wonderful variety of landforms that we see must somehow be caused by differences in these processes. To model erosion from first principles would be a daunting task, because one would have to understand the physics of fluid flow, physical and chemical weathering, sediment transport and deposition, etc. However, if one takes a more empirical and less deterministic approach, it is possible to create realistic appearing synthetic topography with surprisingly simple computer algorithms. You will learn much more about this if you take Brad Werner's class on the Physics of Complex Systems, which deals with erosion and other geomorphological processes in detail. Our algorithm will simulate erosion over a 2-D grid in which elevation, z, is specified as a function of x and y coordinates. We will model the change in elevation with time as a function with the following form: dz/dt = -alpha*W*S + beta*Laplacian(z) where W = water flow at point (x,y) S = slope at point (x,y) alpha = adjustable parameter controlling rate of water erosion beta = adjustable parameter controlling rate of diffusion The effect of water erosion is assumed to be always negative in sign, i.e., only erosion and not sediment deposition is modeled. The diffusion term, however, can be either positive or negative, depending upon whether the point in question is a local high or low. Our approach will be to start off with a given topography model with a small amount of random variation. To make the code run reasonably fast, we will use a 100 by 100 grid in x and y. We will approximate the Laplacian by comparing the z value at a grid point with the average value of the 8 surrounding points, i.e., if z[][] is our elevation array, then zavg = (z[i-1][j-1] +z[i-1][j] +z[i-1][j+1] +z[i][j-1] +z[i][j+1] +z[i+1][j-1] +z[i+1][j] +z[i+1][j+1])/8.0; laplacian = zavg - z[i][j]; We will approximate the slope as the difference between the elevation at the target grid point and the lowest of its 8 surrounding neighbors, i.e., i1 = i - 1; /* find lowest of 8 surrounding points */ j1 = j - 1; if (z[i-1][j] < z[i1][j1] ) j1=j; if (z[i-1][j+1] < z[i1][j1] ) j1=j+1; if (z[i][j-1] < z[i1][j1] ) {i1=i; j1=j-1;} if (z[i][j+1] < z[i1][j1] ) {i1=i; j1=j+1;} if (z[i+1][j-1] < z[i1][j1] ) {i1=i+1; j1=j-1;} if (z[i+1][j] < z[i1][j1] ) {i1=i+1; j1=j;} if (z[i+1][j+1] < z[i1][j1] ) {i1=i+1; j1=j+1;} slope = z[i][j] - z[i1][j1]; Determining the water flux is the trickiest part of the program. One approach would be to compute the area of the upstream basin that drains into each point. However, this would be computationally messy to have to compute. We adopt a simpler method based on the "precipiton" model of Chase (1992). At each time step, we assume that there is uniform rainfall over the entire surface with the volume of water to be added given by an adjustable parameter (termed "rain" in the program). We then assume that this water at each point will flow into the lowest surrounding point. The water volume at each point is the sum of the contribution from the rain plus any water that flows into the point from surrounding points. Depending upon the size of the model, a steady state will be reached in a fairly small number of time steps in which the water volume, W, will be directly proportional to the upstream basin area. Without further ado, here is the complete program: (erosion.c) #include #include /* contains srand and rand functions */ #include #include "/net/rock/shearer/CLASS/COMP/C/nrutil.c" #define XMAX 100 #define YMAX 100 int main() { FILE *fp2; float **z, **z1, **w, **w1, xran, slope, rain, alpha, beta, gamma, zavg, laplacian, temp; int i, j, i1, j1, n, model, menu; srand(5); /* initialize generator */ z = matrix(1,XMAX,1,YMAX); z1 = matrix(1,XMAX,1,YMAX); w = matrix(1,XMAX,1,YMAX); w1 = matrix(1,XMAX,1,YMAX); printf("(1) wedge or (2) corner drop? \n"); scanf("%d",&model); for (i=1; i<=XMAX; i++){ for (j=1; j<=YMAX; j++){ xran = 0.1*( (float)rand()/(float)(RAND_MAX+1) - 0.5); if (model == 1) /* z varies from ~1 at x=1 to ~9 at x=XMAX */ z[i][j] = ((float)i)*9.0/(float)XMAX + xran; else z[i][j] = 9.0 + xran; z1[i][j] = z[i][j]; w[i][j] = 0.0; w1[i][j] = 0.0; } } if (model != 1) { z[1][1] = 7.0; z1[1][1] = 7.0; } rain = 0.1; alpha = 0.01; beta = 0.03; n = 0; while (1) { for (i=2; i 0.0) { temp = alpha*w[i][j]; if (temp > 0.3) temp = 0.3; z1[i][j] -= temp*slope; } else { /* local low spot */ z1[i][j] = z[i1][j1]; /* fill in with sediment */ } zavg = (z[i-1][j-1] +z[i-1][j] +z[i-1][j+1] +z[i][j-1] +z[i][j+1] +z[i+1][j-1] +z[i+1][j] +z[i+1][j+1])/8.0; laplacian = zavg - z[i][j]; z1[i][j] += beta * laplacian; if (z1[i][j] < 0) z1[i][j] = 0.0; if (z1[i][j] < 0.0 || z1[i][j] > 9.99) { printf("Problem with z1 \n"); printf("n = %d \n",n); printf("i,j,i1,j1 = %d %d %d %d \n",i,j,i1,j1); printf("%f %f %f \n",z1[i][j],z[i][j],z[i1][j1]); printf("%f \n",slope); return 0; } } /* end loop on j */ } /* end loop on i */ /* now set edge cells to elevation of adjacent cell */ for (i=1; i<=XMAX; i++){ for (j=1; j<=YMAX; j++){ i1 = i; j1 = j; if (i == 1) i1=2; if (i == XMAX) i1=XMAX-1; if (j == 1) j1=2; if (j == YMAX) j1=YMAX-1; z[i][j] = z1[i1][j1]; w[i][j] = w1[i1][j1]; w1[i][j] = 0.0; } } if (model != 1) { z[1][1] = 7.0; z1[1][1] = 7.0; } ++n; if (n % 100 == 0 ) { printf("n = %d \n",n); while (1) { printf("(1) continue, (2) print topo, (3) save, (4) print flow, (5) stop \n"); scanf("%d",&menu); if (menu == 1) break; else if (menu == 2) { for (j=1; j<=YMAX; j++){ for (i=1; i<=XMAX; i++){ printf("%d",(int)z[i][j]); } printf("\n"); } } else if (menu == 3) { fp2 = fopen("out.erosion","w"); if (fp2 == NULL) { printf("***Can't open file: %s\n", "out.erosion"); return 1; } fprintf(fp2,"%d %d %d\n",XMAX,YMAX,n); fprintf(fp2,"%f %f %f\n",rain,alpha,beta); for (j=1; j<=YMAX; j++){ for (i=1; i<=XMAX; i++){ fprintf(fp2,"%6.3f",z[i][j]); } fprintf(fp2,"\n"); } fclose(fp2); } else if (menu == 4) { for (j=1; j<=YMAX; j++){ for (i=1; i<=XMAX; i++){ temp = sqrt(w[i][j]); /* if (i*j == 2500) printf("test: %d %d %f \n",i,j,temp); */ if (temp < 0.5) printf(" "); else if (temp > 9.4) printf("*"); else printf("%d",(int)temp); } printf("\n"); } } else return 0; } } } } A few comments: (1) The user gets to choose between a starting model in which the topography is wedge shaped with elevations from 1 to 9 and a model of a "plateau" of height 9 in which the corner elevation suddenly drops to 7. (2) For each model a small random elevation perturbation of 0.1 is added to each point. This will determine the shape of the topography that results. To get different topography, use a different value in srand(5) which seeds the random number generator. (3) The values of rain = 0.1; alpha = 0.01; beta = 0.03; seem to give OK results. I have not had time to fully explore what happens when different values are tried. (4) It is possible that the initial random perturbations to the topography will cause local low spots with no outlets. This will result in a negative slope, in which case the program fills in the grid point to the level of the lowest of the surrounding points. In a sense, this simulates sediment filling in the basin. (5) Notice that the main code does not operate on the edges of the model. Obviously the edges cause problems because one cannot examine all of the surrounding cells. Instead of writing special code to handle this, the simple boundary condition is later imposed that the elevation at these points is always the same as the adjacent point. (6) Notice how the temporary arrays z1 and w1 are used to store the new elevations and water volumes for the next time step. Otherwise the calculation would be biased by the fact that some of the z values would have changed before the adjacent cells were computed. Only when all calulations are finished at each time step are the z and w values replaced with their z1 and w1 counterparts. (7) After every 100 iteration, the program allows the user to print a crude screen version of the topography or flow volume, save the current model to a file, and/or to continue the computation. The output file format is ascii with two header lines that give the values of various parameters in the program. The topography models can be plotted using one of the following Matlab scripts: (plotsurf.m) clf reset clear [fid,message] = fopen('out.erosion') [head1,count] = fscanf(fid,'%d %d %d',3); nx = head1(1); ny = head1(2); nit = head1(3); [head2,count] = fscanf(fid,'%f %f %f',3); rain = head2(1); alpha = head2(2); beta = head2(3); [z,count] = fscanf(fid,' %f',[nx,ny]); axis([0 100 0 100 0 10]); hold on; colormap(gray); surfl(z); shading flat; title (['nit,rain,alpha, beta = ', int2str(nit), ' ', num2str(rain),' ', num2str(alpha), ' ', num2str(beta)]); (plotmesh.m) clf reset clear [fid,message] = fopen('out.erosion') [head1,count] = fscanf(fid,'%d %d %d',3); nx = head1(1); ny = head1(2); nit = head1(3); [head2,count] = fscanf(fid,'%f %f %f',3); rain = head2(1); alpha = head2(2); beta = head2(3); [z,count] = fscanf(fid,' %f',[nx,ny]); axis([0 100 0 100 0 10]); hold on; colormap([0 0 0]); meshc(z); title (['nit,rain,alpha, beta = ', int2str(nit), ' ', num2str(rain),' ', num2str(alpha), ' ', num2str(beta)]); (plotcont.m) clf reset clear [fid,message] = fopen('out.erosion') [head1,count] = fscanf(fid,'%d %d %d',3); nx = head1(1); ny = head1(2); nit = head1(3); [head2,count] = fscanf(fid,'%f %f %f',3); rain = head2(1); alpha = head2(2); beta = head2(3); [z,count] = fscanf(fid,' %f',[nx,ny]); colormap([0 0 0]); contour(z,20); title (['nit,rain,alpha, beta = ', int2str(nit), ' ', num2str(rain),' ', num2str(alpha), ' ', num2str(beta)]); (plotmap.m) clf reset clear [fid,message] = fopen('out.erosion') [head1,count] = fscanf(fid,'%d %d %d',3); nx = head1(1); ny = head1(2); nit = head1(3); [head2,count] = fscanf(fid,'%f %f %f',3); rain = head2(1); alpha = head2(2); beta = head2(3); [z,count] = fscanf(fid,' %f',[nx,ny]); z = flipud(z); imagesc([0,1],[0,1],z); colorbar; title (['nit,rain,alpha, beta = ', int2str(nit), ' ', num2str(rain),' ', num2str(alpha), ' ', num2str(beta)]); This erosion simulation is one of the simplest possible algorithms that generates realistic looking topography. More sophisticated approaches might include: (1) Specifying the erosion rate as a more complicated function of water flow and slope. One example would be to include adjustable exponents for W and S. (2) Distinguishing between river channel vs. non-channel flow and erosion. (3) Permitting different erosion rates for different rock types. For example, bedrock and alluvium could be handled separately. (4) Including sediment transport and deposition in areas of decreasing slope. (5) Permitting sub-surface as well as surface water transport. (6) Including an uplift rate function, U(x,y), to more accurately reproduce steady state erosion. (7) Experiment with non-Cartesian grids (random?) to avoid the tendency for streams to occur preferentially at azimuth increments of 45 and 90 degrees. Further Reading: Chase, C., Fluvial landsculpting and the fractal dimension of topography, Geomorphology 5, 39-57, 1992. Rosenbloom, N.A. and R.S. Anderson, Hillslope and channel evolution in a marine terraced landscape, Santa Cruz, California, J. Geophys. Res., 99, 14013-14029, 1994. Tucker, G.E. and R.L. Bras, Hillslope processes, drainage density, and landscape morphology, Water Resources Research 34, 2751- 2764, 1998. Tucker, G.E. and R. Slingerland, Drainage basin responses to climate change, Water Resources Research 33, 2031-2047, 1997 Willgoose et al., A coupled channel network growth and hillslope evolution model 1. Theory, Water Resources Research 27, 1671- 1684, 1991. EROSION SPECIAL PROJECT Experiment with modifications to the erosion program to include two or more of the following: (1) A more resistant layer over a restricted depth interval in the model. (2) An imposed steady uplift rate in the middle of the model to simulate mountain building tectonic activity. (3) Long-term steady strike slip motion across a "fault" in the model. This can be done by shifting the topography on one side of the fault by one grid cell after every 50 to 200 time steps (?, experiment!). The idea is to generate offset stream channels that are diagnostic of strike slip faults. Please turn in both your source code and representative hard copies of images of the topography. Have fun!