ASEG Extended Abstracts - 25th International Conference and Exhibition – Interpreting the Past, Discovering the Future, 2016
25th International Conference and Exhibition – Interpreting the Past, Discovering the Future, 2016
- Articles
-
-
-
Heat flow: The neglected potential field for mineral exploration
More LessAuthors Graeme BeardsmoreThis paper argues that heat flow is a valid ‘potential field method’ for mineral exploration. Certain ore deposits, most notably iron oxide copper gold uranium (IOCG-U) deposits, demonstrably affect the local heat flow and ground temperature conditions. The physics of steady state conductive heat flow is mathematically the same as gravitational acceleration, with buried heat sources analogous to buried masses. Detecting the surface signature of a buried heat source can therefore yield direct evidence of an ore body. The physics is robust, appropriate tools exist for the collection of heat flow data for exploration, and some simple strategies could yield valuable additional information about the subsurface for small marginal cost.
-
-
-
-
Back to basics on broadband seismic amplitudes, phase and resolution
More LessAuthors Andrew LongMarine broadband seismic results typically have observed amplitude spectra that are inconsistent with the earth reflectivity spectra measured from spatially coincident well data. The most notable inconsistency is a visual bias towards stronger ultra-low frequency amplitudes in the 0-10 Hz range on seismic data. Whilst useful in principle for more accurate elastic seismic inversion with less dependence upon low frequency model (LFM) building and for more stable full waveform inversion (FWI), the low frequency bias will degrade seismic resolution if not balanced properly during processing. Furthermore, the observed seismic amplitude spectra on pre-stack data vary with increasing offset and the phase can strongly vary in a frequency-dependent manner. These variations are attributed to three first-order processes: 1. Progressive attenuation of higher frequencies with longer travel paths, 2. NMO stretch, and 3. Offset-dependent tuning, and may easily be compounded by inappropriate parameterization of various processing and imaging steps. A key component of elastic seismic inversion is the extraction and scaling of several angle-dependent and depth-dependent wavelets. These wavelets are ideally independent of the aforementioned processes but inevitably they are not, and the angle ranges used in the inversion impact the significance of NMO stretch and offset-dependent tuning. Even if these various considerations can all be accommodated during processing, the resolution of seismic images is limited by the resolution of the velocity model, the scattering assumptions of the imaging algorithms used, and the a priori information used to constrain imaging and inversion.
After reviewing the principles of seismic wavefield propagation in the contemporary context of broadband seismic methods, largely pursued with the ambition of removing sea surface-related ghost effects, I discuss the additional uncertainties introduced into seismic signal processing by broader bandwidth data-notably at the low frequency end. I consider two rather different ‘broadband seismic’ perspectives: 1. The industry must progress towards higher fidelity and resource-consuming measurements of every source event and every receiver measurement in order to effectively deconvolve the ‘system response’ from the data, or, 2. High-end imaging solutions can automatically eliminate the ‘source term’ without requiring detailed source information-assuming that wavefield separation has robustly accounted for dynamic variations in receiver geometry. A longer term consideration of imaging suggests that the classical sequential processing paradigm is in fact dead, that the definition of ‘noise’ is changing, and that advances in hardware are enabling solutions to long-standing challenges with cross-talk artifacts and irregular illumination. The addition of appropriate a priori information into joint migration and inversion allows historical assumptions about the unknown velocity model having a smooth background to be dismissed, and step changes in velocity model resolution may be achievable. I conclude by discussing how higher resolution velocity models will translate to less non-physical efforts in processing that can corrupt pre-stack amplitude, prestack phase and (elastic) image resolution.
-
-
-
Gravity Gridding in South Australia
More LessAuthors Philip Heath and Laszlo KatonaCreating an image of the acceleration due to gravity over a large area (in this case a large portion of South Australia) is not a trivial process. This paper examines some of the issues involved with creating such an image, and presents some examples. Treating the state gravity as a single dataset highlights outliers and results in ‘dimples’. These are short wavelength features around single points. Treating the data as a compilation of grids and then levelling the grids results in linear artefacts where the survey boundaries meet. Two new approaches have been implemented, involving removing selected data points (based on proximity to adjacent points) and implementing variable density gridding techniques. The resulting grids still have artefacts (notable when viewing a first vertical derivative of the grid), but are smoother and more geologically plausible.
-
-
-
Estimating cover thickness using seismic refraction in the southern Thomson Orogen - An UNCOVER application
More LessAuthors James Goodwin, Tony Meixner, Sarlae McAlpine and Malcolm NicollThe Southern Thomson Project was established to develop a better understanding of the geology and mineral potential of the southern Thomson Orogen. One way in which the Southern Thomson Project is improving this understanding is through the collection of seismic refraction data at 16 greenfields sites to assess the cover thickness (i.e. the amount of regolith and sedimentary basin cover overlying the basement geology). Seismic refraction data were collected using a standard linear array with 48 geophones and a 40 kg propelled weight drop as the energy source. An estimate of the cover thickness was produced from the refraction data using the time-term inversion method. This resulted in the creation of a three-layer model for each site, which accounts for the layers associated with the regolith, sedimentary basin cover and the basement geology.
-
-
-
A Bayesian inference tool for geophysical joint inversions
More LessGeophysical joint inversions seek to exploit the statistical fact that a model that simultaneously satisfies two or more independent data sets is more likely to represent geological ‘reality’ than a model that only satisfies a single data set. Interpreting geophysical data directly rapidly exceeds the capacity of a human as more data are added, so some form of machine assistance is usually required. Conventional inversion techniques can produce a ‘best fit’ model but this might only be one of a large range of possible models that fit the data. Bayesian inference provides a tool to evaluate the relative probability of all possible geological models in a given set, thereby quantifying the amount of information the data is actually providing.
Over 2012-2014, National ICT Australia (NICTA; now Data 61) worked with a number of university, government and industry partners, with support from the Australian Renewable Energy Agency, to build a Bayesian inference software tool for geophysical joint inversions. The tool was initially directed at geothermal energy exploration but is equally applicable to investigating other geological problems. For one geothermal exploration problem, Bayesian inference allowed us to jointly invert gravity, magnetics, magnetotelluric soundings and borehole temperature records to map in three dimensions the probability of encountering granite >270°C beneath the Moomba region of South Australia. The results correlated well with an independent deterministic inversion carried out by Geoscience Australia, but provided a much richer interpretation in probability space.
NICTA released the software tools as open source code on the GitHib platform.
-
-
-
Time Slicing the Cooper Basin
More LessAuthors Witold Seweryn, Dave Cockshell, Peter Hough and Steve FabjancicAn efficient and practical method of generating composite time slices using all available, non-confidential 3D seismic data recorded in the South Australian sector of the Cooper Basin has been developed. Twenty georeferenced tiff images of time slices between 1000 ms and 2900 ms at 100 ms intervals have been prepared and are intended to be made publicly available through the Department of State Development SARIG website. 2D seismic data has also been used in this project in an attempt to ‘fill in’ some of the gaps between 3D datasets. Forty eight 3D and 5076 2D seg-y data files were utilised in this project. All 3D and 2D segy files were scaled to normalise amplitude levels for gridding and imaging.
The ability to time slice 3D volumes of seismic seg-y data is a basic method employed in seismic interpretation. This method is usually applied to a single survey as commonly used software packages are not designed to view and analyse multiple 3D surveys at the same time or such functionality is not yet well implemented. For large scale regional interpretation, such constraint limits interpreting potential, especially in terms of basin wide correlation of main structural features.
Whilst the methodology developed in this project retains the quality and resolution inherent in the variety of surveys used, it provides a cost effective method of generating basin wide time slice images, compared to costly reprocessing required in merging of 3D data sets into one 3D volume.
The time slice data suite prepared complements other value-added regional datasets prepared by the Department of State Development to stimulate petroleum exploration in South Australia.
-
-
-
A robust gradient for long wavelength FWI updates
More LessWe introduce a robust method to produce long wavelength updates in gradient-based Full Waveform Inversion (FWI). The solution introduces dynamic weights in the velocity sensitivity kernel derived from impedance and velocity parameterization of the classical objective function. The new kernel implementation effectively eliminates the migration isochrones produced by the specular reflections, enhances the low wavenumber components in the gradient in heterogeneous media, and is able to deliver velocity updates beyond the penetration depth of diving waves. We use synthetic examples to illustrate how this dynamic weighted FWI gradient successfully recovers the velocity from pre-critical reflections. We also show with dual-sensor streamer data from deep-water Gulf of Mexico how the dynamic weighted FWI gradient can combine both transmitted and reflected energy in a global FWI scheme.
-
-
-
Taming uncertainty in geophysical inversion
More LessAuthors Malcolm Sambridge, Rhys Hawkins and Jan DettmerThe concept of uncertainty in geophysical inversion is often confined to quantification of errors in parameters estimated from some data. A broader definition is to include uncertainty arising from the assumptions made in posing the inverse problem in the first place. These may include assumptions about the physics of the relationship between observations and unknowns, the class of parameterisation assumed for the unknowns, and guesses about the statistical character of random noise contaminating the data. Typically these assumptions are required to arrive at a tractable mathematical problem to solve using geophysical inversion methods. In this paper we outline an inversion approach that allows a broader definition of uncertainty which includes each of these classes of assumption. Including uncertainty in the model parametrization or in the nature of the noise statistics can lead to more realistic inversion results, but not always with increased error bars on the model parameters. For example, relaxing rigid assumptions in the nature of the parametrisation, even in simple problems can result in smaller and more realistic model error estimation. Using the data to decide between different classes of parameterization, physical assumptions or observational noise process is often called ‘model choice’ in statistics, an area that is often overlooked in the geosciences. Over the past 5 years the trans-dimensional inversion approach has increasingly found applications across a variety of inference problems in the geosciences, with new applications appearing regularly.
-
-
-
Airborne IP: Drybones kimberlite VTEM data Cole-Cole inversion
More LessAuthors Domenico Di Massa, Vlad Kaminski and Andrea ViezzoliA VTEM survey was flown over the Drybones kimberlite in 2005, followed by a ZTEM survey in 2009. These data sets were inverted on multiple previous occasions using various 1D, 2D, 3D and plate modelling algorithms. VTEM data showed AIP effects, manifested as negative voltages and otherwise skewed transients. This created artefacts in conventional inversions of VTEM data, which showed some inconsistencies with ZTEM inversions, as well as with the known geology. In 2015 the VTEM data were transferred to Aarhus Geophysics, reprocessed and reinverted using the modified “AarhusINV” code with Cole-Cole modelling. The results are presented in current abstract, they appear to be more interpretable and provide better data fit, than previous inversion attempts.
-
-
-
Potential Field Data Guided Seismic Forward Modelling of Basement Structures: a Case Study from Offshore Nile Delta Basin
More LessAuthors Said Hanafy, Sharaf Eldin Mahmoud and Shastri L NimmagaddaThe offshore Nile Delta region has huge exploration potential, keeping in view the existing commercial oil and gas fields in and around this basin. The majority of the offshore Nile Delta basin is underexplored. Sparsity of geophysical data coverage and ambiguity in the interpretation of existing geophysical data have held back the exploration for many decades. In order to assess and evaluate the exploration potential, it is important to estimate the depth to the basement, the nature of basement and its topography. Identifying the basement reflector based on the current seismic data is challenging, because of limitations of the resolution and record length of the seismic data. Moreover, an estimation of sedimentary column is paramount for further detailed geological and geophysical investigations. The objective of the research is to perform seismic forward modelling of the basement reflector, guided by gravity and magnetic field data, and minimize the ambiguity involved in a single geophysical method of interpretation. The geophysical properties and their contrasts derived from seismic, gravity and magnetic methods of exploration are various criteria used in the integration and modelling process. Seismic stacking velocity data and depth values interpreted from gravity and magnetic field data are used in the modelling process. The uncertainty in modelling the basement character is assessed based on the stacking velocities and seismo-geological cross-sections interpreted from 2D seismic data, instead of cross-sections obtained from unreliable seismically derived depth data. The analysis of errors observed in the seismic data due to younger layer-cake reflections interfering with the basement reflector are analysed and the observed errors in the velocity models are accordingly reconciled based on the potential field data. The potential field data interpretation guiding the seismic forward modelling facilitates the redesign process of seismic data acquisition and processing in the hugely unexplored offshore Nile Delta basin in Egypt.
-
-
-
Gravity gradient data filtering using translation invariant wavelet
More LessAuthors Dailei Zhang, Danian Huang, Junwei Lu and Boyuan ZhuFull tensor gradient (FTG) data is highly useful in hydrocarbon exploration and the detection of some geological targets with small size as its higher detailed information abundance and finer resolution. At the same time, there are some high-frequency Gaussian white noise mixed in the target signal and which has closer frequency range than the conventional gravity data. Thus, one key step before inversion is to remove as much Gaussian white noise as possible and reserve the subtle details. For this pre-processing step, several effective methods have been used, including low-pass filters, least square fitting methods based on Laplace equation and wavelet filtering methods. In this paper, we would utilize the translation invariant wavelet for the reason that it can suppress Gaussian white noise through multi-resolution analysis and at the same time can avoid pseudo-Gibbs phenomenon. The other point different from wavelet method used before is that we applied a mixed threshold constructed according to the curve of both soft threshold and hard threshold. Compared to soft and hard threshold, mixed threshold can keep more details and remove more noise respectively in terms of the energy distribution of signal and noise. Then we process wavelet coefficients with mixed threshold and do inverse transform to recover the data. The results demonstrate that translation invariant wavelet can not only remove the major part of Gaussian white noise, but also reserves high-frequency detailed information of FTG data. Obviously, translation invariant wavelet with mixed thresholding has preferable application effect in filtering FTG data.
-
-
-
Enhancing coal quality estimation through multiple geophysical log analysis
More LessAuthors Binzhong Zhou and Graham O’BrienCoal quality information such as ash content, density, volatile matter and specific energy are important to the coal mining industry for mine planning, design, extraction, beneficiation and utilisation. These parameters are traditionally obtained through laboratory analyses conducted on drill-core samples from exploration drill holes. This process is expensive and time consuming. In this paper, we use a multi-variable data analysis algorithm based on the Radial Basis Function (RBF) neural network methods to estimate coal quality parameters from routinely-acquired multiple geophysical logs such as density, gamma ray and sonic logs. The performance of this RBF-based approach was demonstrated using both self-controlled training data sets and an independent data set from a mine. It was observed that although the density logs play a key role in coal parameter estimation, the use of multiple types of geophysical logs, including logs with different resolutions such as short spaced density log DENB and long spaced density log DENL, improves the estimation accuracy. It is therefore expected that the use of additional geophysical logs such as photoelectric factor (PEF), SIROLOG and PGNAA, which provide data of geochemical constituents, should improve estimates of coal quality parameters.
-
-
-
Integrating core and wireline log datasets - a pathway to permeability from AvO seismic?
More LessAuthors Lahra Lanigan and Jarrod DunneImproved reservoir modelling and simulation is sought by assigning permeability using AvO seismic. The integration of sidewall-core porosity and permeability with wireline logs represents one possible pathway for differentiating the seismic amplitude responses of different reservoir flow-unit facies.
Two sand facies with differing porosity-permeability trends were defined in the reservoir interval of a Paleocene oil field in a shallow offshore clastic setting. Available sidewall-core data from several wells were added to their respective log datasets as discrete, depth-matched points and used to pinpoint the extraction of relevant log values corresponding to each point. The extracted Vp, Vs and density values enabled calculation of AI and Vp/Vs for each sample. Fluid substitutions using the Gassmann equations were then performed on points within the two sand facies.
A statistical comparison of brine and oil porefill cases for each sand facies in this study showed a clear shift due to the hydrocarbon effect. However, only minor differences were observed in an AI-Vp/Vs cross-plot when the two sand facies with common porefill were overlain. Compositional similarity between the sand facies appears the most likely cause for the lack of significant difference in AvO response. Thin carbonates with high AI are also present in the reservoir interval and lateral variation in carbonate volume of as little as 5-10% would overprint the small differences observed between the sand facies. This method should be revisited for sand facies possessing greater compositional differences or where a larger porosity-permeability distinction between facies exists.
-
-
-
Airborne IP detects only fine-grained minerals when compared to conventional IP
More LessAuthors James MacnaeUsing a thin-sheet model, it is possible to predict Cole-Cole parameters of polarizable materials in the near-surface from airborne EM data. With the high frequency of AEM systems, typically more than 100 times higher than ground IP systems, most conventional IP targets will often not show an AIP response. Very fine-grained minerals, around 0.1 mm in average dimension, are however good sources for AIP responses. In 6 examples from Tasmania and NSW comparing AIP with ground IP, 5 have AIP responses that are not coincident with ground responses, but may detect finer grained minerals in the periphery of the alteration system associated with a mineral deposit.
-
-
-
3D aeromagnetic imaging of Iwate volcano, northeast Japan
More LessAuthors Shigeo Okuma, Tadashi Nakatsuka and Ryota KageuraIwate volcano, northeast Japan is an active Quaternary volcano and is comprised of two parts: West-Iwate and East-Iwate. These bodies are underlain by early-middle Pleistocene volcanic rocks. In 1999, fumarolic areas were newly emerged along the ridge between Ubakura and Kurokura Mountains in West-Iwate and Iwate volcano was thought that an eruption was impending in 2000. However, the fumarolic activity has decreased since its peak in July 2001, and the disaster seems to have passed.
In late 2000, a helicopter-borne EM and magnetic survey was conducted over Iwate volcano to better understand the subsurface structure of the volcano related to the ongoing volcanic activity. Recently we have conducted three-dimensional (3D) imaging of Iwate volcano to constrain its subsurface structure. Our model indicates that magnetization highs occupy the main edifice of East-Iwate, which reflects the surface and/or subsurface distribution of basaltic lavas. Meanwhile, magnetization lows are dominant inside the summit caldera of West-Iwate except for a magnetic high over the Onashiro lava flow. Magnetization highs are also distributed on the northern and southern slopes of West-Iwate but local magnetization lows lie on the heads of narrow valleys, corresponding to hydrothermal altered areas. These hydrothermal altered areas are also characterized by resistivity lows observed by the Airborne EM survey.
Although the imaging improved our understanding of the surface and subsurface distribution of volcanic rocks in Iwate volcano, some limitations exist. No information about the magmas which should have intruded during the recent eruptive crisis was obtained by the imaging. The small magnetic contrast between the intruded magmas and their host rocks is the most probable reason.
-
-
-
Large scale magnetotelluric sounding at the periphery of the Songliao Basin, NE China
More LessAuthors Weijun Zhao, Heng Zhu and Qiuhong DingRecently, China Geological Survey (CGS) has launched 3D geological mapping programs from regional to local scales. The project Deep geological survey at the periphery of the Songliao Basin funded by CGS was implemented from 2012 to 2014. Its main goals are to reveal the tectonic framework of the Jarud Basin (JB) as well as to identify the strata distribution of Permian Linxi Formation by integrating new electromagnetic data with existing geophysical and geological data since black mudstones in the Linxi Formation have shown the potential of shale gas. The study area is situated in Jarud Banner and Ar Horqin Banner, Inner Mongolia, China. It is covered dominantly with Cretaceous-Jurassic igneous rocks with exception of the small southeastern part. It is tectonically located in the southern Great Khingan Range, western margin of the Songliao Basin, north of Xar Moron Fault.
Over years of 2012 to 2014, a magnetotelluric survey was carried out at the west of the Songliao Basin. A total of 1559 stations including existing MT data on eight NW and five NE profiles were obtained, covering area that exceeds 10, 000 km2. After dimensionality analysis and static shift removal, the nonlinear conjugate algorithm was used to conduct 2D inversion for TM and TE modes. Inversion results revealed numerous large faults, some of which constitute the boundaries of the Jarud Basin, and modified the tectonic framework. Integrated with well logging and geological data, two Paleozoic fault depressions (Gadasu and Hunnitu) and one Meoszoic depression was discovered. Two Paleozoic sag were inferred to be presence of Linxi Fms. One Mesozoic sag was inferred to be hydrocarbon potential. Attention should be paid to Gadasu sag with area of around 500 km2 since it contains reasonably thick conductive sediments exceeding 4 km in depth which might be black mudstones pertaining to shale gas.
-
-
-
Application of vertical electrical sounding method to identify distribution of hot groundwater around the hotsprings in geothermal prospect area
More LessAuthors Mariyanto, Has Priahadena and Wahyudi W. ParnadiField survey with geoelectric method has been conducted in the Songgoriti geothermal prospect area, Malang, Indonesia. This area is located between Mount Arjuno and Mount Welirang. The aim of this study is to identify the distribution of hot groundwater. Subsurface resistivity data acquisition is done by using Vertical Electrical Sounding (VES) in the four sounding points around the hotsprings with a maximum path length of 160 meters. Data from this measurement is the apparent resistivity that be a response model of the subsurface rock model parameters at each depth. The true resistivity of the subsurface model parameters is determined by inversion modeling. Result of data processing generates a resistivity model of each layer of rock at depth. This study successfully estimate the hot groundwater aquifer layer in the study area. The layer of hot groundwater aquifer is identified by low resistivity in the VES-1 point, VES-2 point and VES-3 point with different depth and thickness. Resistivity of hot groundwater layer is about 19.543.1 ohmmeter with the largest thickness in VES-3 point with the direction of orientation from Mount Welirang to southeast.
-
-
-
Processing of Airborne Gamma-Ray Spectra: Extracting Photopeaks
More LessAuthors Eugene DrukerReceiving information from airborne gamma-ray spectra is based on the ability to estimate the photopeak areas in the regular spectra of natural and other sources. In the airborne gamma-ray spectrometry, extracting the photopeaks of radionuclides from regular one-second spectra is a complicated problem. In the region of higher energies, e.g., above 1.6 MeV, the difficulties are associated with low count rates, while in the region of lower energies, difficulties are due to a significant background level and its statistical noise. In this article a new procedure is proposed to process the measured spectra up to extracting evident photopeaks. The procedure consists of decreasing noises in energy channels along the flight lines, transforming spectra to equal resolution spectra, removing baselines from each spectrum, sharpening details, and transforming spectra back to original channel scale. The resulting spectra are better suited for examining and using the photopeaks. No assumptions regarding the number, positions and magnitudes of photopeaks are needed. Non-negativity of photopeak areas is ensured by the procedure. The proposed technique is likely to contribute to studies of environmental issues, soil characterization and other near surface geophysical methods.
-
-
-
Processing of Airborne Gamma-Ray Spectrometry using Inversions
More LessAuthors Eugene DrukerStandard processing of Airborne Gamma-Ray Spectrometry data generally gives good results when the geological situation is uniform and the conditions of measurements are quite constant within footprint area with possible exception of flight height variations in a small range. Any violation of these conditions leads to certain problems. In reality, violations such as large changes of flight height and/or rugged terrain are not that rare as well as sharp changes in composition of surface rocks. This article proposes an approach where the solutions of inverse problems are used for data processing. The approach is quite natural in the processing of field data measured along the flight lines: it explicitly takes into account one-dimensional models of survey and flight parameters -from topography to sources distribution on the surface. Also, it clearly demonstrates that the inverse problem of Airborne Gamma-Ray Spectrometry data does not have a unique solution. This feature can be used in accordance with the geological problem in hands because various formulations of inverse problems can lead to various geological solutions. The use of the approach is illustrated by several examples given for both flight lines and survey areas.
-
-
-
Magnetotelluric Imaging of a Carbonatite Terrane in the Southeast Mojave Desert, California and Nevada
More LessAuthors Jared R. Peacock, Kevin M. Denton and David A. PonceThe southeast Mojave Desert hosts one of the world’s largest rare earth element (REE) deposits at Mountain Pass, California. Although surface geology has been studied, a full understanding of the carbonatite and associated intrusive suite complex requires subsurface geophysical characterization. In this study, a combination of geophysical methods, including magnetotelluric (MT), magnetics, and gravity are used to create a two-dimensional (2D) geophysical model to a depth of about 10 km. An electrically conductive body is found 2-3 km below and west of the deposit that is associated with a magnetic high that could be connected to a deeper (10 km) conductive body related to possible intrusions or hydrothermal systems. The carbonatite body coincides with a steep magnetic gradient and a bench or terrace in the gravity data that may reflect relative lower-density intrusive rocks. Although carbonatite rocks are typically magnetic, the carbonatite rocks, associated intrusive suite, and host rocks in this area are essentially non-magnetic. Combined geophysical data indicate that the enriched REE deposit may be related to a regional extensive hydrothermal alteration event.
-
Volumes & issues
-
Volume 2019 (2019)
-
Volume 2018 (2018)
-
Volume 2016 (2016)
-
Volume 2015 (2015)
-
Volume 2013 (2013)
-
Volume 2012 (2012)
-
Volume 2010 (2010)
-
Volume 2009 (2009)
-
Volume 2007 (2007)
-
Volume 2006 (2006)
-
Volume 2004 (2004)
-
Volume 2003 (2003)
-
Volume 2001 (2001)
-
Volume 1999 (1999)
-
Volume 1994 (1994)
-
Volume 1987 (1987)
Most Read This Month