ASEG Extended Abstracts - ASEG2012 - 22nd Geophysical Conference, 2012
ASEG2012 - 22nd Geophysical Conference, 2012
- Articles
-
-
-
Time-lapse wave-equation migration velocity analysis
More LessAuthors Jeffrey Shragge and David LumleySummaryTime-lapse analysis of seismic data acquired at different stages of hydrocarbon production or fluid/gas injection has been very successful at capturing detailed reservoir changes (e.g., pressure, saturation, fluid flow). Conventional 4D analysis is performed in the time domain assuming a constant baseline model; however, this procedure becomes difficult when the subsurface is significantly altered by production/injection and large time anomalies and complex 4D coda are recorded. We argue that a more robust 4D analysis procedure in these situations requires iterative wave-equation depth imaging and time-lapse velocity analysis.
Wave-equation depth migration requires accurate knowledge of the velocity field usually obtained by one of two ways. First, data-space methods are where recorded data are matched to those calculated through a background velocity model. Differences between the two datasets are used in tomographic backprojections to generate velocity model updates. Alternatively, imagespace methods are where discrepancies between migrated images (non-flat gathers) are backprojected to estimate velocity model updates. These types of approaches are termed migration velocity analysis (MVA).
This abstract focuses on extending 3D wave-equation MVA (WEMVA) approaches to time-lapse velocity analysis. We discuss the differences between 3D and 4D WEMVA inversion goals, and how we leverage the locality of 4D image perturbations to provide highresolution velocity model updates. We demonstrate the utility of 4D WEMVA analysis in a synthetic CO2 geosequestration experiment by successfully inverting for a velocity perturbation corresponding to a thin layer (<20m) of injected CO2 in a typical North Sea reservoir.
-
-
-
-
Feasibility analysis of drill bit tracking using seismic while drilling technique
More LessAuthors Baichun Sun, Binzhong Zhou, Andrej Bona and Roman PevznerSummaryCheck-shot survey measures the first arrival time with a known depth receiver in borehole to assess formation velocity. This information can be used in correlation with sonic log and surface seismic products for adjustment of interpretation. Check-shot survey can also be implemented with seismic-while-drilling using drill bit noise as the source. This differs from usual check-shot survey as source is in the borehole. It provides a real time, cost saving, and safe measurement.
Check-shot survey needs a known receiver depth, thus velocity can be obtained by fixed wave travel path and the measured first arrival time. However, in seismicwhile-drilling (SWD), drill bit position can vary a lot from vertical drilling to deviated drilling. To address this issue, we present a method that finds the location of the source and estimates the velocity of the formation at the same time. Using a synthetic model, with medium receiver offsets, this method shows good estimation of the drill bit depth location and formation velocity in a layered Earth model.
-
-
-
Architecture and evolution of the West Musgrave Province, and implications for mineral prospectivity
More LessAuthors Aurore Joly, Alan Aitken, Mike Dentith, TC McCuaig, Alok Porwal, Hugh Smithies, Ian Tyler and Shane EvansSummaryThe West Musgrave Province preserves a geological history spanning much of the Proterozoic, including two major Grenville-aged orogenic events, ca. 1345-1120 Ma. These were followed by the intraplate Giles Event, ca. 1080-1050 Ma, which is characterised by 1- the voluminous mafic intrusions of the Giles Complex and 2-the deposition of a thick sequence of bimodal volcanic rocks and sedimentary rocks (the Bentley Supergroup). The architecture resulting from these events was subsequently overprinted by later events, the ca. 600-520 Ma Petermann Orogeny, and the ca. 450-350 Ma Alice Springs Orogeny. Here we use aeromagnetic, gravity and magnetotelluric data, constrained by geological mapping and petrophysical data, to characterise the 3D architecture of the region, and to unscramble its structural evolution. Early deformation events are well preserved in places, although they do not permit a robust interpretation of the regional architecture at that time. The architecture of the Giles Event is extremely well preserved, and a complex polyphase history is conserved. A series of discrete to semi-continuous deformation events is recorded, with a dominant ESE-WNW trend supplemented by additional NE-SW and N-S trends, each of which was active at several times during the event. Later events are manifested primarily as reactivations of earlier structures. Using this architectural framework along with a variety of other spatial datasets (geology, geochemistry etc) GIS-based prospectivity analysis was applied using a mineral systems approach, targeting Ni, Cu, PGEs and orogenic gold. The conceptual method used compares spatial distributions of various targeting criteria, represented by predictor maps. Each predictor map is then related to metal source, fluid pathways, and chemical and physical trap zones. The output prospectivity maps are used as decision-support tools for exploration and reinforce the idea that the area is highly prospective for magmatic nickel-copper and PGE’s deposits.
-
-
-
Application of Airborne EM to Bowen Basin Coal Projects
More LessAuthors Kate E. Godber, James Reid and Guy LeBlanc SmithSummaryIn 2011, several Bowen Basin coal exploration projects in central Queensland were surveyed with the SkyTEM airborne electromagnetic system. The aim of these surveys was to trial the effectiveness of the method in for coal applications, with specific focus on mapping the key aspects of the weathering, overburden and hydrology that are most important for effective development of a coal project.
The first survey, which was over an area near Moranbah with 150m of Tertiary overburden with extensive basalt flows, mapped crucial hydrological and geological features important to efficient further development of the project: These included a) a convoluted network of basalt-filled palaeochannels, some of which extend through and into the coal measures b) location, extent and thickness of major saline aquifers, c) depth to base of Tertiary.
Subsequent surveys over different parts of the Bowen basin, showed that SkyTEM was a very effective tool for mapping the thickness of weathering, a parameter very important for drilling out coal resources. The results are surprising in that they not only provide a detailed map of depth of weathering (including volumes of lateritic sediments requiring excavation), but also map more subtle geological features such as synclines (previously only loosely inferred from surface mapping), major faults, previously unknown dykes, dip and strike of the main stratigraphic units, fresh water aquifers and even possibly also various coal seams (generally because they are saline aquifers).
These surveys represent, at least to the authors’ knowledge, the first example of modern helicopter-borne transient EM in the Australian coal fields. The results from this project and subsequent projects are surprising in that they have had unexpected applications above and beyond those which were originally forecast. Thus, with informed and thorough use, airborne EM promises to deliver valuable geological and geotechnical information to aid in the further development of Australia’s coal resources.
-
-
-
Telfer Region FALCON Airborne Survey, Western Australia
More LessAuthors Campbell MackeySummaryThe Telfer region in the Paterson Province consists of Proterozoic metasediments of the Yeneena Group, of lower-upper Greenschist facies. Sediment lithologies vary from siltstones to quartzites and carbonates. These are intruded by intrusive suites ranging from felsic to mafic, plutonic to hyperbyssal.
Known world-class assets include the Telfer Gold operation, and the O’Callaghan’s Tungsten resource (both 100% Newcrest owned).
The sedimentary and intrusive lithologies have considerable density contrast. These factors, combined with steeply dipping on-edge stratigraphy make the Telfer region ideal for gravity surveying.
Previous gravity coverage in the region has been a combination of small prospect scale surveys, and 1980s vintage 2 x 4 km regional helicopter surveys with barometric levelling control.
In order to get systematic gravity coverage of good resolution, a FALCON survey of 250m line spacing, 60m flight height, was flown in 2010.
The data show a number of blind felsic plutons, some intersected in drilling, some interpreted – especially associated with the O’Callaghan’s tungsten resource, and the Trotman’s Stockwork tungsten prospect. The data also outline much more clearly the NE over SW thrust faults in the region, and the late north-south transfer faults, the most well known of which is the Graben Fault through Telfer Dome.
Three dimensional gravity inversions performed over the O’Callaghan’s and Trotman’s Stockwork tungsten areas in Gocad/VPMG software have defined extra felsic intrusive targets for tungsten and gold mineralisation. Further diamond drilling will increase inversion constraints.
-
-
-
Geophysics at the Hawsons Iron Project, NSW – Eastern Australia’s new magnetite resource
More LessAuthors John Donohue, Quentin Hill and Doug BrewsterSummaryA JORC inferred resource of 1.4 billion tonnes at a magnetite Davis Tube Recovered Grade of 15.5% was defined at The Hawsons Iron Project near Broken Hill in western New South Wales making it the largest magnetite iron ore resource identified to date in eastern Australia.
Reconnaissance surface exploration in 2009 identified the area’s potential for magnetite iron ore fifty years after it was last considered an iron ore play prior to the discovery of Mt Tom Price.
Neoproterozoic Braemar Ironstone Facies sediments host the magnetite mineralisation which generates aeromagnetic anomaly amplitudes greater than 6000 nT. Five target areas for magnetite iron ore were identified from the open-file aeromagnetics.
A realistic Exploration Target was estimated from selective 2.5D magnetic modelling and image processing of the airborne aeromagnetic data. The inferred resource was established within 18 months of commencing the exploration program drilling 72 holes into the five aeromagnetic targets.
The monotonous sedimentary host sequence of siltstones and sandstones with few marker beds makes detailed geological interpretation extremely difficult. A magnetic stratigraphy was constructed from recognizable patterns in wireline magnetic susceptibility traces and correlations of the mineralisation could then be established down dip and across strike.
-
-
-
Structurally constrained lithology characterization using magnetic and gravity gradient data over an iron ore formation
More LessAuthors Cericia Martinez, Yaoguo Li, Richard Krahenbuhl and Marco Antonio BragaSummaryIt is often desirable to extract meaningful lithologic information directly from geophysical data. Here, we describe a multi-faceted approach that combines the known geologic sections of a site from borehole data with recovered density and magnetic susceptibility distributions from 3D inversion of airborne gravity gradient and magnetic data. The technique can produce end-member solutions based on average property values extracted from literature and well records; or, it can be implemented using the geologic sections to extract an appropriate range of property values to address the model smoothing common to modern 3D generalized inversions. To demonstrate our approach, we present results utilizing magnetic and gravity gradient data over part of the Gandarela Syncline iron formation in the Quadrilátero Ferrífero, Brazil.
-
-
-
An iterative approach to optimising depth to magnetic source using the spectral method
More LessAuthors Anthony (Tony) Meixner and Stephen JohnstonSummaryKnowledge of the depth of cover is poor across large areas of Australia. The spectral method is an efficient method of producing reliable depth to magnetic basement estimates across large regions of the continent.
A semi-automated work-flow has been created that enables the generation of depth to magnetic source estimates from windowed magnetic data using the Spector and Grant method. The work-flow allows for the correction of the power spectra, prior to the picking of straight-line segments, to account for the fractal distribution of magnetic sources. The fractal parameter (β) varies with depth and was determined by picking multiple depth estimates in regions of outcropping magnetic basement which have been upward continued to different levels in order to simulate different amounts of burial beneath non-magnetic sediments. A power law function best approximates the decay of β with depth.
An iterative schema used to determine the optimum β where the depths of magnetic sources are unknown, has been incorporated into the workflow. Preliminary testing in a region of known magnetic basement depth has produced encouraging results, although further testing is required. The decrease of β with increasing depth suggests that the fractal distribution of magnetisation becomes less correlated, or fractal, over larger volumes of observation.
-
-
-
From two-way time to depth and pressure for interpretation of seismic geometries and velocities offshore
More LessAuthors Alexey GoncharovSummaryThe effect of water and rock loading on seismic velocities and consequently on interpreted geometries is often underestimated in offshore studies. Direct comparative analysis of interval velocity patterns between areas of significantly different water depth and thickness of rock overburden requires various pressure related changes in velocity to be accounted for. Presentation of velocity models as a function of pressure rather than two-way time, or depth, emerges as a possible solution. An accurate velocity model is essential for meaningful timeto-depth conversion of interpreted seismic horizons. Ideally, it should be based on integration of seismic velocities from well log measurements, refraction seismic surveys and from stacking of multi-channel marine reflection data. In some cases velocities derived from stacking of high quality long streamer marine reflection seismic data correlate reasonably with well log measurements and velocities derived from refraction seismic studies, and provide clues to reasonable depth conversion and lithology interpretation.
-
-
-
Quantitative regularity analysis of offset-vector sampling for seismic acquisition geometry
More LessSummarySymmetric acquisition geometry consisting of identical sampling of shots and receivers, can maintain the spatial continuity of the wavefield automatically, according to symmetric sampling theory. However, asymmetric geometry is often adopted in practical seismic exploration applications. Such geometry can cause uneven sampling and is necessary to be assessed for its sampling performance prior to acquisition. In conventional survey design, based on the common mid-point (CMP) analysis for a horizontally layered earth or common reflection point (CRP) analysis for a complex subsurface structure, the quality of acquisition geometry is generally judged by such bin properties as effective fold, offset scalar and azimuth distributions. However, these conventional approaches are limited by an incomplete understanding of the offset-vector sampling. Therefore, we propose a new method for quantitatively evaluating the continuity of offset-vector sampling including four spatial coordinates of shot and receiver. On the basis of physical potential energy and force-balance principle, it analyzes the regularity coefficient of offset-vector sampling as a whole using potential function model and takes into account fold, offset-scalar and azimuth distribution factors. The combination of regularity coefficients of every bin can produce spatial continuity distribution of offset-vector sampling. Similar to symmetric sampling, this approach emphasize the spatial relationships between adjacent bins rather than single bin attribute, since it aims to maintain the spatial continuity of the wavefield which allows the faithful reconstruction of the underlying continuous wavefield. Using this method, we can quantitatively compare spatial continuity distribution for different seismic acquisition geometries, and then choose the better acquisition scheme.
-
-
-
Investigation of the weathering layer using seismic refraction and highresolution seismic reflection methods, NE of Riyadh city
More LessAuthors Ghunaim T. Al-Anezi, Abdullah M. Al-Amri and Haider ZamanSummaryFive seismic refraction and five high-resolution seismic reflection (HRSR) profiles were carried out in northeastern part of Riyadh city to investigate depth of the weathering layer. Results obtained from seismic refraction survey reveal the depths of weathering layer at 12, 25, 17, 12, and 16 m, respectively. On the other hand, HRSR stack sections illustrate the depths of weathering layer at 14, 28, 20, 13, and 18 m, respectively. The weathering layer is composed of alluvial sediments and gravel, which is underlain by a sequence of limestone and dolomite layer. Seismic results from site no. 2 have been found to be in good agreement with lithological information reported from the adjacent water well. The HRSR data generally reveal better signal-to-noise ratio and enhanced resolution compared to the refraction data. Although, the HRSR data failed in achieving high-quality common midpoint (CMP) stacking profile at site no. 3, it provide an improved image of the subsurface features than the refraction data, recognizing it as a potential seismic technique.
-
-
-
Case study on the application of geophysics to a port expansion project
More LessAuthors Ian James and Greg TurnerSummaryThe first stage of a major port expansion project called for an approach to rapidly and cost effectively screen a broad area to assist with determining the most appropriate location for additional berths, turning pockets and associated channels. A combination of marine geophysical techniques were used to achieve this goal without the need for drilling which, although essential for detailed assessment and planning of the final location, would have resulted in significant delays and extra costs at this stage.
Seismic refraction and continuous seismic profiling (CSP) were used together to characterise the upper 20 m below the seafloor. Sidescan sonar, swath bathymetry and magnetometry were used to check for possible seafloor hazards.
The results showed a region with a deeper sequence of low velocity material and a high seismic velocity ridge. Since the high velocity materials represent a significant impediment to dredging the surveys highlighted the areas more and less suitable for use in the expansion.
In these surveys no unknown seafloor hazards were identified.
The results at the site demonstrate the applicability of marine geophysical surveys to efficiently provide an indication of subsurface conditions over a broad area. This ability makes them a valuable tool for port development investigations.
-
-
-
Mapping of Bedrock Using the High-Resolution Seismic Reflection Technique at Wadi Al Dawasir Region, Saudi Arabia
More LessAuthors Ghunaim T. Al-Anezi, Majed AlMalki and Tariq AlkhalifaSummaryThis paper critically evaluates the utility of seismic data to assist in the interpretation of bedrock. Four highresolution seismic reflection profiles were carried out to provide estimates of the depth to the bedrock and to detect any geological faults present in the area, which may affect the hydro-geological system at Wadi Al Dawasir region, 690 km south of Riyadh city. The bedrock reflection is clearly present in all sites. The depth of the bedrock at site 1 was approximately 750 m, 800 m at site 2, 700 m at site 3, and 950 m at site 4. The bedrock depth obtained from the seismic agreed with the well information. We developed a clear understanding of the bedrock in the study area and mapped its depth as well as mapped some of the clear fault locations.
-
-
-
Full 3D Acquisition and Modelling with the Quantec 3D System - The Hidden Hill Deposit Case Study
More LessAuthors Mehran Gharibi, Kevin Killin, Darcy McGill, William B. Henderson and Trent RetallickSummaryThis abstract will detail Quantec’s 3D system, a proprietary 3D DC/IP and MT acquisition developed by Quantec Geoscience Ltd., Canada by means of a case study over the Hidden Hill Deposit in central Nevada, USA.
An integrated data acquisition, modelling and interpretation of the DC, IP, and MT surveys over a complex geological setting embedding a known epithermal gold and silver-bearing deposit are presented. The survey and interpretation are based on a true 3D acquisition and inversion of the geophysical data.
The results show a very close correlation between the geophysical resistivity and chargeability signatures and the location of the known mineralization. The extremely large data volume collected in an omni-directional current source proved to be essential for a detailed and reliable imaging of the subsurface.
-
-
-
An investigation of statics correction methods for 3D PS-wave seismic reflection
More LessAuthors Shaun Strong and Steve HearnSummaryIn the last decade converted-wave (PS-wave) seismic reflection studies have successfully demonstrated that a more complete geological interpretation can be obtained by integrated interpretation of P-wave and S-wave information, at both the petroleum and coal scales. Full 3D implementation of PS reflection presents particular challenges at the coal scale, because relative offsets and dominant frequencies are both large compared to petroleum-scale reflection.
One of the most difficult steps in the PS processing sequence is estimation of the S-wave receiver statics. Statics are time delays caused by variations in the weathering layer, and changes in source and receiver elevation. These time errors can significantly degrade CMP stack quality, and the final geological interpretation of seismic images. Static errors tend to be much more significant in PS surveys since the S wave travels more slowly through the weathering layer and is therefore more likely to be affected by variations within this layer.
In this presentation we evaluate of a number of different approaches for estimating 3D PS statics solutions. These include a surface-consistent inversion algorithm (analogous to the residual-statics method used in conventional P-wave processing), a so-called 'robuststatistical' method, and PPS refraction analysis. The methods are evaluated using a coal-scale 3D-3C survey acquired in the Bowen Basin. This has been used to examine the comparative performance, and the influence of various algorithmic and geological factors.
The results indicate that the surface-consistent inversion method can fail under some weathering conditions. When this occurs the refraction based method or a robust statistical method are preferred.
-
-
-
Investigation of azimuthal anisotropy in high-fold 3D multicomponent seismic reflection
More LessAuthors Steve Hearn and Shaun StrongSummary3D P-wave seismic surveys can exhibit significant azimuthal variation in stacking velocities, and failure to allow for such variations can introduce smearing into the stacked volume. The problem is likely to be worse in the case of converted-wave (PS) reflection. Firstly, S-waves have lower velocities such that time variations are amplified. Secondly, PS rays are asymmetrical, such that anomalous features may be traversed by different wavetypes (P or S), depending on the direction of travel.
Recently, a coal-scale 3D-PS trial was recorded in the Bowen Basin with the aim of providing a detailed investigation of such azimuthal variation. The survey was designed with extremely high fold (>500). This allowed good quality images to be constructed for data subsets having restricted ray azimuths. Target structures interpreted using different ray-azimuths exhibit significant timing variations (up to 30ms). The observed azimuthal variations may not necessarily indicate true azimuthal anisotropy. They can result from poor statics solutions, and PS imagery is notorious for difficult statics.
On the other hand, the observed variations may be indicative of true geological anisotropy. Based on shearwave splitting models, our PS velocity variations have been modelled in terms of elliptical variation with azimuth, and this approach predicts the orientation of the horizontal stress field. It is interesting that a majority of our data zones indicate a consistent stress orientation. Furthermore, the interpreted horizontal-stress orientation is consistent with observed reverse faulting in the area.
-
-
-
Assessing the presence of hard rock along a gas pipeline alignment with airborne EM
More LessAuthors Niels B. Christensen and James E. ReidSummaryOver the past decades, airborne electromagnetic (AEM) surveys have mostly been used in connection with mineral exploration and a variety of issues in hydrogeophysical mapping. However, increasingly, AEM is used for a wide range of geotechnical purposes: pollution mapping, geotechnical assessment on road and freeway alignments, bathymetry, and depth to bedrock.
We present an investigation using a helicopterborne transient electromagnetic system along the planned trace of a gas pipeline. Oil and gas pipelines are often buried at a depth of a few meters and the cost of construction depends critically on whether the subsurface is composed of soft sediments that can be easily excavated or hard rock formations that require much heavier equipment and possibly have to be blasted. The aim of the AEM survey was to distinguish between the soft, relatively conductive sediments and the hard, relatively resistive bedrock in the upper few meters of the subsurface.
Data were collected with a rather small transmitter moment, but a high repetition frequency that simultaneously allowed high acquisition speed, and reliable data quality. Measurements were inverted with 1D models with both vertical and lateral constraints to produce model sections along flight lines. A novel method of statistical analysis of the set of equivalent models for each inverted model, calibrated against boreholes, improved the estimates of the presence of hard rock along the flight lines.
-
-
-
Inversion methodology for determination of near-surface geology from seismic refraction amplitudes
More LessAuthors Alan MeulenbroekSummaryThe amplitude of a seismic refraction event is determined by the properties of rocks through which the seismic waves travel, the amplitude of the shot and the offset at which the refraction is recorded. A surface-consistent, non-linear inversion scheme, which uses the Levenberg- Marquardt algorithm, has been developed which aims to extract near-surface rock properties from the measured refraction amplitudes.
Perhaps the most important challenge to extracting a unique, geologically plausible solution is in determining initial values of control parameters which dictate how the solution is allowed to progress. If these control parameters are not tailored specifically to the problem in question, convergence to a solution can be very slow and even fail to get off the ground in some cases.
Comparison between results obtained using default control parameters, compared to those obtained using control parameters specifically tailored to the problem, shows a reduction in the error between the true observations and the model-generated observations, while retaining fidelity in the solution. A reformulation of the problem significantly reduces the error but fidelity is apparently compromised. Comments are made on how to achieve an optimal middle ground.
-
-
-
Interpretation of magnetic gradient tensor for automatic locating a dipole source
More LessAuthors Hyoungrea Rim, Young-Sue Park and Hyen Key JungSummaryIn this paper, I propose the algorithm that the location of a magnetic dipole can be detected from the magnetic gradient tensor. I derive the location vector of a vertically magnetizated dipole from magnetic gradient tensor. Deficit of magnetic moment of magnetic dipole makes the induced location information incomplete. However if the observation of magnetic gradient tensor would be collected on one more points, the algorithm is able to point the location of magnetic dipole by clustering the solution of the proposed method. For example, I show that magnetic gradient tensor can be converted as the source location successively by picking common solution area in synthetic case of borehole observation.
-
-
-
“Unsmooth” 1D inversion of frequency domain marine controlled source EM data
More LessAuthors Neil Godber and Peter FullagarSummaryThe goal of geophysical inversion of electromagnetic (EM) data is to recover a model of the geoelectrical properties of the sub-surface. The standard practice for 1D inversion of marine controlled source EM (CSEM) data is to generate smooth conductivity models using least squares (L2-norm) methods. However, sedimentary geology is stratified and piece-wise continuous. As such, smooth resistivity models cannot represent this character. In response to this inconsistency, a means was sought to generate more geologically plausible, piece-wise continuous models.
The common approach in the literature when generating piece-wise continuous inversion models is to regularize L2-norm methods in such a manner as to induce blocky behaviour. Although effective, these techniques are selfconflictive; forcing non-smooth behaviour from an implicitly smooth algorithm. In contrast, L1-norm inversion inherently produces piece-wise continuous models. To investigate the possible utility of this approach, a L1-norm inversion algorithm has been developed and tested on synthetic and real datasets. The L1-norm results were compared with those generated using an industry standard L2-norm algorithm.
The synthetic inversions focused on previously published examples. The real data inversions focused on electric and magnetic field measurements recorded over the main reservoir sand of the Pluto gas field in block WA350-P, North West Shelf, WA.
The L1-norm inversions recovered, to within the resolution limits of the CSEM method, the depth, thickness and resistivity of the synthetic geological models and the Pluto-1 resistivity well log, whilst fitting the input data to within noise. When compared against the L2-norm profiles, the L1-norm inversion more closely represented the stratified character of the sedimentary sequence. It was therefore concluded that L1-norm inversion is an attractive alternative to smooth L2-norm methods when blocky inversion models are desired.
-
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