\documentclass[agupp]{aguplus}
%\documentstyle[12pt,agupp]{article}
\lefthead{LEVIN, MENKE, PARK}
\righthead{MULTILAYERED ANISOTROPY}
\received{ }
\revised{ }
\accepted{ }
\paperid{}
\cpright{AGU}{1998}
\ccc{}
\authoraddress{V. Levin, J. Park, 
Dept. of Geology and Geophysics, 
Box 208109, Yale University, 
New Haven, CT 06520
(email: vadim@hess.geology.yale.edu;
park@hess.geology.yale.edu)}
\authoraddress{W. Menke, LDEO, Palisades, NY
email: menke@ldeo.columbia.edu}
\slugcomment{article submitted to JGR August 1, 1998 }
\setcounter{secnumdepth}{0}

%\input{update.tex}

\begin{document}

\title{Shear-Wave Splitting in the Appalachians and the Urals: A Case for
Multilayered Anisotropy}

      
\author{Vadim Levin}
\affil{Yale University, New Haven, Connecticut}
\author{William Menke}
\affil{LDEO, Palisades, NY}
\author{ Jeffrey Park}
\affil{Yale University, New Haven, Connecticut}
\date{ \today }


% Levin and Park 1997a - GRL
% Levin and Park 1997b - GJI
% Levin and Park 1997c - PAGEOPH

\begin{abstract}
%\baselineskip=20pt

Observations of shear-wave splitting in the Northeastern US
Appalachians and in the foredeep of the Urals vary significantly with
the back azimuth and incidence angle of the incoming phase.  
These variations suggest two or more layers within the upper mantle
with different anisotropic properties.
Synthetic seismograms for simple multilayered
anisotropic structures show that shear-wave splitting parameters tend
to vary substantially with the direction of approach.  Relying on a
subset of back-azimuth and incidence angle may strongly bias the model
inferred, especially if the observations are averaged. On the other
hand, the azimuthal splitting pattern provides additional constraints
on vertical or lateral variation of anisotropic properties in the
Earth.
Using a new error-estimation technique for splitting, we find that
individual measurements from broadband data have
errors on the order of $\delta \phi = 3\deg - 7\deg$ for the fast
direction, and $0.1-0.2$s for the delay of split shear waves.
The azimuthal variation of splitting parameters is broadly
consistent throughout the Appalachian terranes in the Northeastern US,
especially for two long-running stations in the Northeast US -
HRV and PAL. Observations can be separated into two distinct populations, with mean fast-axis azimuths of $N60\deg E \pm 4\deg$
and $N119\deg E \pm 2\deg$. Delay values within each population range from near-zero to $\sim 1$s.
Azimuthal splitting variation  for 
ARU in the foredeep of Uralian mountains is characterized by sharp
transitions between different groups of observations.
Using synthetic seismograms in simple structures, we
develop one-dimensional anisotropic models  under
stations HRV and ARU.  The model for HRV includes two layers of
anisotropic material under an isotropic crust, with fast-axis azimuths
$N53\deg E$ and $N115\deg E$ for the bottom and the top layers,
respectively.  The model for the upper mantle under ARU includes a
layer with a fast-axis at $N50\deg E$ atop a layer with fast-axis
azimuth $N90\deg E$. Our modeling confirms the need for a layer of strong
anisotropy with a slow axis of symmetry in the lower crust under ARU,
reported by {\it Levin and Park} [1997a].
Our results suggest that both Urals and Appalachians possess a relict
anisotropy in the tectosphere, associated with past continental
collision and accretion, underlain by anisotropy with orientation
similar to the local absolute plate motion, suggesting an
asthenospheric component to the signal.


\end{abstract}

\section {Introduction}

Seismic anisotropy, the dependence of seismic velocity on the direction
of propagation, is a common feature of modern global seismic velocity
models [{\it Dziewonski and Anderson}, 1981] and a well documented property of
olivine, the mineral that composes the bulk of the upper mantle
[{\it Crampin et al.} 1984]. A prominent effect of seismic anisotropy is the
so-called splitting of shear waves. Unlike the two shear modes in
isotropic materials, which have equal phase velocities, the two
corresponding modes in anisotropic materials have different
velocities.  As it enters an anisotropic region, a single linearly
polarized shear-wave pulse will split into two pulses of different
velocity. The polarizations of these pulses are related to the
projection of their propagation direction onto the axes of
the anisotropic elastic tensor [{\it Aki and Richards,} 1981].

Shear-wave splitting studies measure and describe the
anisotropy of the Earth.  One selects shear-wave phases that are known
to be linearly polarized prior to entering the anisotropic study
region, and measures their particle motion after traversing the
region.  One seeks a coordinate rotation that
separates the particle motion into distinct "fast" and "slow" pulses,
each of identical shape and linearly polarized in
mutually-perpendicular "fast" and "slow" directions.
The delay between the two pulses is proportional to the strength of the
anisotropic effect, which depends both on the intensity of seismic
anisotropy and the length of the path within the anisotropic material.
The axes of the rotated coordinate system provide information on the
symmetry and orientation of the anisotropic elastic tensor.

In some simple cases,
such as when the tensor has hexagonal symmetry with a horizontal symmetry axis,
the orientation of the fast splitting axis is approximately parallel
to the symmetry axis, regardless of propagation direction.  
This simplification is
often assumed in studies of the upper mantle. Estimates of splitting
time and fast axis direction from many shear waves 
at a given station have been averaged to estimate anisotropic strength
(i.e. delay time) and symmetry-axis azimuth for the mantle beneath that station [e.g. {\it Vinnik et al.}, 1995; {\it Barruol et al}, 1997; {\it Wolfe and Silver}, 1998; {\it Fouch et al.}, 1999]. 
Such "station means" are useful in tectonic settings where uniformity
of the fabric in the lithosphere is likely, e.g. on the ocean floor
[{\it Wolfe and Solomon}, 1998] or in the wake of a hot spot [{\it
Schutt et al.}, 1998].

Station means may be quite misleading, however, in cases where both the
delay time and the fast direction vary significantly with the
propagation direction.
Such variation can occur when an anisotropy tensor is inclined from the vertical, or has a more complicated symmetry, or both (Figure 1). {\it Babuska et al.}, [1993] discusses possible scenarios that involve an inclined orientation of hexagonal and orthorhombic anisotropic tensors. 
This possibility has also been considered by {\it Plomerova et al},
[1996], {\it Levin et al.} [1996], {\it Hirn et al.}, [1998] and
others.
Also, a combination of two or more layers of anisotropy with hexagonal
symmetry and horizontal symmetry axes leads to a systematic variation
of the splitting parameters with the polarization of the incoming shear
wave [{\it Silver and Savage}, 1994; {\it Vinnik et al}, 1995].
For vertically incident shear waves a simple analytic expression
describes the variation of splitting parameters in a simple 2-layer
model, predicting a ${\pi / 2}$ periodicity.
2-layer models with horizontal-axis anisotropy has been invoked to explain
back-azimuth variations in splitting parameters [e.g., {\it Ozalaybey and Savage}, 1994; {\it Russo and Silver}, 1994; {\it Granet et al.}, 1998].

In this paper, we demonstrate that the fast direction and delay
associated with shear waves that sample the Earth's upper mantle
beneath two long-lived mountain belts vary strongly with shear-wave
propagation direction.
We interpret these variations in terms of multilayered anisotropy,
implying either complex deformation in a past collisional event
or, more plausibly, a mix of active and fossil deformations.

The parameters of a split shear wave can be estimated by a grid search
over possible time-delays and fast-axis directions, and using some kind
of goodness-of-fit criteria to select the "best" set of values.
Two criteria have been used: 1) maximal similarity in the pulse shapes
of the two rotated seismogram components, as quantified by
cross-correlation [e.g., {\it Bowman and Ando}, 1984; {\it Iidaka and
Niu}, 1998]; and 2) that the re-assembled "original" pulse has maximal
rectilinearity, as quantified by the  ratio of the rectilinear and
elliptical motion [{\it Kosarev et al.}, 1984; {\em Silver and Chan},
1991].
These two methods give the same results when tested on nearly noise-free data.
Owing to their different treatment of noise, and to complications
induced by multilayered anisotropy, "cross-correlation" and
"rectilinearity" measures can give substantially different results when
applied to noisy data.

