PyRaysum : Software for Modeling Ray-theoretical Plane Body-wave Propagation in Dipping Anisotropic Media

This article introduces PyRaysum , a Python software for modeling ray-theoretical body-wave propagation in dipping and/or anisotropic layered media based on the popular Fortran code Raysum . We improve and expand upon Raysum in several ways: 1) we signiﬁcantly reduce the overhead by avoiding input/output operations; 2) we implement automatic phase labeling to facilitate the interpretation of complex seismograms; 3) we provide the means to correct inaccuracies in the calculated amplitude of free surface reverberations. We take advantage of the modern, object-oriented Python environment to offer various classes and methods to perform receiver function calculation, ﬁltering and plotting. PyRaysum is fully backward compatible with legacy Raysum ﬁles and integrates well with NumPy and ObsPy , two standard libraries for numerical computing and seismology. PyRaysum is built in Python version 3 and requires a Fortran compiler, but otherwise runs on all platforms. The software offers a high-level, ease-of-use user interface and is equipped with complete documentation and testing as well as tutorials to reproduce published examples from the literature. Time-optimized post-processing functions allow for the straightforward and efficient incorporation of PyRaysum synthetic data into optimization or probabilistic parametric search approaches. Non-technical summary We introduce PyRaysum , an overhaul of the popular Fortran computer program Raysum for modeling receiver functions, wrapped in the python programming language. PyRaysum computes synthetic seismograms for subsurface models that consist of a few layers with constant elastic properties. The layer properties may be anisotropic and the layers may be dipping.


Non-technical summary
We introduce PyRaysum, an overhaul of the popular Fortran computer program Raysum for modeling receiver functions, wrapped in the python programming language. PyRaysum computes synthetic seismograms for subsurface models that consist of a few layers with constant elastic properties. The layer properties may be anisotropic and the layers may be dipping. Compared to the original code, PyRaysum is faster, more intuitive to use, and can easily be combined with other programs written in the Python programming language. These enhancements facilitate its usage in the estimation of subsurface properties from seismograms.

Motivation
Modeling teleseismic body-wave propagation in complex media is an important component of passive seismological approaches that aim to decipher upper mantle and crustal seismic velocity structure on the receiver side (e.g., receiver functions and teleseismic shearwave splitting analyses). Modeling wave propagation in highly heterogeneous and anisotropic D media can be performed using full-waveform approaches (e.g. spectral element or nite-di erence methods). Most often, however, a simpler parameterization of the subsurface velocity structure is desirable, because Earth structure is dominantly D or D at the scale of the upper mantle and lithosphere, and faster modeling methods allow searching for a wider range of model parameters that t the data. Among those approaches,

Figure 1
PyRaysum computes three-component synthetic seismograms for a plane wave that travels through a stack of dipping, anisotropic layers. This example shows a 2-layer over a half-space model, where the top and bottom layers are isotropic, and the middle layer is elastically anisotropic (represented by the anisotropic ellipsoid ( ani ) and orientation (plunge and trend ) of the symmetry axis). The top interface of the half-space is inclined, as parameterized by the strike and dip angles. The ray geometry is defined by the back-azimuth and horizontal slowness of the incident plane wave ray vector and the location of the station relative to the origin. seismic velocity structure within subduction zones (e.g., Audet and Bürgmann, ; Nicholson et al., ), orogenic belts (e.g., Schulte-Pelkum et al., ; Sodoudi et al., ), and collisional settings (Schneider et al., ). The original Raysum so ware is written in native Fortran for fast and e cient computations. The input and output to Raysum consist of formatted text les, and users have to write their own scripts for reading, writing, processing and visualizing the synthetic data. This makes it challenging and cumbersome for beginners to quickly produce synthetic data, explore the parameter space e ciently, and combine synthetic with observed data in optimization problems. Furthermore, reading and writing les to disk represent a substantial performance bottleneck for repeated and automated program execution. PyRaysum was developed to remedy these shortcomings, by wrapping Raysum in a modern Python environment with classes and modules to facilitate its usage for beginners and to streamline the modeling approach in optimization or probabilistic search approaches. It employs ObsPy (Krischer et al., ) for handling seismic data, NumPy (Harris et al., ) for data processing, and Matplotlib (Hunter, ) for visualization. The installation, testing, full user interface and Jupyter notebook tutorials are described in the online documentation at https://paudetseis. Listing 1 General structure of a PyRaysum setup. The parameters chosen here reproduce the synthetic data shown in Figure 2. github.io/PyRaysum/index.html. Here we give an overview of the user interface and provide examples that validate the results against published and synthetic examples. We provide timing benchmarks for various program execution options and suggest new applications that arise from the improved performance and transparency of results.

