ASEG Extended Abstracts - ASEG2004 - 17th Geophysical Conference, 2004
ASEG2004 - 17th Geophysical Conference, 2004
- Articles
-
-
-
Practical, 3D surface-related multiple prediction (SMP)
More LessAuthors Ian MoorePrediction of surface multiples via 2D algorithms (eg, SRME) is now routine in data processing, and is often effective in removing those multiples when combined with an adaptive subtraction. There are, however, many situations in which the 2D assumptions made on the geology and the acquisition geometry are invalid to the extent that the resultant errors in the predicted multiples are too large to allow those multiples to be effectively subtracted.
3D algorithms have the potential to overcome these problems, and are therefore very attractive. However, there are many issues to be overcome when implementing a 3D algorithm for use on conventional data. These include inadequacies in the crossline sampling, the limited maximum crossline offset (aperture) of the recorded data, and irregularities in the acquisition geometry. This paper attempts to address these issues through the development of a methodology that predicts the required crossline aperture and the timing errors associated with a given acquisition geometry for a given mode of multiple. Knowledge of the aperture and timing errors is very useful in determining where a 3D prediction is appropriate, and in determining optimum parameters for that prediction.
The results demonstrate the accuracy and value of the predicted errors, and the benefits of a 3D prediction over a 2D prediction for conventional marine data. 3D surface-related multiple prediction is practical for such datasets. The error analysis can also be used to optimise the acquisition geometry and hence the level of success of both 2D and 3D algorithms.
-
-
-
-
Hybrid Velocity Model Prestack Imaging
More LessAuthors Bjørn Müller and Matt LamontDeep-water marine seismic data often suffer from non-hyperbolic moveout distortions generated by highly variable seafloor topography.
The non-hyperbolic moveout is a result of lateral velocity variations across the rugose/dipping water bottom. Time migration assumes no lateral velocity variations and hence is unable to deal with these variations. The result of using a time migration algorithm is both a substandard image and perhaps more importantly, significant lateral movement of events. To overcome this limitation a hybrid depth migration coupled with a time velocity model building procedure is proposed. With this procedure the non-hyperbolic moveout distortions caused by highly variable seafloors are accounted for. The only cost is in the added expense of running a Pre-Stack Depth algorithm rather than a Prestack Time algorithm. The model building effort is comparable.
-
-
-
Some observations on the sedimentary framework of the Loxton-Bookpurnong region, South Australia as defined by borehole, ground and airborne geophysical data - implications for informing the development of groundwater interception schemes
More LessAuthors Tim Munday, Tony Hill, Ben Hopkins, Tania Wilson, Andy Green and Andy TelferThe construction of salt interception schemes (SIS) in the Riverland region of South Australia forms an integral part of a broader strategy to manage saline groundwater intrusion into the River Murray. Results from the inversion of airborne electromagnetic data provided some insight into the distribution and variability of the Loxton sands aquifer. These data indicated regional facies variations associated with the main barrier systems of this prograded strandline sedimentary sequence. The relevance of this information is being followed up in an examination of options for the Bookpurnong highland borefield. A program involving the acquisition and interpretation of neutron, gamma and inductive conductivity borehole logs, NanoTEM ground TEM traverses, and the analysis of HEM data was undertaken to better define the sedimentary and hydrogeological model of the area. This approach has been critical in explaining local scale changes in the sedimentary environment and elucidating reasons for the variable aquifer yield in the Loxton sands aquifer.
-
-
-
3D Full Tensor Gradiometry: a high resolution gravity measuring instrument resolving ambiguous geological interpretations
More LessAuthors Colm A. Murphy and Gary R. MumawThe high-resolution, high-precision gravity gradient measuring technology, Bell Geospace’s 3D FTG (Full Tensor Gradiometry), is presented as an exploration tool to enhance and solve ambiguous geological interpretations from conventional methods. Its high frequency character produces a Total Gravity Field that identifies subtle density contrasts within section, allowing it to be used in detailed hydrocarbon and mineral exploration projects, both offshore and onshore.
This paper presents two case histories describing the results from an airborne and a marine survey. Both studies demonstrate the technology’s ability not only in determining target shape but also in mapping target prospects across large sedimentary basins. The first example presented in this paper images a salt dome onshore Louisiana, USA, making inferences about its emplacement through direct mapping of controlling structure. The second identifies and maps a series of low-density sedimentary deposits on the flanks of the Judd Basin offshore NW Europe.
The implications for exploration initiatives are significant as FTG data reduce risk in geological interpretation, thus facilitating a rapid decision making process.
-
-
-
The Waterloo and Amorac nickel deposits, Western Australia: A geophysical case history
More LessAuthors Paul Mutton and Bill PetersThe Waterloo and Amorac nickel sulphide deposits are significant discoveries located in the north eastern goldfields region of Western Australia. They were found as a result of the systematic application of geophysical techniques during a green fields exploration program. The deposits are 300m from each other, blind, and are overlain by 80m of deeply weathered basement rocks. Despite their close proximity, the deposits are quite separate and different in style.
Interpretation of detailed aeromagnetic data was first used to delineate komatiite-basalt contacts of interest for nickel sulphide exploration. Then these contacts were surveyed with a moving loop TEM survey (Inloop and Slingram receivers) which delineated nine basement conductors. Finally, fixed loop TEM surveys were used to better define the targets for drill testing.
Drill hole TEM surveys have been used since the discovery to guide drilling and better define the mineralisation size, location, orientation, and character. Complicated geometry and multiple conductive sources make careful interpretation and modelling essential.
-
-
-
The in situ stress field of the west tuna area, gippsland basin: implications for natural fracture-enhanced permeability and wellbore stability
More LessAuthors Emma J. Nelson, Richard R. Hillis, Scott D. Mildren and Jeremy J. MeyerThe in situ stress field and natural fracture occurrence in the West Tuna area of the Gippsland Basin were evaluated in order to assess the potential for natural fracture-enhanced permeability in the deep intra-Latrobe group and Golden Beach Subgroup reservoirs, and to investigate wellbore stability issues in the area.
Borehole breakout and drilling-induced tensile fractures (DITFs) interpreted on six image logs from the West Tuna area constrain the maximum horizontal stress orientation to ~138°N. Leak-off test data suggest the upper bound to the minimum horizontal stress is ~20 MPa/km. The vertical stress was derived from density and sonic log data and ranges from 20 MPa/km at 1km to 22 MPa/km at 3km depth. The maximum horizontal stress magnitude was constrained to -40 MPa/km using occurrence of DITFs. The in situ stress regime in the West Tuna area is therefore interpreted to lie on the boundary of strike-slip and reverse (σHmax >σv σv ≈ σhmin).
Natural fractures and wellbore failure (breakout and DITFs) were observed to form preferentially in the cemented sandstone units. Finite element methods were utilised to investigate the far-field and near-wellbore stress distribution between horizontal, interbedded sands and shales. Preliminary modelling indicates that a higher Poisson’s ratio for the shale drives it towards a more isotropic far-field stress state. This decreases the propensity for wellbore failure in the shale layers.
Fracture susceptibility analysis of interpreted fracture sets in the sandstone units suggests that electrically conductive fractures are also optimally oriented to be hydraulically conductive in the far-field. Fractures in the shales are slightly less likely to be open and hydraulically conductive in the far-field due to the transition to a more isotropic in situ stress regime.
-
-
-
Integrated use of Seismic, Ground and Airborne Gravity/Gravity Gradiometer, and Ground Geological Mapping Methods in the Eastern Papuan Basin, PNG
More LessAuthors Andrew Nelson, David Holland, O’karo Yogi, Rodrigo Heidorn and Desmond LeechInterOil has complemented the use of seismic data in its PNG petroleum exploration licences with both ground and airborne gravity/gravity-gradient and aeromagnetic data. Potential field methods may be used to significantly reduce exploration costs in jungle-covered areas due to the relatively high effort and helicopter support required to deploy seismic equipment.
Airborne gravity gradient and aeromagnetic data are used to interpolate and extrapolate reprocessed seismic data in InterOil’s PPL 236. They confirm and extend fault correlations on widely spaced seismic lines.
Seismic data acquired in PPLs 237 and 238 show that faults do not, as suspected by some, sole out at shallow depth (eg 1 to 2 km). Instead, steeply dipping faults continue to several kilometers Surface geological mapping, together with Landsat TM and digital terrain model data, greatly assist in planning, acquisition and interpretation of seismic data. Ground gravity data acquired along seismic lines shows good character which can be correlated with the seismic data interpretation.
Relatively low effort (shallow holes, small charges), optimally located seismic data acquisition is shown to be effective in substantially resolving ambiguous structure and extending interpretation from surface dip and lithology data. Potential field data, once calibrated with seismic data, is demonstrated as a method capable of substantially extending the seismic interpretation.
-
-
-
Enhancing electrical signals with sensor arrays
More LessAuthors M. Norvill and A. KepicIncreasingly, electrical surveys employ multiple sensors to make data collection more efficient. The use of multiple sensors can also be used to improve signal fidelity from each sensor, leading to more accurate geological models and greater depth of investigation. New procedures to remove noise in data collected using an array of sensors have been developed and tested. These algorithms attempt to remove two types of noise: harmonic noise associated with the mains power grid and atmospheric transients due to distant lightning strikes and electrical storms. Field tests demonstrate improvements in signal-to-noise by two orders of magnitude. The algorithms are applicable to any type of electrical or EM survey conducted with a sensor array.
-
-
-
Implementation of volume interpretation in revealing upside potential in a mature field, the sangatta oilfield: A case study
More LessAuthors Susanto B. Nugroho, Bambang S. Murti and Budianto M TohaThe Sangatta field is a shallow, medium sized oilfield, located in prolific Kutai Basin and was discovered in 1939. Field production begun at 1972. Primary production derives from Mid Miocene fluviodeltaic reservoir. Three deep exploratory wells drilled on crestal position in mid 70’ s failed to encounter reservoir quality rock and paleontologically exhibited deeper marine environment. Following conventional structural interpretation based on 3D shot in 1993, a comprehensive 3D volume interpretation was conducted recently. This new method facilitating combination between conventional 3D interpretation with the power of visualization, which enable interpreting of simultaneous multiple seismic attributes as single object that improve degree of confidence. This method eventually has revealed subtle seismic signature that can be interpreted as volumetrically significant geologic body that possess reservoir quality and becoming an attractive exploration target. Further to this study, 4 deep exploratory wells are scheduled to be drilled during this year drilling campaign.
-
-
-
Influence of capillary fringe on the groundwater survey using ground-penetrating radar
More LessAuthors Kyosuke Onishi, Shuichi Rokugawa, Yoshibumi Katoh and Tomochika TokunagaGround-penetrating radar (GPR) is well known as one of the best tools to detect groundwater surface. However, very few people take into account the effect of a capillary fringe over water table with groundwater survey. Water contents vary from 100 percent to several percent in this capillary fringe. The behavior of contaminants is highly influenced by this transition zone. This also has an impact on the design of detailed exploration such as a time-lapse survey. In this study, we took a stepwise approach to examine the influence of the capillary fringe. First, the response of electromagnetic waves reflected from water table was examined in a small scale model in our laboratory. The result indicates that multiple reflections are sometimes detected around the water table and the capillary fringe. It also proved that the wave event which is considered to be reflected on the water table was actually emerged from the transition zone of water content inside the capillary fringe. The waveform of this reflection is not sharp compared with other reflections. This is because the dielectric constant gradually increases in the capillary fringe. In the second step, a field survey was applied to the edge of Kurobe alluvial fan along the seashore. Our instrument having real-time kinematic GPS system enabled us the simple and rapid survey in the wild field with complicated surface roughness. Observed reflections are unshapely detected as the similar shape shown in the scale model experiment. Finally it is concluded that the influence of the capillary fringe should be taken into account the analysis of the detailed groundwater survey.
-
-
-
Seismic delineation of near-field exploration opportunities in the Papuan thrust belt: examples from the SE Gobe Area
More LessAuthors M. ParishHydrocarbon exploration within the Papuan thrust belt has traditionally relied on non-seismic techniques to resolve the complex subsurface structure. Subsequent appraisal drilling for most PNG discoveries has defined a complex pattern of fault-bounded compartments which show differing fluid contacts and pressure isolation.
The Gobe Anticline forms a frontal hangingwall fold at the leading edge of the Papuan thrust belt, with the SE Gobe field forming the easternmost compartment. Experimental seismic data were first acquired over the Gobe area in 1996 with the recording of dip line PN96-204 and in 1998 by acquisition of the SE Gobe experimental dip line PN98-206. This resulted in the successful appraisal of the SE Gobe 7 well, which targeted the southern extent and cutoff of the Iagifu reservoir at depth. Following this success, further seismic acquisition and reprocessing were undertaken during 1999 and 2001 over selected areas of the Papuan thrust belt resulting in over 130 km of 2D seismic data along strike of the SE Gobe trend to define near- field exploration opportunities.
This high-quality seismic dataset significantly enhanced the ability to accurately determine the subsurface structure and identified additional low-risk hangingwall plays at Saunders and Bilip and defined an additional high-risk sub -thrust structure within a footwall play underlying Bilip. These prospects were drilled in 2001 by Chevron (operator for PDL4) and Santos (operator for PPL190) in 2002. Both hangingwall targets proved successful in encountering hydrocarbons. Saunders-1 resulted in a discovery sharing a common OWC at 1260 mss with SE Gobe and extended the known southeast limit of the field. Bilip-1 encountered a deeper down-plunge structural compartment with a GOC established at 1722 mss and an OWC at 1735 mss. The footwall Iagifu sub-thrust target at Bilip was encountered at 2527 mss and was water wet. Geophysical modelling indicates remaining updip attic potential. These results highlight the increasing value of seismic definition in accurately imaging the subsurface and defining additional exploration opportunities within the Papuan thrust belt.
-
-
-
Imaging beneath the Taranaki fault, New Zealand – A feasibility study for wide-angle seismic surveys
More LessAuthors Ingo A. Pecher, Guy Maslen, Vaughan Stagpoole, Derek Woodward and Andrew R. GormanThe Taranaki overthrust fault at the the eastern boundary of New Zealand’s most significant hydrocarbon province, the Taranaki Basin, presents a challenge to the exploration industry. Seismic imaging of potential hydrocarbon traps beneath the fault is difficult using conventional streamers with lengths of up to 6 km, mainly because the high velocities of overthrust greywacke deviate rays to large offsets. We present a modelling study of the feasibility of wide-angle experiments to undershoot the fault. We performed ray-tracing to predict ray coverage at horizons beneath the fault as a proxy for expected image quality. Our results suggest that by increasing the length of seismic streamers from 6 to 12 km, which is possible with modern streamers, ray coverage across potential hydrocarbon traps would increase significantly. Deeper targets, which would be of interest for studying the tectonic evolution of the Taranaki Basin and the source rock distribution, will require “true” wide-angle experiments with ocean bottom seismometers. We also predict significant coverage with P-to-S mode-converted waves. Elastic full-wavefield modelling using a simplified velocity model suggests strong mode conversion on the greywacke/sediment interfaces.
-
-
-
Contact mapping from gridded magnetic data - a comparison of techniques
More LessAuthors Mark Pilkington and Pierre KeatingDelineating the edges of magnetised bodies is a fundamental application of magnetic data to geological mapping in areas of limited exposure. Especially in Precambrian shield-like regions, locating lateral changes in magnetisation of the outcropping crystalline rocks provides spatial information that is crucial in extending mapped geology into sparsely exposed or completely covered areas. Although not all magnetic contacts correspond to lithological contacts, the former provide key information on structural regimes, deformation styles and trends, and magnetic texture.
Many techniques for contact mapping have been developed, some originally based on profile (2-D) data and others designed specifically for grid-based (3-D) data sets. Here, we evaluate five methods applied to gridded data. The first three are based on finding maxima of the horizontal gradient magnitude of the total field (TF-hgm), tilt (TI-hgm) and pseudogravity (PSG-hgm). The fourth and fifth methods rely on locating maxima of the analytic signal (AS) and the 3-D local wavenumber (LW).
Method TF-hgm produces theoretically correct contact locations only when the data are reduced to the pole, and even then may produce false or secondary solutions mimicking contact trends. Method TI-hgm is less sensitive to field direction but also suffers from secondary maxima. Method PSG-hgm is perhaps the most established approach of those mentioned, and in the case of vertical contacts produces reliable maxima. However, knowledge of remanent magnetisation direction is required. Methods AS and LW theoretically produce maxima directly over contacts and are insensitive to magnetisation direction but are more sensitive to noise than the former, which limits their application to higher quality datasets.
-
-
-
UXO Location using Total Field Magnetics in SE Asia
More LessAuthors Timothy Pippett and Stephen LeeUnexploded Ordnance (UXO) is a major problem in a large number of areas throughout the world and SE Asia is no exception. As most of the older ordnance is of ferrous composition, total field magnetics is an eminently suitable tool for its location.
During the past few years work has been carried out on two major infrastructure development projects in SE Asia, one being in Hong Kong and the other in Taiwan. The projects have involved the use of tightly spaced total field magnetic traverses to locate buried UXO. The located items were then removed to allow for safe operations on the sites.
Both sites were intended for public access and thus required a high level of confidence that no UXO remained buried on site.
Both clients have adopted the results of the survey as an approved methodology for the location and removal of buried UXO.
-
-
-
Industrial evolution of depth imaging model building techniques - a Timor Sea case study
More LessAuthors Pierre Plasterie, Llew Vincent, Patrice Guillaume and Volker DirksThe evolution of depth imaging technologies in the past ten years has seen significant changes in the way we built velocity-depth macro models used by pre-stack depth migration algorithms. Both the observation of the data used for tomography inversion of macro models and the performance of inversion themselves have undergone significant changes in the recent years. This has lead to new model building methodologies.
Depth imaging has been (and still is) regarded as a state of the art – yet long and expensive – solution to improve the quality of the interpretation material. How can new trends in model building methodologies provide solutions to faster, yet high quality depth imaging products? We attempt to answer this question.
A case study is presented where a comparison is made between two pre-stack depth migration results obtained with two different depth imaging model building techniques performed on the same dataset. At an industrial scale through this example we try to identify the key elements of the methodology that can make differences to the quality of the final product and the speed at which we obtain it.
For obvious economic reasons, rare are the opportunities to look at new technologies on the same field dataset where the work has been performed both before and after the evolution. However, the opportunity to examine such evolution has been made possible for depth imaging model building techniques with this Timor Sea case study.
-
-
-
An improved pseudo-gravity magnetic transform technique for investigation of deep magnetic source rocks
More LessAuthors David A. Pratt and Zhiqun ShiThe pseudo-gravity transform is one of many possible FFT techniques that can be applied to aeromagnetic data. It enhances the anomalies associated with deep magnetic sources at the expense of the dominating shallow magnetic sources. This transform is an excellent interpretation tool for the detection of deep, magnetic igneous plutons and volcanic piles and the transformed data can be modelled using conventional gravity modelling tools. It is a suitable tool for interpreting deep-seated mineral plumbing systems associated with known, shallow mineral occurrences.
The pseudo-gravity transform is derived by an integration of the total magnetic intensity grid data using conventional FFT tools. Padding of the grid around the region covered by the survey introduces long wavelength artefacts into the transformed grids. These long wavelength artefacts can obscure the targets that are the object of investigation. A variety of regional residual separation procedures is applied to the transformed grid to minimise the impact of the long wavelength artifacts.
The improved pseudo-gravity transform is applied to the Goulburn 1:250 000 survey in the Lachlan Fold Belt of New South Wales to demonstrate the clear separation of deep sources that are difficult to detect or understand in the context of conventional magnetic image analysis. The results are contrasted with other filter techniques. Interpretive modelling of both the magnetic and the gravity transform data show how to derive more relevant geological information from magnetic surveys. By comparing the pseudo-gravity transform results with lower resolution ground gravity data, it is possible to obtain additional geological information by analysing the correlations.
-
-
-
In-situ formation stress determination using sonic logs and In- Borehole-Tomographic-Reconstruction of the near-borehole shear wave distribution
More LessAuthors Carsten Pretzschner, Harald Lindner and Hendrik RohlerFor advanced borehole logging interpretation of cross-dipole shear wave sondes with ‘In-Borehole-Tomographic-Reconstruction’ techniques an effective method for the determination of the mechanical stress distribution in the formation was developed.
The method includes the reconstruction of the 3D-shear wave velocity distribution in the immediate surrounding area of the borehole formation with tomographic algorithms.
The composition of the stress-induced shear wave velocity and comparison to the modelling results allow a quantitative interpretation. This log of the local formation stress anisotropy is excluded from both caliper-effects and the influences of shear wave velocity anisotropies caused by sedimentary effects.
As examples, logs of ratios of the local minimal/maximal formation-stress amplitudes Dsh/SH(z) were generated in a range of 0.86-0.99 for two gas production boreholes.
-
-
-
Practical 3D airborne EM inversion in complex terranes
More LessAuthors Art RaicheMany AEM inversion techniques have been applied successfully to locating isolated, highly conducting targets in uniform hosts. What happens when these targets occur in complex hosts overlain by non-uniform regolith with topographic features? In principle, a full 3D inversion is possible due to the availability of forward modelling programs capable of accounting for the full geoelectric complexity of any terrane. In practice, the inherent non-uniqueness coupled with modelling and data errors plus the requirement of substantial computation times render this approach unappealing for all but the most academic of inverters. On the other end of the complexity scale, simple imaging methods such as CDIs can fail badly for dipping targets and Born style approximations produce wrong models for moderate to high conductivity contrasts. This can be demonstrated by using accurate forward modelling algorithms to compute the EM responses of models produced by these methods.
Practical AEM data interpretation requires a technique that is capable of recognising drill targets in a variety of terranes but one that doesn’t require excess complexity or massive computing times. One such method, LeroiAir, is based on 3D thin sheet models in a layered, conducting host. The implicit forward model is fast, easy to set up and takes account of vortex and current gathering for any conductivity contrast. It is also capable of modelling the effects of other conductors such as faults and palaochannels. As expected, the thin-sheet inversion works well for simple isolated targets. As complexity increases, the inversion quality degrades but still yields useful drill target information in many cases.
-
-
-
Post-processing calibration of frequency-domain electromagnetic data for sea ice thickness measurements
More LessAuthors James Reid, John Bishop, Angus Munro, Andi Pfaffling, Kazu Tateyama and Tony WorbySea ice thickness measurements using electromagnetic (EM) instruments require accurate data. Calibration of sea ice thickness data acquired using a low induction number (LIN) EM sensor can be performed by conducting a geometric sounding at a range of heights over level sea ice of known thickness, and by comparing the observed data with the expected layered-Earth response. Calibration corrections for scaling, phase-mixing and zero-offset errors can be derived using least-squares inversion to minimise the misfit between the observed data and the theoretical response, and can be incorporated in modelling algorithms used to determine sea ice thickness.
This paper presents a case history illustrating identification and correction of calibration errors in low induction number EM data for Antarctic sea ice thickness measurements. Comparison of coincident EM measurements made using three identical LIN instruments showed that measured apparent conductivities disagreed by up to around 100mS/m, resulting in errors in the estimated sea ice thickness of up to 60%. Separate calibration corrections were determined for each instrument by analysis of geometrical sounding data acquired over level sea ice. Sea ice thickness at the calibration site was determined by making a large number of drilled thicknesses over the footprint of the EM instrument, and seawater and sea ice conductivities were determined using independent measurements. After application of the calibration corrections, sea ice thicknesses derived from the three instruments agreed closely with each other and with drilling results.
-
-
-
Footprints of airborne electromagnetic systems over onedimensional Earths
More LessAuthors James Reid and Julian VrbancichWe have used an inductive-limit model to compare footprint sizes for a variety of common airborne electromagnetic survey geometries. The model incorporates the flight height and orientation of the transmitter, and accounts for electromagnetic coupling between the induced current system and the receiver. Horizontal magnetic dipole transmitters are shown to have a smaller footprint than vertical magnetic dipole sources. Given typical survey heights for helicopter and fixed wing airborne electromagnetic systems, the helicopter VCX geometry has the smallest inductive-limit footprint (40 m), and the fixed-wing, towed bird system the largest (550 m for a system measuring the vertical component of B-field).
We also present preliminary calculations of frequency-domain airborne electromagnetic footprint sizes for the case of finite frequency or half-space conductivity. The original definition of the footprint is extended to be the side length of the cubic volume, centred below the transmitter, which contains the induced currents responsible for 90% of the secondary field measured at the receiver. Inphase footprint size for a horizontal coplanar helicopter frequency-domain system is shown to increase from around 3.7 times the flight height at the inductive limit to >9 times the flight height for induction numbers <0.7. The analysis also shows that the quadrature footprint is approximately two-thirds that of the inphase footprint, suggesting a higher spatial resolution for this component.
-
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