If the anisotropic material is homogeneous, an
observed split shear wave is exactly the sum of two pulses of different
polarization, one delayed with respect to the other.  
If, in contrast, the anisotropic material consists of several layers
(or, more generally, 3-D domains) of different anisotropy, then the
observed seismogram has a more complicated form, with a sequence of
pulses corresponding to mode-conversions from the various layer
interfaces.
In general, no rotation exists in which one component of the seismogram
is exactly a delayed version of the other (Figure 2).
Given a high-quality broadband waveform with no interfering seismic
signals, these conversions could perhaps be individually identified and
modeled.
Unfortunately, most $SKS$ data in studies of the upper mantle is
low-passed, and resolving closely-spaced sequence of
pulses is problematic.

Our approach is to retain the two parameter (fast direction and delay)
description for shear-wave propagation in anisotropic media, but to
recognize that this is an "apparent" measurement, without exact
correspondence to an underlying physical process.
This approach was introduced by {\it Silver and Savage}, [1994]
 for the case of two anisotropic layers with horizontal symmetry axes,
and more recently expanded to the case of a smoothly varying medium by
[{\it Rumpker and Silver}, 1998].
The apparent splitting parameters (Figure 3) contain significant
information about the anisotropic medium.
Most importantly, the apparent splitting parameters are
different from what one would expect for a homogeneous medium with the
same "mean" anisotropy, so that some information on the
depth-dependence of the anisotropy is preserved.
Unfortunately, owing to interference between the mode-conversions, the
measured values of the apparent parameters are somewhat sensitive
to the frequency band of the seismic data (Figure 4).
{\it Rumpker and Silver}, [1998] show
that even for a case of vertical incidence in a flat-layered model with horizontal axes of hexagonal symmetry apparent splitting parameters exhibit strong dependence on the ratio between the cumulative splitting effect of the medium and the frequency content of the shear wave.
This sensitivity does not
present any fundamental problem when modeling the apparent
splitting.  One simply compares observed apparent values with
predicted ones that have been computed from synthetic seismograms with
the same frequency content.  However, it  makes difficult the
comparison of data collected by different authors using different
processing schemes.

In this paper we measure apparent splitting parameters at a few
locations within two Paleozoic mountain belts - the Appalachians and
the Urals.  We focus on the anisotropic structure of the upper mantle,
and so use mostly SKS and SKKS phases.  We mainly use
stations with long duration of operation, and thus we are able to obtain
measurements from a wide variety of azimuths and angles of incidence.
As we show below, the strong directional behavior of these data
indicate that the upper mantle beneath these two mountain belts
possesses multiple layers of different anisotropy.

\section {Seismic Anisotropy in the Continental Lithosphere}


{\it Silver} [1996] summarizes arguments, developed over the previous
decade, that anisotropy is primarily a feature of the uppermost few
hundred kilometers of the Earth. This view has been recently challenged
by observations of seismic anisotropy near the core-mantle boundary [{\it
Kendall and Silver}, 1996; {\it Garnero and Lay}, 1997] and within the
transition zone [{\it Vinnik and Montagner}; 1996].  Nevertheless, the
presence of seismic anisotropy within the lithosphere is
well-documented.

The variety of mechanisms that produce seismic anisotropy in the lithosphere
centers on a handful of scenarios.  In the upper crust the strongest
influence is believed to be that of aligned cracks and/or pore spaces
[{\em Babuska and Pros} 1984], for which slower velocities are found for
waves that propagate normal to the average crack plane.  The aspect
ratio of pore/cracks and type of fluid determine the extent
and proportion of anisotropy [{\em Hudson}, 1981; {\em Crampin},
1991].  Alternating thin isotropic layers of higher and lower velocity
can also produce an overall anisotropic effect [{\em Backus}, 1962;
{\em Helbig}, 1994], with the velocities slower normal to the bedding than
along it.  In the lower crust and the uppermost mantle, cracks are
assumed to close in response to increasing overburden pressure [{\em
Babuska and Pros} 1984; {\em Kern et al}, 1993], though field exposures
of (formerly) deep-crustal fluid-filled cracks can be found [{\em
Ague}, 1995].  In the absence of cracks and inclusions, the
lattice-preferred orientation (LPO) of mineral crystals is taken as the
main cause of seismic anisotropy.  Most minerals composing the bulk of
the crust are anisotropic to some degree [{\em Babuska and Cara}, 1991],
as are the olivine and orthopyroxene that predominate in the
upper mantle anisotropy.  Different deformation mechanisms
can lead to the alignment of either the slow or the fast
crystallographic direction in olivine grains [{\em Nicolas et al.},
1973; {\em Ribe}, 1992], but LPO caused by dislocation creep in the
shallow mantle is commonly believed to lead to preferred alignment of
the fast axis [{\em Zhang and Karato}, 1995].  

It is natural to expect that strain-induced seismic anisotropy would be
particularly prominent in plate boundary regions, where deformations
are concentrated. World-wide observations of shear-wave splitting
support this notion [{\it Silver}, 1996]. Present-day regions of active
compression commonly have fast axis of seismic anisotropy aligned
sub-parallel to the strike of the orogen. One explanation for such
orientation is the preferred alignment of "slow" axes of olivine along
the direction of compression [{\it Nicolas et al.,} 1973]. 
{\it Vauchez and Nicolas} [1991] propose an alternative mechanism:
preferred alignment of the olivine "fast" axes along the orogen as a
result of concurrent strike-slip deformation commonly observed during
mountain building.

Some stable continental interiors have anisotropic intensity equal, if
not superior, to actively deforming regions, perhaps because many
now-stable continental regions have experienced plate-boundary
deformation in the past, and have retained a fossil deformation.
Patterns of seismic anisotropy within stable continental masses may
therefore record the tectonic history of these regions.
Seismic stations examined in this work lie within two Paleozoic
mountain belts, the Appalachians and the Uralides.
Both have been loci of continent-building accretionary episodes in
early Paleozoic time.
Both regions are presently embedded within stable continental region,
the North American and Eurasian plates, respectively.

A number of shear-wave splitting studies have examined the Northeastern
US, using both permanent and temporary stations [e.g., {\it Silver and
Chan}, 1991; {\it Fouch and Fischer}, 1995;  {\it Levin et al.}, 1996;
{\it Barruol et al.}, 1997; {\it Fouch et al.}, 1999].
In most cases a single "average" set of parameters was reported for
selected sites (Figure 5).
The average of all values shown in the map may be treated as a
"regional average", and comes out as delay $\tau \sim 0.9$s and 
a fast direction $ \phi \sim 85\deg$ CW from N.
{\it Helffrich et al.} [1994] reported average splitting parameters of
$\phi = N68 \deg E$ and $\tau \sim 1$ s for station ARU, the station we
examine in the foredeep of the Urals.
It is interesting to note that, in both the northeastern Appalachians and
the central Uralides, the fast axes of seismic anisotropy are often not
parallel to the strike of the orogenic belt, and thus do not follow the
pattern reported in active-tectonic regions like the Pyrenees [{\it
Vauchez and Nicolas}, 1991] .


\section {Splitting Parameters and Their Variance}  

Single seismogram estimates are needed to investigate the variation of
splitting parameters with S-wave propagation direction.  We use the
cross-correlation method to find the parameters that best fit the model
that the S-wave is composed of two pulses of identical shape but
orthogonal polarization, one delayed with respect to the other.

Let us suppose that the S-wave polarization lies with the $UV$ plane of
a cartesian $UVW$ coordinate system.  For weakly anisotropic material,
such as the Earth's mantle, the propagation direction will then be
nearly parallel to $W$.  In general, one might need to establish the
relationship between this coordinate system and the usual
$north-east-vertical$ coordinate system used to collect seismic data.
However, the SKS and SKKS core phases typically used in studies of
upper-mantle anisotropy have steep incidence angles, so that the
$UV$ plane is nearly horizontal (typically within $5\deg$ at the surface).
The vertical component of the seismogram is often
contaminated with compressional wave reverberations, so we believe it is
best to use only the horizontal component data.  The resulting
measurement of splitting direction will be slightly biased by this
approach.  However the effect of this bias on earth models can be
avoided simply by comparing these data to synthetic data oriented in the
same fashion.

We seek to find a rotation $\phi$ in the $UV$ plane, and a delay $\tau$ that
maximizes the cross-correlation:
$$
\begin{array}{c}
	C(\phi,\tau ) = {\bf u}^T {\bf v} / ( N s_u s_v) \\
	u_i = u(\phi, t_i) = \cos(\phi) U(t_i) - \sin(\phi) V(t_i) \\
	v_i = v(\phi, t_i) = \sin(\phi) U(t_i - \tau) + \cos(\phi) V(t_i - \tau)
