Research ArticleGEOPHYSICS

Scattered wave imaging of the oceanic plate in Cascadia

See allHide authors and affiliations

Science Advances  14 Feb 2018:
Vol. 4, no. 2, eaao1908
DOI: 10.1126/sciadv.aao1908

Abstract

Fifty years after plate tectonic theory was developed, the defining mechanism of the plate is still widely debated. The relatively short, simple history of young ocean lithosphere makes it an ideal place to determine the property that defines a plate, yet the remoteness and harshness of the seafloor have made precise imaging challenging. We use S-to-P receiver functions to image discontinuities beneath newly formed lithosphere at the Juan de Fuca and Gorda Ridges. We image a strong negative discontinuity at the base of the plate increasing from 20 to 45 km depth beneath the 0- to 10-million-year-old seafloor and a positive discontinuity at the onset of melting at 90 to 130 km depth. Comparison with geodynamic models and experimental constraints indicates that the observed discontinuities cannot easily be reconciled with subsolidus mechanisms. Instead, partial melt may be required, which would decrease mantle viscosity and define the young oceanic plate.

INTRODUCTION

Plate tectonics relies on the transition from a rigid lithospheric plate to a weaker, deeper asthenosphere at the lithosphere-asthenosphere boundary (LAB). The LAB is classically defined thermally, in which case a gradual LAB transition is predicted. Recent observations of sharp seismic discontinuities challenge the classic thermal definition of plates, arguing for a melt- or hydration-defined lithosphere-asthenosphere transition (13). It has also been suggested that the discontinuities are unrelated to the base of the plate, representing frozen-in anisotropy (46), an amplified effect of hydration on seismic waves owing to grain boundary sliding (7) or an enhanced effect of near-solidus conditions on seismic waves (8). Distinguishing the defining mechanism of the plate has proved challenging, perhaps because well-resolved discontinuity results are intermittent and mostly from older lithosphere (912) with more complicated tectonic histories. Observations from young seafloor are lacking.

The plate definition beneath young ocean lithosphere is particularly important for understanding melt dynamics and transport and the formation of the plates. In this region, the mantle is predicted to rise over a broad zone (hundreds of kilometers wide), but melt somehow focuses to linear volcanic chains only tens of kilometers wide at the mid-ocean ridge. The focusing is not easily explained by simple viscous flow models assuming a thermally defined plate (13). Melt could travel along a permeability boundary at the base of a plate (14). In general, tighter seismic constraints beneath the ridges are required to evaluate melt dynamics.

The large Cascadia Initiative Amphibious Array on the west coast and offshore of North America covering the Juan de Fuca and Gorda Ridges is an excellent opportunity to investigate the young oceanic plate in greater detail (Fig. 1) (15). In Cascadia, surface wave tomography images a seismically fast lithospheric plate above slower asthenospheric velocities (1618), either increasing in thickness with age (18) or at a relatively constant depth (17). Teleseismic S-wave tomography finds faster velocities with age, consistent with a thickening lithosphere (19). Constraints on sharp discontinuities from receiver functions are complicated by sediment and water reverberations (20). P-to-S receiver functions from one station on the Juan de Fuca Ridge suggest a conversion from the base of a 10- to 15-km-thick anisotropic layer (21). However, the exact thickness of the plate and the sharpness of the LAB transition have yet to be established beneath the ridge or the plate as it ages out to ~10 million years (My).

Fig. 1 Map of the study region.

Background colors shows bathymetry/topography. Star indicates Axial Seamount, which is the surface expression of the Cobb Hotspot. Inverted triangles show seismometer locations on the ocean bottom (red) and on land (blue). Thin black lines show the locations of cross sections in Fig. 2. Black circles along the cross sections correspond to a spacing of 100 km. Thick black lines show plate boundaries, and serrated line shows the trench. FZ, fracture zone.

Here, we present S-to-P receiver function imaging of seismic discontinuity structure beneath the Juan de Fuca and Gorda Ridges and Plates before they subduct beneath North America at the Cascadia subduction zone and also beneath western North America. Ocean data can be complicated and contaminated by several factors including instrument tilt, compliance, sediment and water column conversions and reverberations, and inaccurate sediment property assumptions during rotation and migration. We perform a series of calculations and corrections to both account for these and ensure that we are not interpreting artifacts. We rotate the data into theoretical P- and S-wave components, perform an extended time multitaper deconvolution, and migrate to depth in three dimensions (see Materials and Methods and the Supplementary Materials) (22, 23). We also stack waveforms according to seafloor age in 2-My bins (24).

