- About this Journal ·
- Abstracting and Indexing ·
- Advance Access ·
- Aims and Scope ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Advances in Civil Engineering
Volume 2011 (2011), Article ID 292904, 10 pages
Simulation of Vertical Plane Turbulent Jet in Shallow Water
1Department of Civil Engineering, Clemson University, 110 Lowry Hall, Clemson, SC 29634-0911, USA
2Department of Civil Engineering, Clemson University, 218 Lowry Hall, Clemson, SC 29634-0911, USA
Received 11 July 2011; Revised 25 August 2011; Accepted 1 September 2011
Academic Editor: Jianqiao Ye
Copyright © 2011 T. N. Aziz and A. A. Khan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
A plane, turbulent, nonbuoyant, vertical jet in shallow water is simulated numerically using a three-dimensional computation model employing standard and renormalized group turbulent closure schemes. Existing data of mean and turbulent flow quantities, measured using laser Doppler velocimeter, are used to assess the two turbulent closure schemes. Comparisons between the measured and simulated flow field data are made in the free jet region, within the zone of surface impingement, and in the zone of horizontal jets at the surface. The results show that the standard scheme performs equally well and in some areas better than the more complicated renormalized group scheme in simulating the mean and turbulent flow quantities in this case.
Turbulent jets have been used extensively for discharging effluent into rivers, lakes, and coastal environments. Detailed theoretical and experimental studies of free, turbulent jets (circular or plane jets in an infinite extent of ambient fluid) have been reported by Albertson et al. , Abramovich , Heskestad , and Rajaratnam . The behavior of free, turbulent jets is studied comprehensively using physical models, computational models, and analytical approaches. However, the behavior of jets in a confined environment due to the presence of a free surface and/or solid boundary has been studied only experimentally or numerically.
Gutmark et al.  conducted an experimental study to investigate the behavior of a plane turbulent jet impinging on a solid boundary. The boundary was located at a distance of 100-times the nozzle width. The mean and turbulent flow properties were measured using a hot wire anemometer. The structure of vortices was also reported. More recently, Maurel and Solliec  experimentally studied impingement of plane turbulent jet on a flat plate located at varying distance from the nozzle. Laser Doppler velocimetry and particle image velocimetry were used to investigate flow properties.
Kuang et al.  made a detailed measurement of mean and turbulent flow quantities using laser Doppler velocimeter in a plane jet issuing vertically upward in a still ambient fluid of finite depth. Based on the measured data, four distinct zones were identified and included zone of flow establishment (ZFE), zone of established flow (ZEF), zone of surface impingement (ZSI), and zone of horizontal jet (ZHJ). The configuration and coordinate directions adopted by Kuang et al.  along with different zones are shown in Figure 1. The nozzle width is and the uniform velocity at the nozzle is given by . As the shear layer develops the core of the uniform velocity diminishes in width. At a distance of from the nozzle, the shear layer penetrates up to the centerline of the jet . This distance is called the potential core. After the potential core the centerline velocity decays as the jet continues to grow. Outside the potential core, the velocity at the centerline of the jet is denoted by , and is the vertical component of the velocity at any point in the domain. The horizontal component of the velocity at any point is given by .
The free jet regions (ZFE and ZEF) in cases of plane jet impingement on solid surface (Case I) and plane jet issuing vertically upward in a still ambient fluid of finite depth (Case II) are similar. However in ZEF for Case II, Kuang et al.  found that the centerline velocities were lower than the free jet and the difference was attributed to the proximity of the water surface. While Beltaos and Rajaratnam  showed that centerline velocity for Case I was as given by the free jet. Using the centerline velocity data of Kuang et al.  and Maurel and Solliec  in ZSI for Case II and Case I, respectively, it was found that although the extent of this zone was the same, the centerline velocities were lower in Case I. The stagnation pressure in Case I manifests itself in the form of surface swell in Case II. The stagnation pressure or the surface swell helps turn the flow from vertical to horizontal direction. The maximum horizontal velocity increases within the zone of surface impingement and then decreases as the horizontal jet expands within ZHJ. In addition to the obvious difference in the velocity profile shape in ZHJ for the two cases (due to the presence of a wall in Case I), the local maximum velocities in ZHJ were found to be lower for Case I (measured by Wygnanski et al. ) than Case II (measured by Kuang et al. ). Analysis of the Reynolds stress data from the two cases also showed differences. Thus, the various zones that form in these two cases, although similar, have subtle differences in the flow characteristics and a detail investigation of jet issuing vertically in shallow depth is warranted.
Due to the difficulty of obtaining an analytical solution, especially in the zones of surface impingement and horizontal jets, the problem of a jet issuing vertically upward in shallow water is studied using physical models. Numerical models have been adopted increasingly to study complex flow problems. To numerically model turbulent flows, the selection of an appropriate turbulence scheme is essential. The two most common turbulence closure schemes are the and the renormalized group (RNG). For free plane turbulent jets the two schemes provide similar results in the zone of established flow with RNG providing better results for the turbulent kinetic energy along the center of the jet . The aim of this study is to assess the applicability of these two turbulence closure schemes in predicting the mean and turbulent flow characteristics of a jet issuing vertically upward in shallow water. The zones of surface impingement (ZSI) and horizontal jets (ZHJ), where both mean and turbulent flow properties are changing rapidly, are of main interest. The computed growth rates, mean flow velocity profiles, and turbulent kinetic energy profiles in different zones of the jet issuing vertically in shallow water are compared with the measured data acquired by Kuang et al. .
2. Mathematical Details
The mass and momentum equations, which constitute the Reynolds-averaged Navier-Stokes (RANS) equations, are given by where represents coordinate directions ( to 3 for directions, resp.), is the time-averaged velocity component ( in directions, resp.), represents time, is the fluid density, is the piezometric pressure, is the kinematic viscosity of the fluid, are turbulent normal and shear stresses, and and are indices which range from 1 to 3.
A constitutive relationship is required to relate the turbulent normal and shear stresses to mean flow velocities and is achieved through the use of turbulent eddy viscosity. The relationship is given by where is the turbulent eddy viscosity, is the Kronecker delta, and is the turbulent kinetic energy per unit mass. The turbulent eddy viscosity is computed using turbulent kinetic energy per unit mass and the dissipation rate of turbulent energy per unit mass () as follows where is an empirical coefficient. To determine the turbulent eddy viscosity and hence the turbulent normal and shear stresses, equations for determining the turbulent kinetic energy () and dissipation rate () per unit fluid mass are required.
In the case of the standard scheme, the and are determined from the following two equations where , , , and are empirical coefficients. The standard values of the coefficients , , , , and used in the scheme are 0.09, 1.44, 1.92, 1.0, and 1.3, respectively.
In the scheme the eddy viscosity is determined from a single turbulence length scale. This means that the turbulent diffusion is restricted to only that length scale for which eddy viscosity is determined. Eddies in a turbulent flow contain a large range of length scales and the turbulent diffusion occurs at all scales. The energy cascades down from larger to smaller size eddies where it is dissipated by viscosity. At some point in this size scale, the energy dissipation is equal to the energy production and is called the inertial scale. This scale is of an order that can be efficiently handled by current computer technology. The renormalization technique uses the inertial scale to describe all the other scales and produces a model that is statistically equivalent to the original Navier-Stokes equations. The RNG scheme uses different equations for determining the and . These equations, as given by Bischof et al. , are where and are Prandtl numbers for and , respectively, represents the mean rate of strain, and is given by The term on the right-hand side of (6) constitutes a major difference between the RNG and schemes. The term is a source term that is a function of , , and . This term is given by where and are constants having standard values of 4.38 and 0.012, respectively, and . The standard values of constants , , and used in the RNG schemes are 0.0845, 1.42, and 1.68, respectively.
The and RNG schemes have been used extensively with the RNG model used for flow involving rapid strain (strong local changes) and low Reynolds number flows. The RNG scheme is found to be more sensitive to the rate of strain because of the presence of the source term . Since renormalized group scheme is favored for locally varying flows with rapid strain, a situation that exists within the zone of surface impingent (ZSI) and zone of horizontal jets (ZHJ), it will be of interest to compare the results from the two turbulence schemes within these regions.
3. Experimental Setup
In this study, the experimental data collected by Kuang et al.  was used to evaluate the accuracy of the two turbulence closure schemes. The layout of a jet issuing vertically in shallow water is shown in Figure 1. The tank used for experimental setup was 6.5 m long (-direction), 0.4 m wide (-direction), and 0.4 m high (-direction). The nozzle was cut across the width of the tank (along -direction). The width of the nozzle () was 2 mm, and depth of water above the nozzle () was maintained at 26 cm (with value of 130). The resulting jet was two-dimensional in the - plane, that is, the variations of the flow field in the -direction were negligible. Three experiments were conducted at nozzle velocity () of 0.90, 1.33, and 1.75 m/s corresponding to the nozzle Reynolds number of 1,560, 2,340, and 3,120, respectively. The measurements of mean and turbulent flow characteristics were made using laser Doppler velocimeter. By choosing appropriate scales for nondimensionalizing, the measurement results were found to be independent of Reynolds number. Complete details regarding experimental setup and results were given by Kuang et al. . The case with the nozzle velocity of 1.33 m/s was considered for comparison with the numerical model results and evaluation of the two turbulent schemes.
4. Numerical Model
A computational fluid dynamics model, FLOW-3D, developed by Flow Science, Inc., is used to numerically solve the RANS equations. The Flow-3D software provides several options for the turbulent closure scheme including Prandtl mixing length model, two-equation models, and large eddy simulation model. In this study, and RNG turbulent closure schemes, with standard coefficients as described earlier, are used. In the Flow-3D model, the Reynolds-averaged Navier-Stokes equations along with the turbulent closure scheme (if needed) are solved using a finite volume/finite difference method in an Eulerian rectangular or cylindrical grid (orthogonal mesh). The model can solve the equations explicitly or implicitly and the momentum advection terms can be approximated using first-order or second-order approximation. In this study, the explicit scheme with second-order accuracy is adopted. The “saw-tooth” representation of the boundaries is eliminated through the use of the fractional area/volume obstacle representation (FAVOR) method . In the FAVOR method, the boundaries are determined independent of the grid generation process. The position of the free surface is predicted using a modified volume of fluid method as part of the solution.
The geometry of the experimental setup, as given by Kuang et al. , was adopted for numerical simulation. A boundary condition in terms of velocity of 1.33 m/s was specified at the nozzle. The water depth was maintained at 26 cm above the nozzle through the use of a continuative boundary condition at the tank boundaries in the -direction. The continuative boundary condition kept the water level the same while allowing flow out of the tank by specifying zero velocity gradients at the tank boundaries.
To model the scenario as two-dimensional flow in the FLOW-3D model, only one cell needed to be specified in the -direction. The specification of one cell in the -direction allowed for 2D computation in the vertical plane. Full 3D computation was also performed using a course mesh and the results were compared with the corresponding 2D computation. The results from the 2D and 3D computations were the same. Since staggered grid was used for computation, where velocities were computed on the cell faces and all other quantities were calculated at the center of the cell, minimum of two cells were needed to determine velocity at the center of the jet. A symmetric mesh was utilized in the -direction on either side of the nozzle. The mesh for the right half of the domain is shown in Figure 2. In the -direction, starting with the cell size of 0.1 cm within the nozzle, the cell size was increased geometrically by a factor of 1.2 up to a distance of 35 cm from the center of the nozzle. After which a constant cell size of 5 cm was used up to a distance of 145 cm. The last segment (between 145 cm and 325 cm) had a cell size of 10 cm. It should be mentioned that the region of interest in terms of comparing computational results with the measured data was within 15 cm on either side of the nozzle. The computational mesh extent in the vertical direction was set to 35 cm (the water depth was set to 26 cm) to account for the bulge at the water surface. The size of the first cell in the vertical direction was set as 0.1 cm. The cell size was increased geometrically by a factor of 1.2 up to a distance of 11 cm. From then on, a constant cell size of 1.5 cm was used. The total numbers of cells in the - and -directions were 146 and 43, respectively. Several other mesh sizes were tried to ascertain that the computational results would be mesh independent. It was found that the results obtained using the current mesh and that obtained by doubling the number of cells had difference of less than 1% (based on the maximum velocity in the ZEF, ZSI, and ZHJ). The computational time step is automatically selected by the model based on stability criterion.
5. Comparison of Results
The computational results obtained using the and RNG turbulent schemes are compared to the measured data  for the nozzle Reynolds number of 2,340. Kuang et al.  determined, by analyzing the variation of longitudinal turbulent intensity along the centerline of the vertical jet, that the zone of flow establishment (ZFE) extended up to of 25 and the zone of established flow (ZEF) ranged from of 25 to 100 (analysis would be shown later). The virtual origin of the jet and the growth rate within the zone of established flow (ZEF) from measured and simulated data are compared in Table 1. The measured and simulated results show a linear growth of the jet. The growth rates are obtained using the longitudinal velocity () profiles across the jet at various locations and determining the position, , from the centerline where the velocity is half the centerline velocity. By fitting a straight line of the form to the data, the growth rate, , and the virtual origin, , are determined. The growth rate of the jet using the scheme compares well with the measured data while the RNG scheme shows a higher growth rate for the jet. The scheme provides virtual origin that is in front of the jet, similar to the measured data. The virtual origin of the jet simulated using the RNG scheme is behind the nozzle. The higher growth rate and the location of the virtual origin in the case of the RNG scheme are the result of wider shear layer predicted by the RNG scheme. This is the direct result of a higher sensitivity of the eddy viscosity to the strain rate.
The variation of longitudinal velocity () and turbulent kinetic energy per unit mass () along the centerline of the jet are shown in Figures 3 and 4, respectively. Although both schemes estimate higher growth rate for the jet in ZEF, the longitudinal velocity along the centerline of the jet is overpredicted by the two schemes when compared to the measured data. The centerline velocity predicted by the RNG scheme is more accurate than that predicted by the scheme. The turbulent velocity fluctuations along the center line of the jet, as shown in Figure 4, show that both schemes predict the ZEF within the range of . However, within this zone, the turbulent velocity fluctuations are overpredicted by both schemes (the scheme providing relatively better accuracy). As shown by the measured data and the two schemes, the turbulent velocity fluctuations increase rapidly beyond of 100 (the zone of surface impingement).
The variations of longitudinal () velocity profiles across the width of the jet at different locations along the jet are shown in Figures 5–7 (the data for the RNG scheme is offset by 0.5). In each case, the predicted velocity profiles are compared with the theoretical curve and/or with a typical measured velocity profile within that zone. The simulated longitudinal velocity profiles in the ZFE are compared to the measured data and a theoretical curve for ZEF in Figure 5. The theoretical curve that fits the measured data  in the zone of established flow is given by where , as determined by the measured data, is used. The simulated velocity profiles using the scheme agree well with a typical measured profile in the ZFE and are below the theoretical line. The velocities are overpredicted by the RNG scheme within the ZFE when compared to the typical measured velocity profile and approach the theoretical curve with the increasing longitudinal distance from the nozzle. The simulated velocity profiles within the ZEF, using the and RNG schemes, are shown in Figure 6. The measured velocity profiles in this region, as shown by Kuang et al. , agree well with (9). This equation is shown in the figure for comparison with the simulated data. The velocity profiles, predicted using the scheme, compare well with the theoretical curve while the velocities are overpredicted by the RNG scheme (the error increasing with increasing distance). The higher shear rate and higher turbulence energy predicted by the RNG scheme result in an overprediction of the vertical velocities and growth rate of the jet in the ZFE and ZEF. In the ZSI, as shown in Figure 7, velocities are underpredicted by both schemes. However, in the case of the RNG scheme, the shape of the velocity profile and the magnitude of velocities near the surface () are not accurately predicted.
The simulated results of the lateral velocity profiles across the width of the jet are compared to the measured and/or theoretical velocity profiles in Figures 8–10 (the RNG results in ZEF and ZFE are offset by 0.2, and by 2.0 in ZSI). Since the velocity profiles are symmetrical about the center of the jet, the following description considers only the left side of the jet. In the ZFE (Figure 8), both schemes underpredict the velocity; however, the accuracy gets better with the increase in the longitudinal distance from the nozzle. The lateral velocity at the center of the jet should be zero based on the symmetry of the jet. However, the measured data show some discrepancy near the center of the jet. In the figures, the theoretical curve is given by  where is used. In the ZEF zone (Figure 9), the lateral velocity profiles predicted by the scheme agree well with (10), while in the case of the RNG scheme, velocity profiles are progressively underpredicted with increasing longitudinal distance.
The zone of surface impingement is where the flow turning and acceleration take place. Starting with zero horizontal velocity, the velocity increases within this zone with horizontal distance. At the same time the vertical velocity diminishes within this zone. This area is associated with high strain rate and, as shown in Figure 10, the horizontal velocity profiles within this area are predicted more accurately by the RNG scheme, especially for of 110 and 120.
The predicted turbulent kinetic energy (TKE) profiles across the jet at different locations along the jet are shown in Figures 11–13 (the RNG results in the ZFE and ZEF are offset by 0.05). In the ZFE (Figure 11), the scheme accurately predicts the shape and magnitude of the TKE; however, the magnitude of the TKE is overpredicted outside the zone of , especially close to the nozzle (). The RNG scheme overpredicts by 30% the TKE near the center of the jet for . In the ZEF, as shown in Figure 12, the TKE profiles predicted by the scheme are self-similar, a property shown by the measured data of Kuang et al.  as well. The TKE predicted by the scheme compares well with the measured data in the center part of the jet (), and the TKE is underpredicted outside of this zone. The TKE profiles predicted by the RNG scheme do not show similarity. The TKE is overpredicted near the center of the jet for low values of . This may be the result of higher sensitivity of the TKE production on the velocity gradients in the case of the RNG scheme. In the zone of surface impingement, as shown in Figure 13 for , the TKE predicted by both schemes agrees well with the measured data near the center of the jet; however, the values are underpredicted in the outer part of the jet.
The variations of lateral velocity in the zone of horizontal jets (ZHJ) predicted by the and RNG schemes are shown in Figures 14 and 15, respectively. Both schemes predict maximum velocity at the surface and similarity of the velocity profiles in the positive flow region. The results obtained using the scheme show a better agreement with the measured data and a Gaussian fit. Kuang at al.  provided a Gaussian fit for approximating the measured data in the positive flow region and is given by where is the maximum lateral velocity. The recirculation zone underneath the surface jets is poorly resolved by the RNG scheme and may be the results of high growth rate of the jet in the ZEF predicted by the scheme. The large width of the jet will interfere with the recirculation flow beneath the surface jets. Kuang et al.  reported that the recirculation zone covered the whole depth and its width was about to . The recirculation zone predicted using the two schemes extended the full depth of the tank and had a length of approximately . The recirculation zone predicted using the RNG scheme, is shown in Figure 16. A bulge in the water surface, due to the interaction of the jet with the free surface, is clearly visible in the figure. The variations of the maximum lateral velocity () for both schemes are shown in Figure 17 along with the fit used by Kuang et al.  to describe the measured data. The fit is given by The predicted profiles of maximum lateral velocity show oscillations across the linear trend line and are similar to oscillations observed in the measured data.
A plane turbulent jet issuing vertically upward in shallow water depth has been simulated using a three-dimensional computational model employing the and RNG turbulent closure schemes. To evaluate the performance of these two schemes, the simulated results from these two schemes using standard coefficients are compared with the data measured by Kuang et al. . The centerline vertical velocities are better predicted by the RNG scheme. The turbulent kinetic energy in ZFE, ZEF, and ZSI is overpredicted by both schemes, with scheme providing better agreement with the measured data. The growth rate in the ZEF predicted using the RNG scheme is higher and the location of the virtual origin is behind the nozzle. This is the result of a wider shear layer predicted by the RNG scheme. The growth rate is predicted accurately by the scheme and the location of the virtual origin is in front of the nozzle as indicted by measured data.
In the ZFE and ZEF, the vertical velocities are overpredicted by the RNG scheme. The lateral velocities are underpredicted by both schemes in the ZFE with RNG scheme providing a better fit. While in the ZEF, the RNG scheme progressively underpredicts the lateral velocity with an increase in distance from the nozzle. The TKE is overpredicted (by about 30%) near the center of the jet by the RNG scheme in ZFE. The TKE profiles predicted using the RNG scheme do not preserve self-similarity shown by the measured data and the values are overpredicted near the center. The scheme overpredicts TKE at the outer edges of the jet in ZFE. It preserves the similarity of TKE profiles in ZEF; however, the values are underpredicted at the edges of the jet.
The onset of ZSI, based on the centerline variation of TKE, is accurately predicted by both schemes. The vertical velocities are underpredicted by both schemes in this zone, with RNG unable to accurately predict the shape of the profile near the surface. The lateral velocities are predicted more accurately by the RNG scheme in ZSI. The TKE is predicted accurately near the center of the jet and underpredicted at the outer edges by both schemes in this zone.
In the zone of horizontal jets, the lateral velocity in the forward direction and the recirculation zone underneath the surface jet are better predicted by the scheme. The larger jet size predicted by the RNG scheme in ZEF may be the reason for poor agreement of the velocity data in the recirculation zone. The decay of the maximum velocity is predicted accurately by both schemes. The scheme, albeit simple, provides more accurate results of the flow field for this problem.
|:||A measure of the jet growth|
|:||Growth rate of the jet|
|:||Width of the nozzle|
|:||Water depth in the tank (26 cm)|
|:||Turbulent kinetic energy per unit mass|
|:||Prandtl number for|
|:||Prandtl number for|
|:||Source term in RNG scheme|
|:||Mean rate of strain|
|:||Time-averaged velocity components|
|:||Velocity in the -direction at any point in the jet|
|:||Maximum lateral velocity in the zone of horizontal jets|
|:||Turbulent normal and shear stresses|
|:||Velocity in the -direction at any point|
|:||Centerline velocity in the -direction|
|:||Jet velocity at the nozzle|
|:||Coordinate directions in|
|:||Distance to the virtual origin|
|Kinematic viscosity of the fluid|
|:||Combination of the fluid and turbulent kinematic viscosity|
|:||Turbulent eddy viscosity|
|:||Density of fluid|
|:||Dissipation rate of turbulence per unit mass|
|:||A function of , , and|
- M. L. Albertson, Y. B. Dai, R. A. Jensen, and H. Rouse, “Diffusion of submerged jets,” Transaction American Society of Civil Engineers, vol. 115, pp. 639–697, 1950.
- G. N. Abramovich, The Theory of Turbulent Jets, MIT Press, Cambridge, Mass, USA, 1963.
- G. Heskestad, “Hot-wire measurements in a plane turbulent jet,” Journal of Applied Mechanics, vol. 32, pp. 721–734, 1965.
- N. Rajaratnam, Turbulent Jets, Elsevier, New York, NY, USA, 1976.
- E. Gutmark, M. Wolfshtein, and I. Wygnanski, “The plane turbulent impinging jet,” Journal of Fluid Mechanics, vol. 88, no. 4, pp. 737–756, 1978.
- S. Maurel and C. Solliec, “A turbulent plane jet impinging nearby and far from a flat plate,” Experiments in Fluids, vol. 31, no. 6, pp. 687–696, 2001.
- J. Kuang, C. T. Hsu, and H. Qiu, “Experiments on vertical turbulent plane jets in water of finite depth,” Journal of Engineering Mechanics, vol. 127, no. 1, pp. 18–26, 2001.
- S. Beltaos and N. Rajaratnam, “Plane turbulent impinging jets,” Journal of Hydraulic Research, vol. 11, no. 1, pp. 29–59, 1973.
- I. Wygnanski, Y. Katz, and E. Horev, “On the applicability of various scaling laws to the turbulent wall jet,” The Journal of Fluid Mechanics, vol. 234, pp. 669–690, 1992.
- T. N. Aziz, J. P. Raiford, and A. A. Khan, “Numerical simulation of turbulent jets,” Engineering Applications of Computaitonal Fluid Mechanics, vol. 2, no. 2, pp. 234–243, 2008.
- C. H. Bischof, H. M. Bucker, and A. Rasch, “Sensitivity analysis of turbulence models using automatic differentiation,” Journal on Scientific Computing, vol. 26, no. 2, pp. 510–522, 2005.
- J. F. Rodriguez, F. A. Bombardelli, M. H. Garcia, K. M. Frothingham, B. L. Thoads, and J. D. Abad, “High-resolution numerical simulation of flow through a highly sinuous river reach,” Water Resources Management, vol. 18, no. 3, pp. 177–199, 2003.