User Interface
Simulating seismic waveforms with PyRaysum is done by setting up: ) a subsurface model; ) the ray and station geometry; and ) a suite of run-control parameters. With this setup, synthetic seismograms can either be generated as ObsPy Streams or NumPy arrays. In Listing , the res1 and arr1 objects both contain the same synthetic seismograms. They di er in that res1 is a Result object that contains the synthetics as feature-rich ObsPy Stream s with additional metadata, while arr1 holds the bare synthetic samples as a NumPy Array . PyRaysum o ers common methods for receiver function post-processing for both objects, where the Ob-sPy-based routines focus on exploratory data analysis, while the NumPy-based ones have an emphasis on computational e ciency.
PyRaysum provides two packages: fraysum and pyraysum . fraysum bundles access to the underlying Fortran routines, which are based on the original Raysum code (Frederiksen and Bostock, ).
pyraysum provides the Python interface for computing, post-processing and plotting receiver functions. The object-oriented interface can be imported directly from pyraysum . Functions for advanced users are stored in three modules: prs allows users to create PyRaysum objects from les, including legacy Raysum les; frs exposes functions for fast, NumPy-based post-processing of fraysum output; plot provides plotting functionality.

fraysum
The fraysum package is generated during compilation of the Raysum code through the NumPy f2py interface generator. This package facilitates access to the run_full() and run_bare() functions, which are the low-level calls to the underlying Fortran code. run_bare() only returns the synthetic -component seismograms, while run_full() additionally returns arrays of times, amplitudes and identi ers of the converted and re ected seismic phases, at the cost of a longer run time. This module lacks any convenience and bookkeeping functionality.

pyraysum
The pyrasyum package exposes the primary function run() to compute synthetic seismograms. It de nes the three classes Model , Geometry and Control that organize the required input parameters and the Result class that holds the results.
run() Results are created through a call to run() , which requires an instance of each input parameter classes representing a subsurface model ( Model ), the ray-and station geometry ( Geometry ) and run control parameters ( Control ) as required positional arguments. The optional keyword arguments rf and mode provide switches to automatically compute receiver functions and skip automatic phase labelling.

Model
The subsurface seismic velocity structure is parameterized as a stack of layers ( Figure ), where each layer is described by its vertical thickness ( thickn , in m), density ( rho , in kg/m 3 ) and isotropic P-and S-wave velocities (vp and vs , in m/s), indexed from top to bottom. Optionally, the layer may be inclined (strike and dip angles in degrees with a right-handrule) and/or anisotropic ( ani is percent anisotropy; Figure ). Anisotropy is parameterized in a simpli ed hexagonal symmetry class as where V and V ⊥ are the seismic velocities parallel and perpendicular to the symmetry axis. P-and S-wave anisotropy are equal and pure elliptical anisotropy is assumed (Porter et al., ; Sherrington et al., ; Levin and Park, ). The orientation of the anisotropy axis is de ned by plunge (degrees down from horizontal) and trend (degrees clockwise from north) angles. Positive anisotropy refers to a fast axis of symmetry, whereas negative anisotropy denotes a slow axis. We note that, unlike all other model parameters that apply to the entire layer properties, specifying layer strike and dip angles refers to the orientation of the top interface. To dene a uniform-thickness dipping layer, the same strike and dip angles must be speci ed at the underlying layer.
The model layers can be accessed and manipulated by their index (see Section . for examples). Convenience methods for the manipulation of a Model instance include: parametrizing vp and vs in terms of VP /VS ( vpvs ); changing model attributes interactively using brief command strings ( change() ); adding, splitting, removing and averaging of layers ( + , split_layer() , remove_layer() , average_layer() ); plotting the subsurface model as a staircase diagram or pro le sketch (plot()); saving the model to le, including legacy Raysum model les (save() ).
Geometry The ray and station geometry are parameterized in terms of the ray back-azimuth angle ( baz , in degrees clockwise from north) and horizontal slowness ( slow , in s/km), and station o set in north and east direction from the model origin (dn and de , in m). Ray parameters can be speci ed as either oats to model single-event waveforms, or iterables to simulate multiple event recordings. Each ray can be accessed and manipulated by index, where the ray indices are associated with those of the three-component synthetic waveforms in the Result object generated from a call to run() (see below).