We determine the features required by the data by comparing to synthetic receiver function predictions from geodynamic models. We perform two-dimensional (2D) geodynamic modeling of the thermal structure and mantle flow field beneath the mid-ocean ridge system assuming pressure- and temperature-dependent viscosity (see Materials and Methods and section S7). We test the models with both dry and mildly hydrated mantle and a range of melt retentions and potential temperatures. We test a range of lateral scales of upwelling by allowing the change in spreading direction across the ridge to occur linearly over a specified distance, which effectively imposes a region of constant strain rate at the ridge. We translate the models to velocity (25) and calculate the predicted synthetic seismograms, deconvolving and migrating to depth in the same way that the data is treated.

RESULTS

Beneath the continents, we image the expected structures (Fig. 2). We image the Moho beneath the continent at 30 ± 5 to 40 ± 5 km, in good agreement with previous receiver function studies (26, 27). We also image a negative phase at 60 ± 5 to 80 ± 5 km depth beneath the continents, in good agreement with P-to-S and S-to-P receiver functions (26) and slightly shallower than the discontinuity at 90 km imaged by S-to-P at longer periods (27), which likely represents the LAB (26, 27).

Fig. 2 Receiver functions and comparison to predictions from geodynamic models.

(A to C) Receiver function cross sections as located in Fig. 1. Red phases correspond to velocity increases with depth. Blue phases correspond to velocity decreases with depth. Inverted triangles show the seismometer locations on the ocean bottom (red) and on land (blue) within 0.5° of the transect and plotted at the appropriate bathymetry/topography (Topo) (black line). Black circles at 125 km depth correspond to 100 km spacing along the cross section. JdF Ridge, Juan de Fuca Ridge. (D) Receiver function predictions for the geodynamic model of the Gorda Plate with up to 1% melt volume retention. (E) Receiver function prediction for the geodynamic model of the Gorda Plate with no melt retention. Seafloor age is indicated (24). (F) Geodynamic model. Background color shows the temperature. Contours show retained melt volume fraction.

We image a positive phase at 0 to 4 km depth beneath the sea surface, shallower than the base of a 6- to 7-km-thick oceanic crust (Fig. 2) (28). This is expected owing to constructive interference from conversions at the sediment-crust interface, the mid-crustal discontinuity, and the oceanic Moho. We also image a negative phase, a velocity decrease with depth, at 20 ± 5 to 45 ± 5 km depth beneath the sea surface. It is shallowest beneath the ridges at 20 ± 5 km depth. It increases in depth with distance from the ridge, reaching 40 ± 5 to 45 ± 5 km beneath the 10-My-old lithosphere at the continental margin. The deepest LAB observations are not an artifact of the accretionary wedge (section S5). The age-depth dependence is best resolved in the 3D model as seen in individual transects (Fig. 2) and map view (Fig. 3), and we report the total range from this model. However, an age-depth relationship is also visible in large age bin stacks (Fig. 4A) and averages of the depths in the 3D model (Fig. 4B). The depth variability in large age bins is more muted than the extremes of the 3D model, as expected, owing to lateral depth variability and also waveform sensitivity. Beneath the oceans, we also image a negative discontinuity in the 60 ± 5 to 80 ± 5 km depth range, particularly beneath fracture zones such as Mendocino and parts of Blanco (fig. S1). It is the strongest in amplitude in the south of Mendocino where it migrates to 70 ± 5 to 80 ± 5 km depth. This feature could either represent a discontinuity or an artifact (section S2). Finally, we image a deep positive phase at 90 ± 10 to 130 ± 10 km depth, particularly beneath and/or near the ridges. It is the shallowest by the Gorda Ridge and the deepest beneath the Axial Seamount, which is the current surface realization of the Cobb Hotspot.

Fig. 3 Map view of receiver function discontinuity depths.

