Summary In this paper we study teleseismic shear-wave splitting due to multiple anisotropic layers. The effects can be described in terms of a multilayered splitting operator that may be calculated using a matrix formulation. As a special case, we consider smooth variations of anisotropy with depth. The model is characterized by the (horizontal)fast-axis directions at the top (φt) and bottom (φb) of the anisotropic region, and the total delay time (Δt). Numerical results show that the particle motion at relatively long periods (T/Δtnbb 5) is clearly elliptical and is a function of the polarization direction of the incident wave. At shorter periods (T/Δtnab1); however, the seismogram is characterized by two distinctly separated arrivals corresponding to the fast- and slow-axis directions at the top of the anisotropic region. Relative amplitudes of the fast and slow arrivals depend on the orientation of the initial polarization with respect to the axis directions at the bottom. Only the high-frequency results agree with predictions based on ray theory. At fixed initial polarization, the longer-period results can be described in terms of an apparent fast-polarization direction (φa) and an apparent delay time (δta). The apparent splitting parameters exhibit a π/2-periodicity as functions of initial polarization (or backazimuth), which is similar to the two-layer case considered by Silver & Savage (1994). Assuming that φt−φb < 45 °, we find variations of φa and δta to be smaller than 10 per cent over most ( ≈ 2/3) of the backazimuth range. In this case, φa≃ (φt+φb)/2 and δta≃Δt. The calculated apparent splitting parameters agree well with direct measurements using synthetic waveforms and may thus be applied to the inversion of observations Long-period apparent splitting parameters for different types of depth variations cannot be resolved uniquely without further constraints from the deformation history at a given station. For example, we find that recently observed azimuthal variations of splitting parameters (Ozalaybey & Savage 1994; Brechner et al. 1998) can be modelled equally well in terms of two distinct layers (a four-parameter model) or smooth variations of fast-axis directions (a three-parameter model). Depth-dependent models with randomly chosen fast-axis directions and delay times in each layer display a characteristic decrease of the most probable apparent delay time, δta (Pmax), with increasing number of layers N. We find that δta (Pmax)∼ 1/√N. The fact that most measurements yield values of δta≃ 1 s (Silver 1996) suggests that the majority of observations can be explained in terms of one or two anisotropic layers, provided the fast-axis directions in different layers are independent.