[MUSIC] So let's explore the potential of this two-dimensional finite difference approach with a heterogeneous model with an example taken from geophysics, actually from seismology. I'm a seismologist, so I'm interested in earthquakes. So what is an earthquake? An earthquake is a rupture that happens usually often on a vertical fault, on a vertical plane, an example is the San Andreas Fault in California. And because there are many slips that happen over thousands of years, actually there is damage created around that relatively narrow fault zone, where this slip is taking place. Now, this actually generates damage, and the damage means that the elasticity changes and the velocities become smaller. And that's what you want to explore. What happens if at depth an earthquake happens, in such a so-called low velocity zone? And you see here an example, a picture or representation of the situation. And I would also like to show you a photo of an epicentral area that was taken in Landers, where actually a quite sizable earthquake happened in the 90s in California. And there you can actually see a displaced road, the road that was going straight over that fault zone. And it's actually displaced by a couple of meters. So, the next step is we're trying to set up a problem. We're trying to determine the finite difference parameters that we need to actually make a useful simulation to understand the science of wave propagation inside the fault zone, and that's going to be the next step. So what we do now actually preparing the finite difference calculation is something that is so fundamental that you actually do whenever you start or whenever you plan a certain simulation task. So that's very, very important. So let's start look again at the, let's look again at the model that we would like to simulate. This is a two-dimensional model now. We have a block, and we decide this is 10 by 10 kilometers. And we would like to understand wave propagation through a fault zone at the center of this model with a width of 200 meters. And now we're going to try and see how many grid points we need, what’s the time step, and so forth for such a simulation. Now for that, let's go step by step through the simulation parameters. As I said, our domain, our physical domain is now 10 by 10 kilometers. The velocity in the background material, the velocity of the p wave, let’s say at this point, it doesn’t matter, is 3,000 meters per second. Now, the velocity inside that fault zone that's indicated here is lower. It's 2,250 meters per second, that’s a typical kind of velocity change we observe in nature. We would like to simulate a seismogram length of about 3.5 seconds, we will later see why. And above an earthquake we often record signals with a dominant frequency of around 10 Hz, and let's assume a maximum frequency here of 30 Hz. Now these parameters already basically set the conditions now for the space time discretization of our scheme. Now, the first thing is what's the minimum wavelength that we'll be propagating in that medium. The minimum wavelength can be calculated by the minimum velocity divided by the maximum frequency. So, that's the 2,250 meters per second divided by 30 hertz, and that's actually 75 meters. Now, the source will be at 5 km depth, that's the next step. We will propagate around 20 dominant wavelengths. The dominant wavelengths, now we get by taking the ratio of the velocity, the minimum velocity divided by the dominant frequency. And by doing that, we'll actually see that we will propagate around 20 dominant wavelengths from the source up to the surface, where we would like to record the ground motion. Now what is the space increment that we need? Now let's assume we've talked about this before that we need 20 grid points per wavelengths to get an accurate result. We will later see in the simulation whether that is actually the case. So for the dominant wavelength that's 225 meters, if we divide this by 20, we will actually obtain a space increment of 11.25 meters. So that's often, we call that dx, dx is 11.25 meters. And here, we assume the same grid increment in both directions, so that dz is equal to dx. So basically, now we can already say what's the size of our computational domain because the maximum, the space domain, is 10 kilometers. So if we divide 10 kilometers by dx, we will already see how many grid points we need along one dimension, that's around 900. And that means for our two-dimensional grid, we will have obtained 900 square grid points for one field, for example, the pressure field at a certain time or the distribution of velocities in our grid. And if we have a double precision number, that would be 8 bytes. And let's assume that we will need around 6 space dependent fields, then we actually obtain a memory requirement of about 40 MB. That's not a lot, that's actually something you could probably already simulate on your smartphone. So what is remaining is to see what is the time increment. Our velocity is 3,000 meters per second, and remember the time increment will be determined by the so-called courant-criterion, it's written again here. So the maximum velocity, multiplying dt by dx, must be smaller than epsilon, which is close to 1. Here, we assume an epsilon of 0.7, and that leads to a time increment of 0.0026 seconds. We said at the beginning that we would like to have a seismogram of a simulation time of 3.5 seconds. So now we know how many grid points will be, how many points in time will be needed. That's 3.5 seconds divided by dt, and that's about 1300 points in time. So again, these kinds of considerations are very important before starting a simulation, and there will also be exercises for you to play around with these concepts. Now, what we did not yet discuss in detail are the boundary conditions. Actually, the boundary conditions for this problem are very simple. This is the scalar wave equation. We can assume that we now do a pressure, a wave-field calculation, and we can simply set the pressure to 0 at the boundaries, which means they will be fully reflecting. Also, there's some other stuff that we will discuss when we look at the actual Python code and the simulation, and then it starts getting real fun. We'll see what the wave field does for such a situation.