Control
The Control class controls the computation of the synthetic waveforms. Various options can be speci ed, namely: the number of samples ( npts ) and the sampling interval ( dt in s) of the seismograms; the polarization of the incoming wave-eld ( wvtype ); the order to which free surface reverberations are computed (multiples, mults ); whether only speci c phases should be computed ( set_phaselist ); the time-alignment (align ) and time-shi ( shift ) of the seismograms; the rotation of seismogram components to le -handed geographical (Z-N-E) or righthanded ray (R-T-Z or P-V-H) coordinate systems ( rot ); and the verbosity of the program execution ( verbose ). Default options exist for each of these parameters.

Result
The Result class holds the output synthetic seismograms, which can be accessed by ray index or as a list of -component ObsPy Stream s with the "stream" or "seismogram" keyword. If the seismograms are computed in a ray coordinate system ( rot equal to "RTZ" or "PVH" ), synthetic receiver functions can be computed on the y using the calculate_rfs() method. They are then stored as an additional list of -component (radial/vertical shear and transverse/horizontal shear) ObsPy Stream objects under the "rf" keyword and can be accessed as the second element of the returned tuple when indexing Result . All functionalities of the ObsPy Stream class are readily available. Convenience methods plot() and filter() allow plotting and ltering the streams or rfs attributes in a single command.
The arrival time, amplitude, phase descriptor, abbreviated phase name, and conversion name of each converted or re ected phase arrival are stored within the stats attribute for each trace within the seismogram streams. The phase_descriptors serve as unique phase identi ers throughout PyRaysum. They are strings that consist of paired indices and letters that fully describe the direction (up or down) and type of rays converted or re ected at speci c interfaces. P indicates a Pwave, S a (fast) S-wave and T a (slow) S-wave. Uppercase and lowercase letters indicate up-going and down-going rays, respectively. In the case of an isotropic medium, S and T arrive at the same time and may both carry some energy. Note that S and T do not imply a polarization, but are chosen as synonyms for S and S to avoid ambiguity with the layer indices. For instance, the phase descriptor of the P-to-S converted wave at the bottom of the topmost layer (index ) in a -layer over half space model ( Figure ) would be P P S. The conversion_names attribute abbreviates the phase descriptors by omitting equal wave types in adjacent ray segments, only indicating the layer indices on top of which a conversion has occurred. The conversion name of the example phase would be P S. Lastly, the phase_names attribute provides the shortest, yet ambiguous phase description, listing only the converted phase types, here PS. The unique set of all phases present can be retrieved with the descriptors() method.

prs
The prs module holds the object-oriented interface described above and additional functions to interact with it. Namely, read_geometry() , read_control() , and read_model() read saved input objects from le and allow the direct use of legacy Raysum les. equivalent_phases() returns the seismic phases that arrive at the same time as a given phase (see Section ).

frs
The frs module holds the functions used to interpret fraysum output. Most importantly, make_array() allocates an array suitable for the repeated postprocessing of similar waveform simulations. Postprocessing can be done with filtered_array() and filtered_rf_array() , the NumPy-based, computationally-e cient functions to compute ltered synthetic seismograms and receiver functions from the output of fraysum.run_bare() .

plot
The plot module holds the plotting functions used by the Result.plot() method. The direct function calls expose more options to customize the plots. stream_wiggles() and rf_wiggles() create plots of multiple-event seismograms or receiver functions with either -or -component panels, respectively, ordered by either back-azimuth or slowness. The function seis_wiggles() creates a line plot of single-event seismograms.

Examples
In this section, we present three simple work ows that showcase usage of PyRaysum: In Section . we forward model three-component waveforms through a simple layer-over-half-space model and compare them with observed data at station G.HYB in Hyderabad, India; in Section . we demonstrate how to interactively manipulate models and quickly examine the resulting receiver function signature; and in Section . we reproduce previously published synthetic receiver functions representative for coastal California, USA, where a strongly anisotropic layer underlies a crustal block.

