JOURNAL OF ENVIRONMENTAL GEOGRAPHY Journal of Environmental Geography 6 (3–4), 39–47. DOI: 10.2478/jengeo-2013-0005 ISSN: 2060-467X THE ADVANTAGES OF USING SEQUENTIAL STOCHASTIC SIMULATIONS WHEN MAPPING SMALL-SCALE HETEROGENEITIES OF THE GROUNDWATER LEVEL László Mucsi *1 , János Geiger 2 , Tomislav Malvic 3,4 1 Department of Physical Geography and Geoinformatics, University of Szeged, Egyetem u. 2-6, H-6722 Szeged, Hungary 2 Department of Geology and Palaeontology, University of Szeged, Egyetem u. 2-6, H-6722 Szeged, Hungary 3 Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, Pierottijeva 6, 10000 Zagreb, Croatia 4 INA Plc *Corresponding author, e-mail: mucsi@geo.u-szeged.hu Research article, received 21 June 2013, accepted 15 September 2013 Abstract In the environmental risk assessment of oil fields, a detailed knowledge of the heterogeneity of groundwater surfaces is absolutely indis- pensable. Based on theoretical considerations, in order to analyse small-scale heterogeneities, we decided that the Sequential Gaussian Simulation (SGS) approach seemed to be the most appropriate one. This method gives preference to the reproduction of small-scale heterogeneities at the expense of local accuracy. To test whether this kind of heterogeneity of the groundwater level corresponds to sedimentological variability, a point bar of the River Tisza (South-Hungary) was chosen. In variograms, the longest range was derived from the large-scale sedimentological heterogeneity of the point-bar, the medium range was in accordance with the radius of the meander and its direction coincided with the depositional strike of the meander, while the shortest range corresponded to the lateral heterogeneity of the deposits where the ground water level was measured. The similarities and differences of the realizations of SGS express the uncer- tainty of the map representation of the ground water surface. The E-type estimates of 100 equiprobable realizations resulted in a very detailed surface. The hydraulic gradient map obtained from the E-type estimates can provide us with a better understanding of the local flow characteristics. Keywords: groundwater level, point bar, sequential Gaussian simulation, sedimentology, heterogeneity, geoinformatics INTRODUCTION The environmental risk assessment of oil fields highly requires a detailed knowledge of the heterogeneity of groundwater surfaces (Vaanlocke et al., 2010; Balderac- chi et al., 2013). In many cases the contamination caused by leaks can only be identified after it has appeared in the groundwater or on the surface. Therefore, the spatial- temporal analysis of the surface of chemically aggressive groundwater could help one pinpoint the most probable location of a natural pipeline leak. In order to monitor a particular environmental sys- tem, we chose the kind of methods suitable for data gathering and analysis. The primary goal of any monitor- ing activity is partly to locate the source of contamina- tion, and partly to analyse its lateral and vertical disper- sion. The ground water surface can be measured only in groundwater monitoring well, therefore the database can contain real data for these points. The groundwater level has to be estimated for other point. In the “traditional” gridding or triangulation meth- od, the numerical inter/extrapolations are formulated according to the method of least squares estimation. The aim of such methods is to improve the local accuracy and demonstrate the regional tendencies (Journel, 1987; Deutsch and Journel, 1998). In this case, the traditional estimation approaches act as “low-pass filters”, filtering out the small-scale heterogeneities affecting the principal tendencies. Kriging has become a very popular estima- tion method (Theodossiou, 1999; Theodossiou and La- tinopoulos, 2006), and it has this filtering property. The kriged surfaces (i.e. surfaces with different geological attributes) appear unnaturally smooth, and the resulting flow and transport dispersion are systematically underes- timated (Deutsch, 2002). The stochastic simulations differ from kriging ap- proaches in the following ways (Carr and Myers, 1985): In most interpolation algorithms (e.g. kriging) the goal is to find the “best” local estimate of the variable without specific regard to the resulting spatial statistics of the estimates taken together. In the solution, the re- production of global features and statistics (histogram, variogram) take precedence over the local accuracy. Based on the simulated groundwater surface, the analysis of contaminant dispersion could also be im- proved. Hence knowledge of the actual position of the groundwater level as well as its flow path may reduce the environmental risk associated with the joint (com- bined?) effects of a high groundwater level and other natural corrosives. 40 Mucsi et al. (2013) STUDY AREA The largest mining and industrial area of the sout h- eastern part of Hungary is the Oil and Gas Plant of MOL (Hungarian Oil and Gas Co.). Here almost every aspect of hydrocarbon exploitation, pipeline transpor- tation and oil processing is carried out ( Fig. 1). Fig. 1 The area surveyed on a Landsat TM 453 (RGB) color composite. The oil, gas and thermal water wells labeled by yellow circles (Mucsi et al., 2004) Within the area surveyed, 100 or even 1000 km long oil and gas plant pipeline networks have been installed below the surface. In the neighbourhood of a pipeline leak, the environmental damage caused can be quite serious. For networks that were built a few decades ago, systematic pipeline monitoring should be performed due to frequently occurring natural leaks. In the area surveyed, a fairly dense pipeline ne t- work links the production wells with tank farms. This area is part of the recent valley of the River Tisza. The oil and gas production here, which began in the 1960’s, has significantly altered the former agrarian landscape (Fig. 2). However, agricultural activity has continued in parallel with mining and oil processing (see Fig. 3). Therefore the risk of a pipeline leak is not just real for MOL, but also for farmers. For the latter group, lower harvest yields due to poorer soil and groundwater quality may be two great causes for concern. At the end of the last century and beginning of this century, the repeated occurrence of a recor d high water level of the River Tisza, a high groundwater level (see Fig. 4) and the extremely dry summers may actually be the consequences of a changing global climate, and justify the need for a more co n- scious use (pro-active stewardship?) of our enviro n- ment. Fig. 2 The alteration of the former agrarian landscape of an aerial photography (1950) and a spy satellite (Corona – KH- 4B) image from 1972 (Hungarian Military Archive, USGS) Fig. 3 Ikonos satellite image from 2004 (INTA SPACETURK) On the surface, Holocene sand, silt and clay can be found owing to fluvial sedimentation. Based on sediment profiles obtained from the wells below the upper silty sandy lens-like series, a laterally extended silt layer appears 4-5 metres beneath the surface, which can be regarded as a flood plain deposit. Under this series, Late Holocene fluvial sand dominates. Originally the monitoring system consisted of 25 wells. Since they are located in a cluster, we decided to supplement the original monitoring system with 17 new wells (Fig.5). The advantages of using sequential stochastic simulations when mapping small-scale heterogeneities of the groundwater level 41 Fig. 4 High groundwater level (aerial photography, Mucsi et al, 2004) Fig. 5 Location map with the new and old monitoring wells THEORETICAL CONSIDERATIONS In general, the ground water level is determined by hydro- logical processes like precipitation (more precisely, its annu- al and long-term distribution), evaporation, transpiration, and pedological-sedimentological and topographical conditions. Here, the heterogeneity of pedological or sed i- mentological factors may be the biggest two factors among the processes of this multivariate system. U n- fortunately, this heterogeneity is not properly reflec t- ed in the maps of the groundwater surfaces. The purpose of this study is to demonstrate a mapping process of the groundwater level that is in consistent with small and medium scale sedimentolo g- ical heterogeneity by using Sequential Gaussian Simu- lation (SGS, Deutsch and Journel, 1992). For this purpose a recent point bar sedimentary body of the River Tisza (Hungary) was chosen. The scales of sedimentological heterogeneity in a point bar In the idealized cross section of a point bar, the fo l- lowing four basic units can be identified: (1) channel lag deposits, (2) laterally accreted bar deposits with accretion surfaces, (3) a natural levee, and (4) a flood plain series (Fig. 6). The bar sequence itself can be subdivided into (a) a lower point bar and (b) upper point bar series (Fig. 6). It is well known that the effective porosity and permeability of a sediment are partly determined by sedimentary textures and structures that permeate the pore space distributions. The influence of sedime n- tary structures and textures on permeability is very significant, not just for the absolute values but also for the anisotropy. A quite trivial conseque nce of this fact is that the steady state level of a groundw a- ter surface can vary according to the sedimentary structures and textures of the actual deposits. A c- cording to the basic rule of sedimentology (Walter’s facies law), lateral facies also occur sup erposed upon one another. It means that the vertical heterogeneity depicted in Fig. 5 can also horizontally affect the spatial position of the groundwater level. The lateral scale of this heterogeneity can vary from a few m e- tres to tens or hundreds of metr es (Pryor, 1987; Mial, 1996). However, according to the GPR images, there is a textural variability that has a range of a few hundred metres to a few kilometres. Summarizing the points above: in the point bar sed- imentary bodies, the lateral variability of sedimentary features can manifest itself on two different scales. They range (1) from several metres to tens or hundreds of metres and (2) from several hundred metres to several kilometres. The basic unit associated with this scale (1) is the scale of the principal sedimentary structure, while for scale (2) the basic unit is the distance between two neighbouring accretion surfaces. Fig. 6 Idealized cross-section of a point-bar sequence (mod- ified after Brunner et al, 2001) The ground water level treated as a regionalized variable When viewed mathematically, two statements can be derived from the above description. They are: (1) The ground water level is a spatial property that can be re- garded as a random variable at each data point; (2) There is a spatial relationship associated with these random variables. The first statement follows from the variability of the textural characteristics of the sedimentary body 42 Mucsi et al. (2013) containing the initial water table. The second statement is a consequence of the fact that a point bar is the result of spatially interacting processes. That is, according to the definition given by Cressie (1991), the water table level is a regionalized variable. A map of absolute position of a ground water level fixes the water table in the real spatial position. This kind of map is a plot of a three-variable function, which gives an estimate of the real position, or more precisely one realization of the several existing estimates. By continuing this process, the “smooth” surface in Fig. 7 tells us that the variability of hydrological elements affecting the water table is “not too large”. This then leads us to ask the following question: “Is the smooth surface simply the result of the gridding technique we applied?” Thus, traditional methods like kriging are really appropriate for describing large-scale heterogeneities. For point bars the latter ranges from a few hundred metres to a few kilometres. The small-scale heterogeneities are simply filtered out via these estimations. In fact the traditional gridding techniques honour the large-scale features quite well. For example, in Fig. 7 the lateral accretion surfaces are clearly recognizable. To overcome the above drawbacks, several methods for stochastic simulation have been developed and they have been regularly used in the mapping of reservoir proper- ties in the oil industry since the mid 1980’s, but they have been widely applied in the environmental sciences as well (Webster and Oliver, 2001). THE APPROACH WE ADOPTED SGS in general Let us define a regionalized variable (RV) and call it z(u). Stochastic simulations are methods in which alternative and equally probable high resolution models of spatial distribution of the z(u) are generated. The realizations obtained from them are called stochastic images. If the realizations honour the input data (the data points), the simulation is called conditional. For the local set of data and conditional statistics, kriging is used as an interpolation technique which pro- vides a simple numerical method that is the “best” in the sense that it is accurate locally. Simulation provides sev- eral alternative but equally probable models, all of which are the “best” reflection of reality in a certain global sense. The differences between the realizations offer an oppor- tunity for measuring the spatial uncertainty. The Gaussian simulations honour the covariance model of the data points. This is why they are suitable for modelling processes with an extremely large continuity. Let us define the common distribution of an RV by Zi (i=1,2,..,N). That is, according to the above theory, let us take the data of all available points and consider the conditioning of this N RV for each type of n data sets. The corresponding conditional cumulative distr i- bution function (CCDF) for this is: F(N)(z1, ... ,zn|(n))=P{Zi≤zi, i=1, ... ,N|(n)} If the ground water table is treated as regionalized variable, this variable must be a random variable for any given data point. That is, in the infinitesimally small neighbourhood of any data point it will have a random value. Hence the data-point value is nothing else than a random value of the corresponding RV derived from the distribution function of the infinitesimally small neigh- bourhood. If processes controlling the ground water table are homogeneous in the region being analysed, then the cumulative distribution function (CDF) surrounding any set of data points can be assumed to be the same as that of over the whole given region. By applying a normal score transformation on the data set, this distribution can be represented in analytical form as well. As is well known, the first two moments determine the normal distribution. This fact was incorporated in the implementation of SGS (Carr and Myers, 1985). SGS is a widely used simulation technique for ex- treme continuous data (e.g. porosity). Details of the vari- ous implementations are given in Deutsch and Journel (1998), and in Deutsch (2002). In essence, SGS models the method of iteration, since the estimation of a grid point is done based on the result of the previous step. This view corresponds to the fact that a contour map is unbiased only after the grid estimation method has been fixed (Brooker, 1979). The expected value (more precisely, the “expected surface”) of the series of realizations is the most character- istic spatial distribution that describes “reality”. The difference between the realizations reflects the uncertainty of the mapping of the property being studied - in our case, the surface of the water table. This uncertainty is independent of the accuracy of the measurement of the Fig. 7 Appearance of the lateral accretion surfaces o n the water table contours The advantages of using sequential stochastic simulations when mapping small-scale heterogeneities of the groundwater level 43 water table taken from each well. In fact, this uncertainty depends just on the heterogeneity of hydrological factors and the way in which this heterogeneity can be described using available well geometry. Pre-processing: the steps involved The first phase was an analysis of the lateral distribution and probability distribution of the data (Fig. 8a-b). In the second phase, two things became immediately apparent: (1) the peculiar clustering of data points; (2) the observation that there are many wells in each cluster, and these have significantly different data values (Fig. 8a). However, the probability distribution of data is quite symmetrical (Fig. 8b). This clustering appearance of wells is fairly common around artificial objects, where small distant wells serve the purpose of environmental protection while the larger dis- tance wells were built for water table observations. Howev- er, this may cause problems in mapping, because the clus- tered wells will shadow the nearby grid points with their own values. Hence the contours will not reveal the true situation. At the same time, this geometry will also affect the relative frequency histogram of the water table values. However, applying a special declustering technique can resolve this problem (Deutsch and Journel, 1998). With this method an appropriate grid system is placed over the area where (in such a way that?) the clustered wells can be sepa- Fig. 9 The variogram surface and the variogram model Fig. 8 The first steps of the geostatistical analysis 44 Mucsi et al. (2013) rated into individual wells. Then some weighting factors can be derived which are small for the very close wells and become larger with increasing distance. In this way the shadow effect of closer wells can be eliminated. In our case the declustering method gave 500 m as the most appropriate cell distance (Fig. 8c). The result of the normal score trans- formation is shown in Fig. 8d. Fig.10 The empirical (a) and model variogram surfaces (b) Variography Our analysis of the spatial continuity of the water table was aided by variography. A variogram surface was used for the visual examination of continuity directions (see Fig. 9a). This method identified a principal direction of NE-SW. The smallest continuity may be identified in the NW-SE direc- tion (Fig. 9b). These directions appear on the smooth sur- face of the water table as well (see Fig. 7). As mentioned previously, this feature corresponds to the lateral accretion surfaces of the point bar. The second figure here shows a variogram model based on the experimental data. As you can see, this is a nested structure consisting of three units. The first com- ponent has a range of 250 m. The second one has a range of about 3.800 m, while the third has a range of 4.100 m. Most likely, the latter is indicated by the ‘smooth’ map of the water table surface. The largest range is almost equal to the length of the point bar, so it may be related to the large-scale pro- cess that forms the point bar sedimentary body. The second largest range is quite similar to the previous one and may be related to the effect of a bend in the river channel going in a SE direction. But the smallest range (250 m) may not be caused by hydrological effects at all, but by sedimentological effects. This may be the result of internal heterogeneity of the point bar affecting the surface of the water table. Hence it is the heterogeneity which was really the object of this study. Fig. 10a depicts the spatial continuity outlined above, while the second figure shows the variogram surface obtained from the variogram model (Panattier, 1996). The similarity is obvious here. Realizations and simulation uncertainties From the SGS, 100 equally probable realizations were generated. Fig. 11 shows 3 of the resulting one hundred images. As can be seen, there are regions where the realiza- tions are not too different. This is caused by the locally low uncertainty of the lateral interpolation (Journel, 1993). Fig. 11 Stochastic realizations showing equally probable posi- tion of the water table The advantages of using sequential stochastic simulations when mapping small-scale heterogeneities of the groundwater level 45 Fig. 12 E-type estimation from the two-hundred stochastic realizations of the relative water table position The one hundred realizations provide one hun- dred values estimated for each grid node. This number is quite sufficient to compute the cumulative distrib u- tion around an infinitesimally small neighbourhood of each grid point. That is, there is a possibility of est i- mating not just the expected value, but also the lower and upper bounds of the confidence interval for each grid node. Moreover, we have a way of explicitly a n- swering questions like “How can we calculate the lateral distribution of probabilities associated with the water table level at 72 cm below the surface?” and “Which regions with a water table of 85 cm have a probability value of 0.8?”. Fig. 12 shows the predicted value map of the relative water surface (i.e. its depth below the surface), while Fig. 13 shows the surface of the water table above the Baltic Sea. The ‘smooth surface’ obtained by traditional kriging and the one obtained by the E-type estimates (expected values) of the SGS are compared in Fig. 14a and Fig. 14b, respectively. Fig. 14b highlights the local effects of underground flow systems. The advantages of this approach in the three-phase flow simulation were described earlier (see, for instance, Geiger and Komlosi, 1996). Their observa- tions confirmed that these approaches make the flow sim- ulations more effective (than earlier approaches?). Fig. 13 E-type estimation from the two-hundred stochastic realizations of the absolute water table position. The black arrows are proportional with the hydraulic potential 46 Mucsi et al. (2013) CONCLUSION Within the area surveyed, 100 or even 1000 km long oil and gas plant pipeline networks have been i n- stalled below the surface. The spatial-temporal analy- sis of the surface of chemically aggressive groundwa- ter helped one pinpoint the most probable location of a natural pipeline leak. Originally the monitoring system consisted of 25 wells. Since they are located in a cluster, the original monitoring system was supplied with 17 new wells whose locations were calculated with geostatistical analysis. The ground water level was measured monthly in these wells and Sequential Gaussian Simulation was used to demonstrate a ma p- ping process of the groundwater level that is in co n- sistent with small and medium scale sedimentological heterogeneity. A map of absolute position of a ground water level fixed the water table in the real spatial position. The developed map was a plot of a three - variable function, which gave an estimate of the real position, or more precisely one realization of the sev- eral existing estimates. The analysis of the spatial continuity of water table was aided by variography. A variogram surface was used for the visual examination of continuity directions. This method identified a principal direction of NE-SW. The smallest continuity was identified in the NW-SE direction. These direc- tions appear on the smooth surface of the water table as well In variograms, the longest range was derived Fig. 14 Comparison of the results coming partly from kriging ( a) and partly from E-type estimation of SGS (b) The advantages of using sequential stochastic simulations when mapping small-scale heterogeneities of the groundwater level 47 from the large-scale sedimentological heterogeneity of the point-bar, the medium range was in accordance with the radius of the meander and its direction coi n- cided with the depositional strike of the meander, while the shortest range corresponded to the lateral heterogeneity of the deposits where the ground water level was measured. The similarities and differences of the realizations of SGS expressed the uncertainty of the map representation of the ground water surface. The E-type estimates of 100 equiprobable realizations resulted a very detailed surface. The hydraulic grad i- ent map obtained from the E-type estimates deter- mined the local flow characteristics. The subsurface contamination can spread along these pathways, so it can get far distance from its source. References Balderacchi, M. ,Benoit, P., Cambier, P Ekloc O. M., Gar- ginid, A., Gemitzi, A., Gurelf, M., Kløve, B, Nakich, Z., Predaai E., Ruzicic S., Wachniewj, P., Trevisana M. (2013). Groundwater. Pollution and Quality Monitoring Approaches at the European Level. Critical Reviews in Environmental Science and Technology 43, 323–408. Brooker, P. 1979. Kriging. Engineering and Mining Journal 1980 9, 148–153. Brunner, D.S.Endres, A.L.Sudicky, E.A. 2001. Detailed ground-penetrating radar survey of a point-bar, Whiteman’s Creek, Ontario for use in a new fully inte- grated 3D surface/ subsurface flow model. – Geological Society of America Annual Meeting, November 5-8, 2001. Boston, Massachusetts. Carr, J.R Myers, D.E., 1985. COSIM: A Fortran IV program for Conditional Simulations.—Computer and Geoscienc- es v.11, 675–705. N. Cressie 1991. Statistics for Spatial Data. John Wiley & Sons, Inc. New York, 892p. Deutsch C.V. 2002. Geostatistical Reservoir Modelling, Ox- ford University Press. Deutsch, C.V., Journel, A. 1998. GSLIB. Geostatistical Soft- ware Library and User’s Guide. – Oxford University Press, New York, 369p. Geiger, J. Komlósi, J. 1996. Szedimentológiai geomatematikai 3-D modellező rendszer törmelékes CH-tárolókban (Sed- imentological, geomatematical 3D modelling system in turbiditic CH reservoirs) Kőolaj és Földgáz 2, 53–81. Journel, A. 1987. Geostatistics for the environmental sciences. – EPA Project no. Cr 811893. Technical report. U.S. EPA, EMS Lab., Las Vegas, NV Journel, A, 1993. Modelling uncertainty: Some conceptual thoughts. In: Dimitrakopoulos, R.(ed): Geostatistics for the Next Century, Kluwer, Dordrecht, 30–43. Miall, A.D. 1996. The Geology of Fluvial Deposits. Sedimen- tary Facies, Basin Analysis, and Petroleum Geology. – Springer Verlag, New York, 582 p. Mucsi, L. 2001. Characterisation of oil-industrial contamina- tion using aerial and thermal images - EARSeL Sympo- sium, Drezda. In: Buchroithner (ed:) A Decade of Trans- European Remote Sensing Cooperation, Balkema, Rot- terdam, 373–377. Mucsi, L., Kiss, R., Szatmári, J., Bódis, K., Kántor, Z., Dabis, G., Dzsupin, M. 2004. Felszín alatti vezetékek környezetszennyező hatásainak felmérése távérzékeléses technológiával (Invesgtigation of possible environmental pollution of subsurface pipelines using remote sensing technologies). Geodézia és Kartográfia 56 ( 4) 3–8. Pryor, W.A. 1987. Permeability, porosity patterns and varia- tions in some Holocene sand bodies. In: Beaumont, E.A., Foster, N.H. (ed) Reservoirs II. Sandstones. AAPGD Treatise of Petroleum Geology Reprint Series. No.4, 127–154. Theodossiou, N. 1999. Evaluation of the distribution of hy- draulic head in an aquifer using the Kriging method. – Scientific Journal of Hellenic Hydrotechnical Associa- tion. Hydrotechnika 9, 3–14. Theodossiou, N., Latinopoulos, P. 2006. Evaluation and opti- mization of groundwater observation network using the Kriging methodology. Environmental Modelling & Soft- ware 21, 991–1000. Vanloocke, R., De Borger, R., Voets J. P. . Verstraete W 2010. Soil and groundwater contamination by oil spills; prob- lems and remedies International Journal of Environmen- tal Studies 8 (1-4), 99-111. Webster, R., Oliver, M. A. (2001): Geostatistics for Environ- mental Scientists, John Wiley & Sons. http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Eklo%2C+O+M%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Gargini%2C+A%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Gargini%2C+A%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Gemitzi%2C+A%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Gurel%2C+M%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Kl%C3%B8ve%2C+B%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Nakic%2C+Z%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Predaa%2C+E%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Ruzicic%2C+S%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Wachniew%2C+P%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Trevisan%2C+M%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Vanloocke%2C+R%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28De+Borger%2C+R%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Voets%2C+J+P%29 http://www.tandfonline.com/action/doSearch?action=runSearch&type=advanced&searchType=journal&result=true&prevSearch=%2Bauthorsfield%3A%28Verstraete%2C+W%29