M I
2= O
o
Ol U1
ISSN 1408-7073
MATERIALS and GEOENVIRONMENT
MATERIALI in GEOOKOLJE
RMZ - M&G, Vol. 65, No. 1 pp. 001-058 (2018)
Ljubljana, March 2018
RMZ - Materials and Geoenvironment RMZ - Materiali in geookolje
ISSN 1408-7073
Old title/Star naslov
Mining and Metallurgy Quarterly/Rudarsko-metalurški zbornik ISSN 0035-9645, 1952-1997
Copyright © 2016 RMZ - Materials and Geoenvironment Published by/Izdajatelj
Faculty of Natural Sciences and Engineering, University of Ljubljana/ Naravoslovnotehniška fakulteta, Univerza v Ljubljani
Associated Publisher/Soizdajatelj
Institute for Mining, Geotechnology and Environment, Ljubljana/
Inštitut za rudarstvo, geotehnologijo in okolje
Velenje Coal Mine/Premogovnik Velenje
Slovenian Chamber of Engineers/Inženirska zbornica Slovenije
Editor-in-Chief/Glavni urednik
Boštjan Markoli
Assistant Editor/Pomočnik urednika
Jože Žarn
Editorial Board/Uredniški odbor
Cosovic, Vlasta , University of Zagreb, Croatia
Delijic, Kemal, University of Montenegro, Montenegro
Dobnikar, Meta, Ministry of Education Science and Sport, Slovenia
Falkus, Jan, AGH University of Science and Technology, Poland
Gojic, Mirko, University of Zagreb, Croatia
John Lowe, David, British Geological Survey, United Kingdom
Jovičic, Vojkan, University of Ljubljana, Slovenia/IRGO Consulting d.o.o.,
Slovenia
Kecojevic, Vladislav, West Virginia University, USA Kortnik, Jože, University of Ljubljana, Slovenia Kosec, Borut, University of Ljubljana, Slovenia Kugler, Goran, University of Ljubljana, Slovenia Lajlar, Bojan, Velenje Coal Mine, Slovenia
Malbašic, Vladimir, University of Banja Luka, Bosnia and Herzegovina
Mamuzic, Ilija, University of Zagreb, Croatia
Moser, Peter, University of Leoben, Austria
Mrvar, Primož, University of Ljubljana, Slovenia
Palkowski, Heinz, Clausthal University of Technology, Germany
Peila, Daniele, Polytechnic University of Turin, Italy
Pelizza, Sebastiano, Polytechnic University of Turin, Italy
Ratej, Jože, IRGO Consulting d.o.o., Slovenia
Ristovic, Ivica, University of Belgrade, Serbia
Šaric, Kristina, University of Belgrade, Serbia
Šmuc, Andrej, University of Ljubljana, Slovenia
Terčelj, Milan, University of Ljubljana, Slovenia
Vulic, Milivoj, University of Ljubljana, Slovenia
Zupančič, Nina, University of Ljubljana, Slovenia
Zupanič, Franc, University of Maribor, Slovenia
Editorial Office/Uredništvo
Technical editors/Tehnična urednika Teja Čeru and Jože Žarn Secretary/Tajnica Nives Vukič
Editorial Address/Naslov uredništva
RMZ - Materials and Geoenvironment
Aškerčeva cesta 12, p. p. 312
1001 Ljubljana, Slovenija
Tel.: +386 (0)1 470 46 10
Fax.: +386 (0)1 470 45 60
E-mail: rmz-mg@ntf.uni-lj.si
Published/Izhajanje
4 issues per year/4 številke letno
Partly funded by Ministry of Education, Science and Sport of Republic of Slovenia./Pri financiranju revije sodeluje Ministrstvo za izobraževanje, znanost in šport Republike Slovenije.
Articles published in Journal "RMZ M&G" are indexed in international secondary periodicals and databases:/Članki, objavljeni v periodični publikaciji „RMZ M&G", so indeksirani v mednarodnih sekundarnih virih: CA SEARCH® - Chemical Abstracts®, METADEX®, GeoRef.
The authors themselves are liable for the contents of the papers./ Za mnenja in podatke v posameznih sestavkih so odgovorni avtorji.
Annual subscription for individuals in Slovenia: 20 EUR, for institutions: 30 EUR. Annual subscription for the rest of the world: 30 EUR, for institutions: 50 EUR/Letna naročnina za posameznike v Sloveniji: 20 EUR, za inštitucije: 30 EUR. Letna naročnina za tujino: 30 EUR, inštitucije: 50 EUR
Transaction account/Tekoči račun
Nova Ljubljanska banka, d. d., Ljubljana: UJP 01100-6030708186
VAT identification number/Davčna številka 24405388
Online Journal/Elektronska revija
www.rmz-mg.com www.degruyter.com/view/j/rmzmag
Table of Contents Kazalo
Original scientific article Izvirni znanstveni članki
Thermal conductivity of rocks and geothermal water	1
Toplotna prevodnost kamnin in geotermalne vode Željko Vukelič
Geostatistical modeling of porosity data in 'oba' field, onshore Niger Delta	21
Oluseun Adetola Sanuade, Akindeji Opeyemi Fajana, Abayomi Adesola Olaojo, Kehinde David Oyeyemi and Joel Olayide Amosunl
Application of analytic element method in hydrogeology	35
Uporaba metode analitičnih elementov v hidrogeologiji Dragan Kaluderovič, Eva Koren, Goran Vižintin
Soil corrosivity and aquifer protective capacity of overburden units in Ado-Ekiti, southwestern Nigeria
Preiskava korozivnosti tal in sposobnosti površinskih plasti za varovanje vodonosnikov na	45
ozemlju Ado-Ekiti v jz Nigeriji
Oyedele, A.A., Olayanju, G.M., Talabi, A.O., Ogunyebi, S.N., Ojo, O.F.
Review paper Pregledni članek
Stochastic - advantages and uncertainties for subsurface geological mapping and	9
volumetric or probability calculation
Stohastika - prednosti in negotovosti v podpovršinskem geološkem kartiranju in računanju prostornin ali verjetnosti Tomislav Malvič
Historical Review Zgodovinski pregled
Instructions to Authors Navodila avtorjem
Historical Review
Zgodovinski pregled
More than 90 years have passed since the University Ljubljana in Slovenia was founded in 1919. Technical fields were united in the School of Engineering that included the Geologic and Mining Division, while the Metallurgy Division was established only in 1939. Today, the Departments of Geology, Mining and Geotechnology, Materials and Metallurgy are all part of the Faculty of Natural Sciences and Engineering, University of Ljubljana. Before World War II, the members of the Mining Section together with the Association of Yugoslav Mining and Metallurgy Engineers began to publish the summaries of their research and studies in their technical periodical Rudarski zbornik (Mining Proceedings). Three volumes of Rudarski zbornik (1937, 1938 and 1939) were published. The War interrupted the publication and it was not until 1952 that the first issue of the new journal Rudarsko-metalurski zbornik - RMZ (Mining and Metallurgy Quarterly) was published by the Division of Mining and Metallurgy, University of Ljubljana. Today, the journal is regularly published quarterly. RMZ - M&G is co-issued and co-financed by the Faculty of Natural Sciences and Engineering Ljubljana, the Institute for Mining, Geotechnology and Environment Ljubljana, and the Velenje Coal Mine. In addition, it is partly funded by the Ministry of Education, Science and Sport of Slovenia.
During the meeting of the Advisory and the Editorial Board on May 22, 1998, Rudarsko-metalurski zbornik was renamed into "RMZ - Materials and Geoenvironment (RMZ - Materiali in Geookolje)" or shortly RMZ - M&G. RMZ - M&G is managed by an advisory and international editorial board and is exchanged with other world-known periodicals. All the papers submitted to the RMZ - M&G undergoes the course of the peer-review process. RMZ - M&G is the only scientific and professional periodical in Slovenia which has been published in the same form for 60 years. It incorporates the scientific and professional topics on geology, mining, geotechnology, materials and metallurgy. In the year 2013, the Editorial Board decided to modernize the journal's format.
A wide range of topics on geosciences are welcome to be published in the RMZ - Materials and Geoenvironment. Research results in geology, hydrogeology, mining, geo-technology, materials, metallurgy, natural and anthropogenic pollution of environment, biogeochemistry are the proposed fields of work which the journal will handle.
Že več kot 90 let je minilo od ustanovitve Univerze v Ljubljani leta 1919. Tehnične stroke so se združile v Tehniški visoki šoli, ki sta jo sestavljala oddelka za geologijo in rudarstvo, medtem ko je bil oddelek za metalurgijo ustanovljen leta 1939. Danes oddelki za geologijo, rudarstvo in geotehnologijo ter materiale in metalurgijo delujejo v sklopu Naravoslovnotehniške fakultete Univerze v Ljubljani. Pred 2. svetovno vojno so člani rudarske sekcije skupaj z Združenjem jugoslovanskih inženirjev rudarstva in metalurgije začeli izdajanje povzetkov njihovega raziskovalnega dela v Rudarskem zborniku. Izšli so trije letniki zbornika (1937, 1938 in 1939). Vojna je prekinila izdajanje zbornika vse do leta 1952, ko je izšel prvi letnik nove revije Rudarsko-metalurški zbornik - RMZ v izdaji odsekov za rudarstvo in metalurgijo Univerze v Ljubljani. Danes revija izhaja štirikrat letno. RMZ - M&G izdajajo in financirajo Naravoslovnotehniška fakulteta v Ljubljani, Inštitut za rudarstvo, geotehnologijo in okolje ter Premogovnik Velenje. Prav tako izdajo revije financira Ministrstvo za izobraževanje, znanost in šport.
Na seji izdajateljskega sveta in uredniškega odbora je bilo 22. maja 1998 sklenjeno, da se Rudarsko-metalurški zbornik preimenuje v RMZ - Materiali in geookolje (RMZ -Materials and Geoenvironment) ali skrajšano RMZ - M&G. Revijo RMZ - M&G upravljata izdajateljski svet in mednarodni uredniški odbor Revija je vključena v mednarodno izmenjavo svetovno znanih publikacij. Vsi članki so podvrženi recenzijskemu postopku.
RMZ - M&G je edina strokovno-znanstvena revija v Sloveniji, ki izhaja v nespremenjeni obliki že 60 let. Združuje področja geologije, rudarstva, geotehnologije, materialov in metalurgije. Uredniški odbor je leta 2013 sklenil, da posodobi obliko revije.
Za objavo v reviji RMZ - Materiali in geookolje so dobrodošli tudi prispevki s širokega področja geoznanosti, kot so: geologija, hidrologija, rudarstvo, geotehnologija, materiali, metalurgija, onesnaževanje okolja in biokemija.
Glavni urednik
Editor-in-Chief
Original scientific paper	Received: November 14, 2017
Accepted: November 29, 2017 DOI: 10.2478/rmzmag-2018-0009
Thermal conductivity of rocks and geothermal water Toplotna prevodnost kamnin in geotermalne vode
Željko Vukelič1,*
University of Ljubljana, Faculty of Natural Sciences and Engineering, Aškerčeva cesta 12, 1000 Ljubljana, Slovenia * zeljko.vukelic@ntf.uni-lj.si
Abstract
In this article, the possibilities of use of geothermal energy in relation to the geothermal gradient and aquifer yield are described. Calculations represent information on potential geothermal water reserves that are not affected by cold return water inflow from the reinjection well after a certain period of production time. The calculations apply for continuous production of geother-mal water from the aquifer without significant pumping breaks.
Key words: geothermal gradient, aquifer, reinjection well
Izvleček
V članku obravnavamo možnosti izrabe geotermalne energije na temelju geotermalnega gradienta in izdatnosti vodonosnika. Izračuni predstavljajo informacijo o potencialnih zalogah geotermalne vode, ki po določenem času izkoriščanja, še ne padejo pod vpliv dotekanja hladne vode skozi reinjekcijsko vrtino. Izračuni veljajo za kontinuirano izkoriščanje vodonosnika brez daljših prekinitev črpanja geotermalne vode iz vodonosnika.
Ključne besede: geotermalni gradient, vodonosnik, re-injekcijska vrtina
2
Introduction
Geothermal energy belongs among renewable energy resources. Geothermal energy finds widespread use in the heating of residential and other buildings, as well as in electrical energy production [1]. Use of geothermal energy for the heating of residential and other buildings represents >80% of all heating energy needs in some countries (e.g. Iceland]. It is very important, therefore, to have knowledge about how to use geothermal energy in an optimal manner and not disturb the balance of underground water reserves. In other words, we need to ensure the renewal of underground water reserves and eliminate the effect of aquifer cooling because of reinjection of used, colder return water.

