Open-access Methodology for correcting spectral ratio absorption estimates for source wave travel paths in DST

Abstract

Downhole Seismic Testing (DST) is utilized for estimating low-strain shear damping ratio (ηs). Estimates of low strain ηs values are very important for predicting and assessing ground amplification during earthquakes. The low strain ηs estimates also provide for a reference of laboratory tests such as the resonant column test. The Spectral Ratio Technique (SRT) is one of the most utilized analysis techniques for processing seismic data to obtain ηs. The SRT is applied in the frequency domain where the spectral ratio of two data sets and relative arrival times are calculated to obtain the quality factor. The quality factor is inversely proportional to the damping ratio. One of the limitations of the SRT is that it doesn’t account for the raypaths of the source waves. In DST the raypaths are of particular importance due to the fact the source waves will refract as they travel from the surface down to the seismic receiver. This paper outlines a technique, so-called FMDSM-SRT, which facilitates adjusting estimated SRT absorption DST values so that source wave raypaths are accounted for. It is of paramount importance to first implement newly optimal estimation algorithms on extensive test bed simulations prior to processing real data sets. This paper outlines the results from processing a challenging test bed simulation with the FMDSM-SRT algorithm.

Keywords:
Downhole seismic testing; Site characterization; Absorption estimation; Fermat’s principle; Spectral ratio technique; Iterative forward modeling

1. Introduction

Downhole Seismic Testing (DST) such as Seismic Cone Penetration Testing (SCPT) (ASTM, 2019; Campanella et al., 1986) is a commonly used geotechnical technique for measuring in-situ (<10-5) shear (Vs) and compression (Vp) wave interval velocities. The in-situ interval velocities form the core (Shear Modulus (G), Poisson’s Ratio (µ) and Young’s modulus (E)) of mathematical theorems to describe the elasticity/plasticity of soils and they are used to predict the soil response (settlement, liquefaction or failure) to imposed loads (whether from foundations, heavy equipment, earthquakes or explosions) (Finn, 1984; Andrus & Stokoe II, 2000; Ishikara, 1982). Baziw (2002) outlined that it was of paramount importance to take into account raypath refraction (Snell’s law and Fermat’s principle) when processing DST data sets for estimating interval velocities. Baziw (2002) also gave the mathematical details of a newly developed algorithm, the so-called Forward Modelling Downhill Simplex Method (FMDSM), which considered Fermat’s principle when estimating DST interval velocities. Baziw & Verbeek (2012, 2014, 2018) published results from processing DST data sets using the FMDSM algorithm where special attention was given to identifying critical soil layers and tomographic imaging. Masters et al. (2024) utilized the FMDSM technique when processing an extensive amount of offshore DST datasets. In Masters et al. (2024) paper it was demonstrated that it was necessary to take into account Fermat’s principle when estimating interval velocities within the top 5m to 10m.

Another important geotechnical parameter estimated with DST is the low-strain shear damping ratio (ηs) (Stewart & Campanella, 1991, 1993; Redpath & Lee, 1986). Stewart & Campanella (1993) state that the DST low-strain ηs is very important for predicting and assessing ground amplification during earthquakes and serves as a reference for medium to large strain ηs values derived from laboratory test such as the resonant column test. The Spectral Ratio Technique (SRT) is the most widely used method for determining DST low strain ηs estimates (Stewart & Campanella, 1991, 1993; Redpath & Lee, 1986; Rebollar, 1984). A major limitation in the application of SRT for processing DST data sets is that the raypath of the source waves is not taken into account. This can lead to significant errors in the estimated damping ratio values.

Baziw & Verbeek (2019, 2021) developed the so-called FMDSMAA algorithm for estimating ηs values from DST seismic data. This algorithm accounts for source wave travel paths when estimating ηs values from peak particle accelerations (maximum amplitudes) of DST seismic traces and it also provides for error estimates. Baziw and Verbeek demonstrated the implementation and performance of the FMDSMAA (AA is the acronym for Absorption Analysis) algorithm by considering challenging test bed simulations and processing real data sets. It was shown that the algorithm was able to obtain accurate absorption estimates which demonstrated the algorithm’s correctness. This paper outlines modifications and additions to the FMDSMAA algorithm so that estimates of absorption from the Spectral Ratio Technique (SRT) are corrected for DST source wave travel paths. In general terms, the FMDSMAA is modified where instead of minimizing the difference between synthetic source wave amplitude ratios and recorded values, the difference between synthetic spectral ratios and measured values are minimized. This technique is referred to as the so-called FMDSM-SRT algorithm.

