Using Volumetric Models

Using Volumetric Models#

*** TO DO: This is a mashup of 2 different sources so it needs editing. It is also very long and probably should be sectioned into different files.

The 3D medium simulated by AxiSEM3D has the shape of a body of revolution with the cross section prescribed by the 2D mesh. Using a local-scale simulation as an example: if the 3D input model is given in Cartesian coordinates, i.e. by nature in shape of a cuboid, AxiSEM3D only reads the subset of this model which fits into the cylinder which is being simulated. 3D models which are smaller than the simulated volume may also be used. In this case regions not covered by the 3D model are filled in with the background model given in the 2D mesh.

The medium properties which can be set using volumetric models are:

  • homogeneous: “VS”, “VP”, “RHO”

  • transversely isotropic: “VSV”, “VSH”, “VPV”, “VPH”, “ETA”

  • fully anisotropic: (this is the stiffness tensor in voigt notation; all components must be provided)

    “C11”, “C12”, “C13”, “C14”, “C15”, “C16”,

    “C22”, “C23”, “C24”, “C25”, “C26”,

    “C33”, “C34”, “C35”, “C36”,

    “C44”, “C45”, “C46”,

    “C55”, “C56”,

    “C66”

Once you have created a new 3D model and saved it as a netcdf file, you need to upload it to the input folder. Then, there are two input files that you need to alter in order to get the model to be read in and represented properly: nparam.model.yaml and inparam.nr.yaml.

inparam.model.yaml#

The portion of the inparam.model.yaml file that needs to change from the 1D case is the 3D Models section. Under the header list_of_3D models, you can add as many as you want. Let us start with etopo1, which you can put as the key in the first indented line. If you are doing checks and want to run the code without the model being read in, but don’t want to rewrite the whole input file, just toggle activated to false.

class_name for all 3D geometric models will be StructuredGrid3D.

class_name for a 3D volumetric model on a structured grid will be StructuredGridV3D.

nc_data_file is also just the name of the 3D model file that you have created and saved in the /input folder.

For coordinates, this project will probably use LATITUDE_LONGITUDE and either DEPTH or RADIUS. Just make sure you get the latter correct: etopo1 is defined as an elevation with respect to some mean sea level whereas the Moho file we are using is defined as a depth with respect to sea level.

Again, you can use the ellipticity correction if you want to but it is unlikely to make an enormous difference.

If you’ve used the provided python scripts, you should not need to change the data_rank, length_unit, angle_unit, or factor flags from those we have given ([0,1], m, degree and 1.0 respectively).