E-h
A-(T2- 71) • t where
E = amount of heat passing through a particular area A (in J);
H = strata thickness (in m); A = particular area through which the heat flux is passing (in m2);
T2 - T1 = temperature difference between strata borders (in K);
t = heat flux passing time (in s).
For North-East (NE) Slovenia, we can write down the thermal conductivity in terms of rock density:
X = 0.142 x p2
(3)
Thermal characteristics of rocks and underground waters
The thermal characteristics of rocks and underground waters can be defined using the following equation [1-3]:
A = a • c • p	(1)
where
A = thermal conductivity (in W/m • K); c = specific heat (in J/kg • K); a = temperature conductivity (in m2/s); p = density (in kg/m3).
The thermal characteristics are not directly dependent on temperature, as with a variance of temperature, physicochemical changes occur and thus initiate the change of thermal characteristics.
Thermal conductivity X
Thermal conductivity is the transition of heat from warmer heat sources to cooler ones. In hydrothermal deposits, heat transition occurs through convection in fluids (physical heat transfer in fluids) and conduction in rocks (physical heat transfer in solid substances, e.g. because of radiation).
Thermal conductivity is determined experimentally in laboratories and in the field, but it can be calculated too, as follows:
Specific heat c
Specific heat is a physical quantity that describes the thermal property of the substance and defines how much energy is needed to increase the temperature of 1 kg mass for 1K, at constant pressure. For the rocks, the average value is cm = 835 (± 15%] J/kg • K; and for water, the value is cw = 4187 J/kg K. The calculation of specific heat of the rock in terms of porosity and at water density 1000 kg/m3 and average rock density 2720 kg/m3 is done as follows:
1	2.72
ck = ^-cw--+{l-^)-cm-—-	(4]
g	g
where
F = porosity.
Specific heat of the rocks in NE Slovenia is
ct = 0.602-e-1177'H + 0.898.
k
Volume-specific heat capacity (cp)
This parameter describes a product of specific heat and density and is defined as the amount of energy needed to raise the temperature of 1 m3 of substance by 1 K at constant pressure.
(cp\ = (cp\• F + (cp)m • (1 - F] (in J/m3 • K] (5]
Volume-specific heat capacity is a direct parameter used for the assessment and calculation of
100 200 300 400
Temperature gradient {m°K/m)
Figure 1: Field determination of temperature gradient.
geothermal energy reserves. In hydrothermal deposits, the rocks are saturated with water; thus, the following relation is derived:
(cp\ = (cp)w • F + (cp)m • (1 - F) (in J/m3 • K) (6)
Volume-specific heat capacity (cp)w is affected by the water density (dependent on mineralization and ranging between 1000 and 1400 kg/m3) and temperature. The range of the term "(cp)m" depends on rock density.
Temperature conductivity or diffusivity (a)
Diffusivity is characterized by the rate of temperature equalization under nonstationary heat conduction and is represented as follows:
1
a = — • (in m2/s) cp
(7)
Geothermal gradient (G)
The temperature at the Earth's surface depends mostly on the radiation of the sun. On average, the temperature from a depth of 30 m is independent of the sun's radiation. The mean value of the geothermal gradient for Europe is 0.03 °C/m, but in Slovenia, the newest research results show the value of geothermal gradient as approximately 0.06 °C/m.
Gt=-
T-Ta
H
X
(8)
where
T = temperature measured at specific depth (K]; T0 = average mean ground temperature for NE Slovenia (11.6°C);
H = depth of temperature measurement (m).
The geothermal gradient is directly proportional to the heat flux but inversely proportional to heat conduction, which changes with depth as
the rock density changes too. Temperature increase in NE Slovenia can be described by the followin g equation:
where W = AW^ t . The total amount of the
pc	c
produced geothermal water is calculated as follows:
Tmf —
576 •H3
[1.09 • (e-
1) + 2.72 • H]2
+ 11.6 (9)
A-h- (cp)L W„r =-:—;-(in m3)
"pc
(cp)v
(11)
Based on the geothermal gradient, the projections for the prospect of use of a single area for geothermal energy are made [8].
Amount of geothermal water that can be produced at constant temperature
Production time at constant temperature represents the relation between production reserves and well output or well yield:
A-h- (cp\
xr =--—-— (years)
c Wp-(cp)w (y )
(10)
A = flow area of the deposit (in m2); h = flow height of the deposit (in m); (cp)w = volume-specific heat capacity of water (in J/m3 ■ K);
(cp)L = volume-specific heat capacity of the deposit (in J/m3 ■ K);
Wpc = annual geothermal water production (in m3/year);
Tr >
A-h- (cp), Vw • (cp)w
(12)
qw = projected production of geothermal water (in m3/day).
Figure 3: Interdependency of temperature and porosity in NE Slovenia.
Results and discussion
We are interested in how long the water will retain the constant temperature without any impact of the cold return water inflow (random data) [4-7]:
VL = flow volume of the deposit (80000000 m3); qW = projected extraction amount of geothermal water (7000 m3/day);
F = porosity of the productive strata (16%); cm = specific heat of the massif 882.2 J/kg • K, and density 2700 kg/m3;
cw = specific heat of water, which is defined as follows:
specific heat of water, at hydrostatic pressure of water 247 bars and temperature 175°C, equals 4271 J/kg • K, but as CO2 (26 m3/m3) is present
in the water, the value should be reduced for ~5% and also equals 4058 J/kg • K. pw = water density in the deposit (923.6 kg/m3) Volume-specific heat capacity of water is derived as follows:
(cp)w = 4058 x 923.6 = 3.748 x 106 J/m3 K.
Volume-specific heat capacity of the deposit is derived as follows:
(cp)l = (cp)w • F + (cp)m • (1 - F) = 3.748 x 106 x 0.16 + 882.2 x 2700 x (1 - 0.16);
(cp)L = 2.60 x 106 J/m3 K.
Figure 4: Aquifer exploitation time (in days).
Exploitation of the water with constant temperature from the deposit can last for the following duration:
A-h- (cp)L _ 68000000 x 2.6 x 106 Tc ~ Qw • (cp)w = 7000 x 3.748 x 106 =
= 6739 days « 19 years.
After the period of 19 years, the water temperature will start to decrease because of colder return water inflow from the reinjection well. This finding stands in case the geothermal water is exploited for production of electrical energy and the reinjection well is not sufficiently far away from the production well.
Conclusions
Geothermal energy belongs among renewable energy resources. Geothermal energy is widely used for the heating of residential and other buildings, as well as for electrical energy production. In NE Slovenia, we have the largest potential for exploitation of geothermal energy. Exploitation of geothermal ener-
gy in NE Slovenia is already causing problems due to excessive production from geothermal aquifers. Water table drawdown is observed in aquifers, meaning that not enough used water is returned back into the aquifer. Calculations present information on the potential geother-mal water reserves that are not affected by the cold return water inflow from reinjection wells after a certain period of production time. The calculations apply for continuous production of geothermal water from the aquifer without significant pumping breaks.
References
[1]	DiPippo, R. (2012): Geothermal power plants: principles, applications, case studies, and environmental impact. Third Edition, Amsterdam: Elsevier, 624 p.
[2]	Golub, M., Kurevija, T. Koscak Kolin, S. (2004): Thermodynamic Cycle Optimization in the Geothermal Energy Production. Rudarsko-Geolosko-Naftni Zbornik, 16, pp. 65-69.
[3]	Golub, M., Kurevija, T., Pravica, Z. (2006): Maximum Energy Output of Geopressured Geothermal Reservoirs in Croatia, International Congress Energy and
Environment 2006, Opatija, Croatia, 25-27 October 2006, pp. 121-130.
[4]	Kaluderovic, D., Vižintin, G. (2014): Modflow USG - the next step in mathematical modelling of underground water. RMZ - Materials andgeoenvironment, 61(4), pp. 241-247.
[5]	Vižintin, G., Vukelič, Ž., Vulic, M. (2008): Monitoring the geothermal potential of deep tertiary aquifers in north-east Slovenia using old abandoned oil and gas wells = monitoring geotermalnog potencijala dubokih tercijarnih vodonosnih horizonata u severnoistočnoj Sloveniji uz pomoč starih napuštenih bušotina nafte i gasa. In: Ristovic, I. (ed.). Savremene tendencije u razvoju energetskog rudarstva = Modern Tendencies in the Development of Energy Mining: zbornik radova = proceedings. Beograd: Rudarsko-geolološki fakultet, pp. 39-52.
[6]	Vukelič, Ž., Vulic, M. (2014): Ocena in natančnost ocene 3D-položaja točk v vrtini = Evaluation of 3D positions and the positional accuracy of points within a borehole. Geodetski vestnik, 58(2), pp. 327-341.
[7]	Vukelič, Ž., Kraljic, M., Dervarič, E. (2013): Lendava -the first geothermal city in Slovenia. 5th Jubilee Balkan Mining Congress, Ohrid, Macedonia, 18th-21st September 2013. Skopje: Association of mining and geological engineers of Macedonia, pp. 1-12.
[8]	Rman, N., Lapanje, A., Rajver, D. (2012): Analiza uporabe termalne vode v severovzhodni Sloveniji = Analysis of thermal water utilization in the northeastern Slovenia. Geologija, 55(2), pp. 225-242.
Review paper	Received: Aug 21, 2017
Accepted: Sep 25, 2017 DOI 10.1515/rmzmag-2018-0003
Stochastic - advantages and uncertainties for subsurface geological mapping and volumetric or probability calculation
Stohastika - prednosti in negotovosti v podpovršinskem geološkem kartiranju in računanju prostornin ali verjetnosti
Tomislav Malvic1
1 University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering, Pierottijeva 6, 10000 Zagreb, Croatia * tomislav.malvic@rgn.hr
Abstract
Stochastic, especially simulation, occasionally could be found in different geological calculations, mostly as the most advanced mapping method. Its main attribute is description of uncertainties that are inherent not only to any geological mapping dataset but also to any volumetric or probability calculation. Here are presented uncertainties in all three cases - mapping, volume calculation and probability calculation - and reasons why and when to use stochastic in them. The stochastic, and consequently simulation, is a recommended tool in case of a low number of data (<15 inputs) or large dataset (>40 inputs), but in both cases, the descriptive statistics needs to be known and is reliable. Almost the same could be applied in volumetric calculation, but the success of stochastic in probability calculation depends on large datasets, with 15 or more inputs.
Key words: simulations, number of input data, stochastic mapping, volume calculation, probability, Croatia
Izvleček
Stohastiko, zlasti simulacije, lahko občasno najdemo pri različnih računih, zlasti pri najsodobnejših metodah kartiranja. Njihova poglavitna značilnost je opisovanje negotovosti, ki so lastne sleherni zbirki podatkov geološkega kartiranja, pa tudi slehernemu volumetrijskemu ali verjetnostnemu računanju. Tu predstavljamo negotovosti v vseh teh treh primerih - pri kartiranju, računu prostornin in računu verjetnosti, ter načine, čemu in kdaj uporabljati pri tem stohastiko. Stohastični pristop in posledično simulacijo je priporočljivo uporabljati v primeru majhnega števila podatkov (manj od 15) ali velike datoteke (nad 40 inputov), vendar mora biti v obeh primerih opisna statistika znana in zanesljiva. Skoraj isto velja tudi za volumetrijo, medtem ko lahko uporabljamo stohastiko v računu verjetnosti ugodnega izida samo v primeru večjih datotek s 15 inputi ali več.
Ključne besede: simulacije, število podatkov, stohastika, kartiranje, računanje prostornine, verjetnost, Hrvaška
Introduction
Stochastic simulation or Gaussian simulation (sequential or indicator] is a special geosta-tistical method based on different algorithms compared to deterministic interpolation methods such as Kriging and Cokriging [1-3]. Differences are a result of extensions introduced in the Kriging algorithm that can have advantages or disadvantages, due to introduction of uncertainties in estimations. Consequently, the selection between Kriging- and Gaussian simulation-based algorithms is very important and asks for experienced professionals [4]. The most common property of simulation is calculation of numerous realisations (values] for each cell in grid (excluded are hard data in conditional ones]. The requirement is input dataset characterised using normal distribution. The total set of realisations is characterised with uncertainties, derived from the size of dataset, variogram model and measurement errors. As input dataset, such an error is also characterised using normal distribution. The simulation obviously calculates an enormous number of new grid values (103 times larger than the input dataset]. Sometimes, it is used for artificially increasing of dataset, combining simulated and input values. As a consequence, descriptive statistics and histogram for analysed data can be easily and clearly calculated. In the grid of 50 x 50 cells, in 100 realisations, totally 250000 values are calculated. Hence, the input dataset of usually 10-20 hard data is enlarged in the scale of 104. Moreover, numerous realisations give as outcomes of numerous maps. All of them are equally probable, and some of them can be selected as "representative", but always at least three. Such a selection is done based on the order of calculation, random sampling, calculation of total map cells' values, etc., but selection always needs to be unbiased. On contrary, if intention is given only to single map as an outcome, then Kriging or Cokriging methods are chosen as algorithms made just for such a purpose. Simulations could be conditional and (rarely] unconditional and also Gaussian and (rarely] indicator. The sequential Gaussian simulation (SGS] could be applied to almost all geological variables characterised using Gaussian distri-
bution, naturally or after transformation. Examples of such variables are porosity, depth, thickness or permeability. There are several subsurface structures in the Croatian part of the Pannonian Basin System analysed by stochastic methods. The Klostar structure (and hydrocarbon field] is the most Croatian geological structure analysed by stochastic geostatistical algorithms till now. An example of how to apply SGS in subsurface mapping of porosity, depth and thickness can be given as a set of maps taken from Zelenika and Malvic [5]. The data were collected from a reservoir of Lower Pontian age in the Klostar field, located in the western part of the Sava Depression. Simulation used present-day depth, thickness, and locations of areas with higher porosity to successfully reconstruct [5] paleo-depositional environments and the distribution of different lithotypes of turbidites (Figure 1]. However, sequential indicator simulation (SIS] is more specific and used only when mapping is based on cut-off (threshold] values defined for the selected variable. Such an application is most often used in lithofacies mapping, when cut-off selected for porosity or thickness, directly or indirectly, indicates lateral changes in lithofacies. Such lateral changes are one of the main properties of turbidite arenite lithofacies in Neogene, Croatian part of the Pannonian Basin System, developed due to paleotopography, current directions and different lateral densities of turbidite. As results in all Croatian hydrocarbon fields reservoir sandstones laterally gradually are changed into marly sandstones, sandy marls, clayey marls and, eventually, marls. One of indicator applications in the Croatian Pannonian Basin System (CPBS) Neogene litho-stratigraphic units had been published for Lower Pontian, the Klostar Ivanic Formation in the Klosar Field (Figures 2 and 3]. Porosity and thickness had been mapped using the mapping indicators and interpreted regarding probability to reach at least the selected porosity value. Resulting maps indicated depositional channel of sandstone, transitional lithofacies to marls, transport direction of turbidites and role of regional fault [6,7].
Figure 1: The first (left) and 100th (right) SGS realisation for porosity. The selection is based on the simple order of calculation. Taken from [5].
Figure 2: Direction of material transport during Lower Pontian, Klostar Structure. Probability maps for porosity >19% (left) and >20% (right). Red indicates fault, and black dot indicates well. Taken from [7].
X
Figure 3: Direction of material transport during Lower Pontian, Klostar Structure. Probability maps for thickness >13 m. Red indicates fault, and black dot indicates well. Taken from [5,7].
Basics of stochastic simulation algorithm
The basic condition for performing SGS is (approximately) normal or Gaussian distribution of input data. Transformation into normality, if possible, is also allowed. Normal distribution is characterised using statistical properties such as expectation and standard deviation (m and a2, respectively), which are basic conditions for calculation of uncertainty range in cells and for estimation of errors.
Properties of input data, type of simulation and zero realisation
All data represented with point values are, so-called, "hard data", which means that they are always constant in space whatever outcome is presented. Consequently, such constants value-in-place in simulation define it as "conditional simulation", where are hard-data un-changeable across the plane or space. It means that inputs will not change whenever simulation algorithm is applied. In addition, there is another type of SGS called "uncowditional simulation". "Unconditional" means that input data are not treated as constants. It is a valid view because each grid is defined with numerous cells. Almost always each cell includes one or zero hard data. However, hard data are considered also as point data, i.e. infinity small point, compared with cells that in the 2D area are always much larger than single point. As a consequence, the so-called "hard-data" if moved through the cell area could easily change value. Consequently, it is why they cannot be considered as constants. Eventually, it is allowed to simulate new value into cells with "hard data". Whatever type is chosen, any simulation is based on the variogram model and Kriging interpolation.
The single Kriging map represents deterministic solution for input dataset and is called "zero (or Kriging] realisation". Such a realisation has known mean value (expectation], standard deviation (m, a2], Kriging variance (gk] and interval margins of simulated values (±3a around m, i.e. probability of 99% to include all possibilities]. Using this realisation, it is possible to perform required number of subsequent realisations for each simulated cell in the grid. As basic Kriging algorithm is very often used, Ordinary Kriging technique (Equation 1, taken from [3]], which is used in the Croatian geological mapping, is the most commonly used Kriging technique until now (e.g. [8]].
rlZt	-Zr a	rlZt	-z2 a .	■ rlZt	-zn a	i-
rlZw	-ZD	rlZw	-z2g .	■ rlZw	-zng	i i
Ytfn	-zo	YeZn	-^a .	■	-zn a	1 1
r^rl		"Kxr -	*ai
¿2		Ylx2 -	*a
		H*n -	*a
-M-		-1	-
(1]
Simulation
Normalised data, sometimes after transformation, with known N(m, a] are used to simulate values into grid cells. Simulation is based on, so called, two methods of introducing "randomness", which makes the entire process unbiased and repeatable (i.e. sequential]. The first randomness method is selection of estimated (simulated] cells inside the grid. When cells are selected, the value is calculated using Kriging (or Cokriging] on hard data. Previously estimated cells in each consequent selection are considered as new "hard data", using the same variogram model as in the first simulation, but repeat the "zero realisation". Hence, each estimated cell is characterised with a new value and also with an uncertainty interval, wide ±3a around expectation in that cell. Eventually, entire realisation has its own variance calculation from "zero realization". It is also known as Krig-ing variance.
The second randomness method, i.e. introduction of stochastic into simulation for the second time, includes cell value estimation and is depended on interval ±3a, different for each cell. Random selection of any value from this interval represents the final cell value in this realisation. Each cell is also characterised using its own probability distribution function (PDF]. This is why it is possible to calculate almost infinite number of equally probable realisations in one simulation.
x
Calculation of numerous realisations
Obviously, it is possible to calculate numerous realisations, which all have the same cell val-
Table 1: Qualitative estimation of simulation reliability and recommended number of realisations.
Completely fulfilled Partially fulfilled
Rarely (irregularly) or not at all fulfilled
Could variable be characterised with normal distribution?	Mostly (like porosity)	Sometimes (after transformation, like permeability)	Not sure; test is needed (variables like depth or thickness]
Could variable be ranked by map values?	Always (any realisation can be ranked according to their cumulative value)		
Which number of realisations represent real spatial uncertainties' characteristic to data?	100 or more realisations could be calculated	10-99 realisations could be calculated	Less than 10 realisations could be calculated (only quick insight]
ues only in hard data locations. Of course, it is valid only if conditional simulation has been performed, which is the most common case in geological mapping. However, the main task is always selection of the most appropriate realisations, regarding researching goal, graphical looks and, sometimes, cross-validation. The main purpose of simulation is presentation of numerous possible solutions, but obviously dozens or hundreds of maps are impossible to present as a single outcome. Hence, there are several quick steps to check how reliable simulation results are and which number of realisation could be appropriate. It is summarised in Table 1.
There is a general opinion that 100 realisations in simulation represents a large enough number of solutions that stochastic is representative for quality input dataset. Such a large set of realisations need to be post-processed to reach several realisations that are the most common characteristic for analytical purpose and simulated space. Such selection techniques are called ranking. The ranked variable is usually cumulative summation of all cells in the observed realisation. For example, if porosity had been mapped in some layer or unit, all cell values in one realisation can be summed and give total porosity for the entire realisation. If there are 100 realisations, such summation can be repeated 100 times, one for each realisation. Eventually, an entire set of 100 realisations can be ranked from the lowest total porosity per realisation (P0], through the median realisation
(P50), to the highest ranked realisation (P100). "P" value defines how many realisations have lesser total score than the observed one, i.e. for the P0, there is 0% of lesser solutions, and for the P100, there is 100% realisations with a low total value.
However, such a method of ranking is not always possible to perform, because summation of cells across realisation would not result in a meaningful value. In such cases, especially when only a small number of realisations are performed (e.g. Table 1, the right-end column], some pure statistical techniques could be used. Selection of only the first and last realisations in sequence or each nth realisation or any number of realisations by the "random seed number" algorithm is allowed. Such a selection is usually applied when simulation includes only a small number of realisations (5 or 10], made from small datasets with the purpose of getting quick visual insight into areas with the largest uncertainties.
To summarise, each new realisation needs to fulfil two conditions: (a] order or simulated cells is defined completely randomly, i.e. "random seed number generator" is applied and (b] consequently, number of "hard data" values inside variogram ellipsoid does not need be equal in the same cell and different realisations.
Stochastic derived from deterministic models
Stochastic and deterministic models are closely connected. In fact, any deterministic model, map or equitation, is an accepted approximation of input dataset, characterised with "artificial certainty". Geological deterministic models could be different, but mathematically, they are often purely numerical, such as calculation of probability of success (POS; e.g. [9-11]), volumetric calculation of structure (e.g. [12]) and graphical outcome like different deterministic maps (e.g. [13]). The last one is described in the previous chapter, and in the following are given models of probability and volumetric calculations.
Stochastic in geological risk or probability of success calculation
Deterministic models are just approximation of stochastic systems, where the variable properties are geologically and statistically known. However, increase in data always partially changes deterministic solution, like maps or statistics. It was proven, e.g. in the calculation of POS that a well-known deterministic model can be used for hydrocarbon discovery in plays and prospects in the Croatian offshore [11] or the CPBS [9]. Even such numerically simple multiplication could be partially described and modified with inclusion of stochastic (e.g. [14,15]). This is why because calculation of such values can be a significant uncertainty process, even if it is applied in geological areas covered with published probability tables for geological categories such as existence of trap, reservoir, source rocks, effective migration of fluids and preservation of hydrocarbons in reservoir. Consequently, Equation 1 is used:
POS = p(t) x p(r) x p(s) x p(m) x p(p).	(1)
Here, POS is the geological probability, p(t) the probability for trap existence, p(r) the reservoir existence, p(s) the source rock existence, p(m) the probability for effective migration and p(p) the probability that hydrocarbons are preserved.
Although raised in hydrocarbon geology, this methodology is easily modified in other geological disciplines, like storage of CO2 in the subsurface [16]. However, it was proven [15] that at least three (sub)categories could be represented by stochastic simulation. Porosity subcategory maps make possible to directly calculate values such as minimum, median and maximum realisations. Hydrocarbon shows and quality of cap rocks are descriptive variables, but often they can be expressed as percentage. Hence, descriptive variables can be transferred into indicators or real number. For example, detection of new gas in the layer can be observed if concentration is >10% (indicator 0 or 1). The measured value of 15% can indicate that seal rocks are not completely impermeable, i.e. probability of sealing is 0.75. Such a probability can be lowered for any critical value of new gas (like 40% is 0.50, 60% is 0.25, 80% is 0.05), which is determined experimentally in-site or in laboratory. This showed that several (sub) categories in POS calculation can be described with several possible values, which is stochastic definition of an event.
Stochastic in numerical integration of structures
Stochastic is a property of numerical calculation of volumes of geological structures. This is the most often used method for calculation of (sub) surface structures. It is based on 2D approximation of 3D space. An object with volume V whose boundaries extend from x = a to x = b is defined by the definite integral (Equation 2) as follows:
V = f A(x)dx.	(2)
J a
The area A(x) of the section made by a plane parallel to the YZ plane has to be known at every point. In practice, the integrand A is defined by the table of values of the definite integral. It can be obtained by several formulas, including the areas (isopachs) measured by the mechanical device called planimeter. Two such formulas are the most often used in geology -trapezoidal and Simpson's rule (e.g. [17-19,12]).
Trapezoidal rule derived name from the first approximation where any integral is represented along a straight line, which then with the interval [a,b] forms the trapezium (Figure 4a). Trapezoidal shape is derived from the first approximation, i.e. straight line, which then with interval [a,b] and axis X forms the trapezium. Hence, the general approximation of trapezoidal formula (Equation 3) depends on a number of segments (Figure 4b):
f
J n
f(x)dx = ^ f f(x)dx ¿=1 Xi-1
1 n
1Yj(xi-xi-i) 1=1
(3 3
fb
Simpson's rule is based on integral I f(x)dx
J a
using approximation of y = f(x) by a parabola (Figure 5a), i.e. polynomial of second degree that passes through the points (a, f(a)), (b,f(b)),
(c,f(c)), where c = — (a + b). The final equation (Figure 5b) is the following:
fb ^ h J f(x)dx « -
f(a) + f(b)
(5)
n-1
+ 2^/(x2i)+4^/(x2i-1) ¿=i ¿=1
[f(xi-i)+f(xi)].
This is valid for any number of equally distanced subintervals, i.e. sections or isopachs (i.e. subintervals+1), which have uniform partition a = xn < x. < ... < x . < x = b, i.e. subinterval
0 1	n-1 n '
b — a
h =-. For example, such a formula for
n
5 sections (or 4 subintervals) would be as (Equation 4) follows:
h,
= 2 + 2fll + 2fl2 + 2a3 + (4)
In practice, again for 5 sections (or 4 subintervals), such a formula becomes as given in the following:
h,
= 3 + 4fl1 + 2a2 + 4a3 + a4).	(6)
Application of these rules is not unique, but some strong recommendations are raised from experience. The Simpson's rule starts approximation with a higher polynomial (2nd order vs. 1st order in trapezium), and the approximation, with numerous sections, is always better. However, there is not strong definition of "nu-
merous sections". For 5 or less sections, almost always simpler trapezoidal rule will lead to better results. For more sections, the Simpson's rule has obvious advantages. Moreover, the Simpson's rule has two versions, resulting from practice. One version, for paired number of sub-intervals, has been proven mathematically and is given here. However also, there is version for even number of subintervals used in Croatian reservoir geology practice for decades ([20]). If in Equation 7 one more subinterval is added, then there will be even number of subintervals and the equation will be as follows:
= 3 + + 4fl2 + 2a3
(7)
+ 4a4 + 2a5) .
Discussion, recommendations and conclusions
Stochastic and deterministic models are closely entangled. In fact, any deterministic model, map or equitation, is only an accepted approximation of natural input dataset or "artificial certainty". Geological deterministic models could be different, but mathematically, they could be divided into (a) purely numerical, like calculation of probability of success or volumetric calculation of geological structure or (b) graphical outcome like different deterministic maps.
Regarding mapping, SIS is probably the most advanced simulation technique that uses original and indicator data for variogram calculation and mapping. Moreover, if Gaussian and indicator simulations are compared, indicator maps sometimes represent more uniform distribution, i.e. differences among realisations are not so large as in Gaussian ones (e.g. [5]). It is a result of variance of indicator variables, and consequently, indicator simulation gives more uniform distribution of cell values. Generally, if indicators are used, the larger number of cut-offs results in larger reduction of in-class "noise" [21]. Eventually, the main purpose of Gaussian simulation is mapping of real values, but the main intention of indicator simulations is probability mapping, i.e. mapping assuming that some cells will have values larger than cut-off. In both cases, it would partially remove the so called "bull's-eye" effect, a very strong feature sometimes observed in deterministic maps. Removal could be even stronger if indicators are used.
Stochastic is an inherent property that is also used for other calculations in numerical geology, such as probability of success and volumetric calculation, which are previously described. Introduction of stochastic in such calculations gives some degree of freedom in selection of categories or fine-tuning of geological models. Consequently, if stochastic is applied, the algorithms for POS calculation or volumetric calculation need to be theoretically well known. Any decision to introduce stochastic (or not)
Figure 6: Decision tree for introducing stochastic in geological mapping and numerical calculations, based on the type of method and number of input data
is based on the type of analyses and number of data (Figure 6).
As general recommendation, when to apply stochastic could be outlined (Figure 6) here:
1.	The strongest criterion is number of input data. Although each dataset includes uncertainties, they are the highest in smaller datasets. On contrary, in large datasets, such uncertainties could be easily and precisely calculated.
2.	Consequently, it means that in "medium-sized" datasets, stochastic could be described, but it does not play an important role in the analytical procedure. Such "medium-sized" datasets are still too small that un-
certainties cannot be precisely numerically calculated (almost as constant) and too large that the representative statistics cannot be calculated.
3.	This is why stochastic is recommended for "small" datasets with <15 inputs or for "large" ones with >40 points.
4.	The calculation of probability of success for any geological category deviates from such recommendations, because it is purely a numerical method, where for <15 points, the porosity cannot be stochastically mapped, as well as other subcategories cannot be reliably estimated with several solutions.
General recommendation for any kind of input dataset, regarding each of three analysed approaches, is clearly summarised in Figure 6, which represents "all-purpose" table that could be applied in all research that include stochastic in geology.
Acknowledgement
This study is financially supported by "Mathematical Methods in Geology II" given by the University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering in 2017 (no. M048xSP17).
References
[1]	Dubrule, O. (1998): Geostatistics in Petroleum Geology. American Association of Petroleum Geologists, Volume 38, Tulsa, 210 p.
[2]	Kelkar, M., Perez, G. (2002): Applied Geostatistics for Reservoir Characterization. Society of Petroleum Engineers, Richardson, 264 p.
[3]	Malvic, T. (2008): Primjena geostatistike u analizi geoloških podataka (Application of geostatistics in geological data analysis). University literature, INA Plc., Zagreb, 103 p. (in Croatian).
[4]	Malvic, T. (2008): Kriging, cokriging or stochastical simulations, and the choice between deterministic or sequential approaches. Geologia Croatica, 61(1), pp. 37-47.
[5]	Novak Zelenika, K., Malvic, T. (2011): Stochastic simulations of dependent geological variables in sandstone reservoirs of Neogene age: A case study of Kloštar Field, Sava Depression. Geologia Croatica, 64(2), pp. 173-183.
[6]	Novak Zelenika, K., Malvic, T. (2014): Utvrdivanje sekvencijskim indikatorskim metodama slabopro-pusnih litofacijesa kao vrste nekonvencionalnih ležišta ugljikovodika na primjeru polja Kloštar (Determination of low permeable lithofacies, as type of unconventional hydrocarbon reservoirs, using sequential indicator methods, case study from the Kloštar Field). Rudarsko-geološko-naftni zbornik, 28(1), pp. 23-38. (in Croatian with English abstract).
[7]	Novak Zelenika, K., Velic, J., Malvic, T. (2013): Local sediment sources and palaeoflow directions in Upper Miocene turbidites of the Pannonian Basin System
(Croatian part), based on mapping of reservoir properties. Geological Quarterly, 57(1), pp. 17-30.
[8]	Balic, D., Malvic, T. (2010): Ordinary Kriging as the most Appropriate Interpolation Method for Porosity in the Sava Depression Neogene Sandstone. Nafta-plin, 30(3), pp. 81-90.
[9]	Malvic, T., Rusan, I. (2009): Investment risk assessment of potential hydrocarbon discoveries in a mature basin. Case study from the Bjelovar Sub-Basin, Croatia. Oil, gas - European Magazine, 35(2), pp. 67-72.
[10]	Režic, M. (2016): Opci model za izračun geološke vje-rojatnosti novih otkriča plina na području Sjevernog Jadrana uz primjer plinskog polja Ika (General model for the calculation of geological probability associated with new gas discoveries in the Northern Adriatic with an example of the Ika gas field), Diploma Thesis. Zagreb: University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering; 48 p. (in Croatian with English summary).
[11]	Malvic, T., Velic, J., Režic, M. (2016): Geological probability calculation of new gas discoveries in wider area of Ivana and Ika Gas Fields, Northern Adriatic, Croatia. RMZ - Materials and Geoenvironment, 63(3), pp. 127-137.
[12]	Malvic, T., Rajic, R., Slavinic, P. & Novak Zelenika, K. (2014): Numerical integration in volume calculation of irregular anticlines. Rudarsko-geološko-naftni zbornik, 28(2), pp. 1-8.
[13]	Husanovic, E. & Malvic, T. (2014): Review of deterministic geostatistical mapping of Croatian hydrocarbon reservoirs and advantages of such approach (Pregled dosadašnjih determinističkih geostatističkih kartiranja ležišta ugljikovodika u Republici Hrvatskoj te prednosti takvoga pristupa). Nafta, 65, 1, 57-68.
[14]	Malvic, T. (2009): Stochastical approach in deterministic calculation of geological risk - theory and example (Stohastički pristup u determinističkom izračunu geološkoga rizika - teorija i primjer). Nafta, 60, 12, 651-662.
[15]	Malvic, T., Velic, J. (2015): Stochastically improved methodology for probability of success ('POS') calculation in hydrocarbon systems. RMZ - Materials and geoenvironment, 62(3), pp. 149-155.
[16]	Gaurina-Medimurec, N., Novak-Mavar, K. (2017): Depleted hydrocarbon reservoirs and CO2 injection wells -CO2 leakage assessment. Rudar-sko-geološko-naftnizbornik, 36, pp. 15-27.
[17]	Atkinson, K.E. (1989): An introduction to Numerical Analysis. 2nd ed., John Wiley and Sons: New York; 712 p.
[18]	Kevo, M. (1986): Numericka integracija (Numerical integration - Slovenian, issue in Serbian). Moj mikro, 2,7, pp. 25-28.
[19]	Quarteroni, A., Sacco, R., Saleri, F. (2000): Numerical Mathematics. Springer-Verlag: New York; 654 p.
[20]	Malvic, T. (2015): Upute za uporabu planimetra (Instructions to measure with planimeter). University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering (handbook), Zagreb, 20 p. (in Croatian).
[21] Deutsch, C.V., Journel, A.G.: GSLIB - Geostatistical Software Library and User's Guide.- 2nd edition, Oxford University Press: New York - Oxford, 369 p.
Original scientific paper	Received: July 30, 2017
Accepted: October 08, 2017 DOI 10.2478/rmzmag-2018-0005
Geostatistical modeling of porosity data in 'oba' field, onshore Niger Delta
Oluseun Adetola Sanuade1*, Akindeji Opeyemi Fajana1, Abayomi Adesola Olaojo2, Kehinde David Oyeyemi3 and Joel Olayide Amosun1
1	Department of Geophysics, Federal University Oye-Ekiti, Ekiti State, Nigeria.
2	Earth Sciences Department, Ajayi Crowther University, Oyo, Nigeria
3	Department of Industrial Physics, Covenant University, Ota, Ogun State, Nigeria
* Corresponding Author: Oluseun Adetola Sanuade Department of Geophysics, Federal University Oye-Ekiti sheunsky@gmail.com +2347038958998
Abstract
A geostatistical approach was used to model porosity of OBA field in onshore area of Niger Delta using simulation technique. The objective is to understand the spatial distribution of porosity and characterize the degree of heterogeneity of underlying formation. Porosity data from twenty-two wells were loaded into SGeMS software. Univariate statistical analysis, experimental semivariogram and Sequential Gaussian Simulation (SGS) were applied on the data. The data was close to normal approximation of Gaussian based of the results of univariate statistics. However, to construct and model horizontal and vertical semivariograms, the data was log-normalized to reduce the coefficient of variation and to get good fit of the model. Parametric semivariogram model shows the range of 72-6480 m, nugget effect of 0.006 and sills of 0.0095, 0.0099 and 0.0111. Six realizations were generated using SGS algorithm and the results suggest that any one of the realizations can independently represents the true picture of the subsurface geology within the study area. Ranking of realizations shows realization 6 as the best and realization 2 as the lowest. This model could be used as an initial condition for simulation of flow.
Key words: geostatistics; stochastic; simulation; porosity; semivariogram
Introduction
One of the common problem in earth sciences is the delineation of subsurface properties from sparse and limited number of wells. Such problem was addressed early via contouring which is simple and gives clue about the trend of the variables but the approach is not unique in its result [1]. Later on this problem was addressed through kriging which is robust, and considered powerful in estimation of the weight and value of the regionalized variable on unsam-pled locations, and the outcome of this estimations is a deterministic model [1]. However, due to the complexity in reservoir properties that vary spatially, and sparseness of information, it is not logical to assign only one value for regionalized variable to construct unique deterministic model of such variable for extremely variable properties. As a result of this situation, stochastic representation of regionalized variable is the best, which follows the concept that any regionalized variable exists as a random variable with specific probability density associated to it [2, 3].
Most geological phenomena are extraordinarily complex in their interrelationship and vary widely in geographical extension, and exact description of the geological feature is neither feasible nor economically possible, therefore the results are uncertain in most cases. Such geological feature can be described by stochastic models that give different possible values of realizations with acceptable measurement of error. A simulation method is more sophisticated than a kriging process in that that it allows the user not only to specify statistical anisotropy in terms of semivariogram parameters as kriging does, but also to model heterogeneity by adding a random factor. Honoring spatial distribution and real value of the measured location is corner stone in reducing risk and uncertainty, and therefore stochastic simulation is the best approach to address such cases. With stochastic simulation, the geoscientist is better positioned to evaluate which geological information is relevant [4]. The ability to understand and predict the possible spatial distribution of a property with uncertainty is critical for understanding geological heterogeneity such as grain size, li-thology, texture, porosity, permeability, and
diagenetic processes. Good prediction of these properties is a sine qua non to proper decision making, planning and management strategy. Therefore, in this study, conditional sequential Gaussian simulation (SGS], one of the stochastic methods, is used on porosity data to characterize and measure geological feature with uncertainty.
Description and Geology of the study area
The Niger Delta is situated at the apex of the Gulf of Guinea on the west coast of Africa (Figure 1A]. The 'OBA' field is located in the onshore area of Niger-delta as seen in Figure 1B. The clastic wedge of the Niger Delta occurs along a failed arm of a triple junction systems which was formerly emanated during the breakup that occurred between South American and Africa plates. This process happened in the late Jurassic [6]. Synrift sediments accumulated during the Cretaceous to Tertiary, with the oldest dated sediments of the Albian age. The thickest successions of synrift marine and marginal marine clastics and carbonates were deposited in a series of transgressive and regressive phases [5]. The Niger Delta clastic wedge prograded into the Gulf of Guinea at increased rate that remained steady so as to respond to the evolution of these drainage areas and continuation of basement subsidence. The movement of marine shales of the Akata formation, which are deep-seated, over-pressured and ductile, within the basin, led to the formation of normal faults in the Niger Delta. The Akata shale is believed to have deformed the clastic wedge of the Niger Delta [5]. Majority of those normal faults are syndepositional and produced at the time of progradation of the Delta. The faulting styles produced in the Niger Delta are simple rollover structures with clay filled channel, growth faults, antithetic fault and collapsed crest (Figure 2].
There are three major lithostratigraphic units in the Niger Delta. These are Akata, Agbada and Benin formations (Figure 3]. Their depositional environments include marine, deltaic to fluvial environments [8].
Figure 1: (A) Location of Niger Delta (B) Location of OBA field (modified from [5]).
Figure 2: Niger Delta oil field structures and their associated traps [7].
The thickness of Akata formation is about 6,400 m at the centre of the clastic wedge. The lithology of the formation are dark gray shales and silts with streaks of sand (having origin tur-
bidite flow]. The age of this Akata is from Pa-leocene to Recent. It grades vertically into the overlying Agbada formation [5].
Figure 3: Stratigraphie column showing Formations in the Niger Delta [10].
Agbada Formation is about 3960 m thick having intercalations of sand and shale. These sands and shales are formed in fluvial-deltaic environment. The age of Agbada Formation ranges from Eocene to Pleistocene. Benin Formation is on the top of the clastic wedge of Niger Delta. The top of this formation is made up of the recent subaerially-exposed delta top surface. The shallow part of this formation consist of non-marine sands deposited in either upper coastal plain or alluvial deposi-tional environments [5]. The age is from Oligocene to Recent [9].
The main oil reservoir in the Niger Delta is the Agbada Formation [5]. The source rock in the Niger Delta is the marine shale of Akata Formation and/or marine interbedded shale of the Agbada Formation [5, 11]. The primary seal rocks are the interbedded shale of the Agbada Formation. However, the Agbada Formation while suitable for petroleum accumulation, is too deep to be relevant to groundwater storage. There arises therefore the major difference between the region where the petroleum geol-
ogist is prospecting for oil, that is, the Agbada Formation, and that, where the hydrogeologist I searching for water - the Benin Formation, in the Niger Delta.
Materials and Methodology
The data consists of about 2396 porosity measurements from 22 wells distributed in the southwest (SW)-northeast (NE) direction in the area. Sampling was done based on the flow direction of the aquifer in the area. The spacing between the well is irregular (Figure 4). Porosity varies between 0.10 and 0.30. The porosity data was estimated from sonic logs using equation 1. The depth of the wells varies from 120 to 500 m and penetrated the Benin Formation
At - Atma Atf - Atma
(1)
where: 0 is the porosity, At is the zone transient time (^s/ft), Atf is the transient time of fluid (189 p.s/ft for water), Atma is the average time of sandstone formation (52.6 [is/ft).
298400 ' 300000 ' 301600 ' 303200 ' 304800
o^^Eirïï^^Sii
Om 5000m 10000m 15000m 20000m
Figure 4: Base map showing well locations.
The data format was changed from excel format to ASCII format and loaded to Minitab 17 well by well to generate various graphical displays of the data [12]. These graphs include histogram, probability density function and cumula-
tive density function. Minitab 17 software was also used to compute statistics summary and perform normality tests. Golden Surfer 9 was used to develop base map of the wells in the study area.
Descriptive statistics was used to interpret different graph categories and to compare results. Then the data was displayed as box plot to check for possible outliers. Our data was found to be approximately close to Gaussian approximation of normal distribution. Semivariogram plots were generated by calculating variogram at different lags and azimuths using SGeMS software. Then experimental semivariograms were fitted with theoretical model (spherical model]. Parameters derived from fitted model were used to identify anisotropy and structural interpretation.
The study area was then gridded, and parameters from semivariogram were used to krige the study area. Afterward the SGS algorithm was applied on the data. This algorithm simulates nodes on a grid in random sequence by first estimating the value at the selected node by kriging with a local neighborhood of conditioning data, and then adding a random component from a normal distribution with zero mean and the kriging standard deviation. After simulation, values are added to the conditioning set for use in simulating additional nodes.
Results and Discussion
Descriptive Statistics
The statistical results of the wells show that the porosity distribution ranges from 0.10-0.30 with a mean of 0.25. This porosity range is typical of sand/sandstone or shale. However, the mean value indicates that it could represent all the data i.e. normally distributed porosity is equal to the mean or higher than mean. Analysis of quartiles for the porosity indicates that 25 % of the porosity data falls below 0.24 and 75 % falls below 0.30. The coefficient of variation (CV), standard deviation and variance of the porosity are 0.21, 0.06 and 0.00349 respectively. Coefficients of kurtosis and skewness for the wells are -0.33 and -0.27 respectively. The value of coefficient of kurtosis indicate that the porosity distribution is platykurtic in shape
Empirical CDF of Porosity
Normal
100
80
„ 60 c S u
0.0	0.1	0.2	0.3	0.4	0.5
Porosity
Histogram of Porosity Normal
200 m
0.06 0.12 0.1S 0.24 0.30 0.36 0.42 0.48 Porosity
Figure 5a: CDF and PDF of wells in the study area.
having flat top which results from large variations within observations. This means that most of the porosity values are less clustered around the mean with a fairly layout uniform of data. This may indicate local variation in terms of geology which may be as a result of faulting/ fracturing in the area [3]. Figure 5a shows the Cumulative Distribution Function (CDF) and Probability Density Function (PDF) for the wells. The CDF and PDF reveal normally distributed porosity data. Figure 5b shows the outliers tests performed on the wells while Figure 5c shows the graphical display of two normality tests for the wells. It was revealed from these tests that the porosity is normally distributed with a very minute outliers. The outlier porosity values were filtered out.
Figure 5b: Outlier tests for the wells.
Figure 5c: Anderson-Darling and Kolmogorov-Smirnov normality tests for the wells.
Experimental Semivariogram Semivariogram of the data shows structural variability and continuity at major direction azimuth 45, minor direction azimuth 135, and vertical direction with range of 6480, 1750 and 72 m respectively (Figure 6a-c; Table 1). Based on these range of values, the porosity in the study area exhibits anisotropic geometry in all directions. Also, the structural variances (C) of the horizontal direction (0.0039 and 0.0051) are greater than the vertical direction (0.0035). This suggests that porosity in the horizontal direction shows more variability than vertical direction because the space between two points in the horizontal direction (near 500 m or more) is greater than vertical direction (less than 30 m). This might indicate change in facies, intensity in cementation, or secondary
porosity that could be caused by fracturing in horizontal direction. The small nugget value shows that there is a small scale variation in the subsurface structure that might be a result of higher degree of cementation or fracture in the subsurface.
The range and nugget effect imply that the length of the spatial autocorrelation is longer than the sampling interval of 72 m in vertical direction. Therefore, the sampling design is good for this study and it is expected that a good spatial structure will be shown on the interpolated map in further study.
The differences in range suggests the elliptic shape of the semivariogram structure. The sills of parametric semivariogram model are 0.0095 in vertical direction and 0.0099 and 0.0111 in horizontal directions (theoretically
Figure 6a: Experimental andsemivariogram model at azimuth 45 (Major direction).
Figure 6b: Experimental and semivariogram model at azimuth 135 (Minor direction).
Figure 6c: Vertical experimental and model semivariogram.
Table 1: Semivariogram parameters. Azimuth	Nugget
Sill
Range(m)
Nugget/Sill ratio (%)
Model Type
45 (Major Direction)	0.006	0.0099	6480	60.6	Spherical
135 (Minor Direction)	0.006	0.0111	1750	54.1	Spherical
Vertical	0.006	0.0095	72	63.2	Spherical
Table 2: Spatial Dependence of Variables [14]. Nugget/Sill Ratio (%)	Inference
<25%
Strong Spatial Dependence
25%-75%
Moderate Spatial Dependence
>75%
Low Spatial Dependence
equal to the variance of the data, [13]. This suggests that the spatial structure exhibits geometrical anisotropy in porosity distribution as seen in Figure 6a-c.
As observed in Table 2, the nugget/sill ratio of porosity for all the semivariogram models are moderate based on the classification by Wei et al., 2007. This suggests that the porosity in the study area has a moderate spatial dependence and as a result, local variations within the study area might be captured.
Kriging
The kriging models (Figure 7] show variability in porosity distribution both laterally and vertically which suggests the heterogeneity of underlying features within the study area. It also reveals interbedded horizontal layering of geologic materials. The kriging variance map (Figure 8] shows principally high errors or uncertainty outside the zone of influence. This may be attributed to data quality or lack of data outside zone of influence
Figure 7: Kriging map.
Figure 8: Kriging variance.
Simulation Realizations
Figures 9-11 show six realizations generated using conditional SGS algorithm. Visual observation shows that the models are similar in terms of variability. The similarity of statistics of the models (Table 3) suggests that each of the realizations can independently represent the true picture of the subsurface geology within the study area. The outputs are set of probabilistic models, which can serve as a measure of uncertainty in predicting porosity distribution within the study area.
Generally, they all show variability in porosity distribution both laterally and vertically, which suggests the heterogeneity of underlying features within the study area. The statistical summary (Table 3) is very similar in all respects. The mean porosity of the six realizations ranges from 0.24 to 0.25, which is very close to the mean of the real data (0.25). The coefficient of variation of the models ranges from 0.067 to 0.081, while that of the real data is 0.081.
Q-Q plot (Figure 12a-f) is a confirmatory test for normality of the generated realizations, which clearly shows positive slope, and alignment of the data along a straight line. Based on the summary statistics and Q-Q plots of the realizations, we ranked the realizations (Table 4). Realization 6 among others appears most true with less uncertainty to represent
Figure 11: Realizations 5 and 6.
the subsurface geology of the study area. This model could be incorporated in any data analysis of the field in order to effectively improve the prediction of porosity and reduce uncertainty. Also, the porosity model could be used as an initial condition for simulation of flow in the field.
Table 3: Statistics summary of realizations
Realization No	Maximum	Median	Minimum	Mean	Standard deviation	CV	Variance
1	0.28	0.24	0.22	0.25	0.041	0.0739	0.0017
2	0.29	0.23	0.22	0.24	0.034	0.064	0.0012
3	0.28	0.24	0.22	0.25	0.038	0.068	0.0014
4	0.29	0.25	0.22	0.24	0.036	0.067	0.0013
5	0.28	0.25	0.22	0.25	0.044	0.079	0.0019
6	0.28	0.25	0.22	0.25	0.045	0.081	0.002
Real Data	0.28	0.25	0.22	0.25	0.099	0.081	0.002
Figure 12: Q-Q Plot of (a) realization 1 (b) realization 2.
(c)
(d)
Porosity
Figure 12: Q-Q Plot of (c) realization 3 (d) realization 4 (e) realization 5.
(e)
uiu
m
<u
lo.
o c
,1
">0.15
					■
				m	
					
					
		m			
					
					
m					
(f)
Porosity
Figure 12: Q-Q Plot of (f) realization 6.
Table 4: Ranking of porosity models. Rank
Realization No
Conclusions
with the result of the semivariogram. However, kriging variance model indicates high value of error outside the zone of influence which may be due to insufficient number of the well data. Six realizations were generated and the results indicates that any one of the realizations can independently represent the true picture of the subsurface geology within the study area with realization 6 ranked as the best and realization 2 as the lowest. We therefore recommend that geostatistical estimation and simulation to be incorporated in any geologic data analysis and that integration of more than one variable from multiple sources will effectively improve the prediction and reduce the uncertainty.
We have used a geostatistical technique to model porosity data in 'OBA' field, onshore Niger Delta. The porosity range in the area is 0.10-0.30 after filtering of outliers, which suggests that the lithology could be sand or shale. The standard deviation for the entire field ranges from 0.04-0.08, coefficient of variation ranges from 0.15-0.28. Based on this, linear and parametric geostatistics were employed to process the data in order to build geostatistical model of the subsurface geology. The semivariogram shows the major direction of continuity "azimuth of 45", spherical geometrical anisot-ropy. Kriging map shows clearly vertical and horizontal heterogeneity of the study area and subsurface interbedded layers, which agrees
References
[1]	Ming-Sheng, Y., Yu-Pin, L., Liang-Cheng, C. (2006): Designing an optimal multivariate geostatistical groundwater quality monitoring network using factorial kriging and genetic algorithms. Environ Geol 50, pp.101-121.
[2]	Al-Kuisi, M., Al-Qinna, M, Margane, A, Aljazzar, T. (2009): Spatial assessment of salinity and nitrate pollution in Amman Zarqa basin: a case study. Environ Earth Sci 59, pp.117-129.
[3]	Sahebjalal, E. (2012): Application of geostatistical analysis for evaluating variation in groundwater characteristics. World Appl Sci J 18(1), pp.135-141.
[4]	Karatas, B.S., Camoglu G., Olgen M. (2013): Spatio-temporal analysis of the depth and salinity of the groundwater using geostatistics integrated with GIS, the Menemen Irrigation System, Western Turkey. Int J Environ (Ekoloji) 22(86), pp.36-47.
[5]	Doust, H., Omatsola, M. (1990): Petroleum geology of the Niger delta. Geological Society, London, Special Publications, 50, pp.365-365.
[6]	Whiteman, A (1982): Nigeria: its Petroleum Geology, Resources and Potential. Graham and Trotman, London.
[7]	Stacher, P. (1995): Present understanding of the Niger delta hydrocarbon habitat: Geology of Deltas. AA Balkema, Rotterdam, pp 257-267
[8]	Weber, K.J. (1987): Hydrocarbon distribution pattern in Nigerian growth fault structures controlled by structural style and stratigraphy. J Petrol Sci Eng 1, pp.1-12.
[9]	Short, K., Stauble, A. (1967): Outline of geology of Niger delta. AAPG Bull 51, pp.761-779.
[10]	Tuttle, M.L., Charpentier, R.R., Brownfield, M.E. (1999): The Niger delta petroleum system: Niger delta province, Nigeria, Cameroon, and Equatorial Guinea, Africa: US Department of the Interior, US Geological Survey
[11]	Ekweozor, C.M., Okoye, N.V. (1980): Petroleum source-bed evaluation of Tertiary Niger Delta. AAPG Bull 64, pp.1251-1259.
[12]	https://www.minitab.com/uploadedFiles/ Documents/getting-started/Minitab17_ GettingStarted-en.pdf.
[13]	Yarus, J.M. and Chambers, R.L. (2006): stochastic modeling and geostatistics: Principles, methods and case studies. AAPG Computer Applications in Geology. 3, pp. 27-36.
[14]	Wei, H., Dai, L. and Wang, L. (2007): Spatial distribution and risk assessment of radionuclides in soils around a coal-fired power plants: a case study from the city of Baoji China. Environ Res. 104, pp.201-208.
Original scientific paper	Received: Sep 30, 2017
Accepted: Nov 26, 2017 DOI 10.1515/rmzmag-2018-0002
Application of analytic element method in hydrogeology Uporaba metode analitičnih elementov v hidrogeologiji
Dragan Kaluderovič1,2, Eva Koren3, Goran Vižintin3-*
1	Center for applied scientifical research in water supply and remediation, Serbia
2	Independent consultant, Serbia
3	University of Ljubljana, Faculty of Natural Sciences and Engineering, Slovenia * goran.vizintin@guest.arnes.si
Abstract
The analytic element method (AEM) has been successfully used in practice worldwide for many years. This method provides the possibility of fast preliminary quantitative analysis of the hydrogeological systems or boundary conditions of the numerical models, as it is shown in the case study of groundwater source of the city of Vrbas. The AEM is also applicable for the initial analysis of a hydrogeological system, which is of particular importance in case of excess pollution that cannot be predicted where it could happen. One example of the application of the AEM is presented in this article. The analytical model is calibrated based on the measured data from several drilled monitoring wells, and this was the base for the numerical model of the contaminant transport. In this case, the AEM enabled the quick access to information on the hydrogeological system and effective response to excess pollution.
Key words: analytical modelling, analytic elements, groundwater flow modelling, mass transport modelling
Povzetek
Metoda analitičnih elementov se že več let uspešno uporablja v praksi po vsem svetu. Kot je prikazano v študiji primera podzemne vode mesta Vrbas, omogoča ta metoda možnost hitre predhodne kvantitativne analize hidrogeoloških sistemov ali mejnih pogojev za numerične modele. Metoda analitičnih elementov se uporablja tudi za začetno analizo hidrogeološkega sistema, kar je še posebej pomembno v primeru onesnaženj, ki jih ni mogoče predhodno napovedati. V tem prispevku je predstavljen primer uporabe metode analitičnih elementov. Analitični model je kalibriran na podlagi izmerjenih podatkov iz več piezometrov, kar je osnova za numerični model transporta onesnaževalcev. V tem primeru je metoda analitičnih elementov omogočila hiter dostop do informacij o hidrogeološkem sistemu in učinkovit odziv na prekomerno onesnaženje.
Ključne besede: Analitično modeliranje, analitični elementi, model toka podzemne vode, model masnega transporta
Introduction
For 3D flow, the equation is
For description of groundwater flow and distribution of the pollution, we have used Darcy's equation and continuity equation. In 1856, the French engineer Henry Darcy conducted research in Dijon. Initially, it was the project of water supplying system that used sand filters. Based on this research and column experiments in the laboratory, he established Darcy's law that describes water flow through sands. It has since been generalised to a variety of situations and is in widespread use today. For the one-dimensional flow, Darcy's law is represented in the form of equation as:
d / ah\ d
-r[Kx — )+-r
dx\ ax) dy
d i ah\ + ~dz \y Hz) = ^
K 3h\ y ay)
ah
at
For the homogenous and isotropic porous media (K = K = K), the equation is
a2 h
dX2
+
a2h ay2
+
a2 h
ax2
Ss ah ~K~at ■
(5)
dh
q = -Kx—A Hx	dx
(1] If the flow is stationary, equation (5] transforms to
where qx is the flux in the positive 'x' direction (discharge per unit area with units of length per time, m/s), Kx the hydraulic conductivity, dh/dx the hydraulic gradient and A the cross-sectional area to flow. The three-dimensional (3D) flow could also be described using Darcy's law. Dar-cy's law is valid only for flow in the continuum region where the representative elementary volume (REV) is defined and for laminar flow through the soil. In most of the aquifers, the dimensions of interstices are small and thus the flow is laminar. The other equation that we use to describe groundwater flow is the continuity equation (conservation of mass); the derivation of this equation is available in Fitts [1]. In the 'x' direction, the equation of conservation of mass is
Ah = 0.
(6)
aqx ah ' — Ss
3X	at
(2)
where qx is the specific flux (discharge per unit area, with units of length per time, m/s) in the positive 'x' direction and Ss is the specific yield. Using equation (1) in equation (2), we get the equation of the one-dimensional groundwater flow as
d
ah\
ah
T" ( Ky- ) — S«--
dx\ sx/
at
(3)
Equation (6) is known as Laplace's equation. Groundwater flow is often presented in hy-drogeological models as two dimensional in the horizontal plane, as it spreads wider horizontally than vertically. In this situation also, groundwater flux could be presented as two dimensional, where h varies with x and y axes but not with the z axis. The alternative formulation of the groundwater flow equation may be obtained by invoking the Dupuit-Forchheimer assumption, where it is assumed that heads do not vary in the vertical direction. A horizontal water balance is applied to a long vertical column with area SxSy extending from the aquifer base to the unsaturated surface. Originally, it is invented for aquifers with atmospheric pressure at the groundwater table, but it could also be used in the case of the artesian aquifers. Derivation of the equation could be found in Fitts [1], and the resulting equation is
d I dh\ d I dh\	dh
dx\ x dX) dy\y dy)	s dt
(7)
N represents the addition of water in the vertical direction (recharge). If Tx = Ty, the above-mentioned equation simplifies to the following form:
/g2h\ (a2h\ N Sah
. 8X2
Tat
Without recharge/leakage, Laplace's equation (8) is
dfy ah	ah
dx ax ' dy ay
AO = 0.
For the stationary groundwater flow, equation (8) transforms into a new form, which is Poisson's equation [1]:
N
Ah = — .	(9)
Background and methods
The analytic element method (AEM) [2, 3] uses the technique of superimposing solutions in linear partial differential equations like (9) using the computer, which superimposes a large number of functions. Each analytical element is a mathematical function that is connected with a boundary condition, where an element represents a well, a stream segment or an area of surface recharge. The AEM can be used in one-dimensional, two-dimensional or 3D flow. Equation (7) represents a two-dimensional flow in an isotropic transmissibility. In the AEM, this and other equations are written in form of discharge potential (O]. Discharge potential is defined using parameters and groundwater level, so the equations are correct:
(13)
Elementary solutions may be superimposed, which is the main principle of the AEM. This method is used in programme AquiferWin32 Version 5 (developed by the American company Environmental Simulations Inc.); the equations in this programme are based on complex functions of different analytical elements, and they have this form: Uniform flow:
O(x, y) = -Q0 (xcosau + ysinau) + C.
(14)
Effect of the well:
Y" T^rffcy)].
1'j=1%
(15)
Regional recharge:
1 2
N
a2 + b2
[(a2sin2ar + b2cos2ar)
(x - x) - 2(a2 - b2) (x - xr) (y - y)
sina cosa (a2cos2a + b2sin2a )
r	r ^	r	rJ
(16)
(10) (y -yr)2-a2b2].
Putting this equation into equation (7), we ob-	Source/pit - round element:
tain	Equation in the element:
ah	X-1" 1 r	i
^ = -V+Ss-.	(U)	-^4[(x-xPj)2-R2PjJNPj.	(17)
Using Laplace's transformation, it is possible to create semi-analytical solution of this non-steady state equation.
When the flow of groundwater is in the steady state with recharge/leakage, Laplace's equation is
Equation outside the element:
Zn i
j=i 4
Rpjln
(x-
■xpj)2 + (y + ypj)21
R2
N
Pj'
AO = N.
(12)
(18)
Source/spring - linear element:
n a L
j=i
(19)
— (Zj — 1)In(Zj — 1)2In[^ (Z" — ~0] — 2)
where x is the x coordinate of the point for which we calculate the piezometric level (L); y the x coordinate of the point for which we calculate the piezometric level (L); Q0 the volumetric flow rate per unit time/unit width of the aquifer (L2T-1); au the angle between uniform flow and the x axis, Qa the capacity of the well (L3T-1); r the distance between the well and the point for which we calculate the piezometric level (L); N the recharge and net infiltration (LT-1); a the length of a, axis of the ellipse (circle) of recharge (L); b the length of b, axis of the ellipse (circle) of recharge (L); xr the x coordinate of the centre of the ellipse (circle) of recharge (L); yr the y coordinate of the centre of the ellipse (circle) of recharge (L); ar the angle
between a and b axes; x . the x coordinate of the
' pi
centre of the circle element (L); ypj the y coordinate of the centre of the circle element (L); R .
v. pi
the radius of the circle element (L); N the in-
v. pi
filtration of the circular element (LT-1); cti the volumetric flow rate per unit length of the linear element (L2T-1); L the length of the linear element (L); zx the starting coordinate of the linear element (L) and z2 the final coordinate of the linear element (L).
Discussion and results
The application and use of the AEM is different than the application of the classic numerical method. Here, we have represented a few fundamental differences. This method is analytical, and there are no errors caused by the numerical approximation of the partial differential equation [4]. The solution derived with this method is not sensitive to the scale of the domain [5], which varies from few tens of metres to several tens of kilometres. The boundary conditions do not influence the accuracy of the solution, ex-
cept right at the boundary where the boundary condition approximation depends on the scale of elements and the degrees of freedom in each element. The biggest limitation of the analytic element method as used in AquiferWin32 is that it simulates only steady-state conditions. In the area of Vrbas city and below it, there are two aquifers, one shallow up to 70 m deep and the second deeper at a depth of approximately 120 m. In this model, the shallower aquifer was simulated; the reason for application of analytic element modelling is that the shallow aquifer was covering a huge area, more than 100 km, and therefore it would be very hard to get physical boundaries of the numerical model. In the AquiferWin32, the model was simulated as a single-layer, homogenous aquifer with fully penetrating elements. Set points of the elevation of the bedrock and the aquifer overlying units were averaged to 17 and 65 a.m.s.l., while the starting value of the hydraulic conductivity was 1.6 x 10-4 ms-1.
Measurements of head in July 1996 in ten piezometers were used to calibrate the model. Small precipitation is a characteristic for this month, as well as the month before and after, and it could be assumed that net infiltration is minimal or none. In this period, the first aquifer was tapped by several wells, a well on the pig farm (5 l/s), well B-10 nearby the existing source in Vrbas (12 l/s), a well at the gravel exploitation location (4 l/s) and a source in Kula (~0 l/s), which were represented as analytical elements of wells.
Based on the geotechnical research documentation, it is determined that bottom of the channel DTD is curved in aquifer sands, so the water level in the channel south from the Savino village is represented with two line elements with constant head 79.75 and 79.65 a. m. s. l. AquiferWin32 and many AEM programmes require input of a groundwater level at one point - that is the reference point, a point with a specified water level regardless of other conditions. The reference point was positioned in the piezometer PP-7, and the value of the groundwater level was 79.01 a.m.s.l., taken as the reference level. This piezometer is chosen as it is positioned far from the influence of well exploitation, which is one of the primary requirements for the selection of the reference point. The gen-
Figure 1: "Scatter" diagram of measured and calculated groundwater levels.
eral direction of groundwater flow is taken as 75°, measured in the clockwise direction relative to east, with a gradient of 1 x 10-4. These piezometers were used in model calibration: PP-2, PP-3, PP-5, PP-7, PP-9, PD-2, PD-3, Ekonomija, Precistac and Vojin Salas (Figure 1]. The first step of calibrating the model was setting the azimuth of the regional groundwater flow, and the best result was for the value of 63°. The other parameter that was modified during the calibration was the regional gradient, and the best result was for the gradient value of 2.9 x 10-4. Changing the elevation of the top and bottom of the aquifer did not show a big influence on the model results, so the input average values were kept. Minimal declination from measured values was for the hydraulic conductivity of 1.7 x 10-4 ms-1, and this value was used in the further simulation. Further improvement of model results came by including the evapotranspiration as a negative vertical balance. The groundwater level in the first aquifer is supposed to be very shallow, close to the ground surface (3 m deep], and the hydraulic connection between groundwa-ter and the "non-permeable" layer exists, so it is possible that the aquifer loses the water by evapotranspiration through this layer. These discharge zones are presented with circular area sink elements and determined based on the comparison of terrine elevation maps and groundwater level maps. Zones with shallower groundwater levels are represented with smaller circle elements, which is a way of discretisation. The smaller elements work with more intensity comparing to bigger circle ele-
Figure 2: The map of isolines of groundwater levels -stationary condition (June 1996) 1:200000.
ments. The process of adjusting the model was finished when measured levels in piezometers and calculated levels were acceptably close to each other (Figure 1]. The isolines of groundwater heads of this model are given in Figure 2. The absolute mean residual value is 0.014 m for all 10 calibration targets; the only bigger deviation is registered on the piezometer PD-2, while other piezometers had small residuals. The causes of the deviation on PD-2 are unknown but may be due to incorrect schematisation, the effects of transient flow in the model or measuring mistakes. The deviation values do not differ that much, so it could be said that the model is well calibrated.
The analytic elements in the simulation are as follows:
•	The channel south from the Savino village is presented as two line elements with a constant level, and these analytical elements participate with 11 l/s as a ground-water source.
•	Circular analytic element already represents the evapotranspiration, with a surface area of 1.38 x 108 m2; the discharge flow got by model calibration is 1.5 x 10-10 ms-1, and it participates in the balance with 20.7 l/s.
•	Circular analytic element, smaller but with a bigger capacity (intensity], represents stronger evapotranspiration (groundwater level closer to the surface]; the surface area is 2 x 107 m2, and the calibrated in-
tensity is 6 x 10-10 ms-1; it participates in the balance with 12 l/s.
• Wells, point elements, work with the capacities that were previously mentioned.
In Figure 3, the nature of the method is shown, and it could be seen that the model has no boundaries, so the solution is valid only for the research area.
Analysing sensitivity, for the purpose of uncertainty of the calibrated model, showed that the smallest differences are caused by changing the elevation values of top and bottom of the aquifer. These values were changed by ±50% of the values accepted at the beginning, and this range represents the upper and the lower elevation values, got by drilling in the research area, so this parameter does not influence much to the solution of this system.
Changes in the statistical calibrating parameters, during the change in the hydraulic conductivity, show bigger response of the results and nonlinearity, i.e. smaller values of the hydraulic conductivity, and cause bigger errors than big values. Input data related to the reference point showed different responses. Changing the azimuth values of the groundwater flow caused small influence to piezometer levels, probably because the levels are fixed by other analytical elements. Sensitivity of the model resulting to changes in the gradient values of the reference point was noticeably higher than that in case of other parameters.
Changes in the evapotranspiration intensity (by changing only the values of bigger analytical elements, because changes in the smaller ones have only the local influence] had a relatively weak influence on the piezometer level. Instead of measuring this parameter, like other elements were measured, values obtained for this parameter were empirical, and values of the intensity were obtained by calibration. It could be concluded that values of the evapotranspiration are the most uncertain. It is important to mention that statistical parameters that indicate the success of calibration obtain their minimum values of the representative input data, and we can use these values for the verification of the success of the calibrating process.
Eile Edit View Add Options Caic Model Window Help	-|a| x|
Figure 3: Display of the domain without boundaries in which one part is the research area.
Evaluation of the impact of oil spill on the environment and aquifer, conducted in two phases
Owing to intentional damage of the pipeline caused by an unidentified person at the chainage 14.5 of the Novi Sad (Pancevo's pipeline], crude oil leaked from the pipeline to the environment. The estimated crude oil leakage was ~54,596 kg. The contaminated site is located on the cadastral plot number 4248 of the local municipality of Kovilj, which is the property of "Vojvodinaput - Backaput" enterprise from Novi Sad. This location was used as the sand excavation site for constructing the Kovilj traffic loop on the Belgrade-Novi Sad highway. The depth of the depression was 3-5 m, which made the accident location difficult to access. The surveys were done in order to contour the pollution, and the mathematical model of the contaminant transport was made based on those surveys.
In general, there are two methods to solve the problem of groundwater flow and transport of the pollution in them: analytical method and numerical method. The main limitation of analytical models is their ability to solve only simple problems; obtaining very precise solutions is their main advantage. However, numerical models are used for simulations in complex surroundings. Numerical models make it possible to replace continual forms of partial differential equations with the final number of algebraic equations, and this procedure provides
Table 1: Measured and computed groundwater heads.
Target head (m)	Computed head (m)	Residual (m)	Well name
75.730003	75.786697	0.056694	PZ-1
75.900002	75.872733	0.027269	PZ-2
75.839996	75.876975	0.036979	PZ-3
76.070000	76.068621	0.001379	PZ-4
76.019997	76.040388	0.020392	PZ-5
the simulations of the systems with irregular special and temporal characteristics.
Mathematical modelling using the AEM - the first phase in modelling
As the first step in mathematical modelling, the AquiferWIN32 Version 5 programme, developed by the North American Company Environmental Simulations Inc., was used. This programme enables simulation of groundwater flow using the AEM in the 2D horizontal plane. The basis of this method is superposition of influences of analytical elements of different shapes and function of regional gradient of groundwater. Although limited to steady-state flow, this method provides a fast and efficient analysis of a hydrogeologic system, which is of particular importance in the early stage of investigation when data are scarce. The programme simulates effects of the following analytical elements: wells, uniform recharge of the entire model surface, circular elements, recharging/discharging wells and linear elements that can be either sinks or sources of water. The number of elements in the model is not limited and depends solely on the computer memory.
Entry data
On the basis of lithological profile that was made on six drilled piezometers at top and bottom of the layer, average input points of 77 and 57 mnm were chosen. Although it is clear that top and bottom of the layer vary (the bottom point is not positively determined due to shallow drills, ~5 m), the limitations of this meth-
Figure 4: The isolines of groundwater heads in the research area of the calibrated model.
od demanded the use of constant values. Initial hydraulic conductivity K taken for this model is 1 x 10-4 m/s, and it was obtained from the data in the literature. It is necessary to mention that sieve analysis gave considerably smaller hydraulic conductivity value and that, determined in this way, it is on the border of certainty. Groundwater head measurements in piezometers (5), done throughout November 2013 (Table 1), served for calibration of the model. In AEM, it is necessary to know the value of piezo-metric heads at one point - the reference point, which serves as a starting value for calculation in which the level remains constant throughout the simulation. The reference point is put into the piezometer PZ-4, and the measured piezo-metric level is 76.07 mnm. By changing the values of the regional gradient, groundwater flow direction and the mentioned parameters, the results were obtained with minor residuals between the measured and computed levels of groundwater (Figure 4). In the course of calibration, the model demonstrated the highest level of sensitivity to the
Table 2: Measured and computed groundwater heads, December 2013. Target head (m) Computed head (m)		Residual (m)	Well name
75.66	75.821607	-0.16160	PZ-1
75.82	75.867556	-0.04755	PZ-2
75.74	75.856221	-0.11622	PZ-3
76.00	76.000148	-0.00014	PZ-4
75.95	76.002048	-0.05204	PZ-5
75.62	75.529781	+0.09021	PZ-6
75.82	75.895146	-0.07514	PZ-7
75.70	75.684449	+0.01555	PZ-8
Figure 5: The isolines of the groundwater heads in the research area of the calibrated model, December 2013.
changes in the values of the regional gradient of the groundwater from 0.01 to 0.001, while the correspondent absolute residual mean (ARM] values varied from 0.67 to 0.042. The gradient value of 0.009 was selected at the end of calibration because it produced the smallest error. Besides, the azimuth of the groundwa-ter flow was changed; 40° was selected, and it demonstrated less sensitivity than the regional gradient.
As a confirmation of the analytic model results, another round of measuring was carried out in December 2013 on eight piezometers (Table 2 and Figure 5], and the model was calibrated. The results were slightly different; the gradient was 0.006, while azimuth of the groundwater flow was 50° (Figure 5]. The results are logical and comply to the previous measuring (Table 2].
Comment on the model results
Although it is obvious that this is a simplified simulation, the calibration results are good. All eight piezometers were well calibrated, which is an excellent result. The results obtained in this way will be used as a basis to come to the boundary conditions of the 2D numerical model of contaminant transport where the spreading of pollution within the aquifer is to be researched.
Numerical model of the contaminant transport in the research area
The boundary conditions from the previous model, obtained using AEM with the help of AquiferWin32 programme, were used as a basis for the 2D numerical model of contaminant transport. The model was practically copied, only this time to the numerical model. The programme used for simulation was MT3DMS [6] and graphic interface was Groundwater Vistas, Version 6 [7], developed by the North American Company Environmental Simulations Inc. The formed 2D numerical model had the same parameters as the model done using AEM, but now discretisation was done with 183 x 141 columns and rows (Figure 6]. In the part necessitating higher precision, the cell size was 2 x 2 m and the total number of cells was 25803. For the values of dispersivity, the following values were given: 5 m for longitudi-
t 10!MM*W»<M<aFrM»«Î0«»C!!MÎBW«0«8®»S!T(l 174 no
200 meteres
Figure 6: Presentation and discretisation of the model.
nal dispersivity and 0.5 m for transversal dis-persivity. These values have been taken from the literature. The taken value for porosity was h = 0.2. The simulation (Figure 7) was done using the literature data for benzene. Benzene was taken as a very soluble contaminant difficult to absorb, and therefore, it is a good choice to simulate. In this simulation, benzene adsorption coefficient (Kd) is 1.52 l/mg and volume density is equal to 1600 kg/m3.
The value for the distribution coefficient K,
d
was taken from the literature [8-10]. Half-life biodegradation value was also included in the simulation to simulate biodegradation, and it was t1/2 = 407.5 days, also taken from the literature [1l, 12]. After 10 years of simulation, the benzene spread out at ~300 m (Figure 7) from the source of contamination, at the place where oil was spilt. The well near the threshing barn was not endangered, which was one of the goals of this simulation - to estimate an environmental impact to well from oil spill.
Conclusions
To estimate the spreading of the pollution dissolved in the groundwater, a simple 2D model was constructed. Naturally, more data would be needed for more precise results (such as monitoring of the groundwater level in the course of one hydrologic year) although the results with parameters at the boundary of certainty demonstrate that receptors are not endangered
Figure 7: Surface of the model that is covered with concentrations of up to 0.005 mg/l after 10 years of simulation with advection, dispersion, adsorption and biodegradation. The characteristic contaminant is benzene.
and that the application of the natural remediation of the aquifer is possible. Use and application of the AEMs significantly differ from the application of numerical methods. Several fundamental differences are mentioned here. The method is analytical, and discretisation of the domain of interest is not a necessity; just because of this, there are no errors due to the model resolution, i.e. due to the numerical approximation of the partial differential equation, which is used for the simulation of the groundwater flow. The solution acquired in this method is not sensitive to the scale of the model, which can range between a couple of dozens of metres and a couple of dozens of kilometres. The schematised aquifer simulated in this analytical method is unlimited in plan; boundary conditions cannot influence the accuracy of the solution. Various analytical elements enable simulation of relatively complex systems, and the solution of the entire system is obtained by superposition of all analytical elements. The biggest limitation of the AEM is its capacity to simulate only steady conditions of the groundwater flow. This shortcoming of the method also has some positive sides; for simulation in unsteady conditions, additional input data (specific yield, starting conditions) are needed, and at the same time, such simulations are considerably more complex to perform. The analytic element model is an ideal tool for researching in accidental pollutions when it is necessary to react swiftly and data are scarce.
References
[1]	Fitts, C. (2012): Groundwater Science, Second Edition, Elsevier: Amsterdam, 692 p.
[2]	Haitjema, H.M. (1995): Analytic element modeling of groundwater flow. Academic Press: San Diego, 407 p.
[3]	Strack, O.D.L. (1999): Principles of the analytic element method. Journal of Hydrology, 226(3), pp. 128-138.
[4]	Hsiao, C., Chang, L. (2002): Dynamic Optimal Ground-water Management with Inclusion of Fixed Costs. Journal of Water Resources Planning and Management, 128(1), pp. 57-65.
[5]	David, R.S., Bernard, E.A. (2006): The Synergistic Powers of AEM and GIS Geodatabase Models in Water Resources Studies. Groundwater, 44(1), pp. 56-61.
[6]	Zheng, C., Wang, P. (1999): MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User's Guide. U.S. Army Corps of Engineers, U.S. Army Engineer Research and Development Center; Vicksburg, Mississippi.
[7]	Environmental Simulation incorporated, 2014: Groundwater Vistas 6, USA.
[8]	Hunt, J.R., Nicholas, S., Udell, K.S. (1988): Non aqueous phase liquid transport and cleanup 1. Analysis of Mechanisms. Water Resources Research, 24(8), pp. 1247-1258.
[9]	Ling, M., Rifai, H.S. (2007): Modeling Natural Attenuation with Source Control at a Dry-Cleaner Site. Groundwater Monitoring Review and Remediation, 27(1), pp. 108-121.
[10]	Miller, C.T., Poirier-McNeill, M.M., Mayer, A.S. (1990): Dissolution of trapped nonaqueous phase liquids: Mass transfer characteristics. Water Resources Research, 26(11), pp. 2783-2796.
[11]	Brauner, J.S., Widdowson, M.A. (2001): Numerical simulation of a natural attenuation experiment with pertoleum hydrocarbon NAPL source. Ground Water, 39(6), pp. 939-952.
[12]	Howard, P.H., Boethling, R.S., Jarvis, W.F. Meylan, W.M., Michalenko, E.M. (1991): Handbook of Environmental Degradation Rates. Lewis Publishers Inc.: Chelsea, Michigan, 725 p.
Original scientific paper	Received: Jun 25, 2017
Accepted: Nov 29, 2017 DOI 10.1515/rmzmag-2018-0004
Soil corrosivity and aquifer protective capacity of overburden units in Ado-Ekiti, southwestern Nigeria
Preiskava korozivnosti tal in sposobnosti površinskih plasti za varovanje vodonosnikov na ozemlju Ado-Ekiti v jz Nigeriji
Oyedele, A.A.1*, Olayanju, G.M.2, Talabi, A.O.3, Ogunyebi, S.N.4, Ojo, O.F.3
1	Department of Physics, Ekiti State University, Ado-Ekiti, Nigeria
2	Department of Applied Geophysics, Federal University of Technology, Akure, Nigeria
3	Department of Geology, Ekiti State University, Ado-Ekiti, Nigeria
4	Department of Mathematical Sciences, Ekiti State University, Ado-Ekiti, Nigeria * akinakin02@yahoo.com
Abstract
Vertical electrical sounding, well inventory and phys-icochemical analysis were conducted to evaluate soil corrosivity and aquifer protective capacity of overburden units in the basement complex terrain of Ado-Ek-iti, southwestern Nigeria. The topsoil is composed of slightly corrosive materials at the eastern, southern and northeastern flanks and the central portion with resistivity values ranging from 60 to 180 Qm. Moderately corrosive/slightly corrosive materials (with resistivity values of 10 < p < 60 Qm) constitute the second layer around the eastern, southern and northeastern flanks. Pockets of areas in the northwestern, southeastern, eastern and central parts of the metropolis are practically non-corrosive with resistivity values in excess of 200 Qm. Zones of good, moderate, weak and poor overburden protective capacity were delineated, with longitudinal conductance (S) values of 0.7 < S < 4.9, 0.2 < S < 0.69, 0.1 < S < 0.19 and S < 0.1 mhos, respectively. On a regional consideration, 23.31%, 18.80% and 57.9% of the study area is characterised by overburden materials of poor, weak and moderate protective capacity, respectively. Only 6.02% of the area indicates good overburden protective capacity.
Key words: Ado-Ekiti, Contaminant, Corrosivity, Geoe-lectrical Survey, Protective Capacity
Povzetek
Vertikalno električno sondiranje, preiskavo vodnjakov in fizikalnokemično analizo so opravili z namenom oceniti korozivnost tal in sposobnost površinskih plasti tal za varovanje vodonosnika na ozemlju kamnin podlage Addo-Ekiti v JZ Nigeriji. Vrhnja plast tal na vzhodnih, južnih in severovzhodnih delih in v osrednjem delu ozemlja sestoji iz nizko korozivnih sestavin z električno upornostjo med 60 in 180 Qm. Nizko do zmerno koro-zivne sestavine (z upornostjo 10 < p < 60 Qm) tvorijo drugo plast na vzhodnih, južnih in severovzhodnih delih terena. Krpe ozemlja v severozahodnih, jugovzhodnih in osrednjih delih prestolnice so pa praktično ne-korozivne z vrednostmi upornosti, večjimi od 200 Qm. Površinske cone tal z dobro, zmerno, nizko in neznatno varovalno sposobnostjo so omejili z ozirom na vrednosti vzdolžne prevodnosti (S), in sicer (0.7 < S < 4.9), (0.2 < S < 0.69), (0.1 < S < 0.19) in S < 0.1 mhos. Regionalno gledano pripada 23.31%, 18.80% in 57.9% preiskovanega ozemlja površinskim materialom neznatne, nizke in zmerne varovalne sposobnosti. Samo 6,02% površine ima površinske plasti z nakazano dobro varovalno sposobnostjo.
Ključne besede: Ado-Ekiti, onesnaženje, korozivnost, geoelektrična preiskava, varovalna sposobnost
Introduction
Performance of civil engineering construction works including utility pipes requires knowledge about the corrosivity of soil. Soil supports man-made structures of all kinds; utilities and infrastructure are buried in it. Buried pipes are susceptible to corrosion and subsequent failure if the host soil medium is corrosive and aggressive. Corrosion refers to the degradation of buried metallic materials whereby metallic substrates are converted into oxides, hydroxides and aqueous salts within a cathode-anode system with corresponding loss in strength, ductility and other mechanical properties. Corrosion of cast iron, ductile iron and steel in soils can lead to a wide range of failures that are often accompanied by a high degree of economic and environmental consequences. For instance, leakage or rupture in pipelines could constitute hazards to the environment. Mitigating measures during design and construction as well as an understanding of the corrosive potential in a given soil environment are thus desirable [1-3]. Soil corrosivity depends largely on the composition of the soil and other environmental factors such as the moisture content, presence of oxygen, pH value, content of dissolved salts and porosity (aeration). The presence and abundance of oxygen promote corrosivity. High concentrations of soluble salts, high moisture content and a pH indicative of an acidic medium readily promote corrosivity. These factors control the soil resistivity. They are thus the main diagnostic factors. There exists a good correlation between the soil resistivity and corrosion rate of the buried metallic materials. Soil corrosivity is inversely related to the soil resistivity [2-4]. Contamination of the hydrogeologic system in metropolitan areas is increasing and has thus led to a critical issue in groundwater quality considerations. Installation of facilities, though essential but capable of provoking permanent damage of the underlying aquifers particularly in areas where residents rely mostly on groundwater, mandates an understanding of the aquifer protective capacity of the overburden units of the host soil medium. The overburden encompasses all geomaterials above the presumably fresh bedrock. Urbanisation and industrialisation remain the predominant con-
tributors of contaminants to the hydrological systems. Leachate from dumpsites, mining activities, buried petroleum pipes/tanks and septic tanks and the widespread use of chemical products such as pesticides, herbicides and solvents portend risks to the groundwater quality status [5-7].
The subsurface layers act as a natural filter to imposed surface pollutants. The ability of the geomaterials to retard and filter percolating fluid is a measure of the protective capacity and a function of transmissivity. Estimating these properties from the traditional methods of pumping tests can be very expensive and time consuming. The electrical resistivity method is significant in in situ determination of subsoil characteristics and conditions [7-9]. The traditional purpose of electrical resistivity survey is to determine the resistivity distribution of the subsurface by taking measurements on the ground surface. The electric conduction in the subsurface is essentially electrolytic through interstitial water in pores and fissures. Groundwater filling the pore spaces constitutes a natural electrolyte with a considerable amount of ions. Soil environments have requisite electrolytic properties for the redox (oxidation-reduction) reactions that take place during corrosion. Soil resistivity indicates the ability of a soil environment to carry corrosion currents [10, 11]. It is noted that burial of utilities and underground storage tanks are restricted to shallow depths.
In the present study, we seek to evaluate the corrosivity and protective capacity offered by the overburden units in the study area.
Geology and description of location
The study area, Ado-Ekiti, southwestern Nigeria, lies within latitudes 7° 32$ and 7° 42' N and longitudes 5° 9' and 5° 22' E (Figure 1). The area is underlain by the basement complex of southwestern Nigeria comprising the migma-tite-gneiss-quartzite complex, charnockitic and dioritic rocks, older granites and unmet-amorphosed dolerite dykes. Charnockite and quartzite ridges, which rise abruptly above the surrounding country rocks, are found in the area (Figure 2). Ado-Ekiti experiences a tropi-
MAPOFSTUDYAREA
5°10'30"E	5"120"E	6'13'30"E	S'lS'lTE	6'16,3G"E	S'lffCE	S'lffSG't
Figure 1: Location map showing the study area (adapted from Google Maps).
Figure 2: Geological map of Ado-Ekiti showing the VES points.
cal climate with distinct wet and dry seasons. The major rivers draining the area include Al-amoji, Ireje and Elemi [12, 13].
Materials and methods
Soil corrosivity is not a directly measurable parameter. One of the proven classifications of soil corrosivity towards the buried metallic materials is based on a single parameter, the soil resistivity. The soil resistivity reflects the main soil properties as it depends on porosity, degree of electrolyte saturation or concentration of dissolved salts in soils. It is thus typically indicative of soil corrosivity. The resistivity method has been adopted for this study [2, 14, 15]. Well inventory data were obtained from 107 hand-dug wells across the metropolis including physicochemical analysis to secure information on groundwater presence and composition, which readily influence the resistivity. Parameters such as electrical conductivity (EC), pH, temperature and total dissolved solids (TDS) were measured in situ in the field [16]. The geoelectric layer parameters were acquired using the Schlumberger arrays of vertical electrical sounding (VES) at 133 locations spread across the study area. ABEM Terrameter SAS 300B with ABEM 2000 booster was used for the measurement of resistance. The apparent resistivity values obtained were plotted against (AB/2) or half the spread length at each station on a bilogarithmic (log-log) paper. Partial curve matching of field curves with relevant Schlum-berger-developed master and auxiliary curves were used to obtain the resistivity values and corresponding thicknesses of the layers. The geoelectric parameters from this manual interpretation were improved upon by using the computer iteration algorithm RESIST Version 1.0 [17].
Depth sounding curves were inspected to determine the number and nature of the layering. The total longitudinal unit conductance, S, for the VES points was computed with the layer parameters as inputs [7, 9, 10]. Each VES station was georeferenced. ArcGIS 10.2.2 was used for the spatial distribution and thresholding of the data.
0therdatainputsincludedsoilmap3000/853/9-65 drawn by the Ministry of Agriculture and Natural Resources on 1:500,000 scale, geological map 3300/6/66/3289/0S prepared by the British Government's Ministry of Overseas Development on 1:250,000 scale and topographical sheet 1000/404/6.68 compiled and drawn by the Federal Survey from photo reduction 244 (Ado-Ekiti), 245 (Ikole), 264 (Akure) and 265 (Owo). Soil characteristics in the study area were considered. The major soil associations found in the study area were evaluated on the account of the type of soil and the associated water-holding capacity.
A Garmin 12-channel global positioning system (GPS) was used to obtain the eastings (longitude), northings (latitude) and elevation above the mean sea level of each point of interest during the fieldwork.
Results and discussion
Field measurements and hydrochemical analyses
The inventory of the shallow-dug wells revealed a depth range of 2.40-14.10 m with a mean value of 6.65 ± 2.47 m and a static water level of 1.55-13.12 m with a mean value of 5.89 ± 2.35 m. The occurrence of groundwater is controlled by a number of factors such as type of parent rock, depth, extent and pattern of weathering, the sand/clay ratio and the degree of fracturing, fissuring and jointing [18, 19]. The shallow groundwater system is characterised by low values of TDS (average 33.57 mg/L), EC of 0.79-77 ms/cm (average 13.20 ms/cm), pH values of 3.50-8.20 (average 6.50) and total hardness of 20-420 mg/L (average 112.23 mg/L). Concentrations of major cations are in the order of Ca > Mg > Na > K with average values of 51.76, 51.29, 8.60 and 6.59 mg/L, respectively. Bicarbonate and chloride are the dominant anions with average concentrations of 87.70 and 53.38 mg/L, respectively. High chloride content plays a major role in the corrosivity of buried metallic materials. Dissolution reactions of many metallic materials involve chlorides [3, 5, 14]. The corrosion process is enabled by the electrolyte between the anodic and cathodic sites. The
Figure 3: Soil map of Ado-Ekiti.
moisture content of soil acts as an electrolyte in the form of soluble salts such as chlorides and sulphates. The moisture content in typical soil samples collected at a depth of 1 m at the OdoAdo area of the metropolis was found to be in the range of 12.1-30.4%. In general, clayey and humus soils hold maximum moisture content than sandy and gravelly soils. An acidic pH value indicates a good electrolyte as more hydrogen ions are available to act as electron acceptors. Slightly acidic/alkaline soil pH levels tend to decrease the soil resistivity and promote corrosivity. High soil resistivity slows down the corrosion activities due to less ionic current flow. Resistivity is thus a function of moisture and the concentration of current-carrying soluble ions. Low electrical resistivity is indicative of good electrical conducting path arising from reduced aeration, increased electrolyte saturation or high concentration of dissolved salts in soils [3, 15, 16].
Geoelectric type curves
The characteristics of geoelectric curves varied greatly as typical of the basement complex terrain. They include the A, AA, H, HA, HK, K, KH, Q and QH curve types and combinations with the H-type curve accounting for 18.11%. This is an indication of the degree of weathering and fracturing [1, 7, 20].
Evaluation of soil corrosivity The three major soil units distributed across the study area include the Iwo, Ondo and Okemesi Associations (Figure 3]. The nature of residual soils in the study area is determined by the underlying geology. The rates of infiltration and permeability are directly interrelated to soil characteristics. The Iwo soil type is underlain by coarse-grained granite, gneiss and charnock-ite. The soil is composed of coarse-textured, greyish brown to brown sandy, fairly clayey soils. The soil type is widespread in the study
Figure 4: Resistivity of the topsoil.
area. The sandy nature of this soil promotes infiltration. The Ondo association is found on medium-grained granite and gneiss underlain areas. The soil comprises fine- to medium-textured, orange brown to brownish red, fairly clayey soils overlying orange, brown and red mottled clay. The Okemesi association is located on quartz schists and gneisses. The soil is composed of very coarse-textured, gravelly, pale grey brown to brown, usually sandy soil [21]. The topsoil in the metropolis varied in composition from clay, sandy clay, clayey sand to sand and laterite with the topsoil resistivity values ranging from 18 to 6410 Qm. The highest % frequency occurred in the resistivity values between 100 and 200 Qm. The low resistivity end (p < 1000 Qm] is diagnostic of alluviumsand horizon, while the high resistivity end (p > 1000 Qm] typifies laterites and compact sand. The wide resistivity range is a consequence of the variable composition of this layer,
degree of fluid saturation (or moisture content) and degree of compaction. Moisture content in soil is significant when considering corrosion potential. A dry soil environment is associated with a high resistivity with practically no corrosion potential. The resistivity decreases rapidly, and corrosion is promoted with increases in moisture content until the saturation point is reached [15, 18, 19]. The thickness of the top-soil is commonly around 1.3 m. Tables 1 and 2 give classification of soil corro-sivity in terms of resistivity. The frequency distribution of corrosivity level within the topsoil is presented in Figure 5. Topsoil materials indicating corrosivity levels of practically non-corrosive (PNC), slightly corrosive (SC) and moderately corrosive (MC) had coverage of 48.87%, 39.1% and 12.03%, respectively.
Figure 5: Corrosivity level of the topsoil..
Table 1: Soil electrical resistivity/corrosivity classification (BS - 1377).
Soil resistivity (Qm)
Soil corrosivity
<10	Severe
10-50	Corrosive
50-100	Moderately corrosive
>100	Slightly corrosive
Table 2: Classification of soil resistivity in terms of corrosivity	
[20,22,23].	
Soil resistivity (Qm)	Soil corrosivity
<10	Very strongly
	corrosive (VSC)
10-60	Moderately
	corrosive (MC)
60-180	Slightly
	corrosive (SC)
>180	Practically
	noncorrosive (PNC)
A large portion of the metropolis (% frequency of 48.87%) is practically non-corrosive with resistivity values of p >180 Qm within the topsoil, particularly areas overlain by lateritic hardpan with relatively high resistivity values. Relatively low resistivity values are indicative of high tendency for corrosivity. Slightly corrosive materials with resistivity values of 60 < p < 180 Qm occupy 39.10% of the top-soil and are observed at the eastern, southern, northeastern flanks and the central portion. Moderately corrosive topsoils with resistivity values of 10 < p < 60 Qm are delineated around Eureka/Oke Ureje (Figure 6). The second layer coincides with the regolith of the H- and HA-type curves, which predominates the area with 18% and 14% occurrence, respectively. The layer is characterised by resistivity values ranging from 3.2 to 5200 Qm (Figure 7). The low resistivity end (p < 60 Qm) is diagnostic of silt or clay horizon with little or no sand content. The resistivity values reflect the varying degree of weathering, the bedrock
Figure 6: Corrosivity map of the topsoil.
Figure 7: Resistivity map of the second layer.
Figure 8: Corrosivity level of the second layer.
structure and mineralogy. It is typically clayey with low layer resistivity values (p < 100 Qm) over basic charnockite and sandy/clayey sand (p > 100 Qm) on fine-coarse-grained granitic/ gneissic rocks. Formation of a clayey regolith is enhanced with intense chemical weathering of the parent rocks.
The soil corrosivity within the second layer indicates a frequency distribution (Figure 8) of 1.5%, 30.08%, 31.58% and 36.84% for VSC, MC, SC and PNC levels, respectively. Figure 9 shows moderately corrosive/slightly corrosive materials around the eastern, southern and northeastern flanks. The northwestern, southeastern, eastern and central parts of the metropolis are practically non-corrosive with resistivity values of p >180 Qm, indicating reduced porosity and negligible fluid content and degree of saturation.
Evaluation of the aquifer protective capacity The longitudinal unit conductance varies widely (0.01-4.40 mhos) across the metropolis (Figure 10). The parameter presents the combination of the thickness and resistivity of the geoelectric layers into a single variable. The qualitative use of this parameter is to demarcate changes in the total thickness of low-resistivity materials, hence its utilisation for evaluating the protective capacity of overburden units in an area [7, 20, 24]. A clayey overburden that is highly impervious presupposes relatively high
longitudinal conductance and offers effective protection to the underlying aquifer. Figurell shows the overburden protective capacity distribution of the study area. The protective capacity of the overburden has been zoned into good, moderate, weak and poor (Table 3).
Table 3: Longitudinal conductance/protective capacity rating [7,24].
Longitudinal conductance (mhos)
Protective capacity rating
>10	Excellent
5-10	Very good
0.7-4.9	Good
0.2-0.69	Moderate
0.1-0.19	Weak
<0.1	Poor
Zones where the conductance is greater than 0.7 mhos are considered as zones of good protective capacity. The portions having conductance values ranging from 0.2 to 0.69 mhos are classified as zones of moderate protective capacity, areas with values ranging from 0.1 to 0.1 9 mhos are classified as areas of weak protective capacity and the zones where the
Figure 9: Soil corrosivity map of the second layer.
Figure 10: Conductance map of the longitudinal unit.
Figure 11: Frequency distribution of aquifer protective capacity.
Figure 12: Aquifer protective capacity map of Ado-Ekiti.
conductance values are less than 0.1 mhos are demarcated as regions of poor overburden protective capacity. Poor/weak aquifer protective capacity of the overburden units are observed around the northwestern and southern axes (Figure 12).
On a regional consideration, 23.31%, 18.80% and 57.9% of the study area is characterised by overburden materials of poor, weak and moderate protective capacity, respectively. Only 6.02% of the area indicates good overburden protective capacity. This scenario suggests ap-
preciable caution in safety practices in well completions and general quest for groundwater protection in the metropolis.
Conclusions
Evaluation of the soil corrosivity and aquifer protective capacity of overburden units in Ado-Ekiti, southwestern Nigeria, had been conducted. Corrosion of cast iron, ductile iron and steel in soils can lead to a range of failures especially in pipelines and buried storage tanks. Integrity assessment of subsurface infrastructure, such as buried steel components, pipelines and steel sheet piles, requires an understanding of the local conditions. Assessment of soil corrosivity is thus germane to design of pipe networks as it provides a useful guide in the selection and prescription of the subsurface steel pipes for a given project and perhaps any required treatment to forestall economic waste and varied hazards associated with the rupture of corroded pipes. Poor, weak, moderate and good aquifer protective capacity zones were delineated in the study area. Areas characterised by poor and weak/moderate aquifer protective capacity should be void of potential contaminant load to ensure overall protection of the groundwater resource. Use of corrosion-resistant pipes is recommended according to the corrosivity level and the design specifications.
Acknowledgement
We acknowledge with thanks all sources of ancillary data for this work. We appreciate the constructive criticisms and useful suggestions of the Editor-in-Chief and the anonymous reviewers.
References
[1] Akintorinwa, O.J, Abiola, O. (2011): Sub-soil evaluation for pre-foundation study using geophysical and geotechnical approach. Journal of Emerging Trends in Engineering and Applied Sciences, 2(5), pp. 858-863.
[2]	Saupi, S.R.A., Sulaiman, M.A, Masri, M.N. (2015]: Effects of Soil Properties to Corrosion of Underground Pipelines: A Review. Journal of Tropical Resources and Sustainable Science, 3, pp. 14-18.
[3]	Roberge, PR. (2008]: Corrosion Electrochemistry. In Corrosion Engineering: Principles and Practice, Roberge, PR. (ed.], McGraw-Hill Education, pp. 35-47.
[4]	Adesida, A., Faleye, E.T, Fatoba, J. (2002]: Electrical Resistivity Survey for Corrosive Soils at WAPCO, Ewe-koro Factory, Ogun State, Nigeria. Journal of Science and Technology Research, 1(1], pp. 22-32.
[5]	Gopal, M. (2010]: Corrosion Potential Assessment, The Geology of part of South-western Nigeria. Geological Survey of Nigeria, pp. 31-87.
[6]	Braga, A.C.O., Filhow, W.M, Dourado, J.C. (2006]: Resistivity (DC] Method Applied To Aquifer Protection Studies. Brazilian Journal of Geophysics, 24(4], pp. 573-581.
[7]	Oladapo, M.I., Mohammed, M.Z., Adeoye, O.O., Adetola, B.A. (2004]: Geoelectrical Investigation of the Ondo State Housing Corporation Estate, Ijapo Akure, Southwestern Nigeria. Journal of Mining and Geology, 40(1], pp. 41-48.
[8]	Ehirim, C.N., Nwankwo, C.N. (2010]: Evaluation of aquifer characteristics and groundwater quality using geoelectric method in Choba, Port Harcourt. Archives of Applied Science Research, 2(2], pp. 396-403.
[9]	Tsepav, M.T., Adamu, Y., Umar, M. A. (2015]: Evaluation of Aquifer Protective Capacity and Soil Cor-rosivity Using Geoelectrical Method. International Scholarly and Scientific Research & Innovation, 9(11], pp. 662-671.
[10]	Mallam, A, Emenike, E.A. (2008]: Preliminary Findings of Subsurface Characteristics From Direct Current Resistivity Survey of The Federal Capital Territory (FCT], Nigeria. International Journal of Pure and Applied Sciences, 2(2], pp. 68-76.
[11]	Ojo, E.O., Adelowo, A., Abdulkarim, H.M and Dauda, A.K. (2015]: A Probe into the Corrosivity Level and Aquifer Protective Capacity of the Main Campus of the University of Abuja, Nigeria: Using Resistivity Method, Physics Journal, 1(2], pp. 172-178.
[12]	Olajuyigbe, A.E. (2010]: Sustainable Water Service Delivery: An Assessment of a Water Agency in a Rapidly Urbanizing City in Nigeria. Journal of Sustainable Development, 3(4], pp. 210-219.
[13]	Rahaman, M. A. (1988]: Recent advances in the study of the Basement Complex of Nigeria. In Oluyide et.al. (eds.], Precambrian Geology of Nigeria. Geological Survey of Nigeria: Kaduna, pp. 157-163.
[14]	Kleiner, Y., Rajani, B, Krys, D. (2010): Impact of Soil Properties on Pipe Corrosion: Re-examination of Traditional Conventions. Water Distribution System Analysis- WDSA 2010, Tucson, AT., USA, pp. 12-15.
[15]	Elarabi, H, Elkhawad, T. (2014): Evaluation of Subsoil Corrosivity Condition around Baracaia Area using the Electrical Resistivity Method, A Case Study from the Muglad Basin, Southwestern Sudan. Journal of Earth Science and Engineering, 4, pp. 663-667.
[16]	Tijani, M.N, Oke, S.A, Olowookere, A.T. (2014): Hydro-geochemical characterization of a shallow groundwater system in the weathered basement aquifer of Ilesha area, southwestern Nigeria. Evolving Water Resources Systems, 364, 475-489.
[17]	Vander Velpen, B.P.A. (1988): Resist Version 1.0. M.Sc. Research Project. ITC: Delft, Netherlands.
[18]	Jayeoba, A., Oladunjoye, M.A. (2013): Hydro-geophysical evaluation of groundwater potential in hard rock terrain of southwestern Nigeria, RMZ - Materials & Geoenvironment, 60, pp. 271-285.
[19]	Teikeu, W.A., Ndougsa-Mbarga, T., Njandjock, P.N, Tabod, T.C. (2012): Geoelectric Investigation for Groundwater Exploration in Yaounde Area, Cameroon. International Journal of Geosciences, 3, pp. 640-649.
[20]	Adeniji, A.E., Omonona, O.V, Obiora, D.N, Chukudebe-lu, J.U. (2014): Evaluation of soil corrosivity and aquifer protective capacity using geo-electrical investigation in Bwari basement area; Abuja. Journal of Earth System Science, 123(3), pp. 491-502.
[21]	Ojo, J.S., Olorunfemi, M.O, Akinluyi, F.O., Bayode, S, Akintorinwa, O.J, Omosuyi, G.O. (2015): Evaluating Soil Erosion Risk in the Basement Complex Terrain of Akure Metropolis, Southwestern Nigeria. Journal of Geography and Geology, 7(1), pp. 56-64.
[22]	Agunloye, O. (1984): Soil aggressivity along steel pipeline route at Ajaokuta, southwestern Nigeria. Journal of Mining and Geology, 21, pp. 97-101.
[23]	Oki, O. A., Egai, A.O, Akana, T.S. (2016): Soil Corrosivity Assessment in the Pre-Design of Sub-Surface Water Pipe Distributary Network in Yenagoa, South-South Nigeria Using Electrical Resistivity. Geosciences, 6(1), pp. 13-20.
[24]	Obiora, D.N, Ajala, A.E, Ibuot, J.C. (2015): Evaluation of aquifer protective capacity of overburden unit and soil corrosivity in Makurdi, Benue state, Nigeria, using electrical resistivity method. Journal of Earth System Science, 124(1), pp. 125-135.
Instructions for Authors
RMZ - MATERIALS & GEOENVIRONMENT (RMZ - Materiali in geookolje) is a periodical publication with four issues per year (established in 1952 and renamed to RMZ - M&G in 1998). The main topics are Mining and Geotechnology, Metallurgy and Materials, Geology and Geoenvironment.
RMZ - M&G publishes original scientific articles, review papers, preliminary notes, and professional papers in English. Only professional papers will exceptionally be published in Slovene. In addition, evaluations of other publications (books, monographs, etc.), in memoriam, presentation of a scientific or a professional event, short communications, professional remarks and reviews published in RMZ - M&G can be written in English or Slovene. These contributions should be short and clear. Authors are responsible for the originality of the presented data, ideas and conclusions, as well as for the correct citation of the data adopted from other sources. The publication in RMZ - M&G obligates the authors not to publish the article anywhere else in the same form.
Specification of the Contributions
Optimal number of pages is 7 to 15; longer articles should be discussed with the Editor-in-Chief prior to submission. All contributions should be written using the ISO 80000.
—	Original scientific papers represent unpublished results of original research.
—	Review papers summarize previously published scientific, research and/or expertise articles on a new scientific level and can contain other cited sources which are not mainly the result of the author(s).
—	Preliminary notes represent preliminary research findings, which should be published rapidly (up to 7 pages).
—	Professional papers are the result of technological research achievements, application research results and information on achievements in practice and industry.
—	Publication notes contain the author's opinion on newly published books, monographs, textbooks, etc. (up to 2 pages). A figure of the cover page is expected, as well as a short citation of basic data.
—	In memoriam (up to 2 pages), a photo is expected.
—	Discussion of papers (Comments) where only professional disagreements of the articles published in previous issues of RMZ - M&G can be discussed. Normally the source author(s) reply to the remarks in the same issue.
—	Event notes in which descriptions of a scientific or a professional event are given (up to 2 pages).
Review Process
All manuscripts will be supervised shall undergo a review process. The reviewers evaluate the manuscripts and can ask the authors to change particular segments, and propose to the Editor-in-Chief the acceptability of the submitted articles. Authors are requested to identify three reviewers and may also exclude specific individuals from reviewing their manuscript. The Editor-in-Chief has the right to choose other reviewers. The name of the reviewer remains anonymous. The technical corrections will also be done and the authors can be asked to correct the missing items. The final decision on the publication of the manuscript is made by the Editor-in-Chief.
Form of the Manuscript
All papers must be submitted via the online system.
The original file of the Template is available on RMZ - Materials and Geoenvironment Home page address: www.rmz-mg.com.
The contribution should be submitted in Microsoft Word. Manuscript should be written in Arial font and 12 pt font with 1.5 line spacing and should contain all figures, tables and formulas. Headings should be written in Arial bold font and shuold not be numbered. Subheadings should be written in Arial italic font. The electronic version should be simple, without complex formatting, hyphenation, and underlining. For highlighting, only bold and italic types should be used.
Annexes
Annexes are images, spreadsheets, tables, and mathematical and chemical formulas. Math formulas should be included in article as editable text and not as images.
Annexes should be included in the text at the appropriate place, and they should also be submitted as a separate document, i.e. separated from the text in the article.
Annexes should be originals, made in an electronic form (Microsoft Excel, Adobe Illustrator, Inkscape, AutoCad, etc.) and in .eps, .tif or .jpg format with a resolution of at least 300 dpi.
The width of the annex should be at least 152 mm. They should be named the same as in the article (Figure 1, Table 1). The text in the annexes should be written in typeface Arial Regular (6 pt).
The title of the image (also schemes, charts and graphs) should be indicated in the description of the image. When formatting spreadsheets and tables in text editors, tabs, and not spaces, should be used to separate columns. Each formula should have its number written in round brackets on its right side.
References of the annexes in the text should be as follows: "Figure 1..." and not "as shown below:". This is due to the fact that for technical reasons the annex cannot always be placed at the exact specific place in the article.
Composition of the Manuscript
Title
The title of the article should be precise, informative and not longer that 100 characters. The author should also indicate the short version of the title. The title should be written in English as well as in Slovene..
Information on the Authors
Information on the authors should include the first and last name of the authors, the address of the institution and the e-mail address of the corresponding author.
Abstract
The abstract presenting the purpose of the article and the main results and conclusions should contain no more than 180 words. It should be written in Slovene and English.
Key words
A list of up to 5 key words (3 to 5) that will be useful for indexing or searching. They should be written in Slovene and English.
Introduction
Materials and methods
Results and discussion
Conclusions
Acknowledgements
References
The references should be cited in the same order as they appear in the article. Where possible the DOI for the reference should be included at the end of the reference. They should be numbered in square brackets. References should be cited according SIST ISO 690:1996 standards.
Book:
[1]	Reynolds, J.M. (2011). An introduction to Applied and Environmental Geophysics. New York: Wiley, 710 p. Unpublished Master thesis or PhD dissertation:
[2]	Trcek, B. (2001): Solute transport monitoring in the unsaturated zone of the karst aquifer by natural tracers. Ph. D. Thesis. Ljubljana: University of Ljubljana 2001; 125 p.
Chapter in an edited book:
[3]	Blindow, N., Eisenburger, D., Illich, B., Petzold, H., Richer, T. (2007): Ground Penetrating Radar. In: Environmental Geology: Handbook of Field Methods and Case Studies, Knödel, K., Lange, G., Voigt, H.J. (eds.). Springer: Berlin; pp. 283-335.
Journal article : Journal title should be complete and not abbreviated.
[4]	Higashitani, K., Iseri, H., Okuhara, K., Hatade, S. (1995): Magnetic Effects on Zeta Potential and Diffusivity of Nonmagnetic Particles. Journal of Colloid and Interface Science, 172, pp. 383-388.
[5]	Mcmechan, G.A, Loucks, R.G, Mescher, P.A, Xiaoxian, Z. (2002): Characterization of a coalesced, collapsed pale-ocave reservoir analog using GPR and well-core data. Geophysics, 67, pp. 1148-1158. doi: 10.1190/1.1500376
Proceedings Paper:
[6]	Benac, C., Grzancic, Z., Sisic, S., Ruzic, I. (2008): Submerged Karst Phenomena in the Kvarner Area. In: Proceedings of the 5th International ProGEO Symposium on Conservation of the Geological Heritage, Rab, Croatia, Mar-janac, T (ed.). Pro GEO Croatia: Zagreb; pp. 12-13.
Electronic source:
[7]	CASREACT - Chemical reactions database [online]. Chemical Abstracts Service, 2000, renewed 2/15/2000 [cited 2/25/2000]. Available on: <http://www.cas.org/casreact.html>.
Scientific articles, review papers, preliminary notes and professional papers are published in English. Only professional papers will exceptionally be published in Slovene.
Units
SI System should be used for units. Physical quantities should be written in Italics (e.g. m, l, v, T). Symbols for units should be in plain text with spaces (e.g. 10 m, 5.2 kg/s, 2 s-1, 50 kPa). All abreviations should be spelt out in full on first appearance.
Manuscript Submission
Please submit your article via RMZ-M&G Editorial Manager System. You can find it on the address: http://edmgr.editool.com/rmzmag/default.htm
Log in as an author and submit your article.
You can follow the status of your submission in the system manager and your e-mail.
Information on RMZ - M&G
—	Editor-in-Chief
Assoc. Prof. Dervaric Evgen
Telephone: +386 1 47 04 616
E-mail address: evgen.dervaric@ntf.uni-lj.si
—	Secretary Vukic Nivesc
Telephone: +386 01 47 04 610 E-mail address: nives.vukic@ntf.uni-lj.si
These instructions are valid from April 2017.