MWH Global

Enviroscan Ukrainian Institute of Speleology and Karstology

Community news

Speleology in Kazakhstan

Shakalov on 04 Jul, 2018
Hello everyone!   I pleased to invite you to the official site of Central Asian Karstic-Speleological commission ("Kaspeko")   There, we regularly publish reports about our expeditions, articles and reports on speleotopics, lecture course for instructors, photos etc. ...

New publications on hypogene speleogenesis

Klimchouk on 26 Mar, 2012
Dear Colleagues, This is to draw your attention to several recent publications added to KarstBase, relevant to hypogenic karst/speleogenesis: Corrosion of limestone tablets in sulfidic ground-water: measurements and speleogenetic implications Galdenzi,

The deepest terrestrial animal

Klimchouk on 23 Feb, 2012
A recent publication of Spanish researchers describes the biology of Krubera Cave, including the deepest terrestrial animal ever found: Jordana, Rafael; Baquero, Enrique; Reboleira, Sofía and Sendra, Alberto. ...

Caves - landscapes without light

akop on 05 Feb, 2012
Exhibition dedicated to caves is taking place in the Vienna Natural History Museum   The exhibition at the Natural History Museum presents the surprising variety of caves and cave formations such as stalactites and various crystals. ...

Did you know?

That yield is the quantity of water discharged from an aquifer [16] (e.g spring or well.) see also well yield.?

Checkout all 2699 terms in the KarstBase Glossary of Karst and Cave Terms

What is Karstbase?

Search KARSTBASE:

keyword
author

Browse Speleogenesis Issues:

KarstBase a bibliography database in karst and cave science.

Featured articles from Cave & Karst Science Journals
Chemistry and Karst, White, William B.
See all featured articles
Featured articles from other Geoscience Journals
Karst environment, Culver D.C.
Mushroom Speleothems: Stromatolites That Formed in the Absence of Phototrophs, Bontognali, Tomaso R.R.; D’Angeli Ilenia M.; Tisato, Nicola; Vasconcelos, Crisogono; Bernasconi, Stefano M.; Gonzales, Esteban R. G.; De Waele, Jo
Calculating flux to predict future cave radon concentrations, Rowberry, Matt; Marti, Xavi; Frontera, Carlos; Van De Wiel, Marco; Briestensky, Milos
Microbial mediation of complex subterranean mineral structures, Tirato, Nicola; Torriano, Stefano F.F;, Monteux, Sylvain; Sauro, Francesco; De Waele, Jo; Lavagna, Maria Luisa; D’Angeli, Ilenia Maria; Chailloux, Daniel; Renda, Michel; Eglinton, Timothy I.; Bontognali, Tomaso Renzo Rezio
Evidence of a plate-wide tectonic pressure pulse provided by extensometric monitoring in the Balkan Mountains (Bulgaria), Briestensky, Milos; Rowberry, Matt; Stemberk, Josef; Stefanov, Petar; Vozar, Jozef; Sebela, Stanka; Petro, Lubomir; Bella, Pavel; Gaal, Ludovit; Ormukov, Cholponbek;
See all featured articles from other geoscience journals

Search in KarstBase