\end{array}
\eqno(1) 
$$
where $s_u$ and $s_v$ are the root mean squared amplitudes of horizontal component seismograms
{\bf u} and {\bf v},
respectively, and $t_i = i \Delta t$, where $\Delta t$ is the sampling interval of
the seismogram. 
After grouping the splitting parameters into a vector ${\bf m} = [\phi,
\tau]^T$, we denote the cross-correlation as $C({\bf m})$.

We estimate the best-fitting vector ${\bf m}^{est}$ using a coarse grid
search followed by refinement with an interpolation algorithm.
We use a grid spacing of $\Delta t$ in  $\tau$ and $5\deg$ in $\phi$.    
We fit (using least-squares) the cross-correlation at the maximal node
and its nearest neighbors with a bi-quadratic function:
$$
C( \phi, \tau ) =  A + B \phi + C \phi^2 + D \tau + E \tau^2 + F \tau \phi
\eqno(2)$$
%(Eqn 2)
where $A, B, C, ... F$ are constants.  
The maximum of $C( \phi, \tau )$ occurs at 
${{\partial C} \over {\partial \phi}} =
{{\partial C} \over {\partial \tau}} = 0$, or:
$$
{\bf m}^{est} = {{[2BE-DF, 2CD-BF]^T} \over {(F^2 - 4CE)}}
%(Eqn 3)
\eqno(3)$$
Numerical tests (not shown) indicate that this technique locates the
maximum cross-correlation within $\sim 1$\% of the grid spacing, at
least for seismograms with the spectral characteristics of typical
$SKS$ phases.

The uncertainty of the splitting parameters can be
calculated by comparing this problem to the linearized inverse problem
${\bf G\cdot m = d}$. Here the goal is to estimate model parameters {\bf m}
and their covariance ${\bf C}_m$ from a data vector {\bf d}.  
A simple least-squares estimate of the model parameters, ${\bf m}^{est}$,
minimizes the misfit function
$$
E({\bf m}) =  { {1 \over N} {({\bf d} - {\bf G}\cdot  {\bf m})^T\cdot  ({\bf d} - {\bf G\cdot  m})} \over {{s_d}^2} }
\eqno(4)
$$
Here ${s_d}^2 = ({\bf d}^T\cdot  {\bf d} )/N$ is the mean-squared amplitude of the data.
If the data have uncorrelated error with variance ${\sigma_d}^2$, the
variance of the estimated model parameters is related to the curvature of the
misfit function:
$$
{\bf C}_m = \left. {{{\sigma_d^2} \over {N {s_d}^2 }} {{\left [{1 \over 2} \nabla_m \nabla_m E \right ]}}^{-1}}\right |_{{\bf m}={\bf m}^{est}}
\eqno(5)
$$
(see {\it Menke} [1989], eqn 3.52).  
This equation quantifies the notion that narrow minima are associated
with precise estimates.
The ratio $s_d / \sigma_d $ can be interpreted as the signal-to-noise
ratio $r_d$.

We apply this formula to maximize the coherence between two
seismograms {\bf u} and {\bf v}, each of length $N$.  
The misfit function is
$$
E({\bf m}) = {{1 \over 2N} \left( {{\bf u} \over s_u} - {{\bf v} \over s_v}\right)^T\cdot\left({{\bf u} \over s_u} - {{\bf v} \over s_v}  \right)}
\eqno(6)
$$
where $s_u$ and $s_v$ are the root mean squared amplitudes of {\bf u} and {\bf v}, respectively.
This definition of misfit is algebraically equivalent to $E=1-C$, where
$C = {\bf u}^T {\bf v} / (N s_u s_v)$ is the cross-correlation between the
two data series. 
The variance of {\bf m} is  given by:
$$
{\bf C}_m = 
\left. {-  {{1 \over { N r^2} } {\left [ {1 \over 2} \nabla_m \nabla_m C \right ]}}^{-1}}\right |_{{\bf m} = {\bf m}^{est} } 
\eqno(7)
$$
where $r$ is the signal-to-noise ratio (assumed the same on the two seismograms).
We estimate the second derivative matrix by differentiating the
quadratic interpolant (2) of the cross-correlation $C(\phi,\tau )$.

The signal-to-noise ratio can be estimated from the value of $C$
itself, by assuming that deviations from perfect correlation are caused
entirely by stochastic noise in the seismograms.
Denoting this noise as ${\bf n}_u$ and ${\bf n}_v$, and ${\bf u}_o$ and ${\bf v}_o$ as noise-free seismic signals, we have:
$$
C = {{{({\bf u}_o + {\bf n}_u)}^T ({\bf v}_o + {\bf n}_v)}
 \over
{\vert {\bf u}_o + {\bf n}_u \vert}^{1 \over 2} {\vert {\bf v}_o + {\bf n}_v \vert}^{1 \over 2}}
\eqno(8)
$$
We assume that the noise-free ${\bf u}_o$ and ${\bf v}_o$ are scaled
versions of each other.
The expected value of the cross-correlation is 
$$
< C > = 
{
{< {\bf u}_o^T {\bf v}_o >} 
\over 
{ \sqrt{<{\bf u}_o^T {\bf u}_o>} \sqrt{< {\bf v}_o^T {\bf v}_o>}[1 + r^{-2}]} 
} 
= {1 \over {[1 + r^{-2}]}}
\eqno(9)
$$
The estimated signal-to-noise ratio is thus $r=(C^{-1} -1)^{-{1 \over 2}} $.

Most time series are oversampled, and thus have correlated noise.  
For this case N in (3) must be replaced with the degrees of freedom
$\nu <N$:
$$
{\bf C}_m = \left. 
- \left({1 \over { \nu r^2}}\right) 
{ {\left [{1 \over 2} \nabla_m \nabla_m C \right ]}^{-1}}\right |_{{\bf m}={\bf m}^{est}}
\eqno(10)
$$
The degrees of freedom can be estimated by computing the autocorrelation
of the normalized misfit
$$
{\bf e}({\bf m}) = 
   {1 \over {\sqrt 2} }
     { 
     \left (
      { {\bf u}({\bf m}^{est}) \over s_u } - 
      { {\bf v}({\bf m}^{est}) \over s_v }
      \right )
     }
\eqno(11)
$$
The ratio $ N / \nu $ is approximately the width of the main peak in the autocorrelation of ${\bf e}({\bf m}^{est})$.

We have tested this method of computing covariance against Monte Carlo
simulations using synthetic SKS phases that have prescribed
signal-to-noise ratios in the 1:1 to 100:1 range.  The results (not
shown) indicate that the above methodology yields accurate
estimates of variance.  Standard deviations generally agree with 
Monte Carlo results to within 20-30\%. 
Standard deviations are typically $3-7\deg$ for the fast axis azimuth,
and $\sigma_{\tau}\sim 0.1-0.2$ s for data used in this study (Figure 6).


\section { Shear-Wave Splitting in NE Appalachians }

Figure 7 shows the pattern of SKS shear-wave splitting for two
earthquakes with different back azimuths (westerly and northwesterly),
observed at all available seismic stations in the NE Appalachian region.  
The observed shear-wave fast directions
differ for the two back azimuths.  The event with 
westerly back azimuth has a northeasterly fast direction, while the
event with the northwesterly back azimuth has an easterly fast direction.
The fast direction appears to be a rapidly-varying function of
back azimuth. Variability is also seen across the region for each
event, indicating that some lateral heterogeneity is  present. This heterogeneity was characterized by {\it Levin et al.} [1996] in terms of differences between "Appalachian" and "Grenvillian" provinces.

We  compiled SKS splitting data for the two longest-running of
these stations, HRV (Harvard, MA) and PAL (Palisades, NY) (Figure 8).
We used observations of $SKS$, $SKKS$ and $PKS$ phases, as well as a
few $S_{diff}$ and $S$ phases from deep-focus events. Core phases
($SKS$ and the like) are $SV$-polarized by the P-S conversion
at the core-mantle boundary, and are useful for the study of the
seismic anisotropy in the upper mantle and lithosphere. In the interest
of broadening the azimuthal coverage of our dataset we also used $S$
and $S_{diff}$ phases from medium-sized events with hypocenters deeper
then 500 km. We assume that these phases encounter anisotropy only in
the "receiver-side" upper mantle and the lithosphere. We also note that
splitting parameters obtained for $S_{diff}$ phases (two observations for HRV,
three observations for PAL) closely match splitting parameters obtained for $SKS$ phases from same events. Thus potential contamination of the $S_{diff}$
signal by the $D''$ anisotropy [{\it Garnero and Lay}, 1997] does not seem strong along these particular paths.

 The data for these two stations are quite similar (Figure 9).  Fast
