ENVS397 Simulating environmental systems Stratigraphic Modelling ENVS397 Simulating environmental systems • Your aim in the practical this week is to construct a *very* simple stratigraphic forward model • Example output from the model is shown on the right • Key elements in the model are: • Definition of an initial topography • Calculation of diffusional erosion and deposition across the topography through multiple forward-in-time steps • Representation of the strata produced by this erosion and deposition, included eroded older layers, and a layer property such as grain size • Plot the results as a 3D final model surfaces, cross-section and chronostratigraphic diagram • You have some of the source code to do this already but you need to write the rest… Stratigraphic Model Outline ENVS397 Simulating environmental systems • The code for the main forward-in-time calculation loop should follow the structure and layout shown in this flow chart … Define and initialise model parameters, matrices and vectors Check the stability of the finite difference solutions Calculate the initial topography are iterations calculated <= total iterations? Calculate and store erosion and deposition in layers Increment iteration count Analyse the strata produced Plot results No Yes Stratigraphic Model Flowchart Calculate diffusion of the topography ENVS397 Simulating environmental systems • To see how diffusional erosion and deposition works, we want to start with a sharp-edged topography that diffusion will smooth… • You have suitable code already from the geomorphic modelling exercise, so that is a good place to start … Initial Topography ENVS397 Simulating environmental systems Calculation of diffusion for x=2:maxX-1 for y = 2:maxY-1 elevSum = newTopog(y,x-1) + newTopog(y,x+1) + newTopog(y-1,x) + newTopog(y+1,x); % 2d forward finite difference solution for diffusion deltaTopog(y,x) = -kappa * dt * (elevSum - (4 * newTopog(y,x))) / (dx * dx); end end smoothedTopog = newTopog - deltaTopog; • Diffusion function applies a simple 4-neighbour central- weighted finite difference approximation of the diffusion equation to the topography h at each grid point x,y • Note well – what might seem like complex maths is typically much simpler when you see it described in a diagram and as code, not just as a formula • In very general terms, think what kind of mathematical operation this is – hint: diffusion smooths the topography ∆,= −∆ −1, + +1, + ,+1 + ,−1 − 4, ∆2 x,yx-1,y x+1,y x,y-1 x,y+1 ENVS397 Simulating environmental systems • The finite difference approximation for diffusion of the topographic surface works well, but do you remember this from the geomorphic modelling solution code? diffusionStabilityCheck = (diffCoeff * dt) / (dx * dx); if (diffusionStabilityCheck >= 0.5) fprintf('Diffusion coefficient %1.4f too big for timestep and gridspacing - UNSTABLE!!!\n', diffusionStabilityCheck); end • Only certain combinations of model time step and grid cell size will generate a numerically stable solution • For example, if the time step dt is too big, so that then the solution is unstable and numerical artefacts will result . Which is bad! • A good way to avoid this is small time steps, but this has the disadvantage that it might require storage of lots of layers in the model • So how can we have a small time step but not too many stored layers? Think about it… x,yx-1,y x+1,y x,y-1 x,y+1 ∆ ∆2 > 1 2 Calculation of diffusion ENVS397 Simulating environmental systems • But there should also be sub loop to calculate the erosion and deposition due to diffusion, within the main model time loop… Are sub- iterations <= total sub- iterationsare iterations calculated <= total iterations? Calculate and store erosion and deposition in layers Increment iteration count No Yes Calculate diffusion of the topography No Yes Calculate diffusion from current topographic surface Increment iteration count Update topography • Why is this a useful way to construct the model? • Hint – what does NOT happen in the sub loop…? Stratigraphic Model Flowchart ENVS397 Simulating environmental systems Strata Produced by Diffusion • At each grid point, deposition or erosion (Δh) calculated according to equations of sediment transport Δh = k(x,y) * s(x) * U(x) • This is actually the diffusion equation that you have seen before • The most basic version of this equation is ℎ = 2ℎ 2 • It can be solved analytically or numerically using a finite difference solution method • In 2D a central weighted FD solution to the diffusion equation is ℎ , = ℎ+1,−ℎ−1, 2 • And you already have code that does this in 3D… X y Δh sealevel We can calculate Δh (same as ℎ , ) for point I,j, using the f.d. equation on the right… Initial topographic surface, time 1 Depositional surface, time 1 Depositional surface, time 2 ENVS397 Simulating environmental systems Strata Produced by Diffusion • At each grid point, deposition or erosion (Δh) calculated according to finite difference solution to the diffusion equation • Once we can calculate topography surfaces in 3D forward through time, how can we store these surfaces as layers/strata? • It’s effectively 4D information, so 3D layers with xyz coordinates, each layer with an age or position in time • Let’s try now, using the Matlab library function, to define a matrix structure that could represent these modelled stratal surfaces in 3D …. X y Δh sealevel Initial topographic surface, time 1 Depositional surface, time 1 Depositional surface, time 2 ENVS397 Simulating environmental systems Representing Deposition & Erosion • What happens between time 3 and time 6 in this case? • If you have studied any Earth Science, this pattern might look very familiar to you, but the question is, how do we code this process and the resulting layer geometry? • Hint – you can see that, for example, depositional surface 2 probably extended higher and has been eroded and truncated along depositional surface 6 • How do we represent that process as an algorithm that creates this kind of erosional and depositional layering structure? X y Initial topographic surface, time 1 Depositional surface, time 1 Depositional surface, time 2 Depositional surface, time 3 Depositional surface, time 6 Depositional surface, time 7 ENVS397 Simulating environmental systems • The diagram suggests that the key thing that happens to produce erosional truncation are • Depositional surfaces produced, with grid nodes increasing in elevation each time step • Erosion surface produced, with one or more nodes at elevations lower than previous depositional surface or surfaces at that node • Higher elevation nodes in previous layers are eroded and removed. What happens to their elevation values? • This is a key algorithm to devise and code that you need to write… X y Initial topographic surface, time 1 Depositional surface, time 1 Depositional surface, time 2 Depositional surface, time 3 Depositional surface, time 6 Previous extent of depositional surface, time 3, before erosion during interval time 4 to time 6 Eroded grid node Erosional grid node Depositional grid node Representing Deposition & Erosion ENVS397 Simulating environmental systems Representing Layer Properties • And finally, since we are modelling layers, it would be useful to represent in the model what these layers are made of • This is often referred to as property modelling, where a grid of cells in a model have one or more physical property associated with them • We already talked about a matrix structure that could represent strata in 3D • What if, instead of storing elevations of the layer top/base surface at each grid node we stored something else, like grain size? • Think about how to store layer properties… X y Δh Sealevel Rock type 1, sandstone Rock type 2, mudstone ENVS397 Simulating environmental systems Representing Layer Properties • A good layer property to calculate to get started is grain size • But how? • Grain size decreases with transport distance • So a very simple algorithm could calculate the distance from the source for the sediment and use this distance as a proxy for grain size • The simplest possible way to model the distance to the source is to calculate the distance to the highest point on the grid, in this simple model the central point • A better way would be the distanced down-slope from the up- slope start of deposition in any layer…. ENVS397 Simulating environmental systems Plotting: 3D View • This looks complicated to do, but it should be easy… • …either with the Matlab library routine surface, or with the patch command • You already have code that does something like this in the geomorphic model from the previous practical exercise • A more sophisticated version would draw a slice through the strata also on the edges of this plot… ENVS397 Simulating environmental systems Plotting: 2D Cross Section • You can draw this cross section with the Matlab graphic primitives library routines line and patch • Use examples from the geomorphic modelling to see how these work • The lines are the elevation of the stored strata • The colours are drawn with patch, using four grid points from two strata layers for each, and the stored rock property for the colour ENVS397 Simulating environmental systems Plotting: Chronostrat Diagram • You can also draw the chronostratigraphic section with the Matlab graphic primitives library routines line and patch • Use examples from the geomorphic modelling to see how these work • But the lines are the age of the stored strata, not the elevation • The colours are drawn with patch, using four grid points from two strata layers for each, and the store rock property for the colour • For example, here is a depth cross-section and an equivalent chronostratigraphic section from CarboCAT Horizontal distance 20 km El ev at io n 9 0 0 m G e o lo gi ca l t im e 3 M y ENVS397 Simulating environmental systems Plotting: Chronostrat Diagram X y Δh Sealevel Rock type 1, sandstone Rock type 2, mudstone X t Δt Initial topographic surface, time 1 Depositional surface, time 1 Depositional surface, time 2 Initial topographic surface, time 1 Depositional surface, time 1 Depositional surface, time 2 • A chronostratigraphic diagram is like a cross section but instead of plotting the elevation of the layers, it plots the geological age of the layers • Since the layes in the model are isochronous timelines, of equal age at all points along their length, they will plot as horizontal lines ENVS397 Simulating environmental systems Stratigraphic Model Flowchart • Follow the instructions for the practical • You will need to design some algorithms, using these slides, especially the flowcharts as a guide, and then code, debug and test those algorithms • Once complete, use the code to answer the questions posted in the practical instructions • When complete submit the Matlab code .m file on Canvas along with the PDF containing the short report on how the model behaves with different k values. Both these will be marked • NB Your submissions will be compared with submissions from all other students and examples with very similar text will be flagged as potential plagiarism etc



































































































































































































































































































































































































学霸联盟