ACTA CARSOLOGICA 33/2 7 115-149 LJUBLJANA 2004 COBISS: 1.01 FORECASTING VERSUS PREDICTING SOLUTE TRANSPORT IN SOLUTION CONDUITS FOR ESTIMATING DRINKING-WATER RISKS PRIMERJAVA RAZLIČNIH METOD NAPOVEDOVANJA PRENOSA RAZTOPIN V KRASU PRI OCENI OGROŽENOSTI PITNE VODE MALCOLM S. FIELD1 1 U.S. Environmental Protection Agency, National Center for Environmental Assessment (8623D) 1200 Pennsylvania, Ave., N.W., Washington, D.C. 20460, phone: (202) 564-3279, fax: (202) 565-0079 e-mail: field.malcom@epa.gov 2 Disclaimer: The views expressed in this paper are solely those of the author and do not necessarily reflect the views or policies of the U.S. Environmental Protection Agency. Abstract UDC: 628.1:551.44 Malcolm S. Field: Forecasting Versus Predicting Solute Transport in Solution Conduits for Estimating Drinking-Water Risks Abstract: Contaminant releases in karstic terranes can cause rapid and devastating affects on drinking-water supplies. Because future contaminant releases are likely it is necessary that local water managers develop release scenarios so as to be prepared prior to an actual contaminant release occurring. Release scenarios may be forecasted using appropriate historical data or they may be predicted using selected measured parameters. Forecasting contaminant releases to drinking-water supplies in karstic terranes is best accomplished by conducting numerous tracer tests from each potential source location to each exposure point so that acceptable solute-transport parameters for each solution conduit may be estimated from analyses of the breakthrough curves. Compositing the numerous breakthrough curves and fitting a quintic spline allows development of a single representative breakthrough curve that may then be used to forecast the effects of a release. Predicting contaminant releases is accomplished by combining basic measured field parameters for selected solution conduits in functional relationships for application in solute-transport models. The resulting breakthrough curve and solute-transport parameters can be used to predict the effects of a release. The forecasting and prediction methodologies were tested using a hypothetical release into a solution conduit developed in a karstic aquifer. Both methods were shown to produce reasonably acceptable results. The prediction methodology produced better time-of-travel results and better mass recovery and exposure concentration results than did the forecasting methodology. Key words: Forecasting, Prediction, Solute Transport, Solution Conduits, Drinking-water risks. Izvleček UDK: 628.1:551.44 Malcolm S. Field: Primerjava različnih metod napovedovanja prenosa raztopin v krasu pri oceni ogroženosti pitne vode Vnos polutantov v kras ima lahko hitre in uničujoče posledice za pitno vodo. Ker je nevarnost izlitja vedno navzoča, je prav, da upravniki vodnih virov predvidijo scenarije onesnaženja še preden se to v resnici zgodi. Različne scenarije vnosa polutantov v kras lahko predvidimo z uporabo že znanih podatkov in meritev v vodon-osniku. Napoved širjenja polutantov je najbolj zanesljiva, če temelji na rezultatih številnih sledenj med različnimi vhodnimi in izhodnimi točkami. Na ta način lahko primerne tranportne parametre določimo z analizo krivulj koncentracij. S kombiniranjem različnih krivulj in uporabo ustrezne aproksimacijske metode, lahko razultate združimo v eno samo krivuljo koncentracij, ki jo uporabimo za predvidenje širjenja polutantov. Širjenje polutantov lahko napovemo tudi z uporabo osnovnih tranportnih parametrov, ki veljajo za kraške prevodnike in ustreznih transportnih enačb. Rezultat so krivulje koncentracij, ki jih tudi uporabimo za napoved posledic možnega vnosa polutantov v kras. Obe metodi smo preizkusili na namišljenem vnosu polutantov v kraški kanal in pri obeh dobili sprejemljive rezultate. Metoda, ki temelji na uporabi transportnih enačb, se je izkazala za boljšo pri oceni potovalnega časa in koncentracije, ki smo ji izpostavljeni. Poleg tega je pri tej metodi delež povrnjenega sledila večji. Ključne besede: Napovedovanje, transport polutantov, kraški vodonosnik, ogroženost pitne vode 1. INTRODUCTION Contamination of air, soil, water and the consequent human health risks associated with the contamination are serious ongoing problems. It is generally accepted that for any contaminant release in any media, three practical questions need to be addressed which are most succinctly stated as where, when, and how much. Brouyere (2003, p. 10) addresses the where, when, and how much in slightly different terms as (Figure 1): 1. How long will it take for the contaminant to reach the downgradient receptor; 2. At what concentration level(s) will the receptor be contaminated; and 3. For how long will the contamination persist? It should be noted in Figure 1 that concentration C, time for first arrival f |, and time for last arrival ^^ are somewhat arbitrary in that the level of importance depends greatly on perception. The concern for C depends on what is considered hazardous for the particular pollutant while the values for t\ and f J depend on a selected concentration and instrument sensitivity A hypothetical example for an inorganic compound may be found in Mull et al. (1988, pp. 77-78) while hypothetical examples for an organic compound and a microbial pathogen may be found in Field (2004). Determining when, where, and at what concentration a pollutant will arrive is at the heart of all contamination investigations. The effect of accidental and deliberate releases of toxic substances to drinking-water supplies need to be forecasted or predicted if water managers are to initiate appropriate actions should a release occur. While somewhat difficult in terms is infinitely more difficult in ground water and especially so for karstic aquifers. The most effective method for evaluating pollutant releases in hydrologic systems is to proactively conduct of tracer tests. Quantitative hydrologic tracer testing is routinely used to establish solute-transport trajectories, to define hydraulic characteristics, and to define solute-transport parameters. Tracer testing also provides a sense of possible consequences resulting from future pollutant releases. It has been documented that conducting repeated tracer tests over an appropriate time period between the same input and output locations makes it possible to reasonably forecast the arrival of pollutants for individual solution conduits. Forecasting in general may be defined as the analysis of historical data and simple extrapolation of components of these data in order to estimate future values. Repeating tracer tests provides the necessary historical data (Kilpatrick and Taylor, 1986; Taylor et al., 1986; Mull and Smoot, 1986; Smoot et al., 1987; Mull et al., 1988, pp. 66-74). However, a sufficient number of tracer tests over an appropriate amount of time necessary for adequate pollutant-transport forecasting remains an unknown quantity. Prediction involves the incorporation of subjective factors into the estimate. Rather than basing pollutant transport on the results of one or more tracer tests previously conducted, predicting pollutant transport allows for some basic facts about the hydrologic system be known or estimated. These known or estimated facts generally entail typical hydrologic and geometric field parameters which may then be incorporated into a solute-transport model for individual solution conduits. By combining the forecasted or predicted results with risk-assessment analyses for acute exposures, water managers can develop a set of alternatives as part of an overall strategy for protecting human health. This set of alternatives could range from no action (i.e., no significant concern), changes in water treatments and water management, disconnecting the water-supply system, announcing a no contact warning, or arranging for the supply of an alternative water source. The effect of accidental and deliberate releases of toxic substances to drinking-water supplies need to be predicted if water managers are to initiate appropriate actions should a release occur. The purpose of this paper is to illustrate how reliable pollutant-transport forecasting and prediction may be possible and necessary for estimating future human health risks. The forecast and prediction models described mainly apply to a Type I karst network but may be applied to more complicated karst networks albeit with increasing uncertainty as complexity increases (Figure 2). The Forecasting methodology is based entirely on a set of breakthrough curves developed at each discharge point that specifically considers temporal variations in discharge (J. Similarly, the prediction methodology was specifically designed to account for increases in discharge (J, diverging pathways, and storage systems. Neither method is applicable to a Type V network, but logically only the prediction methodology could mistakenly be applied to a Type V network because the forecasting methodology cannot be applied to instances where BTCs have not been previously obtained. In addition, the forecasting and prediction methodology will be ineffectual on a basin-wide basis. Forecasting or predicting pollutant behavior in a karstic aquifer on a basin-wide basis requires a distributional numerical model which is not addressed by either of the methods developed in this paper. 2. FORECASTING AND PREDICTING METHODOLOGY The differences between forecasting and predicting methodologies are profound. Although both rely on numerical analyses of selected data, neither is even remotely similar to the other except in terms of discharge. Spring or well discharge is the overwhelming dominating factor controlling transport rate and downstream concentration. Figure 2. Seven simple karst network types that describe tracer migration in karst conduits. Any of these networks may significantly influence tracer tests between the point of inflow (IN) and the point of outflow (OUT) in a karst system. Discharge into the conduit is ti, discharge out of the conduit is tracer mass injected into the conduit is rrj,, and tracer mass recovered is J"^. Note: Any one of these network types may be interconnected with any of the others. Modifled from Atkinson et al. (1973) and Gaspar (1987, p. 64). Figure 3. Relationship between discharge (J and mean transport velocity 1). Figure 5. Relationship between discharge (J and mean traveltime standard deviation C, ■ Figure 7. Relationship between discharge (J and dispersivity (1( • Figure 9. Relationship between discharge (J and peak normalized concentration C„ 2.1. Significance of Discharge As originally documented by Mull et al. (1988, p. 66-71) discharge is the controlling factor that characterizes the hydrologic conditions associated with solute transport. This becomes evident when calculated transport parameters and solute concentrations are plotted as a function of discharge and the relationship between the two established. Figures 3-6 and Figures 9-10 were developed in similar fashion to those originally developed in Mull et al. (1988, p. 66-71) although second degree polynomials were developed here for Figures 4, 5, and 10 rather than first degree polynomials as was originally done by Mull et al. Figures 7, 8, and 11 have been developed here using the same logic as Mull et al. Figures 3-8 and 11 illustrate the relationship that exists between calculated transport parameters and discharge and Figures 9 and 10 illustrate the relationship that exists between concentration values and discharge. An examination of Figures 3-6, 9, and 10 allows for the following determinations (Mull et al., 1988, p. 67): 1. Mean flow velocity increases and mean traveltime decreases as discharge inceases: 2. The standard deviation of mean traveltime decreases as discharge increases because of a decreased travletime for solute dispersal; 3. Longitudinal dispersion increases as discharge increases as a result of increased turbulence and mixing caused by the increase in flow velocity; 4. Peak normalized concentration decreases as discharge increases due to increased dilution; 5. Peak normalized load increases as discharge increases because traveltime is reduced. Figures 7, 8, and 11, allows some additional specific determinations: 1. Dispersivity increases as discharge increases as a result of increased turbulence and mixing caused by the increase in flow velocity; 2. Peclet number decreases as discharge increases because longitudinal dispersion increases faster than does mean velocity; 3. Cross-sectional area increases with increasing discharge. 2.2. Forcasting Methodology To forecast pollutant transport, numerous tracer tests need to be conducted from the same input location and recovered at the same discharge location for each solution conduit. These multiple tracer tests will ideally be spread over a significant period of time to reflect varying hydrological conditions associated with different seasons and/or dry and wet years. The results of the tracer tests then need to be equated in a dimensionless model that can then be related to a specific release of a solute mass. The basic steps involved are as follows (Mull et al., 1988, pp. 76-77): 1. Perform multiple quantitative tracer tests in the flow system under varying hydrological conditions; 2. Determine necessary solute-transport parameters (mean traveltime f , traveltime standard deviation C,, and peak normalized concentration C„ ) associated with tracer tests using an appropriate methodology such as that described in QTrXcer2 (Field, 2002c); 3. Evaluate the relationship between discharge (J and and C^, ; 4. Convert each BTC into a standardized dimensionless BTC using t -J f. (1) (2) 5. Create a composite set of dimensionless BTCs; 6. From the composite set of dimensionless BTCs, numerically or manually create one single dimensionless BTC representative of the composite set; 7. For a given discharge and solute mass spilled, determine the parameters T, iT^ . and C^, from the relationships developed in step 3; 8. Plot the forecasted BTC and determine the full range of solute-transport parameters using an appropriate method such as that developed in QTRACER2. 9. a given discharge and solute mass spilled, determine the parameters / ^ + and ffrom the relationships developed in step 3; 0 '' The forecasting process can be considered onerous because of the need to conduct numerous tracer tests over a selected time period. In addition, step six was originally done manually and required considerable trial-and-error adjustment to obtain a desired fit and selected statistics: mean traveltime 7 , traveltime standard deviation tj,, and traveltime skewness Y^. A numerical procedure for developing one smooth standardized curve through all the composite curves has not been reported previously. However, appropriate methods do exist and greatly facilitate the process. 2.2.1. Quintic-Spline Fitting Polynomial interpolation for fitting is a common method for fitting a smooth curve through a set of data because of stability and robustness. However, end points, monotonicity, convexity, and continuity of derivatives can and do influence the results in contradictory ways (Weston, 2002). Cubic splines are generally more successful but, these have their limitations as well. For example, they do not work very well with non-uniform data and may lead to unpredictable results at the end points of data sets (Weston, 2002). Cubic Hermites are considerably better (Kahaner et al., 1989, pp. 94^95, 110) but are still problematic with this kind of data set. The quintic-spline interpolation algorithm develops a "piecewise quintic spline" function. That is, the interpolant is defined in terms of a set of quintic polynomials, each of which is defined between a pair of consecutive data points. The coefficients of these quintic polynomials are chosen so that the interpolant and its first four derivatives are continuous. This is not enough to uniquely determine the interpolant, and the remaining freedom of choice is used to ensure that the interpolant is "visually acceptable" so that monotonicity in the data results in monotonicity in the interpolant and the interpolant avoids "wiggles." A piece-wise interpolating quintic natural spline function ^(r for a set of monotonic data points (f^, ^ J = 1+2, K 1 jV , produces a spline curve that is both smoother and less volatile in response to variations in knot points (Weston, 2002). According to Grenville (1967) as cited by Weston (2002) if N > 2 then there is an interpolating quintic natural spline that has the following properties: • S(z ) is a fifth degree polynomial in each interval (-j. --j +-1) • Jl and its first four derivatives, upto and including "'"(r ) are continuous in -1 ■ - v ] • = /=1, .... .V The spline function can then be represented in the general form: ,V(.v) = V, + + f \r + /)/ + i-:/ + r/ (3) with f — z — for S r < i=l, 2,..., N. The actual solution for ) is accomplished according to Equations 4a-4d. (4a) (4b) (4c) (4d) where ^j; are knot positions representing sampling times. Because the quintic-spline methodology requires strict monotonicity, the several sets of standardized dimensionless curves need to be combined into a single file and sorted from lowest to highest in terms of time values prior to fitting a smooth curve through the data. This is not a difficult or unreasonable burden and is easily accomplished. 2.3. Prediction Methodology Predicting solute transport requires selected parameters be measured and combined into functional relationships for application in a solute-transport model such as the advection dispersion equation (ADE). The Efficient ^ydrologic Tracer-test Design (EHTD) Program (Field, 2002a, 2003) can be used to predict the effects of a toxic-substance release once source-water areas have been established (Field, 2003, 2004). By using the same measured or estimated parameters intended for tracer-mass estimation and specifying a solute mass, solute-transport parameters and downstream arrival concentrations for individual flow channels are predicted. The forecasting methodology previously discussed (Mull et al., 1988) is similar to the prediction methodology but the forecasting method requires considerably more time and effort. EHTD has been found to reliably reproduces known results with minimal information required. 2.3.1. Physical Parameters Predicting solute mass and transport time requires specific hydraulic and geometric parameters measured in the field be combined in simple functional relationships. These measured parameters and functional relationships can then be applied in a simple hypothetical model as a precursor to application of a mathematical model. For flowing streams in surface channels or solution conduits, measurements for discharge (J, cross-sectional area A , and transport distance .V need to be taken. These three parameters are the most important measurable field parameters necessary for establishing basic hydraulic and geometric controls in any hydrologic system. For porous media, additional measurements or estimates for specific parameters must also be taken. These additional parameters are transverse spread JI' and vertical spread H of the tracer plume and effective porosity n^. Moreover, it is necessary that the type of tracer test be identified as either a natural-gradient test or a forced-gradient tracer test in which a radially-symmetric flow field is created by an extraction well. Functional Relationships Discharge, cross-sectional area, and transport distance can each be measured at a downstream location such as a spring if a karstic aquifer or estimated from Darcy's law if a porous-media aquifer. Table 1 lists some of the equations that can be used with selected field parameters to calculate various solute-transport parameters. EHTD does not use all of the equations listed in Table 1 because several of the required field parameters are not readily obtainable. For example, the first two equations for ft^ in a turbulently flowing stream in open channels or closed conduits are not applicable to solution conduits in karstic terranes because the slope of the water surface S generally cannot be measured unless the conduit is accessible; the third D^ equation requires prior knowledge of the shear velocity which is not readily accomplished. To obtain an estimate of , EHTD uses the method developed by Chatwin (1971) because it tends to not overestimate longitudinal dispersion as is generally the case with most published equations. Combination of the measured parameters into functional relationships is accomplished using some established equations (Table 1) and by a manual estimation procedure. The process of combining parameters into functional relationships in EHTD is described in Field (2002a, 2003). 2.3.2. Mathematical Model EHTD predicts the effects of a toxic-substance release by initially predicting solute transport parameters and estimated solute mass as described in Field (2003). The basic form of the ADE used by EHTD appears as The retardation factor can represents retardation in porous media by (Freeze and Cherry, 1979, p. 404) (6a) represents retardation in fractured-media media by (Freeze and Cherry, 1979, p. 411) (6b) and represents retardation in solution conduits by (Field and Pinsky, 2000) The first-order rate coefficient for solute decay in porous media is given as (Toride et al., 1995,p. 3) (7a) in fractured-media may be represented as (7b) and in solution conduits as '' (7c) If solute decay for the liquid phase equals tracer decay for the sorbed phase Jt j, then the combined first-order decay p is equal to (Toride et al., 1995, p. 4) which is a reasonable assumption commonly employed (Toride et al., 1995, p. 35). The zero-order rate coefficient for solute production "f ) for porous media, fractured-rock media, and solution conduits is given as (Toride et al., 1995, p. 3) (8) where n^ may be taken as 1.0 for flow through individual fractures and solution conduits. Solute production is only applicable to exponential growth of bacteria in water with adequate growth conditions (e.g., high levels of turbidity and nutrient loading) an example of which may be found in (Field, 2004). EHTD solves the ADE as a boundary value problem BVP for a third-type inlet condition which conserves mass (Toride et al., 1995, p. 5). EHTD also allows consideration of an initial value problem IVP for uniform background concentration and production value problem PVP for exponential production (Toride et al., 1995, pp. 9-12). In dimensionless form the ADE appears as (Toride et al., 1995, p. 4) ST P,dz- ez ■ ' • ^ ' (9) where C^ represents the reduced volume-average solute concentration. The solution to the ADE for resident concentration and third-type inlet condition is given as (Toride et al., 1995, p. 6) c^l/., + +(10) where the ^ superscript denotes resident concentration, and the P superscripts denote boundary, initial, and production value problems, respectively (Toride et al., 1995, p. 6). All other parameters are described in the Notations section. p. 8) Boundaiy Value Problem The BVP may be solved for impulse release as a Dirac ) function by (Toride et al., 1995, and for a pulse release for the case where JJ ^ = 0 by (Toride et al., 1995, p. 8) and for the case where JJ ^ * 0 by (Toride et al., 1995, p. 8) where the /T superscript denotes dimensionless equilibrium concentration. The auxiliary functions I , I , and r/' are defined in the appendix. Initial Value Problem The IVPmay be solved for uniform initial concentration by (Toride et al., 1995, p. 10) (14) The auxiliary function Fj is defined in the appendix. (15) Production Value Problem The PVP may be solved for solute production that changes exponentially with distance using (Toride et al., 1995, p. 12) which gives the solution as (Toride et al., 1995, p. 12) (16) where the auxiliary functions Fj'' and Pj^' are defined in the appendix. Selecting a solute mass causes the present mean solute concentration C in EHTD to be overridden and a new t" to be predicted. A typical BTC representing the downstream effects of the release is then produced. 3. METHODS VERIFICATION To evaluate the forecasting and prediction methodologies for estimating the effects of a release of a chemical or biological agent, a solution conduit in a karstic aquifer in which a relatively small spring is used for public drinking water was investigated. Tracer tests have established a connection between a distant karst window (depression revealing part of a subterranean river flowing across its floor, or an unroofed portion of a cave (Field, 2002b, p. 110) and the spring establishing the system as a Type I karst network (Figure 2). The spring is near a local transportation route which serves to illustrate its relative vulnerability to an accidental chemical spill. 3.1. Forecasting Process The forecasting methodology can be verified by examining the quality of the standardized curve fit to the combined dimensionless BTCs in relation to the dimensional BTCs. This is achieved by adjusting the standardized curve by use of calculated values for the relevant parameters, 'I'Tj. and peak concentration C^ that are representative of a particular tracer test. Figure 12. Breakthrough curves for seven tracer tests conducted over a period of one year. Figure 13. Composite plot of standardized dimensionless BTCs for seven tracer tests conducted over a one year period. Parameter Equation Description Volume r y - At Open channel or closed conduit volume Porous media volume |' with regional gradient Porous media volume |' with forced gradient Mean Velocity U Open channel or closed conduit flow velocity ** Porous media flow velocity U with regional gradient Porous media flow velocity U with regional gradient'' Porous media flow velocity U with forced gradient' u = Porous media flow velocity U with tracer recirculation'' Mean Traveltime t Obtained directly from velocity U estimate Requires cross-sectional area A be approximated Traveltime Duration Ij rj_.r / ■ Duration based on type of tracer release' Dispersion Coef ^Jjr Laminar flow in open channels or closed conduits' S K =:10' Cross-sec. area A, m^ 2C1OXI0" Mean velocity ", m h"' 4.WXI0" Dispersion coef. , m^ h^' Retardation I.OOKIÜ*^ I.OOkIÜ*^ Nonequilibrium coef. ^ £.50>:l0-' Mass transfer coef. 9.00x10-' Equilirbrium decay coef , h"' Noequilibrium decay coef. , h"' Equilibrium production coef. T^t Nonequilibrium production coef Table 2. Contaminant release design parameters for the synthetic, forecasted, and predicted data files. Figure 12 is a plot of all seven tracer tests in their originally recorded units. Such a plot shows the range of traveltimes that occurred over a period of one year. It also illustrates the range of concentrations that were measured. Figure 13 is a plot of all seven standardized dimensionless curves developed using Equations 1 and 2 and depicting the similarities between each. From this plot it may be subjectively determined that all dimensionless BTCs appear similar and may be used in combination to develop a single standardized dimensionless curve. Figure 14 is a plot of the quintic-spline fit to all the BTCs. Figure 14 appears to be visually acceptable and is supported by the coefficient of determination f^' = O.^iS 57 • 3.1.1. Forecasting Comparison Results In order to determine the quality of the forecasting methodology, the quintic-spline fit shown in Figure 14 needs to be plotted against one or more of the real data sets. By comparing the fitted model with one or more real data sets, it is possible to determine if the forecasting methodology is a viable method to use. Figure 15 illustrates the quality of the quintic-spline fit to the May 30, 1985 BTC (see Figure 12) and reflects the quality of the regression fits to the traveltime, the traveltime standard deviation, and the peak normalized concentration (Figures 4, 5, and 9). Solute traveltime and spread are a good fit but peak concentration for the quintic spline greatly exceeds the measured BTC. Discrepancies in the transport parameters and concentration parameters reflect the nonperfect regression fits to the measured parameters. For example, the chosen measured data set (May 30,1985) for comparison has a peak normalized concentration (f= [. BT^ rng L"' kg-i) that is slightly less than the regression fit (=2 266 mg L"' kg-') ((J = 116. m' h-') Figure 14. Plot of the quintic-spline fit to the seven composite standardized dimensionless BTCs. (Figures 9 and 15). This discrepancy results from the much greater peak normalized concentration determined by the quintic-spline fit relative to the chosen measured BTC. Using the February 26, 1986 BTC in which a peak normahzed concentration (t.'^ = 2.0 ] 0 mg L"' kg"') is slightly greater than the regression fit in (f = L. S47 mg L"' kg"') (^ = m' h"') (Figures 9 and 16) has the opposite effect. In general, Figures 15 and 16 suggest that the forecasting methodology may not be capable of matching normalized concentrations but may be able to reasonably approximate traveltimes. 3.2. Prediction Process The prediction methodology requires that specific field parameters be measured and applied in the program EHTD. The measured parameters will then be combined in functional relationships (Table 1). However, in the case of longitudinal dispersion , EHTD uses a unique methodology (Chatwin, 1971; Field, 2002a) rather than a published equation which have been found to overestimate (Field, 2002c). 3.2.1. Prediction Comparison Results The prediction methodology may also be evaluated by applying the appropriate measured parameters in the program EHTD. For the May 30,1985 BTC, EHTD used (J - 1 | fi, 2Ä m^ h"', ,V = 9 M m, A = I .'JiS m^ and — 3,57 g all of which conform with original conditions. (Note that A 1 —1—1 ■ 1 ■—■—"—r —1—■—r-1—I—r -■-F-1-1-■-■-■-1-■--- 1—r ■ 1 St ■ 0.0113 Jrt* l" ■ J - h - 1 ] . i.ssol h - : 1 J - 'V 10 20 JO Timf trvn (li) 40 50 Figure 15. Comparison of the quintic-spline fit to the May 30, 1985 BTC where the forecasted peak concentration greatly exceeds the measured peak concentration. Data points represent actual sample-collection times. Identified parameters refer to quintic-spline fit only. was obtained from Figure 11.) For this test, .l/j, = g was fixed so as to force EHTD to use this specific mass (in normal operating mode, EHTD adjusts .l/jj to meet a specific target C' )• No tracer reaction (e.g., retardation), decay, or production was allowed. Applying the measured parameters to EHTD resulted in a visually acceptable fit between the predicted BTC and the actual measured BTC (Figure 17). The larger predicted peak concentration ( = ?. 51 MS L"') produced by EHTD does not appear to be too much greater than the measured peak concentration (C,, =■1.12 MS L"')- For the February 26'! 1986 BTC EHTD used a ^ = 2^3.76 m' h-', ,v = 9 M m, ^ = 2.4^ m^, and 1 /p = 7.1-1 g all of which conform with original conditions. (Note that A was obtained from Figure 11.) For this test, =7.14 g was also fixed so as to force EHTD to use this specific mass. Applying the measured parameters to EHTD resulted in a visually poor fit between the predicted BTC and the actual measured BTC (Figure 18). Specifically, EHTD resulted in a much lower peak concentration (t = B.S I Mg L"') than the measured peak concentration (( ' = 14. H5 l^g L"')-As with the forecasting methodology. Figures 17 and 18 suggest that the prediction methodology may be incapable of matching concentrations but may be able to reasonably approximate traveltimes. Figure 16. Comparison of the quintic-spline fit to the February 26,1986 BTC where the measured peak concentration exceeds the forecasted peak concentration. Data points represent actual sample-collection times. Identified parameters refer to quintic-spline fit only. 4. Hypothetical Contamination Problem Applying the forecasting and prediction methods to a hypothetical pollutant spill serves to illustrate the positive and negative features of the two methods. Consider an accidental spill of 220 L of 5% calcium arsenate solution (Ca3(AsO4)2) into the same karst window as was used in the tracing studies (see Figure 12). This quantity of calcium arsenate solution, a common agricultural chemical, contains 4.14 kg of arsenic. Using a spring discharge equal to 91.80 m3 h-1, the program CXTFIT was run to develop a two-region nonequilibrium synthetic data set for comparison with forecasted and predicted analyses of the release. A two-region nonequilibrium data set was created to impose a relatively long tail on the data set. The degree of nonequilibrium and mass transfer coefficient temporally averaged concentration Cf > and most of the geometric parameters were better represented by the predicted parameter results. From Table 3 and Figures 19 and 20 it would appear that predicting the effects of a release is a slightly more reliable method than the forecasting methodology. In terms of solute transport rates it would definitely seem that the prediction methodology is the better of the two. In addition, the critically important temporally averaged concentration Ci necessary for human health exposure and risk assessments suggests that the prediction methodology may also be the better of the two methods. Temporally integrated exposure Cf , obtained from and temporally average (17) and temporally averaged exposure Cf , obtained from (18) are critically important parameters used to calculate risks because they are more representative of actual exposures than are other measures of contaminant concentration (e.g., C^, C", C/^). Exposure and risk assessments may sometimes use the peak or mean concentration, (C^, C ) or the temporally integrated exposure (Cf ) but all three tend to overestimate the exposure greatly For example, ^ would seem to be a relevant parameter, but the time that an individual might be exposed to the peak concentration for an instantaneous contaminant release becomes vanishingly small. Altematively, the temporally integrated concentration C/^ is really a summation of all concentrations that an exposed individual will contact and/or ingest. Neither C^, C , or Cf^ are reliable parameters for use in an exposure/risk assessment. The temporally average concentration Cf appears to be an overly small concentration (Figures 19 and 20 and Table 3) but this concentration is aVuch more appropriate measure of exposure than Figure 19. Forecasted BTC for the hypothetical release of220 L of a 5% solution of calcium arsenate resulting in 4.14 kg of arsenic. Identified parameters refer to quintic-spline fit only; temporally average concentrations refer to synthetic data (dotted line) and quintic-spline fit (dashed line). Figure 20. Predicted BTC for the hypothetical release of220 L of a 5% solution of calcium arsenate resulting in 4.14 kg of arsenic. Identified parameters refer to predicted fit only; temporally average concentrations C- refer to synthetic data (dotted line) and predicted fit (dashed line). In Parameter Synthetic Forecasted Predicted Solute Transport Parameters First arrival , min 408.00 35.63 421.95 Peak arrival fp, h 16.80 19.19 18.14 Mean traveltime f , h |[(.S2±4.I9 Jü.WlJ.Sö Final arrival , h 59.80 55.70 55.20 Mean velocity u , m h"' 4S.59iHJ.]4 4S.S6±7..'i6- Dispersion coef. , m^ h"' 502.96 345.98 500.29 Peclet Number P 88.34 115.70 89.30 Geometrical Parameters Volume y , m^ 1 I.'JSKIO' Cross-sec. area A , m^ 1.89 2.10 1.88 Surface area A , m^ I.OTX.IO'^ Sorption coef. , m 0.00 0.00 0.00 Conduit diameter fi , m 1.55 1.63 1.55 Hydraulic head loss // ^ , m 1.7-IkIO"'" l-IJ^IO"'" Friction factor 0.44 0.41 0.44 Viscous-flow sublayer J'j, mm 4.20 4.80 4.17 Solute Mass Parameters First arrival , mg L"' 1.0 J K 10"^ I.67KIO"- ^.(i-JKIO"* Peak concentration C^, mg L"' 5.66 9.50 6.59 Mean concentration , mg L"' 4.15 6.53 6.43 Mass recovered I/,, kg 4.14 4.57 4.14 Percent recovery 100.01a 110.44 100.05b Temporally integrated conc , mg L"' 45.10 49.80 45.12 Temporally averaged conc C/ , mg h"' fl.15x10-' a The percent solute mass recovered for the synthetic data exceeded 100% because the actual mass recovered equaled 4.1402 kg. b The percent solute mass recovered for the predicted data exceeded 100% because the actual mass recovered equaled 4.1420 kg.(18) is , C , or • This is so because Cf is scaled by the time of exposure to be a more realistic measure of that part of the BTC that an individual is actually likely to be exposed through contact and/or ingestion of the contaminated water. 5. CONCLUSIONS Forecasting and/or predicting the results of contaminant releases prior to occurrence can be beneficial to water managers of karstic aquifers with relatively simple solution-conduit systems by giving them the opportunity to properly prepare appropriate responses to a release. Forecasting the results requires that water managers take a very proactive stance by thoroughly investigating the transport properties of their aquifer. Results obtained from the investigations may be readily combined and analyzed to provide a measure of what to expect should a real release occur. Predicting the results of contamination is less of a proactive measure, but still represents a means by which water managers may develop a means of protecting their aquifer by examining possible contaminant release outcomes. Results obtained using the prediction methodology can provide water managers with the information necessary to develop appropriate protection strategies prior to a contaminant release occurring. A comparison of the forecasting and prediction methodologies to a synthetic BTC demonstrated the effectiveness of the two methods. In general, the differences between the parameter results for the synthetic data and the forecasting and prediction methods are relatively small and probably does not warrant a specification as to whether the forecasting or prediction methodologies is most reliable. If the results of multiple tracer tests are available or the time to conduct them can be taken then the forecasting method may be an appropriate method to choose provided it is recognized that that this method may produce some significantly erroneous results. In lieu of multiple tracer-test results or the time to acquire them, the prediction methodology may well suffice although with less certainty than might be desired. In most instances, the prediction methodology will likely yield better results provided the basic field measurements have been taken and a Type V karst network is not in evidence. It may also be expected that as aquifer complexity increases, the degree of uncertainty will increase as well for both the forecasting and prediction methodologies. NOTATION A spring cross-sectional area (L2) J. surface area (L2) h one half fracture width (L) P dimensionless partition coefficient for mobile- and immobile-fluid regions for nonequilibrium conditions C solute concentration (M L-3) characteristic concentration for dimensionless parameters (M L-3) first arrival solute concentration (M L-3) C mean solute concentration (M L-3) peak normalized concentration (M L-1 kg-1) peak solute concentration (M L-3) C reduced volume-averaged solute concentration (M L-3) standardized concentration (dimen.) Cf^ temporally averaged concentration (mg L-1 h-1) temporally integrated concentration (mg L-1) fj solution conduit diameter (L) diffusion coefficient (L2 T-1) axial dispersion (L2 T-1) friction factor (dimen.) s, input concentrations for pulse injection; i = Kif ~ 3 ~ function of values such that^^tr ^-j 7 zero-order production term (M L-3 T-1) dimensionless exponential production (growth) constants for the PVP ■ ~ '' ^^ Vj zero-order production term for the liquid phase (M L-3 T-1) v.. zero-order production term for the adsorbed phase (M L-3 T-1) r'- . , , ../'My + pJiVj dimensionless production — ^ ^ V. mean traveltime skewness coefficient (dimen.) r/ auxiliary functions for equilibrium transport [see Appendix] H solute-migration zone thickness (L) fff hydraulic head loss (L) fracture and/or solution conduit distribution coefficient (L) Kf porous media solute distribution coefficient (L3 M-1) I. characteristic distance (L) solution conduit length (L) peak normalized load (mg s-1 kg-1) solute mass released (M) solute mass recovered (M) effective porosity (dimen.) "«J multiplier for estimating tracer test duration (dimen.) P. Peclet number (dimen.) Ph bulk density (M L-1) spring discharge (L3T-1) r solution conduit radius (L) solute retardation (dimen.) R slope of the water surface or energy gradient (dimen.) sinuosity factor (dimen.) mean traveltime standard deviation (T) t time (T) 1 mean residence time (T) solute release time (T) first measured arrival time (T) last measured arrival time (T) peak time of arrival (T) r dimensionless time — f dimensionless mean residence dimensionless pulse time ■= ^ ; = In ^ = öj u mean flow velocity (L T-1) shear velocity (L T-1) p solute decay (T-1) dimensionless equilibrium decay _ L [ i}Ji, + A'jA.), (+ ^ ). (+ "^^Ji, > P/ liquid phase solute decay (T-1) Jt, sorbed phase solute decay (T-1) Pi dimensionless decay coefficient for equilibrium conditions Pi dimensionless decay coefficient for nonequilibrium conditions V volume (L3) viscous-flow sublayer (L) jr width of solute-migration zone (L) dimensionless mass transfer coefficient for nonequilibrium conditions .v distance (L) J: angle between the hydraulic gradient azimuth and direction of either of two joint sets (dimen.) dimensionless distance — L APPENDIX: DEFINITIONS FOR THE AUXILIARY FUNCTIONS r,- (A1) (A2) (A3) I — urft ■J I +— 1 I'T ^aFTJJ^ (A4) (A5) (A6) REFERENCES Atkinson, T. C., Smith, D. I., Lavis, J. J., Whitaker, R. J., 1973. Experiments in tracing underground waters in limestones. J. Hydrol. 19, 323-349. Brouyere, S., 2003. A quatitative point of view of the concept of vulnerability. In: Zwahlen, F. (Ed.), Cost Action 620: Vulnerability and Risk Mapping for the Protection of Carbonate (Karst) Aquifers; Final Report. European Union, Switzerland, pp. 10-15. Chatwin, P C., 1971. On the interpretation of some longitudinal dispersion experiments. J. Fluid Mech. 48(4), 689-702. Eckstein, Y, Xu, M., 1995. Use of weighted least-squares method in evaluation of the relationship between dispersivity and field scale. Ground Water 33(6), 905-908. Field, M. S., 2002a. Efficient hydrologic tracer-test design for tracer-mass estimation and sample-collection frequency. 1. Method development. Environmental Geology 42(7), 827-838. Field, M. S., 2002b. A Lexicon of Cave and Karst Terminology with Special Reference to Environmental Karst Hydrology. Tech. Rep. EPA/600/R-02/003, 214 pp., U. S. Envir. Prot. Agency, Washington, D.C. Field, M. S., 2002c. The QTRACER2 Program for Tracer-Breakthrough Curve Analysis for Hydrologi-cal Tracer Tests. Tech. Rep. EPA/600/R-02/001, 179 pp., U.S. Envir. Prot. Agency, Washington, D.C. Field, M. S., 2003. Tracer-Test Planning Using the Efficient Hydrologic Tracer-Test Design (EHTD) Program. Tech. Rep. EPA/600/R-03/034, 175 pp., U.S. Envir. Prot. Agency, Washington, D.C. Field, M. S., 2004. Assessing the risks to drinking-water supplies from terrorists attacks. In: Mays, L. W. (Ed.), Water Supply Systems Security. McGraw-Hill, New York, pp. 6.1-6.26. Field, M. S., Pinsky, P F., 2000. A two-region nonequilibrium model for solute transport in solution conduits in karstic aquifers. J. Contam. Hydrol. 44, 329-351. Fisher, H. G., List, E. J., Koh, R C. Y, Imberger, J., Brooks, N. H., 1979. Mixing in Inland Coastal Waters, 483 pp. Academic Press, Inc., New York, N.Y. Freeze, R. A., Cherry, J. A., 1979. Groundwater, 604 pp. Prentice-Hall, Engelwood Cliffs, N.J. Gaspar, E., 1987. Modern Trends in Tracer Hydrology, 137 pp. Vol. II. CRC Press, Boca Raton, Fla. Grenville, T. N. E., 1967. Introduction to spline functions. In: Theory and Applications of Spline Functions. Academic Press, New York. Kahaner, D. K., Moler, C., Nash, S., 1989. Numerical Methods and Software, 495 pp. Prentice-Hall, Inc., Englewood Cliffs, N.J. Kashefipour, S. M., Falconer, R. A., 2002. Longitudinal dispersion coefficients in natural channels. Water Research 36, 1596-1608. Käss, W., 1998. Tracing Technique in Geohydrology, 581 pp. A.A. Balkema, Rotterdam, The Netherlands. Kilpatrick, F. A., Taylor, K. R., 1986. Generalization and applications of tracer dispersion data. Water Resour. Bull. 22(4), 537-548. McQuivey, R. S., Keefer, T. N., 1974. Simple method for predicting dispersion in streams. J. Environ. Eng. Div., ASCE 100, 997-1011. Mull, D. S., Liebermann, T. D., Smoot, J. L., Woosley, Jr., L. H., 1988. Application of Dye-Tracing Techniques for Determining Solute-Transport Characteristics of Ground Water in Karst Terranes. Tech. Rep. EPA/904/9-88-001, 103 pp., U.S. Envir. Prot. Agency, Region IV, Atlanta, Ga. Mull, D. S., Smoot, J. L., 1986. Ground water flow characteristics described by quantitative dye tracing in karst terrane in the Elizabethtown area, Kentucky. In: Crawford, N. C., Quinlan, J. F. (Eds.), Proceedings of the Environmental Problems in Karst Terranes and Their Solutions Conference. National Water Well Association, Bowling Green, Ky., pp. 407-421. Parriaux, A., Liskay, M., Müller, I., della Valle, G., 1988. Guide Pratique Pour L'usage des Traceurs Artificiels en Hydrogeolgie, 51 pp. Societe Geologique Suisse, Groupe des Hydrogeologues, GEOLEP EPFL, Lausane, Switzerland. Smoot, J. L., Mull, D. S., Liebermann, T. D., 1987. Quantitative dye tracing techniques for describing the solute transport characteristics of ground-water flow in karst terranes. In: Beck, B. F., Wilson, W. L. (Eds.), Karst Hydrogeology: Engineering and Environmental Applications; Proceedings of the Second Multidisciplinary Conference on Sinkholes and the Environmental Impacts of Karst. Florida Sinkhole Research Institute, College of Engineering, University of Central Florida, Orlando, Fl., Orlando, Fla., pp. 269-275. Sweeting, M. M., 1973. Karst Landforms, 362 pp. Columbia University Press, New York. Taylor, G., 1953. Dispersion of soluble matter in solvent flowing slowly through a tube. Proc. Royal Society of London Ser. A 219, 186-203. Taylor, G., 1954. Conditions under which dispersion of a solute in a stream of solvent can be used to measure molecular diffusion. Proc. Royal Society of London Ser. A 225, 473-477. Taylor, K. R., James, R. W. J., Helinsky, B. M., 1986. Traveltime and Dispersion in the Shenandoah River and its Tributaries, Waynesboro, Virginia to Harpers Ferry, West Virginia. Tech. Rep. WRI 86-4065, 60 pp., U.S. Geological Survey, Towson, Md. Toprak, Z. F., §en, Z., Savi, M. E., 2004. Comment on "longitudinal dispersion coefficients in natural channels". Water Research 38, 3139-3143. Toride, N., Leij, F. J., van Genuchten, M. T., 1995. The CXTFIT code for estimating transport parameters from the laboratory or field tracer experiments; Version 2.0. Tech. Rep. 137, 121 pp., U.S. Salinity Lab., Riverside, Calif. Webster, D. S., Proctor, J. F., Marine, I. W., 1970. Two-Well Tracer Test in Fractured Crystalline Rock, 22 pp. U.S. Geological Survey Water-Supply Paper 1544-I, Washington, D.C. Weston, S., 2002. An introduction to the mathematics and construction of splines. Tech. rep., Ad-dix Software Consultancy Limited, http://www. addix. com/products/spline/ SplineMaths Overview.pdf [accessed 21 May 2004]. Worthington, S. R. H., 1991. Karst Hydrogeology of the Canadian Rocky Mountains. Ph.D. Dissertation, 227 pp., McMaster University, Hamilton, Ontario, Canada. The editors apologize to the author and readers that the outward form of the article is not such as expected, due to technical problems.