direction azimuths tend to fall into the northeasterly and easterly
populations discussed above, resulting in a bimodal distribution of the
azimuthal angle between pairs of measurements (Figure 10).  The
assumption of two populations (with mean azimuths of $N60\deg E \pm
4\deg$ and $N119\deg E \pm 2\deg$) is statistically superior to the
assumption of only one population (with mean azimuth of $N95\deg E \pm
5\deg$) at the 99.9\% significance level (computed via the $F$-test).
The means of these two populations are also different at the 99.9\%
significance level (computed via the $t$-test).  We have examined the
SKS seismograms that were used as input to the splitting parameter
estimation procedure (Figure 11).  No anomalies that might cause
spurious parameter estimates are apparent, giving us confidence that
observed variation of splitting parameters is real.

{\it Barruol et al.} [1997] measured shear-wave splitting at HRV using
several of the earthquakes in our data set.
They report $\phi = N86\deg E (or N89\deg E), \tau = 0.99s (or 0.65s)$
(values given for two different processing techniques).
This result is quite similar to our "single population" mean of $N95
\deg E$.
It is interesting to note that the plot of all HRV data in {\it
Barruol et al.} [1997] (their Figure 5) and the table listing
individual values (electronic supplement table 2) contain members of
both populations of fast directions.
Another analysis of HRV data was done by {\it Fouch and Fischer} [1995]
to compare with data from the MOMA portable array, but included only
those events that occurred while the array was active (1995
through early 1996). 
They reported an average $\phi = N 118\deg E , \tau
= 1.13$s,  similar to our mean over easterly back-azimuths
($N119\deg E$).  In both cases the measured average value is
biased by the event distribution, which is dominated by northwesterly
events from  NW Pacific earthquakes. The discrepancy in reported
sets of splitting parameters most likely reflects differences in
the distribution of data with back azimuth.  
Although of dubious value given
the systematic fluctuations of the data, the mean value of the
splitting direction (calculated as a model result, and discussed below) is
$N89\deg E$, which matches the $\phi$ value reported by {\it
Barruol et al.} [1997].

We  modeled the seismic data using splitting parameters derived from
synthetic seismograms of SKS phases.  We  compared two different
algorithms for computing the synthetic seismograms in vertically
stratified, anisotropic media: a propagator matrix method [{\it Levin
and Park}, 1997b], and a ray method.  Both give near-identical results.  
We use a grid-search over anisotropic parameters to find a best-fitting
earth model, where the goodness-of-fit criteria minimizes the misfit
$$ E = \sum_i {
{({\tau_i}^{observed} - {\tau_i}^{predicted})^2 \over
{\sigma^2}_{\tau}} + {({\phi_i}^{observed} - {\phi_i}^{predicted})^2
\over {\sigma^2}_{\phi}} } $$ 
Here $\sigma_\tau$ and $\sigma_\phi$ are
the standard deviations of the delay time $\tau$ and fast azimuth
$\phi$, respectively.  We have examined two classes of earth models:
one or two anisotropic mantle layers placed between an isotropic
crustal layer and an isotropic mantle halfspace.  We search only for
the thickness of the layers and the orientation of the anisotropic
tensors. The anisotropic medium is constrained to consist of 30\%
orthorhombic olivine and 70\% isotropic olivine, a mixture that is
about 6\% anisotropic for shear waves.  The best-fitting one-layer model has an
anisotropic layer that is 58 km thick, and the two layer model has top
and bottom layers that are 60 and 90 km thick, respectively. 
Parameters for our preferred  hexagonally symmetric model are indicated
in Table 1.
Tensor orientations for our preferred orthorhombic model are indicated
in Table 2 and Figure 12.
Both hexagonal and orthorhombic two-layer models correctly
capture the variation of splitting parameters with back azimuth, while
the one-layer models do not. 
Figure 13 compares results for the 1- and 2-layer orthorhombic models.  
The variance reduction of the two-layer model is roughly three times
greater than the one-layer model, an amount that is statistically
significant to the 99\% level (computed via the $F$-test).  
The two-layer orthorhombic model gives fast-axis azimuths of $N53\deg
E$ and $N115\deg E$ for the bottom and top layers, respectively, which
are close to the means of the two observed azimuthal populations.
The fast-axis strikes for hexagonal symmetry are only slightly
different, at $N50\deg E$ and $N100\deg E$ for the bottom and top
layers, respectively.
However, the symmetry axes are tilted: only $15\deg$ above the
horizontal in the top mantle layer, but $40\deg$ below the
horizontal in the lower layer.

To test whether the symmetry-axis tilts are significant
we compare the observed back-azimuth pattern of the apparent fast
direction with those predicted by our hexagonal and orthorhombic
models, and also by the 2-layer splitting operator of {\it Silver and
Savage} [1994] (Figure 14). 
To construct the operators, we computed the delays $\tau_k$ a
vertically-incident shear phase would experience in each layer of the
orthorhombic and hexagonal models, and used respective fast
directions.
Predicted patterns are quite similar although, as one
would expect, the approximation and synthetics show greatest difference
at discontinuities in the pattern, where waveform complexity is the
greatest. Also, patterns from forward modeling are only approximately
periodic because of the tilts of the anisotropy axes.  These violations
of ${\pi} / 2$ periodicity may possibly serve as diagnostic traits in
choosing the preferred model. 
Our present collection of data is too limited to uniquely resolve the
tilt of the symmetry axes, though they seem to prefer some deviation
from the horizontal.

The broad spatial coherence of shear-wave directions throughout the NE
Appalachians (as evidenced in Figure 5) suggest that there is a strong
vertical stratification of the anisotropy in the upper mantle beneath
that region.  Earth models with two anisotropic upper-mantle layers can
fit the observed SKS splitting data well. More complicated
vertically-stratified models cannot be ruled out, but are not required
by the data.  Some lateral heterogeneity is also present in the
splitting data.  We note, however, that even large heterogeneity in the
splitting data does not necessarily translate into large heterogeneity
in structure.  Given the quick variation of the parameters with
back azimuth, just a few degrees rotation of either the earth model or
the incoming shear wave can lead to widely different values of the
splitting parameters.  A complete characterization of any 3-D "anisotropic
domains" responsible for the observed lateral heterogeneity will
require extensive back-azimuth coverage at many stations.

We interpret the top anisotropic layer to represent the continental
lithosphere associated with the Appalachian orogen, and the bottom
layer to represent the asthenosphere.
A conceptual model of layered anisotropy under HRV is
presented on Figure 15.  Alignment of the fast axes of olivine in the
upper layer is near-normal to the strike of geomorphological features
in New England, where the trend of the Appalachians rotates from northeast to north. The alignment
in the lower layer of the model is more in line with the overall strike
of the Appalachian Orogen, as well as the hypothetical "edge" of
the North American cratonic keel.

\section{Shear Wave Splitting at Station ARU in the Foredeep of Urals}

The GSN station ARU (Arti, Russia) is located on the easternmost edge
of the East European platform (Figure 16), in the Uralian foredeep
[{\em Ivanov et al.} 1975; {\em Zonenshain et al.} 1984] that accumulated
sediments in a passive continental margin setting through most of the
Paleozoic. The continental collision at the final stage of Uralian
orogeny formed an extensive thrust-sheet complex to the east of ARU,
and may also have deformed the sediments underlying ARU.
[{\it Levin and Park}, 1997a] studied crustal anisotropy beneath ARU
using the anisotropic receiver function method.
This method, which utilizes $P-S$ conversions in layered media, is
sensitive to anisotropy adjacent to interfaces where the conversions
occur.
{\it Levin and Park}, [1997a] show that receiver function data are consistent
with a vertically-stratified earth structure in which both the
uppermost and the lower crust are anisotropic, with hexagonal symmetry
and strong tilt in the symmetry axes.  Unlike  typical models of
seismic anisotropy in the mantle, the crustal model for  ARU
contains a layer of anisotropy with a slow symmetry axis. A conceptual
model that would exhibit this type of anisotropy is a layer of
imbricated peridotite and metapelite lenses in the region of the
crust-mantle transition, as described by {\it Quick et al.} [1995] in the
Ivrea deep crustal exposure in the Alps.
Upper mantle anisotropy under ARU has previously been studied by {\it
Helffrich et al}, [1994], who report a mean fast axis strike of $N 68\deg
E$ and a splitting time $\tau\sim 1$ s.