Conversions and multiples in a Seismogram
Station G.HYB is located on a seismically transparent cratonic crust. It yields very clear receiver function data (Saul et al., ) that show the direct P-to-S Moho conversion (PS) and rst-order free-surface reverberations (PpS and PsS) arriving ∼ s, ∼ s and ∼ s a er the direct P-wave, respectively. These data are consistent with a km-thick crust with an average S-wave velocity VS of . km/s and a P-to-S wave velocity ratio VP /VS of . (Saul et al., ). Here we use seismic data recorded at station G.HYB for a magnitude M . earthquake that occurred on January , in the Philippines, as a test case for waveform modeling. This wave front arrives due east at the station, and its direct P-waveform exhibits a relatively simple, Ricker-II-like shape due to the large earthquake focal depth ( km; Figure ). We model these threecomponent waveforms with PyRaysum (Listing ) using the appropriate source-receiver geometry and the seismic velocity model proposed by Saul et al. ( ). In the Control parameters, no time alignment ( align=0 ) is applied and we specify a geographic coordinate system ( rot="ZNE" ). The synthetic waveforms reproduce the main phase arrivals, although the convolved source wavelet distorts this comparison. The phase names, times and amplitudes are the phase_names output by PyRaysum and facilitate the understanding and description of seismograms. A more complete example is part of the online PyRaysum documentation.

Interactive exploration of receiver functions
We next demonstrate how the e ect of changes in the subsurface structure on receiver functions can be explored. In Listing , we specify rays with evenly spaced back-azimuths and a constant horizontal slowness of . s km −1 . We set the rotation of the coordinate axes to the P-V-H system, which are oriented parallel to the P-, SH-, and SV polarization of the incoming wave eld as predicted by the ray back-azimuth, slowness, and isotropic velocities of the topmost model layer (Kennett, ). This rotation ensures that as much converted energy as possible is mapped to the radial and transverse components and removes the constant-amplitude, zerolag signal on the radial component. Receiver function calculation is straightforward with a simple argument rf=True in the call to run() . Figure a shows the receiver function signature of the isotropic crust modeled in the previous example with a horizontal Moho as the only interface. The PS-conversion is clearly visible on the radial component at s.
To explore the e ect of a possible dipping Moho, we next set a • eastward dip of the interface (Listing and Figure b). The e ect is an undulation of the Moho conversion in timing and amplitude with a period of • (so-called -θ variations). Converted waves from the west arrive earlier and with a lower amplitude on the radial component, due to the relatively shorter ray path and lower layer-orthogonal incidence angle. Conversely, conversions from the east arrive later and with a higher amplitude. Energy from northerly and southerly directions gets converted into the dip direction, evident from the positive and negative amplitudes on the transverse component. The P-coordinate vector is not aligned with the actual ray polarization, because the dipping interface results in an apparent slowness of the ray that is di erent from the actual slowness of the layer. Therefore, some energy appears on the radial and transverse components at time s.
Next we rotate the interface by • -so that it dips west-and increase its dip to • (Listing and Figure c), thus producing an interface orthogonal to that in the previous example. The undulation pattern of the conversion is more pronounced between and s. The greater misalignment of the P-coordinate vector with the ray polarization leads to signi cant converted energy on the P-component that gets mapped to the receiver function as an arti cial secondary pulse between and s through deconvolution.
To explore anisotropy in the topmost layer, we set the layer dip back to horizontal and instead vary the anisotropic parameters in the model (Listing and Figure d and e). With a fast symmetry axis (orthogonal to plane of slow axes) trending northward and plunging Interactive exploration of receiver functions. The sections correspond to simple layer-over-half-space subsurface models. An interpretation is given in Section 3.2. The code to produce the images is given in Listing 2.
• , the converted wave eld gets separated into a fast S -wave-arriving earlier, N-S-polarized parallel to the fast anisotropy axis, with a strong move-out patternand a slow S -wave, arriving later, E-W-polarized, without a move-out pattern characterized by • periodicity (i.e., -θ variations). Next we switch from fast-to slow-axis anisotropy and rotate the now unique slow anisotropy axis such that it is aligned with one of the in-clined slow anisotropy axes from the previous example (Listing and Figure e). The polarization of the S and S waves remains the same. Now the S -wave shows the move-out pattern, as the wave encounters azimuthallyvarying velocities.
Listing 2 Exploration of simple layer-over-half-space models with either a dipping interface or anisotropy. Definitions of Listing 1 are assumed.