(A) Shallow negative ocean LAB discontinuity. (B) Positive discontinuity at the base of the melt triangle. Axial Seamount, which is the surface realization of the Cobb Hotspot (yellow star), Juan de Fuca Ridge (JdFR), Gorda Ridge (GR), Blanco Fracture Zone (BFZ), and the Mendocino Fracture Zone (MFZ) are indicated on the map. Thin black lines show the locations of cross sections in Fig. 2. Thick black lines show plate boundaries.

Fig. 4 Averaged age-depth relationship of the LAB.

(A) Colored lines show large (2 My) age bin stacks of S-to-P receiver functions listed at the bin center. Amplitude is normalized to minimize the effects of lateral heterogeneity and emphasize the age-depth trend. (B) Receiver function LAB depths from averages of individual bins in the 3D model (red circles) and from waveform stacks of large age bins over the entire study area (black circles) and over just the Gorda Plate (cyan circles) compared to the predicted receiver functions for the geodynamic model (1350°C) that includes melt retention and upwelling over a broad zone, presented as background shading [red (Moho) and blue (LAB)]. A 1200°C isotherm from a half-space cooling model is also shown (green line). We present Gorda separately because the age bin stacks are dominated by Juan de Fuca.

DISCUSSION

The depth of the negative phase at 20 to 45 km generally coincides with the gradual drop in velocity from the seismically fast lithosphere to the slower asthenosphere in teleseismic surface waves (18) and full-waveform tomography of ambient noise and local earthquakes (17). Our waveforms are only sensitive to changes in velocity rather than in absolute velocity. However, our slowest velocity (3.9 km/s) is similar to that from surface waves [~3.9 km/s (17) and ~4.0 km/s (18)]. The age-depth trend of the negative phase generally agrees with thickening of the seismically fast plate with age imaged with surface waves beneath the Juan de Fuca Plate (1618) and inferred from teleseismic S-wave tomography (19). The deepening of our phase beneath the aging Gorda Plate agrees with subtle thickening in one result (18) and inferred subtle thickening from body waves (19) but agrees less with the constant thickness of another (17). The overall agreement with the gradual drop in velocity from surface waves suggests that the phase represents the LAB at the base of the plate.

Beneath the ridge, our LAB discontinuity at 20 ± 5 km agrees well with a discontinuity in P-to-S modeling, at ~21 to 26 km beneath the sea surface, interpreted as the base of a 10- to 15-km-thick mantle layer (21). The result suggests that a fast lithosphere exists beneath the 0-My-old lithosphere at the ridge axis, which is also consistent with imaging from the full-waveform tomography of ambient noise and local earthquakes (17). It agrees less well with the lack of a fast lid observed by both teleseismic (18) and ambient noise (16) surface waves. One possibility is that a fast lithosphere exists but is too thin to be imaged by the latter methods. A seismically fast lithosphere at the ridge axis is predicted by geodynamic modeling of intermediate-to-slow spreading ridges that incorporates lateral cooling (29).

In the 3D receiver function model, the LAB phase remains constant in strength until just west of the continental margin. This strength is in general agreement with the strong velocity anomaly in body-wave tomography interpreted as the accumulated weak buoyant and partially melted mantle that has traveled up along the slab (30). Approaching the continent, our discontinuity loses amplitude just west (by 1° to 2°) of where the anomaly also tapers off in the body-wave prediction.

We do not image the dipping, negative discontinuity at 50 to 125 km depth at the continental margin near ~42°N using land stations and multimode P-to-S conversion techniques likely related to the descending subducting slab (31, 32). Instead, the LAB signals are significantly reduced or nonexistent in the ocean-continent transition region. This is because our S-to-P receiver function gridding and smoothing scheme favors flat- to gradual-dipping discontinuities. Previous S-to-P and P-to-S results from beneath the continent using a similar methodology to ours did not detect this dipping feature (26, 27).