2. Absorption analysis and the spectral ratio technique

The amplitude of a seismic wave propagating in soils decays due to geometric spreading (due to the change in wave front), apparent attenuation (due to mode conversion, reflection-refraction at an interface) and material losses (intrinsic attenuation or absorption). In general terms, the soil profile acts as both an attenuator and a low pass filter as the seismic wave travels through it (Aki & Richards, 2002; Gibowicz & Kijko, 1994; Sheriff & Geldart, 1982). The signal amplitude A within a homogeneous medium at distance x from the source is related to the amplitude A0 at distance x0 by:

A x = A 0 x 0 / x n e α x x 0 (1)

In Equation (1) it is assumed that the decay is due to only geometric spreading and absorption (α in Equation (1)). The Quality Factor, Q, is related to the absorption coefficient as follows:

Q = π / α λ (2)

In Equation (2)λ is the source wave’s wavelength. The Quality Factor is a desirable term to define the absorption of a medium because it is nondispersive. The logarithmic decrement, δ, and fraction of critical damping or damping ratio, ηs, are expressed as:

δ = ln amplitude amplitude one cycle later = α λ = π / Q (3)
η s = δ 2 π = 1 2 Q (4)

In geotechnical engineering the shear damping ratio is typically defined as ηs= 1/(2Q). The unit of δ is nepers (Np).

Attenuation from refraction at a layer boundary is quantified by the transmission coefficient. The amplitude decay term in Equation (1) (x0/x)n corresponds to geometric spreading where n = 1 based upon the conservation of the energy flux of a traveling seismic wave and spherical divergence. Spherical divergence is given as:

E r = E 0 4 π r 2 (5)

where r is the expanding radius from the source, E0 is the initial source energy and Er is the energy at distance r from the source. The amplitude of the source wave decays proportional to the distance r from the source as outlined as follows:

A r = E 0 4 π r 2 (6)

Some researchers have noted that the amplitude of the seismic wave generally does not decay as 1/r where r is the travel distance, and this is why the exponent n is incorporated into Equation (1).

The apparent attenuation due refraction at an interface is quantified by the transmission coefficient. The transmission coefficient quantifies the loss of energy when transitioning from layer 1 to layer 2 is outlined below in Equation (7) for horizontally polarized shear wave (SH).

T 12 = A 2 A 1 = 2 Z 1 Z 1 + Z 2 = 2 ρ 1 V 1 c o s θ 1 ρ 1 V 1 c o s θ 1 + ρ 2 V 2 c o s θ 2 (7)

In Equation (7)T12 is the transmission coefficient for the source traveling moving from layer 1 to 2, A1 is the amplitude of incident wave, A2 the amplitude of refracted wave, ρi is the medium density of layer i, θ1 denotes the incident angle, θ2 is the refraction angle, V1 is the medium velocity of layer 1, and V2 is the medium velocity layer 2 (note: Zi=ρiVi is the acoustic impedance).

The ability to solve the absorption coefficients demands that several assumptions be made regarding the properties of the medium of propagation. These include the assumptions that the medium is an attenuating, velocity dispersive, causal, and linear system. Some properties of seismic medium absorption is proportional to frequency, absorption is inversely proportional to velocity, absorption is proportional to pore saturation, absorption is proportional to porosity in saturated rocks and soils (e.g., absorption is greater in clean sands than shales), absorption is dependent upon pore fluid type (e.g., absorption is greater for brine-watered saturated rocks than methane saturated rocks), and absorption is inversely proportional to formation pressure.

There are several mathematical techniques which facilitate the derivation of absorption from seismic time series. The most popular methodology is the spectral ratio technique (SRT) (Rebollar, 1984; Campillo et al., 1985; Blakeslee et al., 1989; Kvamme & Havskov, 1989; Stewart & Campanella, 1993; Howie & Amini, 2005).

2.1 Spectral ratio technique equation derivation

The mathematical representations of an elastic and anelastic plane wave in the depth domain (z) and frequency domain (ω=2πf) are outlined as follows:

Elastic Plane Wave:

u z , ω = u 0, ω e i k z (8)

