Functional model doesn't support HTDP

[ccrook]

The HTDP point motion model used in the continental US employs a couple of algorithms not implemented in the current functional model definition. HTDP is distributed as Fortran code in https://geodesy.noaa.gov/TOOLS/Htdp/Htdp.shtml. All the data, including grids, is implemented as Fortran static data.

The algorithm implemented in the code is as follows (TBC):

• Calculate a plate velocity
  • Identify plate/region using a point in polygon routine. Each region either has a velocity grid or X,Y,Z rotational velocities.
  • If a gridded region then interpolate velocity from grid
  • Otherwise calculate velocity from X,Y,Z rotational velocities
• Apply coseismic displacements. For each earthquake the deformation is defined by a geophysical fault plane model and a radius of influence. Check if the point is within the radius and the epoch is after the eq. If so evaluates the model at the point using the Okada (1985) formulae and adds to the total.
• Apply post-seismic displacments. These are defined by gridded models each with a exponential decay time function. For each grid check check if the epoch is after the eq and the point is within the grid extents. If so evaluate interpolate the grid, multiply by the time function, and add to the total displacement

This is a very different implementation from the current DMFM. In particular it adds requirements for

• polygon definition of the extent of a componet, used to determine if it applies at a location using a point in polygon routines to identify
• Definition of a spatial model using the Okada formulae for calculation of the deformation due to a fault model

The Okada formulae define displacement of an infinite uniform elastic half space due to a uniform slip on a rectangular fault plane. They are widely used in geophysical modelling of fault motion. A common approach in geophysical analysis of fault movement is to divide a fault plane into multiple rectangular subfaults each with its own slip vector. This provides an approximation to a varying slip vector across the fault plane. In HTDP the models range from 1 to 911 subfaults.

There are three complexities in adding the Okada formulation to the functional model:

• the interpolation coordinate system must be mapped to a flat earth coordinate system used as the half space in which the model is evaluated. Typically geophysical software may use one of many approximate formulae to do this (eg scaling latitude and longitude offsets from a reference point as in HTDP).
• there is no common convention used to parameterise the fault model. For example the fault location may be defined the centre of the fault, the centre of the top of the fault plane, one corner of the fault plane. The slip vector may be defined as dip slip, strike slip components (and possibly a component perpendicular to the fault plane) or by the length and rake of the slip vector. And so on!
• the Okada formulae are quite complex - much more so than any of the other formulae used - which will be a disincentive to implement the DMFM, and also potentially make it less efficient to evaluate. In particular where there are a large number of subfaults the formulae must be evaluated for each subfault.

The incentive for using the Okada formulae is that it is able to represent a good first order approximation to the fault movement in 8 or 9 parameters, which is very compact compared to a comparable nested grid. A counter argument is that it is generally not a very good approximation, particularly near the fault, and that it is inefficient to calculate.

While it would be possible to build a nested grid structure to represent the HTDP model to an arbitrary level of precision (except if it includes surfaced fault models - TBC) they still would not match the model exactly.

In summary - if we want the DMFM to support the authoritative HTDP definition then it needs these added functions. If it is enough to replicate HTDP to an arbitrary level of accuracy then the current DMFM can do that, possibly with the addition of a polygon extent definition, which may be sensible in any case.

Reference: Okada Y (1985) Surface deformation due to shear and tensile faults in a half-space. Bull Seismol Soc Amer 75:1135-1154

[demiangomez]

Hi @ccrook,

Maybe I'm mistaken, but I think that HTDP is implemented in a similar way as the trajectory prediction model for Argentina. We have an externally-calculated grid which is interpolated (say, linearly) to estimate the jump at a given coordinate. As you mention, in this case there is an issue with discontinuities. I think that incorporating the Okada formulation to the functional model is not necessary. Also, as you noted, Okada uses a halfspace not a sphere as the underlying coordinate system. There are multiple situations when this is ok (if the earthquake is small enough) but it is not appropriate for large earthquakes. So, there could be compelling arguments to use other dislocation models, such as Politz (1996) which uses a sphere. Thus, I think Okada (or any other dislocation model) should not be incorporated in the functional model.

Best,

[ccrook]

From discussion at OGC meeting 22 Feb 2021:

CP: HTDP is to be deprecated CP: Possibly limited utility in exactly replicating deformation. RS: Trans4D? CP: Still under development. PM: HTDP->TRANS4D. 3D+time.
PM: Still in discussion at NGS. Splitting off deformation model from other transformations.

In summary:

HTDP is being replaced with Trans4D. This is the same as HTDP 4.0, and appears to be functionally very similar to previous HTDP apart from 3d support.

NGS are still evaluating options for HTDP and so resolution is to keep a watching brief on HTDP developments but not to consider modifying our functional model specification to support HTDP implementation until NGS have resolved the future path.

This will is therefore closed for now!