The positive phase beneath the ridge at 90 to 130 km is in agreement with the increase in velocity with depth that occurs in surface wave models at >50 km depth, that is, beneath the low-velocity zone interpreted as containing a small amount of partial melt beneath the Juan de Fuca Ridge (1618). The phase occurs over a broad 1.5° × 1.5° area that is offset just west of the Juan de Fuca Ridge in the region of the Cobb Hotspot, in agreement with the location of the strongest anomaly in surface wave (18) and teleseismic S-wave tomography (19). We also image this feature over a narrower 0.5° wide zone beneath the Gorda Ridge, where surface waves find less evidence for a strong melt signature (1618). This could be caused by the broader lateral sensitivity of surface waves in comparison to receiver functions. The depth of the phase beneath the ridges (95 ± 10 to 115 ± 10 km) is consistent with the expected solidus depth for a slightly damp mantle [about 100 parts per million (ppm) water] and a potential temperature of 1350° to 1375°C. On average, the phase is deeper (125 ± 10 km) near the Cobb Hotspot, which is consistent with a higher mantle potential temperature (by ~50°C), additional hydration (100 ppm), or some combination of the two. The broader area of the phase near the Cobb Hotspot is also consistent with the effects of a broader hotspot anomaly. Our result agrees with geochemical estimates of mantle melting, which suggest greater temperatures in the vicinity of the Cobb Hotspot with deeper, more extensive melting and a broader melting column, relative to the adjacent Juan de Fuca Ridge (33, 34).

Comparison with the seismic predictions from geodynamic modeling shows that the strong negative discontinuity at 20 to 45 km depth cannot be explained by a purely thermally defined plate (Fig. 2). The observed age-depth trend with only mild thickening with age beneath young seafloor is inconsistent with half-space cooling (Fig. 4B) and also cannot be explained by passive upwelling with lateral heat conduction. Instead, a broad zone of upwelling is needed in our geodynamic models (Fig. 4B). The broad upwelling and weak age dependence at the youngest ages needed in our models may be related to the hypothesized upwelling from beneath the slab (30). In addition, the predicted LAB velocity gradient in depth from a thermal model is very gradual. The synthetic receiver function cross sections calculated from this model do not match the strength of our observed phase (Fig. 2, A and E). The thermal model also does not predict the age-depth dependence of the discontinuity in the receiver functions, an overall effect of a gradient that is too gradual to give a strong response at depth beneath older ages. Finally, the thermal model does not predict the subtle positive phase at 95 to 115 km beneath the ridges and at 125 km beneath the Cobb Hotspot.

In addition, alternative subsolidus mechanisms cannot unambiguously explain our result. Variability in polarity with azimuth and/or focal mechanism is expected for azimuthal anisotropy. There is no evidence for this in our data. Radial anisotropy caused by aligned olivine and/or compositional layering cannot explain a velocity contrast of this magnitude (35). If an enhanced effect of near-solidus conditions (8) is included, then the predicted negative receiver function phase beneath the 10-My-old seafloor is too shallow (26 km) to explain our observations (45 km). Elastically accommodated grain boundary sliding predicts an increase in the sharpness of the LAB velocity gradient with age (7), which would be realized as larger amplitudes and/or more impulsive phases at older ages. This is not observed (Fig. 2). Finally, none of these mechanisms predict the subtle velocity increase with depth at 95 to 115 km beneath the ridges and at 125 km beneath the Cobb Hotspot. Therefore, melt is our preferred mechanism given the current knowledge, recognizing that near-solidus conditions and elastically accommodated grain boundary sliding are areas of ongoing research that may evolve.

Instead, the strength of the S-to-P LAB phase is more consistent with geodynamic models with a small amount of retained partial melt beneath the plate. The melt gives rise to a strong sharp velocity gradient at the shallow limit of the melt triangle solidus that explains the observed receiver functions (Fig. 2). Retained partial melt beneath the ridge is also supported by slow surface wave velocities (18) and high attenuation (36). Melt retained along this boundary in models that include a broad 50-km-wide mantle upwelling reproduces the observed age-depth trend in the data. Melt also decreases mantle viscosity, and its existence likely defines the plate (2) in its young oceanic form beneath Cascadia.

In addition, geodynamic models with retained melt explain the positive phase at 90 to 130 km depth as conversions from the base of the melt triangle. Shallower observations beneath the Gorda Ridge and deeper observations beneath the Cobb Hotspot are well predicted given the hotter and/or more hydrated hotspot conditions, as previously discussed. The LAB discontinuity at 20 to 45 km depth is stronger than the deeper phase at 90 to 130 km as predicted in the geodynamic models. In the models, this is the result of melt retention over the broad melt triangle beneath the plate, with greater concentrations just beneath the plate. Some mechanism must exist to focus this melt to the ridge axis where most volcanism occurs. One possibility is that melt travels along the base of the plate to the axis of volcanism where it is erupted at the mid-ocean ridge.