Anelastic Medium:

u z , ω = u 0, ω e i k z e α ω z (9)

Parameter k in Equation (8) is the wavenumber (Sheriff & Geldart, 1982). Equation (9) implies that anelastic attenuation is a linear phenomenon. As was previously stated, α denotes the absorption and is defined as:

α ω = ω / 2 V ω Q ω (10)

where V(ω) - dispersive velocity and Q(ω) - dispersive quality factor

Experimental data have shown that Q is almost identically constant with respect to frequency over a large bandwidth (10mHz to 10MHz) (Toksoz et al., 1979). This bandwidth encompasses frequencies encountered in seismic wave propagation. Equation (9) can also be rewritten in terms of a complex wave number as

u z , ω = u z 0 , ω e i k ω z z 0 , k ω = ω / V ω + i α ω (11)

Substituting Equation (10) into (11) and taking real parts (to have physical meaning) gives

ln u z , ω u z 0 , ω = ω z z 0 2 V ω Q = ω T T ω 2 Q (12)

In Equation (12)TT(ω) is the travel time difference between depths z and z0. The arrival time difference is strictly a function of frequency; however if the distance (z-z0) is not to large and V(ω) is not to dispersive, then TT(ω) will be approximately constant with respect to frequency; therefore, Q can be determined from the dominant frequencies of the wave understudy, the relative arrival time and the natural logarithm of the spectral ratio by Equation (12) as follows:

Q = T T 2 × S R S , S R S = ln u z , ω u z 0 , ω / ω (13)

The propagating waves are also affected by non-intrinsic amplitude variation (i.e., frequency independent) such as geometric spreading. The non-intrinsic amplitude variations are lumped into the intercept term of Equation (12) to give

ln N z u z , ω N z 0 u z 0 , ω = ω T ω 2 Q + ln N '
where : N = N z N Z 0 (14)

In general terms, Equation (14) indicates that attenuation due to geometric spreading and apparent attenuation are included in the intercept term and do not affect the intrinsic Q estimation derivation provided they are independent of frequency.

A significant problem with implementing the SRT on DST data sets is that the SRT assumes the source waves travel along the same path as shown in Figure 1. Source wave trajectories adhere to Fermat’s principle, meaning the raypath travels along the trajectory that minimizes the travel time between points. This results in raypath refraction as shown in Figure 2. The schematic illustrated in Figure 2 outlines a DST investigation where seismic source waves are generated at the surface. The source waves travel down into the soil profile where they are recorded by seismic sensors vertically located at increasing depths. Parameter X denotes the source radial offset from the vertically installed seismic sensors. Parameters Di, Vi, and αi denote the depth of the seismic sensor, interval velocity and absorption value for soil layer i. respectively. Finally, in Figure 2dij is the travel path for each source wave j within layer i. As is shown in Figure 2, the source waves refract at the soil interfaces according to Snell’s law and Fermat’s principle.

Figure 1
Assumed source wave travel paths when implementing SRT.
Figure 2
Near surface DST where SRT assumptions are not valid.

3. FMDSM-SRT algorithm

Baziw Consulting Engineers Ltd. (BCE) developed the so-called FMDSM-SRT algorithm for correcting DST SRT absorption estimates for source wave travel paths. The FMDSM-SRT utilizes iterative forward modeling (IFM) for parameter estimation. IFM is based on iteratively adjusting the parameters until a user specified cost function is minimized. The IFM technique implemented by the FMDSM-SRT algorithm is the downhill simplex method (DSM) developed by Nelder & Mead (1965). The DSM searches for the minimum of the costs function by taking a series of steps, each time moving a point in the simplex away from where the cost function is largest. The simplex moves in space by variously reflecting, expanding, contracting, or shrinking. The simplex size is continuously changed and mostly diminished, so that finally it is small enough to contain the minimum with the desired accuracy.

Equations (15) to (18) define the expected intrinsic source amplitudes for traces recorded at depths D1, D2, D3, and D4 for the source raypaths illustrated in Figure 2.

A 1 = e α 1 d 11 x 0 (15)
A 2 = e α 1 d 12 x 0 + α 2 d 22 (16)
A 3 = e α 1 d 13 x 0 + α 2 d 23 + α 3 d 33 (17)
A 4 = e α 1 d 14 x 0 + α 2 d 24 + α 3 d 34 + α 4 d 44 (18)