Reproduction of earlier work
To demonstrate the handling of layer dip and anisotropy in multi-layered models and show additional plotting functionality of PyRaysum, we reproduce gure of Porter et al. ( ). In Listing , the de nition of the dipping lower-crustal layer model dipm implies that the top of the half space and that of the lower layer are dipping • , striking east. The layers of the anisotropic model anim are at. The bottom layer is characterized by % hexagonal anisotropy with a slow symmetry axis ( ani[1] = −20.) trending south, and plunging • down from horizontal. The layer con guration is visualized in Figure a and

Validation
We validate the synthetic seismograms created with PyRaysum by comparing them with those obtained using the matrix propagator method (Kennett, ; Thomson, ), as implemented in the Telewavesim package for Python (Audet et al., ), for the same seismic velocity models. The two test cases are the isotropiclayer case with free surface re ections ( mod1 in Listing ; Figure ) and the anisotropic lower crust case with direct conversions only (anim in Listing ; Figure b and d). Telewavesim seismograms were generated using the code given in Listing . Note that Telewavesim lacks the capacity of synthesizing seismograms for dipping layers and inherently always includes all theoretical phase arrivals within the time window. Additionally, Telewavesim results do not provide a good in nite frequency approximation and therefore need to be ltered for comparison with PyRaysum results.
For the isotropic -layer model, the Telewavesim and PyRaysum traces agree with a cross-correlation coecient of 0.9975, if multiples of the direct P-wave ( rstorder multiples) are considered (mults=2; Figure a). However, the amplitude of the PpS arrival di ers noticeably. This is the case because, with the mults=2 option, no second-order multiple is calculated. For instance, neither re ections of conversion (e.g., PSpP) nor conversions of re ections (e.g., PsP) are considered. As we will discuss below, this behaviour is deliberate, because the implicit inclusion of these phases results in a large overhead for multi-layered models, leading to long run times and eventually segmentation faults. The limitation of mults=2 to rst-order multiples ensures that the phase with the largest amplitude is present in the synthetic seismograms.
In the present case, the free-surface re ection of the P-to-S Moho conversion PSpP contributes significantly to the amplitude of the PpS phase. Such equivalent phases can explicitly be included using the Control.set_phaselist() method. The unique set of phase descriptors present can be retrieved from a Result object and be used to compute the additional equivalent phases (Listing ). With the inclusion of equivalent phases, the cross-correlation coe cient between the Telewavesim and PyRaysum seismograms improves to 0.9996 (Figure a).

