This folder contains results of a surface wave anisotropy study, including the results of a joint inversion with SKS splitting constraints. The files are arranged according to the figure in which those results are presented. All files are .mat format files, and almost all contain MATLAB "structures", with fields that hold the data. Each file is explained below. Capitalised words in filenames below are wildcard stand-ins for the names of the different geographic regions described in the text. Fig1_phi_histogram_REGION.mat - contain histograms of the anisotropic fast azimuths (phi) in each regionat all measured periods. The histograms are generated by taking 10000 random samples of the cosine and sine coefficients of the phase velocity at those periods. Cosine and sine terms are measured with uncertainties (see file Fig2_Fig3_SW_anisotropy_allregions.mat), but due to non-linearity in the computation of fast azimuths from these two terms, the best way to estimate uncertainty in the azimuths is to sample cos/sin and compute phi for each sampled pair) - load('Fig1_phi_histogram_REGION.mat') will load a structure, "anis_obs" with the following fields: - "swperiods" = [15x1] vector of periods of the Rayleigh waves measured - units are in seconds - "phibins" = [37x1] vector of bin edges used to compute the histogram - units are in degrees from North - "phihist" = [15x36] matrix of counts of samples in each bin, at each period. Fig2-3_mean_phase_velocity.mat - contains a [15x1] vector of Rayleigh wave phase velocities averaged across the study region. Units are in km/s. Values correspond to each of the periods specified in Fig2_Fig3_SW_anisotropy_allregions.mat or Fig1_phi_histogram_REGION.mat. Fig2-3_SW_anisotropy_allregions.mat - contains measurements of surface wave anisotropy at each period in each region, and the results of linearized 1D inversions of these measurements for 1D profiles of anisotropy as a function of depth (these inversions are of the surface wave data alone, prior to inclusion of SKS constraints) - load('Fig2-3_SW_anisotropy_allregions.mat') will load a structure array, "SWanisres". This structure array has 10 dimensions, one for each sub-region in the study. - for each sub-region [E.g. for the third region, type "SWanisres(3)"] you will see the following fields - "name", "lat", and "lon" give the name, and centroid locations of each sub-region - "period" is a [15x1] vector of periods of the Rayleigh waves measured - units are in seconds - "cos" and "sin" contain the results for the cosine and sine 2-theta terms. These fields are themselves structures, which contain the following fields: - "Vs_anis" - another structure that contains the surface wave measurements at each period: - "dVs_obs" - [15x1] vector of observed sin/cos terms (a1, a2 in the paper) at each period - units are km/s - "dVs_std" - [15x1] vector of estimated standard deviations for the sin/cos terms - units are km/s - "dVs_pre" - [15x1] vector of predicted sin/cos terms using the final model of the linearised inversion - units are km/s - "dVs_resid" - [15x1] vector of observed-predicted misfits for sin/cos terms - units are km/s - "import" - [15x1] vector of importance, a measure of significance - "dVsmod" - another structure that contains the surface wave 1D inversions: - "z" - [25x1] vector of depths of layer centres at which model is defined - units are km - "Vs" - [25x1] vector of inverted sin/cos terms in each layer - units are km/s - "Vs_std" - [25x1] vector of inverted sin/cos standard deviations - units are km/s - "ranki" - [25x1] vector of rank of each layer - a measure of resolution - "rank" - total rank of the model = a measure of the number of effective model parameters - "smthco" - smoothing coefficient for this inversion - RMS - root mean squared error - std - standard deviation - so, if you wanted to plot the observed (blue circles) and modelled (red line) dispersion curve for cosine terms for the sixth sub-region (Juan de Fuca plate interior), you would enter: >> clf; hold on >> plot(SWanisres(6).period,SWanisres(6).cos.Vs_anis.dVs_obs,'ob') >> plot(SWanisres(6).period,SWanisres(6).cos.Vs_anis.dVs_pre,'-r') - or, if you wanted to plot the inverted profile of sine terms as a function of depth for the sixth sub-region (Juan de Fuca plate interior), you would enter: >> clf; hold on >> plot(SWanisres(6).sin.dVsmod.Vs,SWanisres(6).sin.dVsmod.z,'-k') >> set(gca,'ydir','reverse') Fig4_kernels.mat - contains shear velocity kernels used in 1D surface wave inversion - load('Fig3_kernels.mat') will load a structure array, "K". This structure array has 15 dimensions, one for each period - for each period [E.g. for the fourth period, type "K(4)"] you will see the following fields - "period" gives the period, in seconds - "z" - [25x1] vector of depths of layer centres at which model is defined - units are km - "Kval" - [25x1] vector of sensitivity values relating fractional (i.e. unitless) changes in shear velocity at that depth to fractional changes in phase velocity at that period - units are 1/km Fig4_layered_anisotropy.mat - contains the results of the 1D surface-wave-only inversion but constrained to have uniform structure within three layers. - load('Fig3_layered_anisotropy.mat') will load a structure, "anislay", which contains the layerised model for all ten sub-regions - this structure has the following fields - "names" - a {10x1} cell array of region names - "top" - a [3x10] matrix of (3) layer upper boundaries for each (10) region (all are the same: 11,40,90) - units are km - "bottom" - a [3x10] matrix of (3) layer lower boundaries for each (10) region (all are the same: 40,90,210) - units are km - "lons", "lats" - [10x1] vectors of sub-region centroid locations - "cV2theta"/"sV2theta" - [3x10] matrices of coefficients for the cosine ("cV...") and sine ("sV...") 2-theta terms in each of the three layers, for each of the ten regions - units are km/s - "cV2std"/"sV2std" - as for the ...2theta terms, but these are standard deviations - units are km/s Fig5-6_joint_inv_REGION.mat - contains the results of the joint surface wave and SKS inversion - load('Fig4_joint_inv_REGION.mat') will load the following six structures: - "SKSdat" is a structure containing the SKS data for this region from Bodmer et al. 2015, including the following fields: - "avanis" - another structure containing the mean measured splitting at stations in this regions - "allstaanis" - complete results from Bodmer, with every measured split at every station in the region - "ev2use" - list of the 16 most commonly recorded events, for use in the forward modelling - "SWmod" full surface wave model, with the following fields - "SWanisres" - replica of the contents of Fig2_Fig3_SW_anisotropy_allregions.mat, but only for the one region - "Vmod3d" - 3D model of isotropic velocities for the area. This is actually a 90x90 structure array, each structure is a 1D profile at a single lat,lon point (there is a grid of 90 positions in latitude and 90 in longitude that spans the area) - "anislay" - a replica of the contents of Fig3_layered_anisotropy.mat, but only for the one regions - "covm" - a 6x6 covariance matrix for cosine (upper left 3x3 cov matrix) and sine (lower right 3x3 cov matrix) terms in the three layers for the layerised inversion - "icovm" - as for "covm", but the inverse covariance (computed once at the outset to save time) - "mav" - [6x1] vector comprising [3x1] cosine terms and then [3x1] sine terms that are the average values (in km/s) in each layer. The same as [anislay.cV2theta;anislay.sV2theta] - "par" - a structure that contains all the parameters from this run of the inversion - not particularly useful without the codes - "prior" and "posterior" - ensembles of models in the prior (empirically derived by running the algorithm with no likelihood term) and posterior. These have some of the following fields: - "Niter" - the total number of iterations - "Nstored" - the number of stored models - "logLike" - the log (base e) of the likelihood of each stored model - "logPrior" - the log (base e) of the prior probability of each stored model - "cterm", "sterm" - [Nstored x 3] matrix of cosine and sine terms for all 3 layers, for all stored models (units are km/s) - "anisa", "anisp" - [Nstored x 3] matrix of anisotropy amplitude and fast azimuth (p is for 'phi'), respectively, for all stored models (just computed from sine and cosine terms, above) - "E_T" - [Nstored x 16] matrix of transverse energy calculated for each of the 16 earthquakes used for the forward modelling step - "final_model" - final model from the posterior distribution. This is a structure with the following fields: -"cterm", "sterm", "anisa", "anisp" which are, respectively, cosine, sine, amplitude, and azimuth (p is for 'phi') terms - all terms have the following fields: - "av" - [3x1] vector of average values for each layer - "minmax" - [3x2] matrix of maximum (first column) and minimum (second column) values for each layer - "sig1" - [3x2] matrix of upper (first column) and lower (second column) bounds for 68th percentile (= 1 sigma, if gaussian) - "sig2" - [3x2] matrix of upper (first column) and lower (second column) bounds for 95th percentile (= 2 sigma, if gaussian) Fig7_forwardmodel_REGION.mat - contains the results of sampling the final model 500 times in each region and forward modelling to illustrate predicted splitting. - load('Fig5_forwardmodel_REGION.mat') will load the following variables: - "par" - a structure that contains all the parameters from this run of the inversion - not particularly useful without the codes - "phi_SC" - [500x16] matrix of fast azimuths (in degrees from North) for all 500 samples calculated by forward modelling the splitting and then inverting using the Silver and Chan "minimum energy" method for the 16 earthquakes specified above. - "dT_SC" - [500x16] matrix of splitting times (in seconds) for all 500 samples calculated by forward modelling the splitting and then inverting using the Silver and Chan "minimum energy" method for the 16 earthquakes specified above. - "Emaps" - [91x41x500x16] matrix of the energy surface computed for the forward modelling using the Silver and Chan "minimum energy" method and grid searching through all possible combinations of splitting parameters. The 91 rows correspond to azimuths on the range [-90:2:90], and the 41 columns correspond to splitting time (units of seconds) values on the range [0:0.1:4]. The 500 layers are for the 500 samples of the posterior, and the 16 layers correspond to the 16 earthquake arrivals for which the computation was conducted. - "wts" - a [500x1] vector of weights for each of the samples based on the measured posterior likelihood of that sample.