We have compiled splitting data from SKS and other core-refracted
phases for this station (Figure 17).  Several populations of
fast-axis strikes are evident. For example, events from the
northwest have westerly strike, while those from the east have more
northwesterly strike.  The best-fitting one-layer anisotropic earth
model (Figure 18, top), with fast-axis at $N73\deg E$
azimuth, does not reproduce this pattern well.  A combination of the {\it
Levin and Park} [1997a] crustal model with the mantle model of  {\it
Helffrich et al.} [1994], does better (not shown), but does poorly at
the ENE back azimuths. 
A model with a second mantle layer improves the fit significantly
(Figure 18, bottom).
This best-fitting model (Table 3) has a 58 km  upper mantle
layer with a fast-axis striking $N50\deg E$ atop a 140 km
layer with fast-axis plunging 40 $\deg $ to the east (Figure 19).

Figure 20 illustrates a conceptual model of anisotropic layering under
ARU.  Lower crust and uppermost mantle under the Uralian foredeep are
characterized by the common direction of anisotropy $\sim 50\deg$,
which is significantly oblique to the trend of Urals.  It may be
related to the deformation within the East European platform that predates the formation of the orogen.  The
anisotropy-inducing fabric in the lower part of the lithosphere is
aligned nearly normal to the strike of the Uralian Orogen.

%\newpage
\section {Discussion of Geodynamic Implications}

In both Appalachians and Uralides we find evidence for at least two distinct
layers of seismic anisotropy in the mantle. The exact depth and
thickness of layers is subject to assumptions about the percent of
aligned minerals within the volume. 
Anisotropic intensities in our models follow measurements done on
hand-samples of peridotite from ophiolites [{\it Christensen}, 1984],
and are likely to represent the upper bound in percent of alignment.
If the alignment of minerals is weaker, the required thicknesses of the
anisotropic layers will increase.
The tilt of inferred anisotropy-inducing fabric with respect to the
horizontal is a function of the symmetry system chosen for the
anisotropy. 
For the Appalachian stations, both  orthorhombic anisotropy, with
nearly horizontal fast axes, and hexagonal anisotropy, with tilted fast
symmetry axes, satisfy the data.
We suspect that this tradeoff will be common in studies of this type,
and difficult to resolve without a clear indication of the best
symmetry choice from mineral texture studies.
Significantly, in both modeling exercises the horizontal azimuth of the
fast axis does not depend on the symmetry system and the method of
computing synthetic seismograms.
We thus believe that robust elements of our models are the minimum
number of anisotropic layers, their vertical sequence, the orientation
of fast direction and the cumulative anisotropy within each of them.

The overall regional consistency of observations in the Northeastern
Appalachian region (Figure 7) argues against a "local" character for
shallow-mantle structures revealed by shear-wave splitting. Rather, the structure
modeled using the HRV dataset appears to be common to the area of
Appalachian terranes. On the other hand, some regional variation is
present (e.g., {\it Levin et al}, [1996]). Whether it is caused by the
"true" lateral variation in anisotropic features of the subsurface, or
simply reflects changes in geometry of observation, is an open
question.  

In case of the Uralian foredeep we have a "spot" measurement, and
arguing for its regional extent is harder.  Major geologic structures
of the Ural Orogen are very consistent along strike, making it almost
2-dimensional.
An average splitting direction of $0 \deg$ (sub-parallel to the strike of the Urals) reported by [{\it Vinnik et al.}, 1992] at the station SVE near Yekaterinburg (Figure 16) falls on the opposite side of the main
Uralian Fault Zone. This major suture divides accreted terranes from the East European Platform [{\it Zonenshain et al.}, 1984]. The splitting at SVE most likely reflects a different structure in the crust and the uppermost mantle. 

Given the uncertainties discussed above, the interpretation
is necessarily tentative.
If we consider the anisotropy to be "frozen in" or "fossil," confined
within the continental keel, in both Appalachian and Urals we may infer
at least two distinct past tectonic episodes.
{\it Abbott} [1992] describes a conceptual model for accumulating
material with different deformation fabrics within the body of a
continent, via underplating the continent with oceanic-plate material
during successive episodes of subduction.
A similar model, only with oceanic "slabs" stacked in the lateral
direction, has been advocated for Western Europe [e.g. {\it Babuska et
al.}, 1993, {\it Plomerova et al}, 1996].
The "frozen fabric" explanation appears plausible for the upper
layers in both modeled regions. In Urals the lowermost crust has the same
orientation ($N50\deg E$, from {\it Levin and Park} [1997a]) as the
upper layer of our mantle model. 
Interestingly, the sense of anisotropy
reverses from the crust (slow axis) to the mantle (fast axis). We
believe that crustal anisotropy in the crust is  
imposed by fine layering of materials with contrasting properties,
while in the mantle it is imposed by
preferred alignment of peridotite minerals. The strike of fabric in
lower crust and uppermost mantle under ARU is oblique to the strike of
the orogen, and may reflect a tectonic episode predating its
formation.
Analysis of crustal $P-SH$ conversions [{\it Levin and Park}, 1997b]
under HRV did not reveal strong crustal anisotropy of the kind seen
under ARU. The upper layer of the mantle anisotropy has fabric
orientation that is roughly normal to the geomorphological features and
main tectonic boundaries in the region. Numerous subduction episodes,
both east-verging and west-verging, took place during the formation of
the Appalachians, and the fabric we reconstruct is likely to be a
remnant of one of them.  

An alternative mechanism for anisotropy in the mantle -- active flow in
the asthenosphere [e.g., {\it Vinnik et al}, 1992] -- appears more
suitable for the lower layers in our models.
For the Appalachians the strike
of fabric in the lower layer of our model aligns closely with the
absolute plate motion vector of $\sim N65\deg E$ [{\it Gripp and
Gordon}, 1990]. It also aligns with the hypothetical edge of the
North American continental keel.  
The keel edge would direct the orientation of asthenospheric flow if
one assumes that mantle moves "west" relative to the North American craton, and
around its keel [{\it Fouch et al.}, 1999].
The Urals are located in the middle of the Eurasian
continent which is near-stationary with respect to the hotspot
reference frame. 
The direction of possible motion for Eurasia is approximately
east-west, if one ignores the error bars on the rate of motion in {\it
Gripp and Gordon}, [1990].
This direction, though poorly constrained, aligns with the fast axis of
the lower layer in the anisotropic structure we infer under ARU.

Our results clearly contradict the notion that fast-axis strike for
shear-wave splitting in a region of compressional deformation should
align with the strike of the orogen [{\it Vauchez and Nicolas}, 1991;
{\it Nicolas}, 1993; {\it Vauchez and Barruol}, 1996]. 
In two regions of Paleozoic deformation we find 
patterns of seismic anisotropy consistent with vertically heterogeneous
anisotropic structure.
It is possible that in such complex regions a subset of measurements will
indeed line up with the strike of an orogen.
Our experiments with multilayered anisotropic structures show that
direct mapping of splitting parameters (or their averages) onto
tectonic features can be misleading.



%\newpage
\section{Conclusions}

Observations of shear-wave splitting in the Northeastern US
Appalachians and in the foredeep of the Urals vary significantly with the
back azimuth and incidence angle of the phase.
To analyze these datasets properly we developed a new technique for
estimating uncertainties of splitting parameters.
Using this technique we find that typical errors of the shear-wave
splitting parameters determined from low-passed broadband data from GSN
station HRV are $3\deg - 7\deg$ for the fast direction, and
$0.1-0.2s$ for the delay.

Experiments with synthetic seismograms generated in simple multilayered
anisotropic structures show that  splitting parameters tend
to vary significantly with the back-azimuth of the analyzed shear wave.  
A  restricted subset of back azimuths may strongly bias any
model derived from observations, especially if the observations are
averaged.
On the other hand, the azimuthal variation pattern provides important
constraints on vertical or lateral variation of anisotropic properties
in the Earth.

On the basis of data from well-recorded events with different
back azimuths, splitting parameters appear to
be broadly consistent throughout the Appalachian terranes in the
Northeastern US.  
(This consistency weakens for stations 
west of the Appalachians.) 
A close similarity in back-azimuth  dependence of splitting
parameters is found in data from two long-running stations in the
Northeast US - HRV and PAL.  Good back-azimuth coverage at
these two stations allows us to separate observations into two
statistically significant populations.  Within these populations mean
azimuths are $N60\deg E \pm 4\deg$ and $N119\deg E \pm 2\deg$, and
delay values vary within each population from near-zero to $\sim 1s$.
The exact values of delays, as well as individual estimates of fast
direction, are affected by the filter parameters chosen when low-passing.  
The back-azimuth dependence of splitting parameters for the station ARU
near the Urals is characterized by sharp transitions between different
groups of observations.