Performance
We compared the run times of typical calls to Raysum with comparable calls to PyRaysum, as well as different processing options of PyRaysum, using an AMD EPYC™ GHz Central Processing Unit. We rst tested the run time of the current Rayusm version . with input models consisting of to layers over a half-space, for rays, calculating all rst order mul-Listing 4 Code illustrating the use of equivalent phases. Definitions of Listing 1 are assumed. tiples (mults=2). The model, geometry and parameter les were read from disk; seismic traces, arrival times and phase descriptor les were written to disk. The computation completed in about ms for the single layer model with a factorial increase to about s for the layer model (blue line in Figure ). At layers, Raysum encountered a segmentation fault, because the number of computed phases exceeded the maximum maxph , which is de ned during compilation. Using PyRaysum, the same results can be obtained faster. With NumPy array output and an otherwise identical con guration, the layer model completed in ms and the layer model in . s (teal line in Figure ). These numbers convert to a -to -fold decrease in computation time, primarily because the in-/output overhead is avoided.
We now illustrate the computational cost of the treatment of phase combinations (setting of Control.mults ) and the overhead of handling the metadata-rich ObsPy objects against the bare NumPy array. We use the ObsPy-based receiver function computation and ltering provisions as given in Listing and measure the run time of Listing with di erent options for Control.mults : Obs-M -F-RF Compute all conversions and rst order multiples: ctrl6 = Control(mults=2) Obs-M -F-RF Only compute direct conversions: Obs-M -F-RF Only compute two explicit phases: the direct P wave and one P-to-S conversion: ctrl6.set_phaselist(["1P0P", "1P0S"]) (for the -layer case). This implicitly sets ctrl6.mults=3 .
The pink-shaded lines in Figure  Run times on the order of a second are usually acceptable when executing code a few times manually. For frequent and automatic calls that do not require bookkeeping of metadata, results can be obtained faster by directly calling the Fortran routine run_bare() . We next illustrate the computational cost of three NumPy-based post-processing options, keeping the number of phases constant as in the last case examined above. The post processing steps timed are (Listing ): Num-M -F-RF Compute seismograms. Then compute receiver functions trough spectral division and lter them using filtered_rf_array() .
Num-M Only compute synthetic seismograms.
In Listing , rfarray holds the processed data. It is allocated once before the (possibly subsequent) calls to filtered_rf_array or filtered_array .
For all three cases, run time is constant with respect to the number of model layers (green-shaded lines in Figure ). Compared to the ObsPy-based post processing, time spent for ltering and spectral division (Num-M -F-RF) is approximately halved, with only instead of ms, implying that as much time is required for bookkeeping. Sparing spectral division (Num-M -F) saves about / th of run time, or ms in the present example. Filtering is computationally cheap, taking only about / th of time, or ms (Num-M ).

Outlook and Future Work
With the time-e cient, NumPy-based post-processing, PyRaysum can be used as a forward code in parameter estimation problems. For example, the common problem of nding the optimal crustal thickness and bulk VP /VS ratio from the time and amplitude of the Moho conversion and reverberations in receiver functions (e.g. Zhu and Kanamori, ) can be re-formulated as a minimum mis t problem and be generalized to multilayered models. First tests indicate that, e.g., SciPy's dual_annealing() function can nd a minimum mist model from teleseismic receiver functions for the thickness, VS, and VP /VS of three layers representative for the continental crust and subducting slab of the Cascadia subduction zone within a few hours. Automatic phase labels (Figures and ) can help to identify more multiples in complex receiver functions and facilitate a more thorough understanding of the scattered wave eld. Anisotropy is currently parameterized as a single parameter, ani (Equation ). Internally, Raysum handles anisotropy in its most general form using the Christo el equation. This in principle allows to explore other theoretical crystal symmetry classes (e.g. general hexagonal, orthorhombic), as well as anisotropy predicted and measured for speci c minerals and rocks (e.g. Brownlee et al., ), which can be implemented through a new de nition and internal handling of the Model class.

Conclusion
PyRaysum is a modern and fast incarnation of Raysum (Frederiksen and Bostock, ), the popular algorithm to compute seismograms that result from the plane wave propagation through dipping, anisotropic, layered media. New features include bypassing hard disk readwrite operations, automatic labeling of seismic phases, easy de nition of explicit phase lists, inclusion of equivalent phases, ObsPy integration, e cient receiver function post processing, and interactive manipulation of subsurface models. With these improvements, PyRaysum makes it possible to play with receiver functions in a simple and e cient Python environment, as well as to invert for subsurface properties using state-of-the-art inverse modeling algorithms.

ACKNOWLEDGMENTS
The authors thank Andrew Frederiksen for extensive discussion, and for his continued work on Raysum; and Michael Bostock for discussion. First and foremost, please cite their work (Frederiksen and Bostock, ) when using PyRaysum. The comments of two anonymous reviewers helped to improve the manuscript and the initial major release of PyRaysum. WB was located on the traditional, ancestral, and unceded territories of the Musqueam People while performing this research. PA acknowledges and respects that his workplace stands on unceded Algonquin territory. This work is funded by the German Research Foundation, grant BL

Code availability
PyRaysum is open-source so ware and is under active development on GitHub (https://github.com/ paudetseis/PyRaysum, last accessed January ). A snapshot of the most recent version . . is hosted on zenodo (Audet and Bloch, ). The full documentation is included with the current snapshot and online at https://paudetseis.github.io/PyRaysum/ index.html. Raysum is available at https://home.cc. umanitoba.ca/~frederik/Software/.