The things that you will need to change are as follows. Close and careful attention is needed when changing these, as mistakes can be very challenging to debug. The order that we go through these is not the order that they appear in the input file, but we think it is more logical.

  • nc_variables and nc_var: these are the variable names used in the netcdf file. If you have kept the defaults, the first one should be [latitude, longitude] and the second one should be depth, moho and elevation for etop01 and C11, C12… for seismic anisotropy. Note that these are just the labels that we used in our netcdf file, they can be called anything and just because you set ‘depth’ or ‘radius’ as the name does not mean that the code will treat it as such unless you also set the flags below.

  • vertical should be DEPTH for the Moho and seismic anisotropy, and RADIUS for topography. Make doubly sure this is consistent, if you have files which you think are the wrong way round and do not want to remake them, you can always swap factor from 1.0 to -1.0.

  • undulation_range: interface should be the depth or radius of the interface you are undulating, as defined in the background 1D model. If you are using depth as your variable (as we used for the Moho), then this should be, for example, -24400.0, as the Moho is 24.4 km below the surface in PREM. If you are using radius as your variable, as we do for the surface topography, this should be 6371000.0. Remember that you will define the undulations as elevations above the boundary in consideration. So, at a subduction trench where the Moho is at 7 km below the surface, you would define the gridpoint for the Moho’s depth as +17400 (-7000 + 24400 = 17400 m). Similarly, in Tibet where the Moho is at 80 km depth, the gridpoint would be -55600 (-80000 + 24400 = -55600 m).

  • depth_below_solid_surface: This option is only relevant if there is an ocean or other fluid medium at the free surface and a solid medium below, and if depth is used as a vertical coordinate. In this case, if set to true, depth will be measured as depth below the ocean floor rather than as depth below the ocean surface.

    Lots of models in seismology (e.g., PREM) do not include the ocean, so if you add one on top, you need to either add one right on top (making the total radius 6371 km), or remove 3km of rock first (to keep the radius as 6368 km). This flag allows you to define whether your interface coordinates are with respect to the seafloor (true) or the ocean surface (false). You might also want to think about whether, when you read in a 3D model for subsurface features, the depths should be defined with respect to the solid surface in the 1D case, or the 3D case – as the position of the surface may have moved considerably in between these.

  • data_rank: This option allows for flexible handling of row-major and column-major input data. AxiSEM3D is written in C++ and hence uses column-major data storage by default. If the input data is also in column-major format (e.g. default python output) the data rank option should be set to [0,1,2]. If the data is in row-major format (e.g., default MATLAB and Fortran output) the order needs to be reversed, to [2,1,0]. This option also allows the user to easily swap axes if they need to debug their model setup.

  • undulated_geometry: This option is only relevant if geometric models (including ellipticity) and volumetric models are used at the same time, and it determines the order in which the models are read. As explained in the next section, geometric models effectively deform the mesh. If this option is set to false, the volumetric model is read in first and then deformed together with the mesh. If the option if set to true, the mesh is deformed first and the volumetric model is read afterwards, such that it is unaffected by the geometric model.

  • undulation_range: min_max tells the code in which range it is permissible to stretch or shrink elements to accommodate your deformed boundary. For example, if you set the Moho to be at 24 km and then undulate it, elements on either side need to deform to accommodate the undulated boundary. min and max define the lower and upper edges of the deformation zone, which should be at least twice as large as the undulation range. For example, in Crust 1.0 the Moho depth varies between 75 km and 7 km, so doubling the deepest depth and halving the shallowest one in this range would give you min_max = [3500.0, 150000.0]. Again, just make sure that the signs are consistent and correct: 3500 m is actually a shallower depth (i.e., the minimum deformed depth), rather than the maximum! Either way the interface needs to be contained within the bounds you have set, else the undulation will not work!

  • whole_element_inplane: This option is used to make 3D models adhere to the mesh geometry. If it is set to false, each point inside an element reads the volumetric model individually. If set to true, it is ensured that the model boundaries are placed at the edge of elements. Use this option if a 3D model borders on a solid-fluid boundary or a sharp velocity contrast.

For volumetric models, nc_variables*should be set in the order [lat, lon, depth] or [lat, lon, radius]. The names should match the names in the netcdf file.

Data rank should map to the order of the variables in the netcdf file. For example, if the netcdf file is in the order radius, longitude, latitude, then the data rank would be set as follows: data_rank: [2, 1, 0].

Note on Efficiency.
Before starting the time-loop, AxiSEM3D attempts to locally reduce 3D input models to 1D as well as to decrease the level of anisotropy. Hence, the performance of the simulation is automatically optimised and does not depend on the structure of the input data.

However, large volumetric models are a bottleneck in terms of working memory, since the entire model has to be read once before being distributed to the parallel nodes. If the model can be divided up into separate sections with the meshed background model in between this is highly recommended.

inparam.nr.yaml#

Fourier series introduction#

This input file, which we did not touch upon earlier, is used to set the azimuthal (i.e. variation in longitude) parameters of the simulation. You will recall that AxiSEM3D works by taking a 2D mesh and representing the third dimension through a pseudospectral (Fourier) expansion. The order of this expansion determines the complexity of the entire simulation in the azimuthal direction. Note that this does not just affect the degree to which 3D models are accurately represented, lower Fourier order expansions also reduce the azimuthal complexity of the seismic wavefield that you are going to record.

To understand this another way, the Fourier expansion in AxiSEM3D is akin to the representation of any normal mathematical function in terms of sines and cosines. The more terms you include, the more accurate the function’s approximation is. For an arbitrary function, an infinite number of terms is needed for an exact solution - but we can make do with a smaller (finite) number of terms Nu, Nr if we are happy for the solution to only be accurate to within a given degree. what is the conversion between these two? Not sure it is consistently defined…