Using synthetic seismograms computed in flat-layered media we developed
one-dimensional models of seismic anisotropy distribution under
stations HRV and ARU.  The model for HRV contains two layers of
anisotropic material under an isotropic crust, with fast-axis azimuths
of $N53\deg E$ and $N115\deg E$ for the bottom and the top layers,
respectively. Depending on the choice of symmetry for the elastic
tensors, these axes are tilted (hexagonal symmetry)
or near-horizontal (orthorhombic symmetry). Assuming 30\% orthorhombic
olivine and 70\% isotropic olivine, a mixture that is about 6\%
anisotropic, the vertical dimensions are 60 and 90 km for the top and
bottom layers, respectively.
The model for ARU includes crustal structure that was constrained
using Ps converted phases [{\it Levin and Park}, 1997a].  
Assuming hexagonal symmetry of the upper mantle anisotropy, the model for  ARU
predicts a $\sim $ 60 km layer with a fast-axis at $N50\deg E$ atop a 140 km
layer with a fast-axis plunging 40$\deg$ towards $N90\deg E$.

The analysis performed in this paper was made possible by good
azimuthal coverage of observations. These are generally obtainable
through prolonged observation.  Data from short deployments, even in
stable continental regions, apparently run the risk of bias from
an uneven distribution of seismicity.

\begin{acknowledgements}
We used {\bf GMT} software \markcite {[{\it Wessel \& Smith,} 1991]}
to prepare figures. Comments by Matt Fouch, Ray Russo and 2 anonymous reviewers helped improve this paper. We thank Matt Fouch for providing a copy of {\it Fouch et al.}, [1999] prior to publication.
This work has been supported by NSF grant EAR-9707189.  
\end{acknowledgements}

\begin{references}

\newcommand{\GJI}{{ \em Geophys.\ J. Int.},\ }

\reference
Abbott, D., The case for accretion of the tectosphere by buoyant subduction, \grl {\em 18}, 585--588, 1992.

\reference
Ague, J. J.,  Deep crustal growth of quartz, kyanite and garnet
into large aperture, fluid-filled fractures, northeastern Connecticut,
USA, {\em J. Metamorphic Geology, 13}, 299--314, 1995.

\reference
Aki, K., and P. G. Richards,   {\em Quantitative Seismology},
W. H. Freeman, San Francisco, 1981.

\reference
Babuska, V., and Z. Pros, Velocity anisotropy in granodiorite
and quartzite due to the distribution of microcracks, \gjras {\em 76},
121--127, 1984.

\reference
Babuska, V., and M. Cara, {\em Seismic Anisotropy in the Earth},
Kluwer Academic, Dordrecht, 1991.

\reference
Babuska, V.,  J. Plomerova, and  J. Sileny, Models of seismic
anisotropy in the deep continental lithosphere, \pepi, {\em 78},
167--191, 1993.

\reference
Backus, G. E., Long-wave elastic anisotropy produced by
horizontal layering, \jgr {\em 67}, 4427--4440, 1962.

%\reference
%Barruol, G., and D. Mainprice, A quantitative evaluation of the
%contribution of crustal rocks to the shear-wave splitting of
%teleseismic SKS waves, \pepi, {\em 78}, 281--300, 1993.

\reference
Barruol, G., P. G. Silver and A. Vauchez, Seismic anisotropy in the
Eastern United States; deep structure of a complex continental plate,
\jgr {\em 102}, 8329--8348, 1997.

\reference
Bostock M. G. and J. Cassidy, Variations in shear wave splitting across the
canadian National Seismic Network, \grl {\it 22}, 5-8, 1995.


\reference
Bowman, J. R. and M. A. Ando, Shear-wave splitting in the upper mantle wedge above the Tonga subduction zone, \gjras, {\em 88}, 25--41, 1987.

\reference
Crampin, S., E. M. Chesnokov and R. G. Hipkin, Seismic anisotropy - the
state of the art: II, {\gjras}, {\em 76} , 1--16, 1984.

\reference
Crampin, S., Wave propagation through fluid-filled inclusions of
various shapes: interpretation of extensive-dilatancy anisotropy, \GJI
{\em 104}, 611--623, 1991.

\reference 
Christensen, N. I., The magnitude, symmetry and origin of upper mantle
anisotropy based on fabric analyses of ultramafic tectonites, \gjras
{\em 76}, 89--111, 1984.

\reference
Dziewonski, A. M., and D. Anderson, Preliminary reference earth model,
\pepi, {\it 25}, 297--356, 1981.

\reference
Fouch, M. J., and  K. M. Fischer, A search for split shear waves
beneath the Eastern United States, {\it EOS}, {\em 76}, 412, 1995.

\reference
M.J. Fouch, K.M. Fischer, and E.M. Parmentier, Shear-wave splitting,
continental keels, and patterns of mantle flow, submitted to \jgr, 1999.

\reference
Garnero, E. J., and T. Lay, Lateral variations in lowermost mantle shear
wave anisotropy beneath the North Pacific and Alaska, \jgr {\em 102},
8121--8135, 1997.

\reference
Granet, M., A. Glahn and U. Achauer, Anisotropic measurements in the Rhinegraben
area of Massif Central: Geodynamic implications, {\it  PAGEOPH, 151}, 333-364, 1998.

\reference
Gripp, A. E., and R. G. Gordon, Current plate velocities relative to
the hotspots incorporating the NUVEL-1 global plate motion model, \grl
{\em 17}, 1109--1112, 1990.

\reference
Helbig, K., {\em Foundations of  Anisotropy for Exploration
Seismics}, Pergammon/Elsevier Science, Oxford, 1994.

\reference
Helffrich, G., P. G. Silver, and  H. Given, Shear-wave splitting variation
over short spatial scales on continents, \GJI, {\em 119}, 561--573, 1994.

\reference
Hirn, A., J. Diaz, M. Sapin and J.-L. Veinante, Variation of shear-wave residuals and splitting parameters from array observations in southern Tibet,
{\it PAGEOPH, 151}, 395-405, 1998.

\reference
Hudson, J. A., Wave speeds and attenuation of elastic waves in
material containing cracks, \gjras  {\em 64}, 133--150, 1981.

\reference
Iidaka, T. and F. Niu, Evidence for an anisotropic lower mantle beneath eastern Asia; comparison of shear-wave splitting data of SKS and P660s, \grl {\it 25}, 675-678. 1998.

\reference 
Ivanov, S. N., A. S. Perfiliev, A. A. Efimov, G. A. Smirnov, V. M.
Necheukhin, and G. B. Fershtater, Fundamental features in the
structure and evolution of the Urals, {\em Amer. J. Sci., 275A},
107--130, 1975.

\reference
Kendall J. M. and P. G. Silver,  Constraints from seismic anisotropy on
the nature of the lowermost mantle, {\it Nature}, {\em 381}, 409--412,
1996.

\reference
Kern, H., C. H. Walther, E. R. Fluh, and  M. Marker, Seismic
properties of rocks exposed in the POLAR profile region - constraints
on the interpretation of the refraction data. {\em Precambrian Research
64}, 169--187, 1993.

\reference
Kosarev, G. L., L. I. Makeyeva, and L. P. Vinnik, Anisotropy of the
mantle inferred from observations of $P$ to $S$ converted waves, \gjras
{\em 76}, 209--220, 1984.

\reference
Levin, V.,  W. Menke, and  A. L. Lerner-Lam, Seismic anisotropy
in northeastern US as a source of significant teleseismic $P$ traveltime
anomalies, \GJI, {\em 126}, 593--603, 1996.

\reference
Levin, V., and J. Park, Crustal anisotropy in the Ural Mountains foredeep
from teleseismic receiver functions, {\em Geophys. Res.
Letts., 24}, 1283--1286, 1997a.

\reference
Levin, V., and J. Park, P-SH conversions in a flat-layered medium with
anisotropy of arbitrary orientation,  \GJI, {\em 131}, 253--266, 1997b.

\reference
Menke, W., Geophysical Data Analysis: Discrete Inverse Theory, Revised
Edition, {\it Academic Press}, New York, 1989.

\reference
Nicolas, A., Why fast polarization directions of SKS seismic waves are parallel to mountain belts, \pepi, {\em 78}, 337--342, 1993.