Your search for simulation (Keyword) returned 143 results for the whole karstbase:
Showing 16 to 30 of 143
Mixed transport reaction control of gypsum dissolution kinetics in aqueous solutions and initiation of gypsum karst, 1997, Raines M. A. , Dewers T. A. ,
Experiments with gypsum in aqueous solutions at 25 degrees C, low ionic strengths, and a range of saturation states indicate a mixed surface reaction and diffusional transport control of gypsum dissolution kinetics. Dissolution rates were determined in a mixed flow/rotating disc reactor operating under steady-state conditions, in which polished gypsum discs were rotated at constant speed and reactant solutions were continuously fed into the reactor. Rates increase with velocity of spin under laminar conditions (low rates of spin), but increase asymptotically to a constant rate as turbulent conditions develop with increasing spin velocity, experiencing a small jump in magnitude across the laminar-turbulent transition. A Linear dependence of rates on the square root of spin velocity in the laminar regime is consistent with rates being limited by transport through a hydrodynamic boundary layer. The increase in rate with onset of turbulence accompanies a near discontinuous drop in hydrodynamic boundary layer thickness across the transition. A relative independence of rates on spinning velocity in the turbulent regime plus a nonlinear dependence of rates on saturation state are factors consistent with surface reaction control. Together these behaviors implicate a 'mixed' transport and reaction control of gypsum dissolution kinetics. A rate law which combines both kinetic mechanisms and can reproduce experimental results under laminar flow conditions is proposed as follows: R = k(t) {1 - Omega(b)() zeta [1 - (1 2(1 - Omega(b)())(1/2)]} where k(t) is the rate coefficient for transport control, and Omega(b)() is the mean ionic saturation state of the bulk fluid. The dimensionless parameter zeta(=Dm(eq)()/2 delta k() where m(eq)() = mean ionic molal equilibrium concentration, D is the diffusion coefficient through the hydrodynamic boundary layer, delta equals the boundary layer thickness and k() is the rate constant for surface reaction control) indicates which process, transport or surface reaction, dominates, and is sensitive to the hydrodynamic conditions in the reactor. For the range of conditions used in our experiments, zeta varies from about 1.4 to 4.5. Rates of gypsum dissolution were also determined in situ in a cavern system in the Permian Blaine Formation, southwestern Oklahoma. Although the flow conditions in the caverns were not determinable, there is good agreement between lab- and field-determined rates in that field rate magnitudes lie within a range of rates determined experimentally under zero to low spin velocities A numerical model coupling fluid flow and gypsum reaction in an idealized circular conduit is used to estimate the distance which undersaturated solutions will travel into small incipient conduits before saturation is achieved. Simulations of conduit wall dissolution showed-member behaviors of conduit formation and surface denudation that depend on flow boundary conditions (constant discharge or constant hydraulic gradient and initial conduit radius. Surface-control of dissolution rates. which becomes more influential with higher fluid flow velocity, has the effect that rate decrease more slowly as saturation is approached than otherwise would occur if rates were controlled by transport alone. This has the effect that reactive solutions can penetrate much farther into gypsum-bearing karst conduits than heretofore thought possible, influencing timing and mechanism of karst development as well as stability of engineered structures above karst terrain

Strategic ground water management for the reduction of karst land collapse hazard in Tangshan, China, 1997, Wang H. T. , Li Y. X. , Wang E. Z. , Zhao Z. Z. ,
Karst land collapse is a geological hazard caused mainly by human activities such as the pumping of ground water. The present study is aimed at determining a management strategy to minimize the risk of karst land collapse in the urban area of Tangshan, China. A method of groundwater management for multi-aquifer systems with groundwater confined and unconfined interchange properties has been developed. A groundwater simulation model and a management model are constructed for the studied area. The data used by the simulation and management models include the geometry of the aquifer systems, the parameters of hydrogeology, the distribution of recharge, and discharge in the simulation period and the historical water table elevations. The results of the model calculations show that the total balance between recharge and discharge for both Quaternary and bedrock aquifers is 4.7 x 106 m(3) in the period from November 1989 to October 1990. The optimal annual average pumping rate from the underlying bedrock aquifer is some 25.8 x 10(6) m(3), which is 54.4% of that which was pumped out of the aquifer in the hydrological year of 1990. If the management plan is carried out, the groundwater level in the bedrock aquifer will recover its confined state in most of the areas at risk of land collapse. This is the most resistant state to collapse in terms of groundwater flow.

Risk assessment methodology for karst aquifers .2. Solute-transport modeling, 1997, Field Ms,
Ground-water flow and solute-transport simulation modeling are major components of most exposure and risk assessments of contaminated aquifers. Model simulations provide information on the spatial and temporal distributions of contaminants in subsurface media but are difficult to apply to karst aquifers in which conduit flow is important. Ground-water flow and solute transport in karst conduits typically display rapid-flow velocities, turbulent-flow regimes, concentrated pollutant-mass discharge, and exhibit open-channel or closed-conduit how Conventional groundwater models, dependent on the applicability of Darcy's law, are inappropriate when applied to karst aquifers because of the (1) nonapplicability of Darcian-flow parameters, (2) typically nonlaminar flow regime, and (3) inability to locate the karst conduits through which most flow and contaminant transport occurs. Surface-water flow and solute-transport models conditioned on a set of parameters determined empirically from quantitative ground-water tracing studies may be effectively used to render fate-and-transport values of contaminants in karst conduits. Hydraulic-flow and geometric parameters developed in a companion paper were used in the surface-water model, TOXI5, to simulate hypothetical slug and continuous-source releases of ethylbenzene in a karst conduit. TOXI5 simulation results showed considerable improvement for predicted ethylbenzene-transport rates and concentrations over qualitative tracing and analytical ground-water model results. Ethylbenzene concentrations predicted by TOXI5 simulations were evaluated in exposure and risk assessment models

Numerical versus statistical modelling of natural response of a karst hydrogeological system, 1997, Eisenlohr L, Bouzelboudjen M, Kiraly L, Rossier Y,
Structural and hydrodynamic characteristics of karst aquifers are mostly deduced from studies of global responses of karst springs (hydrographs, chemical or isotopic composition). In this case, global response is often used to make inferences with respect to infiltration and ground water How processes as well as on the hydrodynamic parameters. Obviously, the direct verification of these inferences is very difficult. We have used an indirect method of verification, introducing well defined theoretical karst structures into a finite element model and then analysing the simulated global response according to the currently accepted interpretation schemes. As we know what we are introducing into the numeric model, the consistency of the interpretation may be checked immediately. The results obtained in the hydrogeological study of two karst basins in the Swiss Jura and from 2-D and 3-D numerical simulations show the difficulty of finding structural parameters and hydrodynamic behaviour from statistical methods alone, i.e. correlation analyses discharge-discharge and precipitation-discharge. In effect, our first results show that the form of the correlograms depends on several factors besides the structure of the karst aquifer: (i) on the form of the floods. in other words the contrast between quick Row and base How, (ii) on the frequency of hydrological events during the period analysed and (iii) on the type of infiltration processes, in other words the ratio of diffuse infiltration to concentrated information. Obviously, the variability of a karst hydrograph is a result of a combination of these factors. Distinction between them is not always possible on hydrographs, and therefore on correlations (discharge-discharge and precipitation-discharge). (C) 1997 Elsevier Science B.V

Numerical simulation as a tool for checking the interpretation of karst spring hydrographs, 1997, Eisenlohr L, Kiraly L, Bouzelboudjen M, Rossier Y,
A schematic representation of karst aquifers may be that of a high hydraulic conductivity channel network with kilometre-wide intervals, surrounded by a low hydraulic conductivity fractured limestone volume and connected to a local discharge area, the karst spring, The behaviour of the karst spring (hydrographs, chemical or isotopic composition, etc.) represents the global response of the karst aquifer to input events. The available data an karst aquifer hydraulic parameters are limited, Global response is therefore more easily obtained and is commonly used to make inferences on the recharge and groundwater How processes, as well as on the hydraulic parameter fields. Direct verification of these interpretations is, obviously, very difficult. We have used an indirect method of verification, consisting of introducing well-defined theoretical karst structures into a finite element model and then analysing the simulated global response according to presently accepted interpretation schemes. As we know what we put into the numerical model, the validity of any interpretation may be checked. The first results indicate that some of the generally accepted interpretations are not necessarily true. In particular: (i) separation of simulated recession hydrographs into several components shows that different exponential components do not necessarily correspond to aquifer volumes with different hydraulic conductivities: (ii) non-exponential parts of recession hydrographs do not always give information about the infiltration process: and (iii) the recession coefficient of the baseflow (i.e. the last, nearly exponential part of the recession hydrograph) depends on the global configuration of the whole karst aquifer, not just on the hydraulic properties of the low hydraulic conductivity volumes. (C) 1997 Elsevier Science B.V

Simulation of the evolution of maze caves., 1997, Clemens T. , Hiickinghaus D. , Sauter M. , Liedl R. , Teutsch G.
The development of cave systems in carbonate rocks depends on a variety of boundary and initial conditions. Among the cave systems, two main types of geometries can be distinguished: the dendritic and maze pattern. A numerical model has been developed capable of modeling the genesis of karst systems in complex geological environments. It is applied to simulate the development of the above mentioned two different types of cave geometries. The results confirm that a prerequisite for the development of maze caves is evenly distributed recharge (White 1969). However more important for the development of maze caves is that flow through the system is restricted by an overlying less conductive horizon, e.g. a sandstone caprock. Thus the feed back mechanism of higher flow rates leading to higher dissolution rates and therefore the preferential development of a small number of tubes does not dominate the evolution of the karst aquifer. This hydraulic restriction furthers the development of other conduits also to achieve a significant diameter.

Simulation of the evolution of maze caves, 1997, Clemens T. , Huckinghaus D. , Sauter M. , Liedl R. , Teutsch G.

Evolution of size distributions of natural particles during aggregation: modelling versus field results, 1998, Atteia O,
In this paper a discretized model simulating aggregation of size distributions jointly with sedimentation and transport is presented. A review of the current theory provides some helpful hints about the relative importance of each aggregation process, i.e. Brownian motion, shear flow and differential sedimentation, which are tested by using collision efficiency factors. The novel aspect of the model arises from the use of a varying mean particle diameter in each size class. This allows both non-steady-state and steady-state calculations and free choice of size classes. A comparison with a classical approach shows the exactitude of the results and the improvment obtained for several cases. The simulations gave a family of curves characterized by three parts corresponding to peri-, and orthokinetic aggregation and to sedimentation. The role of collision effciency is crucial in the relative extent of each part of the size distribution. The comparison with a series of data from a karst spring showed that the model was able to fit most of the particle size distributions using significant values of each parameter. This allowed information about particle aggregation and transport within a non-accessible aquifer to be inferred.

Simulation of the development of karst aquifers by using a pipe flow model coupled to a continuum model, 2, Model verification and sensitivity analysis., 1998, Clemens T. , Huckinghaus D. , Sauter M. , Liedl R. , Teutsch G.

Influence of aperture variability on dissolutional growth of fissures in karst formations, 1998, Hanna R. B. , Rajaram H. ,
The influence of aperture variability on dissolutional growth of fissures is investigated on the basis of two-dimensional numerical simulations. The logarithm of the fissure aperture before dissolution begins is modeled as a Gaussian stationary isotropic random field. The initial phase of dissolutional growth is studied up to the time when turbulent flow first occurs at a point within the fissure (the breakthrough time). The breakthrough time in variable aperture fissures is smaller than that in uniform fissures and decreases as the coefficient of variation of the aperture field (sigma/mu) increases. In comparing uniform and variable aperture fissures in limestone, the breakthrough time with sigma/mu = 0.1 is about a factor of 2 smaller than that in a uniform fissure. The breakthrough time is reduced by about an order of magnitude with sigma/mu = 2.0. The mechanism leading to reduced breakthrough times is the focusing of flow into preferential flow channels which are enlarged at a faster rate than the surrounding regions of slower flow. Dissolution channels are narrower and more tortuous as sigma/mu. increases. Investigations of the influence of reaction rate reveal that the influence of aperture variability is more pronounced in rapidly dissolving rock. In uniform fissures in rapidly dissolving minerals, breakthrough times are very long since water becomes saturated with respect to the mineral within a short distance of the entrance to the flow path. However, in variable aperture fissures, breakthrough occurs rapidly because of selective growth along preferential flow channels, which progressively capture larger fractions of the total flow. These results partly explain why conduits develop rapidly in gypsum, although previous one-dimensional studies suggest that conduit growth will not occur

Oxidation of organic matter in a karstic hydrologic unit supplied through stream sinks (Loiret, France), 1998, Alberic P, Lepiller M,
The aim of this paper is to appraise the ability of the oxidation of riverine organic matter in the control of limestone dissolution, in a karst network. Biogeochemical processes during infiltration of river water into an alluvial aquifer have already been described for an average flow velocity of 4-5 m d(-1) (Jacobs, L. A., von Gunten, H. R., Keil, R, and Kuslys, M. (1988) Geochemical changes along a river-groundwater infiltration flow path: Glattfelden, Switzerland. Geochim. Cosmochim. Acta 52, 2693-2706; Von Gunten, H. R., Karametaxas, G., Krahenbuhl, U., Kuslys, M., Giovanoli R., Hoehn E. and Keil R. (1991) Seasonal biogeochemical cycles in riverborne groundwater. Geochim. Cosmochim. Acta 55, 3597-3609; Bourg, A. C. M. and Bertin, C. (1993) Quantitative appraisal of biogeochemical chemical processes during the infiltration of river water into an alluvial aquifer. Environ. Sci. Technol. 27, 661-666). Karstic drainage networks, such as in the River Loire-Val d'Orleans hydrologic system (Fig. 1), make possible flow velocities up to 200 m h(-1 a) and provide convenient access to different water samples several tens of km apart, at both extremities of the hydrologic unit (Chery, J.-L. (1983) Etude hydrochimique d'un aquifere karstique alimente par perte de cours d'eau (la Loire): Le systeme des calcaires de Beauce sous le val d'Orleans. These, Universite d'Orleans; Livrozet, E. (1984) Influence des apports de la Loire sur la qualite bacteriologique et chimique de l'aquifere karstique du val d'Orleans. These, Universite d'Orleans). Recharge of the karstic aquifer occurs principally from influent waters from stream sinks, either through coarse alluvial deposits or directly from outcrops of the regional limestone bedrock (Calcaires de Beauce). Recharge by seepage waters From the local catchment basin is small (Zunino, C., Bonnet, M. and Lelong, F. (1980) Le Val d'Orleans: un exemple d'aquifere a alimentation laterale. C. R. somm. Soc. Geol. Fr. 5, 195-199; Gonzalez R. (1992) Etude de l'organisation et evaluation des echanges entre la Loire moyenne et l'aquifere des calcaires de Beauce. These, Universite d'Orleans) and negligible in summer. This karstic hydrologic: system is the largest in France in terms of flow (tens to hundreds of m(3)/s) and provides the main water resource of the city of Orleans. Chemical compositions of influent waters (River Loire) and effluent waters (spring of the river Loiret) were compared, in particular during floods in summer 1992 and 1993 (Figs 2-4). Variation of chloride in the River Loire during the stream rise can be used as an environmental tracer of the underground flow (Fig. 2). Short transit times of about 3 days are detectable (Fig, 2) which are consistent with earlier estimations obtained with chemical tracers (Ref. in Chery, J.-L. (1983) These, Universite d'Orleans). Depending on the hydrological regime of the river, organic carbon discharge ranges between 3-7 and 2-13 mg/l for dissolved and particulate matter respectively (Fig. 3). Eutrophic characteristics and high algal biomasses are found in the River Loire during low water (Lair, N. and Sargos, D. (1993) A 10 year study at four sites of the middle course of the River Loire. I - Patterns of change in hydrological, physical and chemical variables in relation to algal biomass. Hudroecol. Appl. 5, 1-27) together with more organic carbon rich suspended particulate matter than during floods (30-40 C-org % dry weight versus 5-10%). Amounts of total organic carbon and dissolved oxygen (Fig. 3) dramatically decrease during the underground transport, whereas conversely, dissolved calcium, alkalinity and inorganic carbon increase (Fig. 4). Anoxia of outflows map start in April. Dissolution of calcium carbonates along the influent path outweighs closed system calcite equilibrium of inflow river waters (Table 3). The impact of organic matter oxidation on calcite dissolution may be traced by variations of alkalinity and total carbonates in water. Following, Jacobs, L. A., von Gunten, H. R., Keil, R. and Kuslys, M. (1988) Geochemical changes along a river-groundwater infiltration flow path: Glattfelden, Switzerland. Geochim. Cosmochim. Acta 52, 2693-2706), results are shown graphically (Fig. 5). Extent of reactions is controlled by the consumption of dissolved O-2 and nitrate for organic matter oxidation and by the release of Ca2 for calcite dissolution (Table 2). The karstic network is considered to behave like a biological reactor not exchanging with the atmosphere, with steady inhabitant microbial communities (Mariotti A., Landreau A, and Simon B. (1988) N-15 isotope biogeochemisrry and natural denitrification process in groundwater: Application to the chalk aquifer of northern France. Geochim. Cosmochim. Acta 52, 1869-1878; Gounot, A.-M. (1991) Ecologie microbienne des eaux ei des sediments souterrains. Hydrogeologie, 239-248). Thus, energy requirements only are considered, not carbon assimilation. Moreover, there is no necessity to invoke any delay for nitrification enhancement, as observed elsewhere, after waste water discharge into the river (Chesterikoff, A., Garban, B., Billen, G. and Poulin, M. (1992) Inorganic nitrogen dynamics in the River Seine downstream from Paris (France). Biogeochem. 17, 147-164). Main microbial processes are assumed to be aerobic respiration, nitrification and denitrification. Reactions with iron and manganese, real but not quantitatively important, were neglected. Sulphate reduction and methane formation, certainly not active, were not considered. Denitrification, which is suggested by low nitrate and ammonium concentrations and anoxia in the outflow, is known to be rapid enough to be achieved in a short time (Dupain, S. (1992) Denitrification biologique heterotrophe appliquee au traitement des eaux d'alimentation: Conditions de fonclionnement et mise au point d'un procede. These, Universite Claude Bernard, Lyon). Reaction are somewhat arbitrary but conform to general acceptance (Morel, M. M. and Hering, J. G. (1993) Principles and Applications of Aquatic Chemistry. Wiley, New York). Anaerobic ammonium oxidation (Mulder A., van de Graaf, A. A., Robertson, L: A. and Kuenen, J. G. (1995) Anaerobic ammonium oxidation discovered in a denitrifying fluidized bed reactor. FEMS Microbiol. Ecol. 16, 177-184). although possible, was not considered. In fact, C/N ratio of the reactive organic matter has only mild repercussions on the results; i.e. in the same range as the analytical errors for alkalinity and total carbonates. The objective was simply to roughly confront characteristics of outflowing waters and the calculation. Respective roles of aerobes and denitrifiers, for instance, are not certain. Several periods during low water or floods were selected with various ranges for calcium dissolution or nitrate and oxygen concentrations. The result is that in most cases simulation and data are in reasonable accordance (Fig. 5). Amounts of organic matter in River Loire are generally sufficient to sustain the process (Table 3. Particulate organic matter is probably the most reactive. The balance of oxidation of organic matter indicates that about 65 mu g C-org/l.h are oxidized during the transport without much variation with the river regime or organic discharge. It is concluded that limestone dissolution is directly dependent on organic matter oxidation, but variation occurs (7-29 mg CuCO3/l) with the level of bases that can be neutralized in the River Loire water. (C) 1998 Elsevier Science Ltd. All rights reserved

Flow pattern variability in natural fracture intersections, 1999, Kosakowski G. , Berkowitz B. ,
We use numerical simulations to examine the variability of flow patterns in representative fracture intersection geometries. In contrast to existing studies of perfectly orthogonal intersections, we demonstrate that more realistic geometries lead to a rich spectrum of flow patterns. Moreover, numerical solutions of the Navier-Stokes equations in these fracture intersections indicate that non-linear inertial effects become important for Reynolds numbers as tow as 1-100. Such Reynolds numbers often exist in naturally fractured formations, particularly in karst systems and in the vicinity of wells during pump tests

Determining karst transmissivities with inverse modeling and an equivalent porous media, 1999, Larocque M. , Banton O. , Ackerer P. , Razack M. ,
Flow simulation is difficult to implement in heterogeneous media such as karst aquifers, primarily because the structure of the rock is extremely complex and usually unknown. The aim of this study was to verify the possibility of using inverse modeling and an equivalent porous media to identify transmissivities in a slightly karstified aquifer, the La Rochefoucauld karst (Charente, France), Different simulation scenarios were tested: using two spatial discretizations with different finite-element cell sizes and using measured or interpolated heads. The inverse modeling was performed with the downscaling parameterization procedure, using a finite-element representation of bidimensional ground water flow. The inverse modeling converged satisfactorily with all scenarios: head residuals were small and spring flow rates and the river/aquifer exchanges were adequately stimulated. The scenario using small cells and measured heads generated a highly heterogeneous transmissivity field, indicating an overparameterization of the problem. The calibrated transmissivities and simulated heads of this scenario proved less reliable overall than those of the other scenarios. The use of interpolated heads generated more uniform transmissivities as a result of the head smoothing. A rotation of the initial parameter mesh showed that the scenarios using interpolated heads generate the most stable and reliable results. The scenarios with interpolated heads could therefore be used when head measurements are limited or are unevenly distributed over the aquifer. Overall, the calibrated transmissivities reproduced the entire range of transmissivities measured in the field using different methods. The results indicate that inverse modeling and an equivalent porous media can be used to determine transmissivities in a moderately karstified aquifer

Simulation of the development of karst aquifers: role of the epikarst., 1999, Clemens T. , Huckinghaus D. , Liedl R. , Sauter M.

Simulation of the development of karst aquifers: role of the epikarst, 1999, Clemens T, H_ C, Liedl R, Sauter M,

Results 16 to 30 of 143
You probably didn't submit anything to search for