MATERIALS AND METHODS

S-to-P receiver function imaging

We used events located at epicentral distances 55° to 80° with magnitudes >6 for land stations and >5.8 for oceans. The lower threshold for ocean stations was to ensure that we considered all data that might be usable from the smaller oceanic data set.

Tilt and compliance corrections were applied to the ocean-bottom data (37, 38). We used a 5-hour noise window before each event to calculate the two transfer functions, vertical versus horizontal and also vertical versus pressure. We applied the tilt correction using a smoothed transfer function. The correction was applied to a low-frequency band where coherence was >0.8, as discussed by Bell et al. (38). For compliance, we again used the low-frequency band where coherence is >0.8. The band where the correction was applied was typically <0.05 Hz, depending on water depth.

The data were rotated into the P- and S-wave components using a transformation matrix, either free surface or oceanic depending on whether the station was located on land or ocean (39). For land stations, we assumed that P-wave velocity (VP) = 5.5 km/s, S-wave velocity (VS) = 3.2 km/s, and density = 2900 kg/m3. For ocean stations, we used sediment velocities determined at each station from P-to-S delay times from the sediment-crust conversion (see section S1) and density = 1000 kg/m3. The waveforms were bandpass-filtered from 0.04 to 0.25 Hz.

We visually inspected both the P- and S-wave components near the theoretical S-wave arrival on 18,968 waveforms, 2522 of which were from ocean-bottom stations. We selected waveforms with visible arrivals on the S-wave component (fig. S2). The S-wave amplitude should be larger than amplitudes before or after its arrival. Some energy was predicted on the P-wave component owing to strong conversions from the base of a sediment layer. The S-wave arrival should also be within 10 s of the theoretical arrival (40). We hand-selected a window around these visible S-wave arrivals to use as the source waveform in the deconvolution (fig. S2). After handpicking the data, we were left with 343 waveforms from ocean-bottom stations and 4678 waveforms recorded on land.

The data were deconvolved using an extended time multitaper method (22, 23). The deconvolved signals were multiplied by −1 so that polarity was consistent with that typical for P-to-S receiver functions, where positive phases correspond to velocity increases with depth (fig. S2). The waveforms were then migrated and stacked on a 0.5° × 0.5° grid with a depth spacing of 1 km. The migration model for land stations assumed IASP91. For ocean stations, we used the sediment thickness and velocity at each station based on P-to-S delay times (see section S1) and an expected relationship between sediment thickness and VS (41) and also VP from active source results (42). For delays more than 1.5 s, fixed velocities were assumed (see section S1). We assumed ocean crustal velocities with a linear gradient from 3.0 to 3.7 km/s and VP/VS = 1.75 between 0 and 7 km depth and ocean mantle lithosphere velocities VP = 8.04 and VP/VS = 1.8. We corrected for station elevations in the migration process. We only used bins with more than three waveforms in the stack. We smoothed the bins over the Fresnel zone of the waves. The Fresnel zone was approximated assuming an S-wave period of 7 s and a VS that varies in depth according to the IASP91 model. We implemented a minimum Fresnel zone cutoff of 50 km, which was in effect at the depths of the interpreted LAB.

Error bars represented the maximum variation in depth we expected from extreme migration model variations. For instance, keeping the VS constant, an increase in the assumed mantle VP by 10% would cause depth variations of ~5 km or less for discontinuities in the 20 to 45 km depth range. For deeper discontinuities, the error was on the order of 5 to 10 km, given that velocity anomalies (including very large VP/VS) of this strength are probably not realistic over very large swaths (>100 km) of the mantle or beneath most land stations. We showed discontinuity depths only in cases where the receiver function amplitude exceeded the background noise level determined by histogram evaluations of the amplitude structure, generally these corresponded to contrasts on the order of 5% or greater (Fig. 3 and fig. S1).

Geodynamic modeling