\reference
Nicolas, A.,  F. Boudier, and A. M. Boulier, Mechanism of flow in
naturally and experimentally deformed peridotites, {\em Am. J. Sci.},
{\em 10}, 853--876, 1973.

\reference

Ozalaybey, S., and M. K. Savage, Double-layer anisotropy resolved from S phases,
\GJI {\it 117}, 653-664, 1994.

\reference 
Park, J., 1996. Surface waves in layered anisotropic structures, \GJI
{\em 126}, 173--183.

\reference
Plomerova, J., J. Sileny and V. Babuska, Joint interpretation of
upper-mantle anisotropy based on teleseismic P-travel time delays and
inversion of shear-wave splitting parameters, \pepi, {\it 95},
293--309, 1996.

\reference
Quick, J. E., S. Sinigoi and  A. Mayer, Emplacement of mantle
peridotite in the lower continental crust, Ivrea-Verbano zone,
northwest Italy, {\it Geology, 23}, 739--742, 1995.

\reference
Ribe, N., On the relation between seismic anisotropy and finite
strain, \jgr {\em 97}, 8737--8747, 1992.

\reference
Rumpker, G. and P. G. Silver, Apparent shear-wave splitting parameters in the presence of vertically varying anisotropy, \GJI, 790--800, 1998.

\reference
Russo, R., and P. G. Silver, Trench-parallel flow beneath the Nasca Plate
from seismic anisotropy, {\it Science, 263}, 1105--1111, 1994.

\reference
Silver, P. G., and M. K. Savage, The interpretation of shear-wave 
splitting parameters in the presence of two anisotropic layers, 
\GJI {\it 119}, 949--963, 1994.

\reference
Silver, P. G., and W. W. Chan, Shear-wave splitting and
subcontinental mantle deformation, \jgr {\em 96}, 16429--16454, 1991.

\reference
Silver, P. G.,  Seismic anisotropy beneath the continents:
Probing the depths of geology, {\em Ann. Rev. Earth Planet. Sci., 24},
385--432, 1996.

\reference
Schutt, D., E. D. Humphreys and K. Dueker, Anisotropy of the Yellowstone Hot Spot wake, eastern Snake River Plain, Idaho, {\it PAGEOPH, 151} 443-462, 1998.

\reference
Vauchez, A. and A. Nicolas, Mountain building; strike-parallel motion
and mantle anisotropy, {\it Tectonophysics}, {\em 185}, 183--201,
1991.

\reference
Vauchez A. and G. Barruol, Shear-wave splitting in the Appalachians and
Pyrenees: importance of the inherited tectonic fabric of the
lithosphere, \pepi {\em 95}, 127--138, 1996.

\reference
Vinnik, L. P.,  L.I. Makeyeva, A. Milev, and A. Yu. Usenko,
Global patterns of azimuthal anisotropy and deformations in the
continental mantle, \GJI {\it 111}, 433--447, 1992.

\reference
Vinnik L. P.,  V. G. Krishna, R. Kind, P. Bormann and K. Stammler, Shear-wave splitting in the records of the German Regional Seismic Network, \grl {\it 21}, 457-460. 1994.

\reference
Vinnik, L. P , R. W. E. Green and L. O. Nicolaysen, Recent deformations
of the deep continental root beneath Southern Africa. {\it Nature},
{\em 375}, 50-52, 1995.

\reference
Vinnik, L. P. and  J.-P. Montagner, Shear-wave splitting in the
mantle Ps phases, \grl {\em 23}, 2449--2452, 1996.

\reference
Wessel, P. \& W. H. F. Smith, 1991. Free software helps map and display data,
\eos {\it 72}, 441.

\reference
Wolfe, C. J. and P. G. Silver,  Seismic anisotropy of oceanic upper mantle; shear-wave splitting methodologies and observations, \jgr {\it 103}, 749 - 771, 1998. 

\reference
Wolfe, C. J. and S. C. Solomon, Shear-wave splitting and implications for 
mantle flow beneath the MELT region of the East Pacific Rise, {\it Science, 280}, 1230-1232, 1998.

\reference
Zhang, S. and S. Karato, Lattice preferred orientation of olivine
aggregates deformed in simple shear, {\em Nature, 375}, 774--777, 1995.

\reference
Zonenshain, L. P., V. G. Korinevsky, V. G. Kazmin, D. M. Pechersky, V. V.
Khain, and  V. V. Matveenkov, Plate tectonic model of the
South Urals development, {\it Tectonophysics, 109}, 95--135, 1984.

\end{references}


\newpage
\onecolumn

\begin {table}
\caption{
Anisotropic structure (hexagonal symmetry) consistent with shear-wave
splitting observations at HRV.
``Depth'' indicates the bottom of each layer.
The parameters B and E scale peak-to-peak variations of compressional
and shear velocity, respectively, each with $\cos 2\eta$ azimuthal
dependence \markcite{[{\it Park}, 1996]}.
The angles $\theta$ and $\phi$ define the tilt (from vertical) and
strike (CW from N) of the axis of symmetry within each anisotropic
layer.
}
\begin{tabular}{cccccccc}
\hline 
depth (km) & Vp, km/s & Vs, km/s & $\rho$ ($\mbox{g/cm}^3$) & B (\% of Vp) & E (\% of Vs) &
$\theta$ (\deg)  & $\phi$ (\deg) 
\\
\hline
35 & 6.00 & 3.60 & 2.5 & 0 & 0 & - & - \\

135 & 8.00 & 4.60 & 3.1 & 3 & 3 & 75 & 100 \\

235 & 8.00 & 4.60 & 3.1 & 4 & 4 & 130 & 50 \\
 
$\infty$ & 8.20 & 4.7 & 3.3 & 0 & 0 & - & - \\
 
%\hline  
\end{tabular}
\end{table}

\vspace{20pt}

\newpage

\begin {table}
\caption{
Anisotropic structure (orthorhombic symmetry) consistent with shear
wave splitting observations at HRV. Velocity values in anisotropic
layers are the isotropic averages of respective elastic tensors. 
The anisotropic medium is modeled as 30\% orthorhombic
olivine and 70\% isotropic olivine, a mixture that is $\sim 6$\%
anisotropic.
The angles $\theta$ and $\phi$ define the tilt (from vertical) and
azimuth (CW from N) of the symmetry axes (fast, intermediate, slow)
within each anisotropic layer.}

\begin{tabular}
{cccccccccc}
\hline 
depth (km) & Vp, km/s & Vs, km/s & $\rho$ ($\mbox{g/cm}^3$) &  $\theta_f$ (\deg)  & $\phi_f$  (\deg) & $\theta_i$ (\deg)  & $\phi_i$ (\deg) & $\theta_s$ (\deg)  & $\phi_s$ (\deg) \\
\hline
40 & 6.8 & 3.9 & 2.85 & - & - & - & - & - & \\
100 & 8.3 & 4.8 & 3.3 & 89 & 115 & 80 & 25 & 10 & 210 \\
190 & 8.3 & 4.8 & 3.3 & 90 & 53 & 67 & 323 & 23 & 143 \\
$\infty$ & 8.3 & 4.8 & 3.3 & - & - & - & - & - & \\

%\hline  
\end{tabular}
\end{table}

\vspace{20pt}
\newpage

\begin {table}
\caption{
Anisotropic structure (hexagonal symmetry) consistent with $P$-coda
receiver functions [{\it Levin and Park}, 1997a] and shear-wave
splitting observations at ARU.
Depth indicates the bottom of each layer.
The parameters B and E scale peak-to-peak variations of compressional
and shear velocity, respectively, each with $\cos 2\eta$ azimuthal
dependence \markcite{[{\it Park}, 1996]}.
The angles $\theta$ and $\phi$ define the tilt (from vertical) and
strike (CW from N) of the axis of symmetry within each anisotropic
layer.
}
\begin{tabular}{cccccccc}
\hline 
depth (km) & Vp, km/s & Vs, km/s & $\rho$ ($\mbox{g/cm}^3$) & B (\% of Vp) & E (\% of Vs) &
$\theta$ (\deg)  & $\phi$ (\deg) 
\\
\hline
1 & 5.00 & 2.95 & 2.1 & -15 &  0 & 45 & 345 \\
 
20 & 6.40 & 3.70 & 2.3 & 0 & 0 & - & - \\

23 & 5.8 & 3.30 & 2.3 & 0 & 0 & - & - \\

33 & 6.60 & 3.80 & 2.6 & 0 & 0 & - & - \\
 
