Simul2017: a flexible program for inversion of earthquake data for 3-D velocity and hypocenters or 3-D Q

23-Feb-2024 Please use the updated version Simul2023 at zenodo.org/records/10695070

      The program 'simul2017' provided in this zenodo site is a flexible program for simultaneous inversion of travel-time data for 3-D velocity and hypocenters.  We have used it extensively since the 1980s, and it has evolved over that time with modifications to the original code of Thurber (1983).  Eberhart-Phillips (1990) developed the Simulps code to use P and S travel- times to solve for Vp and Vs.  From Simulps, other groups have evolved simul versions with useful modifications (Haslinger and Kissling, 2001), which are not included in simul2017.  All S travel-time ray-paths are calculated using the 3-D Vs structure in the Eberhart-Phillips (1990) code and all subsequent modified codes such as Eberhart-Phillips and Reyners (1997, 2012) and Eberhart-Phillips and Fry (2017).  
       
       These parameterize velocity on a 3-D grid of nodes with velocity linearly interpolated between nodes, allowing flexible gridding through linking of nodes (Thurber and Eberhart-Phillips, 1999).  The solution is obtained by damped least squares, with no other smoothing applied, to obtain a model that has few artefacts and stays close to the initial model where there is low resolution.  For typical earthquake travel-time data, the Vs model is poorly constrained relative to the Vp model, less representative of crustal structure and difficult to use for interpreting Vp/Vs (Eberhart-Phillips, 1989).  Thus, as described by Eberhart-Phillips and Reyners (1997), it is preferable to solve for Vp and Vp/Vs, when using local earthquake travel-time data.  This parameterization is retained for group velocity, with the Herrmann (2013) Vp kernels related to Vp model inversion parameters and Vs kernels related to partial derivatives of Vp and Vp/Vs model parameters (Eberhart-Phillips and Fry, 2017; Eberhart-Phillips et al., 2022).
       
       Rietbrock (2001) implemented Q inversion in simul, using t* observations which describe path attenuation from spectral decay.  This has been extremely useful in many settings, since Q has ability to characterize many features more distinctly than velocity (eg: Eberhart-Phillips et al., 2020; Eberhart-Phillips, 2016).  It can be used for Qp or Qs, and does not need to solve for hypocenters. The code was extended for joint teleseimic and local inversion for velocity or Q (Eberhart-Phillips and Fry, 2018).
       
       The simul2017 code is a complete fortran code.  The beginning of the code has descriptions of input parameters.  An old users manual for an early version may be of some use (Evans et al., 1994). Our work has taken the approach of adding new features, without eliminating previous options.  There are examples in the examples directory.  For further knowledge of input formats, the input subroutines can be reviewed. The include file has parameters related to array size and the user should consider adjusting these. There should be no need to edit the fortran code.  A notable input parameter is 'nitmax', specifying the number of iterations, relocation only (0), or create synthetic data (-1).  The coordinate system model cartesian coordinates and distances are computed with Transverse Mercator conversion, and earth-flattening transformation is used. The coordinates are defined at the beginning of input file 2 (stations), and older options of short-distance conversion or NZMG49 are allowed.  The model 0-km depth equates to a present sea-level datum, the z= -1 km depth grid is 1 km above sea-level, and velocity is linearly interpolated between gridded nodes. Travel-time ray-tracing includes station elevations.  There are several options of travel-time data format, selected by input parameter kttfor, and used in subroutine input4.