Geophysical Prospecting - Volume 72, Issue 2, 2024
Volume 72, Issue 2, 2024
- ORIGINAL ARTICLES
-
-
-
Elastic properties of unconsolidated sandstones of interest for carbon storage
More LessAuthors Colin M. Sayers and Sagnik DasguptaAbstractUnconsolidated sandstones are attractive targets for underground storage of carbon due to their high porosity and permeability. Monitoring of injection and movement of CO2 in such formations using elastic waves requires an understanding of the acoustic properties of the sandstone. Current approaches often use the so‐called soft‐sand model in which a Hertz–Mindlin model of the acoustic properties at high porosity is mixed with the acoustic properties of the mineral phase to predict the acoustic properties over the entire porosity range. Using well‐log data from two unconsolidated sand formations of interest for CO2 storage, we discuss the limitations of this model and provide an alternative approach in which the mechanical properties of grain contacts are obtained by inversion, and the properties of infill material lying within the pore space are estimated. The formations considered are the Paluxy Formation in Kemper County, Mississippi, and the Frio Formation near Houston, Texas. The ratio of the normal to shear compliance of the grain contacts is found to be significantly less than unity for both formations. This implies that the grain contacts are more compliant in shear than in compression. However, the grain contact compliance is higher and the ratio of the normal to shear compliance is lower for the Frio example than for the Paluxy case, and this may lead to sliding at grain contacts with low shear compliance and transport of grains during fluid flow, particularly if CO2 acts to weaken any cement that may be present at the grain contacts. Such transport was suggested by Al Hosni et al. in explaining why the magnitude of the time‐lapse effect due to the injection of CO2 at the Frio CO2 injection site is greater than predicted using conventional rock physics models. A simple model of the mechanical properties of infill material lying within the pore space suggests that the bulk and shear moduli of infill material in the Paluxy case are significantly higher than the Frio case, consistent with the lower grain contact compliance in the Paluxy case.
-
-
-
-
Shale distribution effects on the joint elastic–electrical properties in reservoir sandstone
More LessAbstractWe investigated the effect of shale distribution on the joint elastic wave and electrical properties of shaly reservoir sandstones using a dataset of laboratory measurements on 75 brine‐saturated (35 g/L salinity) rock samples (63 samples from the literature, 12 newly measured samples). All the data were collected using the ultrasonic (700 kHz) pulse‐echo measurement technique for P‐ and S‐wave velocities (Vp, Vs), attenuations (Qp−1, Qs−1), and a four‐electrode method for resistivity under elevated hydrostatic confining pressures between 10 and 50 MPa (pore fluid pressure 5 MPa). The distribution of volumetric shale content was classified by comparing the calculated dry P‐wave modulus to the modified Upper Hashin–Shtrikman bound for quartz and air mixtures, assuming pore‐filling shale. This scheme in particular allowed us to distinguish between pore‐filling and load‐bearing shale distributions according to idealized definitions, which provides new insight into the joint ultrasonic properties and resistivity behaviour for shaly sandstones. In resistivity–velocity space, the resistivity of load‐bearing shale increases with increasing velocity which form a more distinct trend with steeper gradient compared to those for partial pore‐filling shale and clean sandstones. Moreover, the pore‐filling shale trend straddles the clean sandstone trend and meets the load‐bearing shale trend between 100 and 150 apparent formation factors. In resistivity–attenuation space, the highest attenuations exist when the volumetric shale content is close to the frame porosity (for Qp−1 in particular), at the transition between pore‐filling and load‐bearing shales. The results will inform the development of improved rock physics models to aid reservoir characterization from geophysical remote sensing, particularly for joint seismic and controlled source electromagnetic surveys.
-
-
-
Fluid identification by nonlinear frequency‐dependent amplitude variation with offset inversion based on scattering theory
More LessAuthors Tianjun Lan, Zhaoyun Zong and Weihua JiaAbstractPre‐stack seismic amplitude variation with offset inversion is a crucial technique in seismic exploration, employed to estimate reservoir elastic parameters and thus reservoir fluid properties. However, traditional amplitude variation with offset inversion methods are based on elastic theory and linear approximation, neglecting the inelasticity of medium and nonlinear theory. To overcome this limitation, a quadratic scattering coefficient equation of the viscoelastic fluid factor is derived, which provides the basis for the equation for amplitude variation with offset inversion. Traditional amplitude variation with offset inversion methods typically neglect seismic dispersion and attenuation, failing to account for the influence of the seismic wave velocity attenuation and frequency variation. Quality factors of P‐ and S‐waves represent the degree of attenuation of seismic waves. To comprehensively address the effects of seismic wave dispersion and attenuation, a novel method called pre‐stack seismic nonlinear frequency‐dependent amplitude variation with offset inversion has been developed. This method builds upon the new quadratic scattering coefficient and is utilized for reservoir fluid prediction. The reliability and stability of the method are verified through synthetic and filed data examples. Further analysis reveals that the method is more reasonable and accurate compared to the traditional linear amplitude variation with offset inversion method. The results demonstrate that the proposed pre‐stack seismic nonlinear frequency‐dependent amplitude variation with offset inversion method can effectively identify reservoir fluids, providing a novel solution for reservoir fluid identification.
-
-
-
Permeability evolution in fine‐grained Aji granite during triaxial compression experiments
More LessAuthors Kazumasa Sueyoshi, Ikuo Katayama and Kazuki SawayamaAbstractTriaxial compression experiments were carried out on samples of fine‐grained Aji granite to measure the evolution of permeability during deformation prior to failure under confining pressures of 20 and 40 MPa. During the initial stages of deformation, a small decrease in permeability was observed, due to the closure of pre‐existing microcracks; permeability then increased with increasing differential stress. During deformation, permeability varied by up to two orders of magnitude, and we observed a small pressure dependence, with a larger variation observed at 20 MPa than at 40 MPa. This suggests that more cracks developed during brittle deformation under the lower confining pressure. The observed increase in permeability during our experiments was approximately proportional to inelastic volumetric strain, which corresponded to the volume of dilatant cracks. On the other hand, prior to brittle failure, we observed a further increase in permeability that was greater than the inelastic volumetric strain, suggesting crack aperture opening accelerated at this stress level (>∼80%). The permeability enhancement resulting from the crack dilation affects the gas‐sealing capacity of the fluid‐saturated caprock. Our experimental findings would be beneficial for the safety assessment in applications such as the long‐term storage of various gaseous wastes.
-
-
-
Well‐log integration and seismic‐to‐well tie off George V Land (Antarctica)
More LessAuthors Davide Gei, Giuliano Brancolini, Laura De Santis and Riccardo GelettiAbstractCenozoic sediments from the continental rise off the George V Land consist of silty/clayey facies variably rich in diatom ooze; these sediments hold a record of the glacial history of the Wilkes Subglacial Basin and are important for estimating the contribution of the East Antarctic Ice Sheet to global sea level changes during past transition to warmer climates. This is fundamental to predict future scenarios related to the global warming. The petrophysical properties of Antarctic marine sediments are influenced by the ice sheet dynamics and may affect the amplitude of seismic reflections. Through a seismic‐to‐well tie procedure, we investigate the origin of high amplitude reflections from Miocene‐early Pliocene deposits identified in two seismic lines crossing at the Integrated Ocean Drilling Program Expedition 318 drill site U1359. Downhole and core log measurements are preconditioned and merged to obtain complete velocity and density records from the sea floor to the bottom of the deepest of the four wells drilled at this site. We generate a synthetic trace by convolving the reflectivity series with the seismic wavelet obtained from the sea‐floor reflection and match the synthetic trace to the seismic data with a time variant cross‐correlation procedure. This procedure established a robust time‐depth relationship, not achievable from the available small‐offset seismic data. To our knowledge, this is the first seismic‐to‐well tying in the George V Land area. Based on results from synthetic data, the anomalous high amplitude seismic package can be linked to changes in density of sediments. Such changes are interpreted as representative for the alternation of diatom‐rich (warm climate) and silty‐clay layers with ice‐rafted debris (cold climate) inside the deposits. We suggest that the analysis of the characteristics and the distribution of similar seismic anomaly around Antarctica can give insight into the modality of past Antarctic ice sheet dynamics.
-
-
-
Revised models for evaluating hydrocarbon saturation of tight sandstone formations based on dielectric dispersion logs
More LessAuthors Xiaoyu Wang, Guangzhi Liao, Peiqiang Zhao and Zhiqiang MaoAbstractAccurate estimation of hydrocarbon saturation is significant for log interpreters to evaluate tight sandstone formations. It is difficult to use conventional wireline logs to calculate hydrocarbon saturation in complex reservoirs. Dielectric dispersion logs are processed using dielectric models to obtain accurate saturation assessment, especially in high salinity shaly sand formations. However, the dielectric models consist of complex refractive index method model and shaly sand model ignore the impact of clay effects on complex refractive index method model, and the dielectric dispersion logs generated by them appear several unphysical phenomena in the clay‐bearing formations, such as the permittivity data at the highest frequency are greater than those at the second highest frequency and the conductivity data at the highest frequency are lower than those at the second highest frequency. This paper modifies the dielectric models by introducing a combined conductive phase in complex refractive index method model to correct the unphysical behaviours. Revised models are combined with four permittivity and four conductivity logs acquired by dielectric scanner at 24, 102, 360 and 960 MHz to inverse the water saturation, water salinity, textural index for water‐bearing pore and fraction of clay‐related phase. The complex refractive index method model interprets the logs measured at 960 MHz and shaly‐sand model interprets other dielectric dispersion logs. The particle swarm optimization algorithm is used to search inversion parameters. The results of forward simulation show that the revised models can correct unphysical behaviours of dielectric dispersion logs generated by the original models. Moreover, we apply original models and revised models to laboratory core data. The results of experiments show that the mean relative error values of the original models are 19.38% and 23.60% and of the revised models are 5.54% and 6.28% in two cases. It shows that the revised models are more accurate than the original models to interpret dielectric dispersion logs of tight sandstone reservoirs.
-
-
-
Influencing factors of coal elastic parameters based on the generalized Gassmann equation: A case study in Yuwang colliery, eastern Yunnan province
More LessAuthors Yuyan Che, Yanhai Liu, Guangui Zou, Chaochao Jin, Fei Gong and Sandugash SatibekovaAbstractCoal elastic parameters are an important index reflecting the material composition and porosity, which can be used to guide reservoir evaluation and the safe mining of coal and coalbed methane. In this study, coal samples collected from coal seams #2, #3, #7 + 8 and #9 in the Yuwang colliery, eastern Yunnan, were studied, and the influences of organic matter, ash content and porosity on the elastic parameters of coal samples were investigated. The results show that the elastic parameters of coal are negatively correlated with organic matter content and porosity, and positively correlated with ash content. This means that the bulk modulus and primary‐wave velocity of coal decrease with a gradual increase in organic matter content or porosity and increase with a gradual increase in ash content. The main reason for this is that the ash modulus is typically higher than that of coal. When the ash content gradually increases, the contribution of ash to the equivalent modulus of coal increases, whereas the contribution of organic matter and pores decreases. In this study, the primary‐wave velocity of the target coal seam was fitted based on the generalized Gassmann equation. The calculation results were consistent with the observations and the relative error was within 7%. According to the distribution behaviour of elastic parameters, this provides adequate data support for the prediction of gas content in coal mines and ensures safe and green mining of coal resources.
-
-
-
Detection and picking of shear wave arrival for stiffness evaluation of highly porous chalk
More LessAbstractElastic wave velocities of compressional and shear waves propagating through sedimentary rocks are often coupled with information of bulk density to derive the rock stiffness. Acquiring the transit time of compressional and shear waves often involves manual picking of wave arrival times from wave trains recorded in the laboratory or by well‐logging tools. Picking the compressional wave arrival time is commonly accepted as straightforward. Oppositely, detecting the shear wave arrival and picking its arrival time is often troublesome because the transmitted shear wave partly converts to compressional waves and back to a secondary shear wave, concealing the transmitted shear wave arrival in the wave train. In laboratory settings, we illustrate the difficulty of shear wave detection in wave trains recorded on highly porous chalk plug samples from the Danish North Sea Basin. Wave trains were recorded on plugs dry, Tap‐water or Isopar‐L saturated during uniaxial strain compaction. The recorded shear wave trains showed two distinct features, which could be interpreted as the transmitted shear wave first arrival; we denoted them as early and late arrivals. However, as only one feature can mark the arrival of the transmitted shear wave, we propose a semi‐empirical disclosure strategy combining a graphical representation of stacked wave trains with rock physical modelling. By stacking recorded wave trains in a graphical strain–time–amplitude domain, we demonstrate that an early shear wave feature marks a converted shear to compressional to shear wave and not the transmitted shear wave. We used physical modelling to identify early shear wave features and illustrate the consequences of adopting a falsely interpreted shear wave on stiffness properties.
-
-
-
Seismic fluid identification method based on the joint PP‐ and SH–SH‐wave stochastic inversion
More LessAuthors Ying Lin, Guangzhi Zhang, Baoli Wang, Zhenyu Zhu, Jianhu Gao and Lin LiAbstractPre‐stack seismic inversion is an effective method for elastic parameter inversion using seismic data, which facilitates seismic fluid identification. However, pure PP‐wave inversion has issues of strong multi‐solution and limited prediction accuracy. Therefore, we propose a seismic fluid identification approach based on the joint PP‐ and SH–SH‐wave stochastic inversion. First, the linearized SH–SH‐wave amplitude variation with offset approximation parameterized by shear modulus and density is derived. Numerical simulations demonstrate that the SH–SH‐wave amplitude variation with offset approximation has a good accuracy. Reflection coefficient contribution analysis indicates that the new formulation has better parameter sensitivity to shear modulus and density than the PP wave amplitude variation with offset approximation derived by Russell, which helps one to improve the inversion of shear modulus and density. On this basis, we construct a joint inversion equation of PP and SH–SH waves for a Russell fluid indicator, a shear modulus and density and present a novel joint stochastic inversion method based on the ensemble smoother with multiple data assimilation. Stanford VI‐E model tests reveal that the Russell fluid indicator factor, shear modulus and density obtained from the joint PP‐ and SH–SH‐wave inversion have higher identification accuracy and smaller relative errors than those from pure PP‐wave inversion. Furthermore, field data tests indicate that this method has practical applicability in seismic fluid identification.
-
-
-
Three‐axis borehole gravity monitoring for CO2 storage using machine learning coupled to fluid flow simulator
More LessAuthors Taqi Alyousuf, Yaoguo Li, Richard Krahenbuhl and Dario GranaAbstractThe field of geophysics faces the daunting task of monitoring complex reservoir dynamics and imaging carbon dioxide storage up to several decades into the future. This presents numerous challenges, including sensitivity to parameter changes, resolution of obtained results and the cost of long‐term deployment. To effectively store CO2 subsurface, it is necessary to monitor and account for the injected CO2. The gravity method provides several advantages for CO2 monitoring, as changes in fluid saturation correspond directly and uniquely to observed density changes. Three‐axis borehole gravity has demonstrated significant promise as a next‐generation tool for reliably monitoring reservoir dynamics across a range of depths and sizes. However, the gravity inverse problem is highly ill‐posed, necessitating regularization that incorporates prior knowledge. To address this issue, we propose using a feed‐forward neural network, a machine learning method, to invert time‐lapse three‐axis borehole gravity data and monitor CO2 movement within a reservoir. By training the neural network on models that analyse changes in density and corresponding gravity responses resulting from perturbations made to the reservoir model, we can create scenarios that train the algorithm to identify unexpected CO2 migration in addition to the normal movement of CO2. Our method is demonstrated using reservoir models for the Johansen formation in offshore Norway. We convert reservoir saturation models into density changes and generate their corresponding three‐axis gravity data in a set of boreholes. Our results show that the developed machine learning inversion algorithm has high reliability and resolution for imaging density change associated with CO2 plumes, as demonstrated in the Johansen reservoir models utilized by the simulator. We also investigate machine learning inversion using regularization parameters and show that it is robust, with a strong tolerance for higher levels of noise. Our study demonstrates that the developed machine learning algorithm is a powerful tool for inverting three‐axis borehole gravity data and monitoring the migration and long‐term storage of injected CO2.
-
-
-
Gravity inversion by second order approximation applied to the Styrian Basin (Austria)
More LessAuthors Harald GranserAbstractAn approximative second order gravity inversion scheme by a truncated power series expansion is applied to derive the thickness of the Neogene, mostly clastic, sedimentary section of the Styrian Basins in South–East Austria, which are sub‐basins of the Pannonian Basin System. These sub‐basins with a derived thickness of up to 4 km are of interest for geothermal exploitation because of the increased geothermal gradient and heat flow observed in the Pannonian Basin in general and a geothermal gradient of 4°–5°/100 m measured in some wells in the Styrian Basin. The Styrian Basin also has been an area for hydrocarbon exploration in the past 50 years, with oil and gas show encountered in several exploration wells and one sub‐commercial gas discovery. The Miocene and Plio–Pleistocene volcanism in the Styrian Basin caused by Miocene crustal thinning is discussed in terms of the influence to the gravity inversion taking the aeromagnetic field into account. The volcanism is of relevance for the geothermal prospectivity but poses problems for the single layer‐based gravity inversion scheme. Results are discussed from a computational side comparing observed and calculated gravity fields but also the match with well data is discussed. In terms of gravity inversion methodology, the presented can be viewed as an approximative fast‐track approach.
-
-
-
2D fast Fourier transform analytical solutions in all space for all gravity and magnetic components
More LessAuthors Olivier BoulangerABSTRACTForward modelling of potential field data is an important part of optimization algorithms used to invert large datasets such as those involving rugged terrain or borehole data. Two‐dimensional fast Fourier transform modelling with a prism or a dipole is one of the most efficient methods compared to the forward modelling in the space domain. However, the exact solution of a prismatic source is limited to the case of a half‐space with the computation of data on a horizontal datum above the topography. Starting from the three‐dimensional Fourier forward modelling analytical formulation for a prism, an integration according to the wavenumber w is accomplished which allowed to find a two‐dimensional Fourier exact analytical formulation outside, at the interfaces of, and inside a prism for all potential field components. This new formulation requires the calculation of only four integrals. The gravity and magnetic fields are computed with this two‐dimensional fast Fourier transform formulation in the entire domain and compared with the analytical space domain and the three‐dimensional fast Fourier transform formulations. From the three‐dimensional calculated field, each component can be interpolated with the tri‐linear interpolation method along a borehole or on a drape surface simulating an airborne survey. Based on experiments demonstrated in this work, the two‐dimensional formulation in the Fourier domain gave accurate results with greater speed of execution in comparison to modelling in the space domain. The forward modelling method is tested on real gravity data from the north of Alberta (Canada).
-
-
-
Adopting normalized full gradient method for regional‐scale gravity modelling: A case study for Northwestern Iran
More LessAuthors Ako Alipour, Khalil Motaghi, Zahra Mousavi, Hamideh Cheraghi and Seyed Abdoreza SaadatAbstractThe normalized full gradient was developed to determine anomalous bodies, such as oil and gas fields or simple geological structures studies. We believe that even in complicated geology, normalized full gradient is practical. We introduce data preprocessing and use step‐by‐step simple‐to‐complicated synthetic tests to develop previous researchers’ ideas for regional‐scale gravity modelling. One of the most important steps of the normalized full gradient is the determination of the N optimum value. We found that prevalent methods such as the standard spectral or maxima method are feasible in simple structures only. So, we have suggested the imaging criteria routine for complicated cases. We trace maximum normalized full gradient responses to detect the normalized full gradient responses at the increasing harmonic numbers as the transition of the extensive part of the anomaly to the sharp part of that. With imaging criteria for the determination of N optimum values, the complicated synthetic test results show the success of the normalized full gradient to understand complicated gravity signals. In the real case, we have studied the Northwestern Iran normalized full gradient model of the Bouguer ground gravity data beneath the seismic profile and prepared a P receiver function depth section to uncover the geometry of the Moho boundary and important interfaces in the crust. We suggest the inferred synthetic model from the Bouguer ground gravity anomaly and P receiver function depth section to normalized full gradient trustworthy test in real cases. According to the synthetic test results, we understand the frame of the normalized full gradient responses in the semi‐real case and truthful responses in the real case. Along with this, we study the second ground gravity profile of Northwestern Iran in a good resolution to uncover the deeper structures. The real case results show the possibility of Moho offset and thinning lithosphere beneath the North Tabriz Fault lithospheric boundary, the possible source of Sahand volcanic centre at the west side of the Moho offset beneath North Tabriz Fault, the deep root of the Sabalan volcanic centre in the lower crust and the lithospheric and asthenospheric wedge with the density contrast beneath Sahand–Sabalan volcanic centres. One of the most important results of our study is the lithosphere–asthenosphere boundary offset and stepped Moho possibility beneath the Talesh Mts next to the South Caspian Basin boundary.
-
-
-
Magnetotelluric images of the medium enthalpy Bakreswar geothermal province within a granitic gneissic complex, Eastern Indian Peninsula
More LessAuthors Roshan K. Singh, Ute Weckmann and Shalivahan SrivastavaAbstractThe Bakreswar geothermal province represents a medium enthalpy geothermal system with its Bakreswar and Tantloie hot springs. It lies within the Chotanagpur Granite Gneissic Complex in the eastern part of the Indian Peninsula. The province has a high heat flow and a high geothermal gradient of 90°C/km. Magnetotelluric data from 95 sites in a frequency range of 10 kHz–10 Hz were acquired over the Bakreswar geothermal province to obtain an electrical conductivity model and map the geothermal reservoir with its fluid pathways and related geological structures. Subsurface conductivity models obtained from three‐dimensional inversions of the Magnetotelluric data exhibit several prominent anomalies, which are supplemented by gravity results. The conductivity model maps three features which act as a conduit (a) a northwest–southeast trending feature, (b) an east–west trending feature to the south of the northwest–southeast trending feature (which lies 1 km north of the Oil and Natural Gas Corporation fault marked by previous studies) and (c) shallow conducting features close to Bakreswar hot spring. The northwest–southeast trending feature coincides with the boundary of the high‐density intrusive block. This northwest–southeast trending feature provides the pathway for the meteoric water to reach a maximum depth of 2.7 km, where it gets heated by interacting with deep‐seated structures and then it rises towards the surface. The radiogenic process occurring within the granites of Chotanagpur Granite Gneissic Complex provides the heat responsible for heating the meteoric water. The northwest–southeast and east–west trending features are responsible for the transport of meteoric water to deeper depths and then towards the shallow regions of the Earth. The near surface features close to the Bakreswar hot spring are responsible for carrying the water further towards the hot spring. The resistivity of these structures plotted as a function of salinity and temperatures for saline crustal fluids suggests the involvement of meteoric water. Further, applying Archie's law to this resistivity suggests that the conduit path has a porosity greater than 10%. This study successfully maps the anomalous structures which might foster the migration of geothermal fluid in Bakreswar geothermal province.
-
Volumes & issues
-
Volume 74 (2026)
-
Volume 73 (2024 - 2025)
-
Volume 72 (2023 - 2024)
-
Volume 71 (2022 - 2023)
-
Volume 70 (2021 - 2022)
-
Volume 69 (2021)
-
Volume 68 (2020)
-
Volume 67 (2019)
-
Volume 66 (2018)
-
Volume 65 (2017)
-
Volume 64 (2015 - 2016)
-
Volume 63 (2015)
-
Volume 62 (2014)
-
Volume 61 (2013)
-
Volume 60 (2012)
-
Volume 59 (2011)
-
Volume 58 (2010)
-
Volume 57 (2009)
-
Volume 56 (2008)
-
Volume 55 (2007)
-
Volume 54 (2006)
-
Volume 53 (2005)
-
Volume 52 (2004)
-
Volume 51 (2003)
-
Volume 50 (2002)
-
Volume 49 (2001)
-
Volume 48 (2000)
-
Volume 47 (1999)
-
Volume 46 (1998)
-
Volume 45 (1997)
-
Volume 44 (1996)
-
Volume 43 (1995)
-
Volume 42 (1994)
-
Volume 41 (1993)
-
Volume 40 (1992)
-
Volume 39 (1991)
-
Volume 38 (1990)
-
Volume 37 (1989)
-
Volume 36 (1988)
-
Volume 35 (1987)
-
Volume 34 (1986)
-
Volume 33 (1985)
-
Volume 32 (1984)
-
Volume 31 (1983)
-
Volume 30 (1982)
-
Volume 29 (1981)
-
Volume 28 (1980)
-
Volume 27 (1979)
-
Volume 26 (1978)
-
Volume 25 (1977)
-
Volume 24 (1976)
-
Volume 23 (1975)
-
Volume 22 (1974)
-
Volume 21 (1973)
-
Volume 20 (1972)
-
Volume 19 (1971)
-
Volume 18 (1970)
-
Volume 17 (1969)
-
Volume 16 (1968)
-
Volume 15 (1967)
-
Volume 14 (1966)
-
Volume 13 (1965)
-
Volume 12 (1964)
-
Volume 11 (1963)
-
Volume 10 (1962)
-
Volume 9 (1961)
-
Volume 8 (1960)
-
Volume 7 (1959)
-
Volume 6 (1958)
-
Volume 5 (1957)
-
Volume 4 (1956)
-
Volume 3 (1955)
-
Volume 2 (1954)
-
Volume 1 (1953)
Most Read This Month