In general terms, the intrinsic amplitudes recorded at each subsequent DST depth of acquisition are mathematically expressed as follows:

A i = A 0 e α 1 d 1 i x 0 + j = 2 i α j d j i , j 1, i > 1 (19)

If you then consider the ratio of the amplitudes, whether in absolute terms or globally normalized, the unknown amplitude initial A0 drops out of the set of equations. The equation describing the intrinsic amplitude ratio estimated DST absorption values from the measured seismic traces is outlined below.

A i + 1 m A i m = e α m i D i + 1 D i (20)

where αmi is the estimated SRT absorption value for interval i and Di denotes the total travel distance for raypath i.

SRT absorption estimates αmi are obtained from the real data sets by substituting Equation (10) into Equation (13)

α m i = ln u D i + 1 , ω u D i , ω / D i + 1 D i = S R S * Δ ω / D i + 1 D i (21)

Based on the above the proposed FMDSM-SRT algorithm for estimating absorption coefficients can then be described as follows:

  • Utilizing the standard interval velocity FMDSM technique (Baziw, 2002; Baziw & Verbeek, 2012, 2014), obtain estimates of Vi and dij.

  • For the depth increments under analysis determine the measured SRT absorption estimates (αmi - using Equation 21) from the recorded DST traces each depth increment.

  • Implement FMDSM-SRT algorithm to calculate the synthesized amplitudes with Equation (19) based on assumed absorption coefficients, whereby the absolute difference between the measured and synthesizes amplitude ratios is iteratively minimized.

    minαii=1n1AiAi+1eαmiDi+1Di(22)

where n is the number of layers or absorption coefficients to be estimated.

4. FMDSM-SRT test bed simulation

The test bed simulation outlined in this paper uses the same one outlined by Baziw & Verbeek (2019). However, in this case measured SRT absorption values are utilized instead of PPAs. This was done for consistency and comparison of results with the FMDSMAA algorithm. Table 1 below outlines the FMDSM-SRT test bed simulation working parameters. Columns 1 and 2 of Table 1 are the assumed depths and associated source wave arrival times of the simulated DST. Column 3 of Table 1 contains the estimated interval velocities with the depths and arrival times input into the FMDSM algorithm, which accounts for raypath refraction. Column 4 of Table 1 outlines the true absorption values for each soil interval. Column 5 of Table 1 contains the wavelength of the source wave as it travels through the soil layers of varying interval velocities (λ = V/f). Column 6 of Table 1 outlines the quality factor obtained by implementing Equation 2. Finally, column 7 of Table 1 contains the associated estimated measured SRT absorption values for each depth increment. The SRT measured absorption values were calculated by implementing Equations 19 and 23. As is evident from reviewing column 7, there are numerous nonsensical absorption values such as the negative values at 3m and 5m.

Table 1
FMDSM-SRT test bed example parameters (f = 100 Hz).
α m i = ln A i + 1 A i / D i + 1 D i (23)

Table 2 and Figure 3 outline the output of the FMDSM-SRT algorithm after inputting the values given in Table 1, and the raypaths determined by the FMDSM algorithm which are illustrated in Figure 3. Figure 3 also illustrates the estimated FMDSM-SRT algorithm Q values. Column 4 of Table 2 outlines the very large percent differences between the true absorption values and the DST measured SRT values. For example, depth intervals to 2m to 3m and 4m to 5m have percent differences of 266 and 299, respectively. Column 5 of Table 2 contains the FMDSM-SRT estimate, α/, while column 6 outlines the percentage differences between the true absorption values and the FMDSM-SRT estimate. As is evident from the precent differences outlined in Table 2, the FMDSM-SRT was able to correct the DST SRT values for the source wave raypaths and obtain accurate estimates of the true absorption values. Column 7 of Table 2 outlines the low FMDSM-SRT error residuals.

Table 2
FMDSM-SRT test bed example results.
Figure 3
Output from the FMDSM-SRT algorithm illustrating source wave raypaths and estimated Q values.

5. Conclusion

