|
| 1 | +# 2024 Cape Mendocino Mw 7.0 earthquake - shakeMovie visualization |
| 2 | + |
| 3 | +--- |
| 4 | + |
| 5 | + |
| 6 | + |
| 7 | +This example will setup a model with the in-house mesher to simulate the |
| 8 | +2024 Cape Mendocino earthquake using the USGS finite-fault solution. |
| 9 | + |
| 10 | +About the event: |
| 11 | +* The Mendocino earthquake happened 2024-12-05 18:44:21 UTC with an estimated moment magnitude Mw 7.0. |
| 12 | + The USGS provides the following event page about the earthquake: |
| 13 | + [2024 Mendocino Mw7.0 earthquake](https://earthquake.usgs.gov/earthquakes/eventpage/nc75095651/executive) |
| 14 | + |
| 15 | + We will use the derived finite-fault solution of this event, provided on this detail page: |
| 16 | + [Mendocino Mw7.0 finite-fault solution](https://earthquake.usgs.gov/earthquakes/eventpage/nc75095651/finite-fault?source=us&code=nc75095651_1) |
| 17 | + |
| 18 | + To download this USGS finite-fault solution in a CMT-format for our SPECFEM3D example, you can get it into the `DATA/` folder by: |
| 19 | + ``` |
| 20 | + > cd ./DATA/ |
| 21 | + > wget https://earthquake.usgs.gov/product/finite-fault/nc75095651_1/us/1733529167722/CMTSOLUTION |
| 22 | + ``` |
| 23 | + |
| 24 | +Given the size of this earthquake and its wide-spread affected area, we will use a lon/lat range of about 4 x 4 degrees. |
| 25 | + |
| 26 | +Our region of interest: |
| 27 | + |
| 28 | + Longitude range ~ -126.5 to -122.0 |
| 29 | + Latitude 38.2 to 42.2 |
| 30 | + |
| 31 | + |
| 32 | + |
| 33 | +## Step-by-step |
| 34 | + |
| 35 | +We'll setup first the surface topography of this region, together with an internal interface to build a spectral-element mesh. |
| 36 | +Our intention here is to run a simulation for creating a shakeMovie visualization, thus the focus is on creating a rather simple mesh. |
| 37 | +However, you can use this as a starting point to build your own, more detailed and sophisticated, meshes. |
| 38 | + |
| 39 | +Furthermore, we want to use a (realistic) velocity model for this region. |
| 40 | +We will take one of the EarthScope (IRIS) [EMC models](https://ds.iris.edu/ds/products/emc-earthmodels/). |
| 41 | +For this region, a nice model seems to be [WUS256](https://ds.iris.edu/ds/products/emc-wus256/) by Rodgers et al. (2022). |
| 42 | + |
| 43 | + |
| 44 | + |
| 45 | +1. **Setup model**: |
| 46 | + |
| 47 | + We will first setup the topography surface and tomographic model for meshing our region. |
| 48 | + In this example folder, we provide a bash script `setup_model.sh` to setup these model files. |
| 49 | + |
| 50 | + Just run the setup script with the target region: |
| 51 | + ``` |
| 52 | + > ./setup_model.sh -126.5 -122.0 38.2 42.2 |
| 53 | + ``` |
| 54 | + |
| 55 | + This might take a while to complete... |
| 56 | + |
| 57 | + There are mainly two steps executed: |
| 58 | + * Topography data: it will use the script `utils/scripts/run_get_simulation_topography.py` to download all the needed topography data and create an interface file in folder `./topo_data/`, |
| 59 | + with the main surface topography stored in file `topo_data/ptopo.xyz.1.dat`. |
| 60 | + The script determines in what UTM zone this region lies, |
| 61 | + as well as how the interface definition should look like for the mesher setup. |
| 62 | + It then automatically updates the setup files `DATA/Par_file`, `DATA/meshfem3D_files/Mesh_Par_file` and `DATA/meshfem3D_files/interfaces.dat` to reflect region, UTM zone and interface definition. |
| 63 | + |
| 64 | + * Velocity model: it will download the EMC model into a folder `IRIS_EMC/` |
| 65 | + and create the tomography model for our region using the script `utils/scripts/run_convert_IRIS_EMC_netCDF_2_tomo.py`. The tomography model file `tomography_model.xyz` will be put into folder `DATA/tomo_files/`. |
| 66 | + |
| 67 | + Note that for meshing, we use a second, shifted internal interface to make it easier to use doubling layers. These doubling layers help coarsening the mesh with depth. This way, we will have a mesh that is mostly refined at the top to capture surface waves (which is what we are mostly interested in for visualizing the shaking). The mesher parameter file `DATA/meshfem3D_files/Mesh_Par_file` is setup accordingly to include these doubling layers, working together with the `DATA/meshfem3D_files/interfaces.dat` file which defines the surface and internal interface and number of element layers with depth. |
| 68 | + |
| 69 | + |
| 70 | +2. **Wave simulation**: |
| 71 | + |
| 72 | + After you have completed the mesh setup, you can run the in-house mesher `xmeshfem3D` and `xgenerate_databases` to create the spectral-element mesh of our region. |
| 73 | + The seismic wave propagation solver `xspecfem3D` then creates the needed movie data files for our visualization. |
| 74 | + |
| 75 | + To run the simulation, just type: |
| 76 | + ``` |
| 77 | + > ./run_this_example.sh |
| 78 | + ``` |
| 79 | + |
| 80 | + |
| 81 | + |
| 82 | +3. **Visualization**: |
| 83 | + |
| 84 | + For our simulation here, we turned on the surface movie and shake map outputs in `DATA/Par_file` like: |
| 85 | + ``` |
| 86 | + # save AVS or OpenDX movies |
| 87 | + # MOVIE_TYPE = 1 to show the top surface |
| 88 | + # MOVIE_TYPE = 2 to show all the external faces of the mesh |
| 89 | + CREATE_SHAKEMAP = .true. |
| 90 | + MOVIE_SURFACE = .true. |
| 91 | + MOVIE_TYPE = 2 |
| 92 | + MOVIE_VOLUME = .false. |
| 93 | + SAVE_DISPLACEMENT = .false. |
| 94 | + MOVIE_VOLUME_STRESS = .false. |
| 95 | + USE_HIGHRES_FOR_MOVIES = .true. |
| 96 | + NTSTEP_BETWEEN_FRAMES = 100 |
| 97 | + HDUR_MOVIE = 0.0 |
| 98 | + ``` |
| 99 | + |
| 100 | +  |
| 101 | + |
| 102 | + To visualize the corresponding output data (`OUTPUT_FILES/moviedata***`), we can create movie snapshot files as `OUTPUT_FILES/AVS_*.inp` files: |
| 103 | + ``` |
| 104 | + > ./xcreate_movie_files.sh |
| 105 | + ``` |
| 106 | + |
| 107 | + Similar for the shake map, we can plot the peak-ground velocity (PGV) values by |
| 108 | + ``` |
| 109 | + > ./xcreate_shakemap.sh 2 |
| 110 | + ``` |
| 111 | + |
| 112 | + You can use for example [Paraview](https://www.paraview.org) to look at these *.inp files. |
| 113 | + |
| 114 | + |
| 115 | +That's it, time for you to play! Here below are some more details to try out... |
| 116 | + |
| 117 | +## VS30 surface data (optional) |
| 118 | + |
| 119 | +USGS provides a global [Vs30 model and data](https://earthquake.usgs.gov/data/vs30/) set. |
| 120 | +We can download and extract a corresponding Vs30-interface for our region: |
| 121 | + |
| 122 | +``` |
| 123 | +> ln -s ../../../utils/scripts/run_get_simulation_USGS_Vs30.py |
| 124 | +> mkdir -p USGS_VS30 |
| 125 | +> ./run_get_simulation_USGS_Vs30.py -126.8 38.1 -121.9 42.3 |
| 126 | +``` |
| 127 | +This will create an interface file `./USGS_VS30/interface_vs30.dat`. |
| 128 | + |
| 129 | +Note that we use a slightly extended region to make sure all the mesh points will find corresponding Vs30 interface values. (This is a technical detail: the mesher converts the lon/lat coordinates of the corner points specified in `Mesh_Par_file` to UTM coordinates and determines the minimum/maximum UTM coordinate values. The latter will be used for the mesh grid points creation. Due to the UTM zone distortion, this can lead to grid points that if we would convert their coordinates back to geographic coordinates, are slightly beyond the minimum/maximum longitude and latitude range as specified in `Mesh_Par_file`). |
| 130 | + |
| 131 | +In case we want the meshing to be aware of this Vs30 interface |
| 132 | +and incorporate the Vs values for the top 30m GLL points, we can move it to our `DATA/` folder: |
| 133 | +``` |
| 134 | +> mv ./USGS_VS30/interface_vs30.dat DATA/ |
| 135 | +``` |
| 136 | +The `xgenerate_databases` binary will check if this file is found in `DATA/` and if so, will use the corresponding interface data |
| 137 | +to over-impose the Vs values on the corresponding grid nodes. |
| 138 | +To avoid creating artificial rock properties, it will also scale Vp and density accordingly for crustal rocks (Brocher 2005, *Empirical Relations between Elastic Wavespeeds and Density in the Earth's Crust*, BSSA). |
| 139 | + |
| 140 | +Note: |
| 141 | + |
| 142 | + Using Vs30 velocities requires a very high-resolution mesh at the top, |
| 143 | + otherwise these low velocities will "smear down" to deeper depths because of the |
| 144 | + GLL interpolation within the spectral elements. |
| 145 | + |
| 146 | + |
| 147 | + |
| 148 | +Regarding meshing for high resolutions, here are some further notes to consider: |
| 149 | + |
| 150 | +* USGS Vs30 - with a 30-m thickness layer: |
| 151 | + |
| 152 | + In principle, we would want an element with 30 m thickness at the top for these Vs30 values. |
| 153 | + The element below then would resort to the background model velocities. This would create a "sharp" model interface. |
| 154 | + Unfortunately, this will also create a tiny time step size for the whole simulation (unless we would use the local-time stepping feature `LTS_MODE` in `DATA/Par_file` to circumvent this global time step size restriction). |
| 155 | + |
| 156 | + Could we use a larger top element and have a somewhat "smoother" transition from Vs30 values to background velocities? |
| 157 | + |
| 158 | + The minimum GLL point distance within a spectral element depends on the order of the element, i.e., for a 4th-order element (NGLL == 5) the minimum distance is given by `d_min = dx * element_size` with the factor `dx = 0.1726` for `NGLL == 5`. |
| 159 | + |
| 160 | + For example, having an element size of 500 m leads to a minimum GLL point distance of `dx * 500 m = 86.3 m` for NGLL == 5. That is, the velocity profile sampled within the top element assigns the Vs30 velocity to the uppermost GLL point and then resorts back to the background model velocity for the second GLL point at a depth of 86.3 m. |
| 161 | + |
| 162 | + To create such mesh elements at the top, an interface for element sizes with 500-m thickness can be constructed by scaling down the surface topography interface accordingly: |
| 163 | + ``` |
| 164 | + > awk '{printf("%.6f\n",$1-500.0);}' topo_data/ptopo.xyz.1.dat > ./DATA/meshfem3D_files/ptopo.xyz.500m.dat |
| 165 | + ``` |
| 166 | + This interface can be added to the `DATA/meshfem3D_files/interfaces.dat` file and assign a single element layer between this 500-m interface and the top (topography) surface. |
| 167 | + |
| 168 | + To give it a try, a corresponding `interfaces.vs30.dat` file is provided in folder `DATA/meshfem3D_files/`. You could replace the default `interfaces.dat` file with this one. However, you will also need to decrease the time step size `DT` in the main parameter file `DATA/Par_file` before re-running the simulation for comparison. |
| 169 | + |
| 170 | +## Blender visualization |
| 171 | + |
| 172 | + |
| 173 | + |
| 174 | +Once created, the shakeMap file `OUTPUT_FILES/AVS_shaking_map.inp` can also be visualized with [Blender](https://www.blender.org), a rather professional 3D visualization software. To make things easy, we have some python scripts in folder `utils/Visualization/Blender/python_blender/` to help render images. |
| 175 | + |
| 176 | +With Blender installed on your system, you could try: |
| 177 | +``` |
| 178 | +> ln -s ../../../utils/Visualization/Blender/python_blender/plot_with_blender.py |
| 179 | +> ./plot_with_blender.py --vtk_file=OUTPUT_FILES/AVS_shaking_map.inp --title="Mendocino (PGV)" --color-max=0.25 --transparent-sea-level --background-dark --vertical-exaggeration=2.0 |
| 180 | +``` |
| 181 | + |
| 182 | +## Reference solution |
| 183 | + |
| 184 | +For comparison, we provide a reference solution in folder `REF_SEIS/` with corresponding output files. |
| 185 | +The simulation was run in parallel using 16 MPI processes, each using a single Nvidia A100 GPU. The total simulated time is 140 s with a time-to-solution of ~2 min. |
0 commit comments