How much smaller than infinity the number of terms/Fourier order Nr can be depends on a number of different factors:

  1. The complexity of the 3D model that you add to AxiSEM3D. If you have a very long-wavelength model with relatively smooth variations, like SEMUCB, you need far fewer terms than if you have something that is strongly varying and densely sampled, like ETOP01. For a 1D model, you never need more than Nr = 5.

  2. The complexity of the wavefield that you want to represent: in general, you will only add in complex 3D models if you want to see what impact they have on synthetic seismograms. However, if you wanted for some reason to use a complex 3D model but were not interested in the fine detail of the wavefield, you could get away with a lower Fourier order.

  3. The complexity of the source: as noted before, if you are only using a 1D model you can save a small amount of time by reducing the Fourier order needed to represent the source term. A second rank moment tensor (i.e. something from CMT) requires Nr = 5, a point force requires Nr = 3, an isotropic/monopolar pressure term requires Nr = 1. For most 3D models the value of Nr needed to capture their variations will be far larger than 5, so there is no saving to be made even if you only use a point force or pressure term.

Remember that the overall resolution of the simulation (i.e., the smallest scale that you can resolve meaningfully) is still limited by the mesh resolution, so increasing the Fourier order too much becomes pointless as you are not sensitive to these small spatial scales anyway. For this reason you should leave bound_Nr_by_inplane = true.

Choosing input parameters#

You will probably be wondering by now, ‘how small a Fourier order can I get away with?’ - the answer is of course that this depends on how accurate a simulation you want, and what much computational resources you have available - a higher Fourier order increases both runtime and memory demands. There are four options that you can choose in AxiSEM3D for how the Fourier order is calculated, which may help with this question. These are, type_Nr = :

  • CONSTANT: here, the same Nr is used throughout the mesh. As mentioned previously, for a 1D model with a moment tensor 5 is both necessary and sufficient.

  • ANALYTICAL: this is like constant, but with slightly finer user-control, in that you can choose how Nr varies through the interior of the planet. As an example, if you wish to investigate the effects of crustal structure on surface waves, you might want a high expansion in the mantle, but not in the core. This is thus a finer level of control than in CONSTANT, which can save you both time and memory.

  • STRUCTURED allows you to specify even finer degrees of control, by setting the fourier order on individual elements. You can thus save even more time and resource with respect to ANALYTICAL, but a more detailed knowledge of how to create the appropriate grid of points is needed.

  • POINTWISE allows you to go even further by having the code try and optimise the variation of Nr through the mesh. See below for more details.

Within ANALYTICAL, you can specify the variation of Nr with depth that you wish to use within the depth_dependent_AxiSEM3D_default option. Using the options control_depths and Nr_at_control_depths you can make a pointwise specification of how Nr varies, remembering to give your depths in metres.

For example, control_depths: \[0., 6371.\] and Nr_at_control_depths:\[100, 100\] would set Nr = 100 across the whole of the Earth’s radius. Changing Nr_at_control_depths to [100, 1000] would set Nr to 100 at the surface, and 1000 at the planet’s centre, with linear interpolation in between. You can add as many control points as you want into these arrays. For example, if you wanted to look at high-resolution crustal structure, you might do something like:

control_depths: \[0, 50e3, 200e3, 6371e3\] and Nr_at_control_depths:\[1000,1000, 100, 10\] which would give you Nr = 1000 in the top 50 km, decreasing to 100 in a linear fashion by 200 km depth, and then linearly again from 100 at 200 km to 10 the centre.

Wavefield scanning#

Wavefield scanning is an advanced feature of AxiSEM3D. It involves the code trying to ‘learn’ the most optimal Nr profile it can, using past runs as a guide. In a way, this lets the code do all the optimisation described in the steps above itself, such that the runtime and memory requirements are minimised.

This method works by analysing the frequency spectrum at all the interpolation points in the mesh. For a given point, if the energy above a given frequency is below some threshold value (which you can set), the code determines that frequencies above this point are not contributing to the overall solution, and the Fourier series can thus be truncated. In subsequent runs, the value of Nr at these points can therefore be reduced as they are determined to have a reduced azimuthal wavefield complexity.

If you want to try this, you can set enable_scanning to true and type_Nr = POINTWISE. Follow the instructions in the input file to make sure that you match the filenames and copy the relevant files across.