The low-strain DST ηs estimates are very important for predicting and assessing ground amplification during earthquakes. In addition, DST ηs estimates provide for a reference of laboratory tests such as the resonant column test. RCT can suffer from various disadvantages such as sample disturbance and sample preparation which can shift the estimated ηs values. The RCT results can be adjusted so that the low-strain RCT ηs estimates agree with the low-strain in-situ DST ηs estimates. The most popular methodology for estimating DST is the spectral ratio technique. The major shortcoming of the SRT on DST data sets is that the SRT assumes that the source waves travel along the same path. This assumption can result in very large errors in the ηs estimates. This paper has outlined the mathematical details of a new technique, the so-called FMDSM-SRT algorithm, which corrects the estimated DST SRT ηs estimates for source wave raypaths. The effectiveness of the FMDSM-SRT algorithm was demonstrated by a challenging test bed simulation. It is of paramount importance that optimal estimation algorithms first carry out sufficiently challenging test bed simulations for algorithm verification prior to implementation on real data sets.

List of symbols and abbreviations

k wavenumber

u plane wave

x0 source wave reference distance

A0 source wave reference amplitude

DST downhole seismic testing

G shear modulus

Q quality factor

PPA peak particle acceleration

RCT resonant column test

SRS spectral ratio slope

SRT spectral ratio technique

T transmission coefficient

TT travel time

V source wave velocity

α absorption coefficient

α/ FMDSM-SRT estimated absorption coefficient

αm SRT measured absorption coefficient

δ logarithmic decrement

η damping ratio

λ source wave’s wavelength

Data availability

The datasets generated analyzed during the current study are available from the corresponding author upon request.

  • Discussion open until August 31, 2025.
  • Declaration of use of generative artificial intelligence
    This work was not prepared with the assistance of Generative Artificial Intelligence (GenAI)