42 & 7.60 & 4.40 & 3.0 & -15 & -15 & 65  & 230  \\
 
100 & 8.00 & 4.60 & 3.2 & 6 & 6 & 80  & 50  \\

240 & 8.00 & 4.60 & 3.2 & 4 & 4 & 130  & 90  \\

$\infty$ & 8.20 & 3.3 & 4.7 & 0 & 0 & - & - \\
 
%\hline  
\end{tabular}
\end{table}


\newpage

Figure 1. Shear-wave fast direction and delay for SKS waves received at
a hypothetical station from a variety of directions and apparent
velocities (in km/s). Data are shown as a bar centered on the nominal
back azimuth and apparent velocity, with the bar's orientation parallel
to the azimuth of the the fast direction, and its length proportional
to the delay. Near-zero delays are plotted with open circles.  (Left)
For an earth model with a single hexagonally-anisotropic layer
with horizontal symmetry axis overlying an isotropic halfspace; (Right)
for a model in which the symmetry axis plunges 45 degrees.  Note that the
splitting parameters vary most slowly in the horizontal case.

Figure 2. Synthetic SKS phases for two earth models. (Right) Two 50 km
hexagonally-anisotropic layers overlying an isotropic halfspace.
The symmetry axis is horizontal in both layers, and has an azimuth of
N30E in the top layer and N60E in the bottom layer. (Left) One 100 km
 anisotropic layer overlying an isotropic halfspace. The
anisotropic tensor is the arithmetic mean of the tensors in the
two-layer case.  Radial (top row) and transverse (second row)
horizontal component synthetic seismograms (bold) are computed by
convolving the impulse response functions (solid) with a pulse that has
a bandwidth similar to typical SKS phases.  This SKS pulse ($N270\deg
E$, 20 km/s) (third row) is reconstructed using apparent
splitting parameters computed by the cross-correlation method.  Note
that the two layer case has the larger error (bottom row).

Figure 3. Apparent splitting parameters for SKS waves received
at a hypothetical station. (Right) For an earth model with two 50 km
 hexagonally-anisotropic layers overlying an isotropic halfspace.
The symmetry axis is horizontal in both layers, an has an azimuth of
N30E in the top layer and N60E in the bottom layer. (Left) For an earth
model with one 100 km  anisotropic layer overlying an isotropic
halfspace. The anisotropic tensor is the arithmetic mean of the tensors
in the two-layer case.  Splitting parameters are computed by the
cross-correlation method from synthetic SKS phases.  Note that the two
layer case has the more complex pattern.

Figure 4. Influence of the
filtering applied to the synthetic waveforms on the apparent splitting parameters. Waveforms are simulated in a 2-layer model developed for station HRV (see Table 2). Apparent splitting parameters measured
from waveforms with an upper spectral limits of 0.45Hz (bold), 0.35Hz (solid) and 0.15Hz (dotted)
are shown, with a lower limit of 0.05Hz common for all. For all cases tested
the estimate of the apparent fast direction tends to vary greatly when the delay is the smallest. The lowest bandpass (0.15Hz) gives most unstable results throughout.

Figure 5. Map of NE Appalachian region, with "average" parameters of
seismic anisotropy plotted at points where they were constrained by
various workers. Arrow azimuths correspond to the fast directions
determined for the particular site, and are scaled with estimated
delay.  The compilation is extracted from the "Anisotropy
Resource Page" maintained by Derek Schutt ($
http://darkwing.uoregon.edu/\sim schuttd/aniso\_source.html$). Individual
values are from {\it Silver and Chan}, [1991]; {\it Bostock and Cassidy}, [1995]; {\it Barruol et al.}, [1997] and
{\it Fouch et al.}, [1999].

Figure 6. Histograms of single seismogram splitting parameter errors for HRV data.  {\it (top) }
One $\sigma$ error in delay time. {\it (bottom) } One $\sigma$ error in
fast direction azimuth.

Figure 7. Map of NE Appalachian region showing shear-wave splitting
data for two earthquakes with different propagation directions (large one-sided arrows) observed in 1995.  Splitting azimuth and delay are shown at
each station as two-sided arrows aligned with the fast direction and scaled with
delay. Symbols are color-coded by event, black for one and white for the other. Note that splitting directions for the two events is quite different, yet
is fairly consistent across the region for each event.


Figure 8. Shear-wave splitting data for PAL and HRV.  Note that the
pattern is similar for the two stations, but varies rapidly with
back azimuth.  Splitting directions for $S_{diff}$ phases (two data points from back azimuth $335\deg$ for HRV,
three data points from back azimuths $257-263\deg$ for PAL) closely match those obtained for $SKS$ phases from the same events. Thus contamination of the $S_{diff}$ signal by the $D''$ anisotropy [{\it Garnero and Lay}, 1997] does not seem to occur along these particular paths. $S$ phases from South American earthquakes with hypocenters deeper then 500 km are included to provide coverage from the south.

Figure 9. Shear-wave splitting data for the subset of PAL and HRV data
that fall within $3 \deg$ of each other, superimposed on one another
(black, HRV; grey, PAL).  Note overall similarity of pattern.

Figure 10. Histogram of the angle between pairs of fast azimuth
measurements for the station HRV.  The distribution is bimodal,
suggesting that two distinct fast azimuths are present.

Figure 11. SKS waves observed at HRV for 4 different backaziumths, A)
SSE, B) WNW, C) NNW and D) NNE. (Left) Observed radial (top row) and
transverse (3rd row) horizontal component seismograms.  The significant
energy observed on the transverse component is an effect of the seismic
anisotropy.  Corrected radial (2nd row) and transverse (bottom row)
component  seismograms, where the effect of propagation through the
anisotropic medium has been removed.  (Right) Particle motion diagrams
{\it (top)} before and {\it (bottom)} after correction. Note that the
energy on the transverse component has been greatly reduced, and the
particle motion made significantly more linear, indicating that the
splitting parameters have been correctly calculated.

Figure 12. Anisotropic tensor orientations for (top) 1 layer and
(bottom) two layer models for HRV.  The fast, intermediate and slow
axes are denoted F, I and S, respectively.  Top and Bottom layers are
denoted T and B, respectively. 

Figure 13. (Left) Shear-wave splitting patterns for best-fitting
one-layer (top) and two-layer (bottom) HRV models. (Right) Observed
(gray) and predicted (black) shear-wave splitting parameters.


Figure 14. Observed and predicted variation of the apparent fast direction
at HRV. Observations are shown by triangles with error bars. A subset
of "robust" data points (circled) was chosen so that $\sigma_\tau <
\tau / 3 $. Clearly, "robust" and "poor" data points follow the same
pattern. Crosses show values of fast direction reported for HRV by
{\it Barruol et al.} [1997] (as given in electronic supplement table
2). Thick lines show patterns predicted by our orthorhombic (solid)
and hexagonal (dotted) models, thin lines show predictions for
equivalent 2-layer splitting operators [{\it Silver and Savage}, 1994].
While all models capture the periodicity of the pattern, the spread of
values and the deviations from the $2{\pi}$ pattern are not matched by
the splitting operator predictions.

Figure 15. A schematic representation of the model for seismic anisotropy distribution under HRV. 

Figure 16. Regional setting of GSN station ARU (open triangle) and the
geometry of inferred anisotropy in the lithosphere. Grey lines -
topography contours at 330 and 560 m.  Wide grey line:  the main
Uralian Fault Zone.  Open arrow: the symmetry axis of
lower-crustal anisotropy. Grey bars: inferred fast axes of anisotropy
for the two layers in the mantle. 

Figure 17. Shear-wave splitting data for ARU (Arti, Russia). See Figure 1 for plotting conventions.

Figure 18.  (Left) Shear-wave splitting patterns for best-fitting one-layer (top)
and three-layer (bottom) ARU models. (Right) Observed (gray) and predicted (black)
shear-wave splitting parameters.

Figure 19. Symmetry axes of best-fitting hexagonal tensors for ARU.  Crust and mantle are labeled
C and M, respectively. Top layer and bottom layers in the mantle are labeled T and B,
respectively. In the bottom layer of anisotropy the symmetry
axis plunges $40 \deg$ to the east, which is equivalent to an upward tilt of $40 \deg$ to the west.

Figure 20.  A schematic representation of the model for seismic anisotropy distribution under ARU (marked by a flag), with the lowermost crust
anisotropy (dark arrow marked "slow") from the model of {\it Levin and Park}, [1997a].


\end{document}

\end{article}