We modeled mantle flow, temperature, melting, and depletion (43) in 2D for a spreading ridge with an imposed half spreading rate of 25 mm/year (see section S7). In the model shown (Fig. 2), we imposed a region of constant strain rate within 50 km of the spreading ridge, which effectively distributed deformation over this region and allowed for a fairly broad mantle upwelling. Our model domain was 800 km wide by 400 km deep, with variably spaced elements, with a minimum element spacing 2 km wide by 1 km deep near the ridge axis. We tested a variety of mantle potential temperatures between 1350° and 1450°C and used a mantle adiabat of 0.3°C/km with 0°C at the top of the model domain. In the model shown, we assumed 100 ppm water in the background mantle and used the melting model of Katz et al. (44). Melt retention was allowed, and melt buoyancy effects were included. We chose a permeability of 1.5 × 10−14 that minimized the effects of melt retention. The models were run to steady state.

We converted the thermal structure to seismic velocities (25), assuming a primarily olivine mantle. The method accounted for attenuation effects based on laboratory scalings, and we assumed a grain size of 20 mm, which was predicted by geodynamic modeling of grain size evolution beneath the oceanic lithosphere (45). We also discuss other choices for grain size (section S6). The seismic velocity structure was calculated for each element in the model. We then calculated the predicted receiver functions from the seismic velocities using a 1D reflectivity code (46) assuming a 6-km-thick crust and a 3-km-deep water column, with the seismometer located on the seafloor.

For purely thermal models (that is, assuming no melt retention), the predicted LAB phases were too weak and also lacked the depth dependence to match our observations (Fig. 2E). The thermal model also could not explain the observed positive phases at 90 to 130 km depth. Including enhanced effects of near-solidus conditions on seismic velocity (8) similarly has an age-depth trend that is more muted (20 to 26 km) than our observation (20 to 45 km) for the 0- to 10-My-old lithosphere. Finally, we also explored the effects of melt on seismic velocity, assuming a linear relationship, where up to 1% retained melt volume yields a 7.9% reduction in shear velocity (47), recognizing that this relationship could be much different depending on melt geometry. We allowed a maximum of 1% melt to be retained in the mantle for the calculation of seismic velocity reduction to minimize the effects of melt. These models show much better agreement with our result (Fig. 2D and fig. S6).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/2/eaao1908/DC1

section S1. Sediment property calculation

section S2. Water column/sediment artifacts

section S3. Alternate migration models

section S4. Dipping layers

section S5. Model reliability

section S6. Comparison between geodynamic models and receiver functions

section S7. Additional geodynamic modeling details

table S1. P-to-S delays measured from Cascadia Initiative Array stations.

fig. S1. Map view of negative receiver function discontinuity at 60 to 80 km.

fig. S2. Example of good-quality ocean-bottom seismograms and receiver functions.

fig. S3. Average delay times at each station between P-wave arrival and P-to-S conversion from the base of the sediment layer.

fig. S4. Sediment VP and VS used in the waveform rotation and migration derived from P-to-S sediment conversion delay times at each station.

fig. S5. Comparison of sediment thickness estimate to global sediment thickness (50).

fig. S6. Example of receiver functions recorded at single stations.

fig. S7. Example of receiver function and geodynamic models.

fig. S8. Comparison of geodynamic models and seismic predictions.

References (4859)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: The seismic data used in this study are available online at the IRIS (Incorporated Research Institutions for Seismology) data management center (http://ds.iris.edu/ds/nodes/dmc) and were downloaded using the SOD (Standing Order for Data) package (http://seis.sc.edu/sod). We thank G. Abers, J. Gaherty, and D. Forsyth for helpful discussions and three anonymous reviewers. Funding: We acknowledge funding from the Natural Environment Research Council (NE/M003507/1 and NE/K010654/1) and the European Research Council (GA 638665). Author contributions: C.A.R. conceived the project, picked the ocean-bottom seismic data, processed the data, performed the synthetics tests, and wrote the paper. N.H. performed tilt and compliance corrections, ran the geodynamic models, translated the models to velocity, calculated expected receiver functions, and contributed to the figures, manuscript, and interpretation. C.A.R. and N.H. characterized the sediment properties. S.T. picked the land data. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article