References

  • Aki, K., & Richards, P.G. (2002). Quantitative seismology (2nd ed.). University Science Books.
  • Andrus, R.D., & Stokoe II, K.H. (2000). Liquefaction resistance of soils from shear-wave velocity. Journal of Geotechnical and Geoenvironmental Engineering, 126(11), 1015-1025. http://dx.doi.org/10.1061/(asce)1090-0241(2000)126:11(1015)
    » http://dx.doi.org/10.1061/(asce)1090-0241(2000)126:11(1015)
  • ASTM D7400/D7400M-19. (2019). Standard Test Methods for Downhole Seismic Testing ASTM International. Retrieved in September 23, 2024, from https://www.astm.org/d7400_d7400m-19.html
    » https://www.astm.org/d7400_d7400m-19.html
  • Baziw, E., & Verbeek, G. (2012). Deriving interval velocities from downhole seismic data. In Proceedings of the 4th International Conference on Geotechnical Site Characterization (ISC-4) (pp. 1019-1024). Leiden: Balkema.
  • Baziw, E., & Verbeek, G. (2014). Identifying critical layers using SCPT and seismic source moveout. In Proceedings of the 3rd International Symposium on Cone Penetration Testing (pp. 357-364). Leiden: Balkema.
  • Baziw, E., & Verbeek, G. (2018). NMO-SCTT: A unique SCPT tomographic imaging algorithm. In Proceedings of the CPT’18 Conference (pp. 137–142). Leiden: Balkema. https://doi.org/10.1201/9780429505980-14
    » https://doi.org/10.1201/9780429505980-14
  • Baziw, E., & Verbeek, G. (2019). Unique algorithm which takes into account source wave travel paths when estimating absorption values from DST data sets. In Proceedings of the DFI 44th Annual Conference on Deep Foundations (pp. 468-479). New York: Taylor & Francis Group.
  • Baziw, E., & Verbeek, G. (2021). Implementation of the FMDSMAA algorithm. In Proceedings of the 6th International Conference on Geotechnical Site Characterization (ISC-6) (pp. 139-145). Leiden: Balkema.
  • Baziw, E.J. (2002). Derivation of seismic cone interval velocities utilizing forward modeling and the downhill simplex method. Canadian Geotechnical Journal, 39(5), 1181-1192. http://dx.doi.org/10.1139/t02-061
    » http://dx.doi.org/10.1139/t02-061
  • Blakeslee, S., Malin, P., & Alvarez, M. (1989). Fault-zone attenuation of high-frequency seismic waves. Geophysical Research Letters, 16(11), 1321-1324. http://dx.doi.org/10.1029/GL016i011p01321
    » http://dx.doi.org/10.1029/GL016i011p01321
  • Campanella, R.G., Robertson, P.K., Gillespie, D., & Clemence, S.P. (1986). Seismic cone penetration test. Geotechnical Special Publication, 116-130. Retrieved in September 23, 2024, from https://cir.nii.ac.jp/crid/1574231874283762432
    » https://cir.nii.ac.jp/crid/1574231874283762432
  • Campillo, M., Plantet, J., & Bouchon, M. (1985). Frequency-dependent attenuation in the crust beneath central France from Lg waves: data analysis and numerical modeling. Bulletin of the Seismological Society of America, 75(5), 1395-1411. http://dx.doi.org/10.1785/BSSA0750051395.
    » https://doi.org/ http://dx.doi.org/10.1785/BSSA0750051395
  • Finn, W.L. (1984). Dynamic response analysis of soils in engineering practice. In C.S. Desai & R.H. Gallagher (Eds.), Mechanics of Engineering Materials (pp. 253-275). John Wiley & Sons Ltd..
  • Gibowicz, S.J., & Kijko, A. (1994). An introduction to mining seismology Academic Press Inc.
  • Howie, J.A., & Amini, A. (2005). Numerical simulation of seismic cone signals. Canadian Geotechnical Journal, 42(2), 574-586. http://dx.doi.org/10.1139/T04-120
    » http://dx.doi.org/10.1139/T04-120
  • Ishikara, K. (1982). Evaluation of soil properties for use in earthquake response analysis. In International Symposium on Numerical Models in Geomechanics (pp. 237-259). Londres: Springer-Verlag.
  • Kvamme, L.B., & Havskov, J. (1989). Q in Southern Norway. Bulletin of the Seismological Society of America, 79(5), 1575-1588.
  • Masters, T., Cazarez, I., Tejani, N., & Verbeek, G. (2024). Evaluation of the challenges present in obtaining, processing, and interpreting useful data from offshore seismic cone penetration testing. In Proceedings of the 7th International Conference on Geotechnical Site Characterization (ISC-7) (pp. 2176-2183). CIMNE.
  • Nelder, J.A., & Mead, R. (1965). A simplex method for function optimization. The Computer Journal, 7, 308-313. http://dx.doi.org/10.1093/comjnl/7.4.308
    » http://dx.doi.org/10.1093/comjnl/7.4.308
  • Rebollar, C.J. (1984). Calculation of Qβ using spectral ratio method in Northern Baja California, Mexico. Bulletin of the Seismological Society of America, 74(1), 1371-1382. http://dx.doi.org/10.1785/BSSA0740010091
    » http://dx.doi.org/10.1785/BSSA0740010091
  • Redpath, B.B., & Lee, R.C. (1986). In-situ measurements of shear wave attenuation at a strong motion recording site Report Prepared for USGS Contract No. 14-08-001-21823. Retrieved in September 23, 2024, from https://www.usgs.gov/
    » https://www.usgs.gov/
  • Sheriff, R.E., & Geldart, L.P. (1982). Exploration Seismology (2nd ed.). Cambridge University Press.
  • Stewart, W.P., & Campanella, R.G. (1991). In situ measurement of damping of soils. In Proceedings of the 2nd International Conference on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics (pp. 83-92). The University of Missouri-Rolla.
  • Stewart, W.P., & Campanella, R.G. (1993). Practical aspects of in-situ measurements of material damping with the seismic cone penetration test. Canadian Geotechnical Journal, 30(2), 211-219. http://dx.doi.org/10.1139/t93-018
    » http://dx.doi.org/10.1139/t93-018
  • Toksoz, M.N., Johnston, D.H., & Timur, A. (1979). Attenuation of seismic waves in dry and saturated rocks; I, laboratory measurements. Geophysics, 44(4), 681-690. http://dx.doi.org/10.1190/1.1440969
    » http://dx.doi.org/10.1190/1.1440969

Edited by

Publication Dates

  • Publication in this collection
    17 Mar 2025
  • Date of issue
    2025

History

  • Received
    23 Sept 2024
  • Accepted
    23 Jan 2025
location_on
Associação Brasileira de Mecânica dos Solos Av. Queiroz Filho, 1700 - Torre A, Sala 106, Cep: 05319-000, Tel: (11) 3833-0023 - São Paulo - SP - Brazil
E-mail: secretariat@soilsandrocks.com
rss_feed Acompanhe os números deste periódico no seu leitor de RSS
Reportar erro