Hydrodynamic dispersion characteristics of lateral inflow into a river tested by a laboratory model

Groundwater and river-water have a different composition and interact in and below the riverbed. The riverbed-aquifer flux interactions have received growing interest because of their role in the exchange and transformation of nutrients and pollutants between rivers and the aquifer. In this research our main purpose is to identify the physical processes and characteristics needed for a numerical transport model, which includes the unsaturated recharge zone, the aquifer and the riverbed. In order to investigate such lateral groundwater inflow process, a laboratory J-shaped column experiment was designed. This study determined the transport parameters of the J-shaped column by fitting an analytical solution of the convective-dispersion equation for every flux on individual segments to the observed breakthrough curves of the resident concentration, and by inverse modelling for every flux simultaneously over the entire flow domain. The obtained transport-parameter relation was tested by numerical simulation using HYDRUS 2-D/3-D. Four steady-state flux conditions (i.e. 0.5 cm hr −1, 1 cm hr−1, 1.5 cm hr −1 and 2 cm hr −1) were applied, transport parameters including pore water velocity and dispersivity were determined for both unsaturated and saturated sections along the column. Results showed that under saturated conditions the dispersivity was fairly constant and independent of the flux. In contrast, dispersivity under unsaturated conditions was flux dependent and increased at lower flux. For our porous medium the dispersion coefficient related best to the quotient of the pore water velocity divided by the water content. A simulation model of riverbed-aquifer flux interaction should take this into account. Correspondence to: G. Wyseure (Guido.Wyseure@ees.kuleuven.be)


Introduction
Knowledge of the aquifer and river water interaction is important for understanding the continuum of groundwater and surface water hydrology.The significance of groundwatersurface water interaction is however difficult to quantify (Valett et al., 1994) and is commonly ignored in watermanagement considerations or policies.
Groundwater has different dissolved minerals, contains less oxygen, and has a more constant temperature as compared to river water.On the one hand riverbed-aquifer flux interactions result in specific dissolved minerals from the aquifer moving into the riverbed, on the other hand down penetrating flow from the river moves substances like oxygen and organic matter into the riverbed and the aquifer below.The part of the riverbed subject to exchange of fluxes is called the hyporheic zone, which also acts as an important heat source and sink that affects river water temperatures (Brown et al., 2005;Moore et al., 2005;Cozzetto et al., 2006) and the solubility of oxygen (Ricci and Balsamo, 2000;Hahn, 2006).
Many studies analyzed the river-subsurface interaction by comparing the difference of tracer concentrations between river water and the subsurface flowpaths as reviewed by Marion et al. (2003) and Zaramella et al. (2006).Conceptual models of the river solute advection dispersion model, such as the Transient Storage Model (TSM) (Bencala and Walters, 1983), are widely adopted to analyze and predict the solute exchange between river water and bed sediment in longitudinal sections of rivers (e.g.Gooseff et al., 2003;Lin and Medina, 2003;Jonsson et al., 2003;Cozzetto et al., 2006;Ge and Boufadel, 2006;Kazezyilmaz-Alhan and Medina, 2006;Zaramella et al., 2006).The transient storage and exchange in the bed sediment is assumed to be governed by Published by Copernicus Publications on behalf of the European Geosciences Union.P. Y. Chou and G. Wyseure: Dispersion of lateral inflow into a river flow-induced pressure differences over the riverbed (i.e.advective pumping) (Wörman et al., 2002;Harvey et al., 2003;Rehg et al., 2005).The flow-induced pressure differences can be due to riverbed irregularities and associated waves or flood hydrographs.The waves in river flow generally have small amplitude and high frequency, while flood hydrographs peaks have a rapid increase associated with high pressures and slower recession.Many studies have investigated the stream-subsurface interactions resulting from topographical features of the riverbed (Harvey and Bencala, 1993;Wroblicky et al., 1998;Anderson et al., 2005;Gooseff et al., 2005;Wondzell, 2006;Boano et al., 2007).
In the TSM the physical mechanics of transient storage are conceptually lumped into a single storage zone (Packman and Bencala 2000;Runkel et al., 2003).Kazezyilmaz-Alhan and Medina (2006) assumed in their version of the TSM that the solute concentration in the river and in the storage zone varies only along the longitudinal direction of river.The study of Zaramella et al. (2003) concluded that TSM did not explicitly represent the subsurface.From their perspective of hyporheic exchange of metals Zaramella et al. (2006) pointed out that there is a need to assess local transport within the hyporheic zone.One major limitation of the TSM approach is that the spatial and temporal dynamics of lateral groundwater inflow and the wide range of solute residence times are not considered.
Several studies show that in ephemeral rivers the groundwater flow dominates the river-subsurface flow regimes.Storey et al. (2003) used a three-dimensional groundwater flow model (MODFLOW).They concluded from their model-study that exchange flows are up to twice as strong, but more variable in time, at the sides of the stream than near the centre.They also revealed that vertical flow paths beneath the channel are more persistent under the range of conditions modelled than lateral flow paths into the banks.This study reinforced the need to represent the subsurface in a spatially distributed way.The study of Malcolm et al. (2004) indicated that rapid changes in groundwater-river interactions occurred during hydrological events.They concluded that the differences in concentrations between river and subsurface water increased with depth into the hyporheic zone and that during low flows groundwater lateral flow dominates the groundwater-river interactions.The study of Wondzell (2006) found that hyporheic exchange was little affected by river discharge, but was rather influenced by the hydraulic gradients between the river and the floodplain.Consequently, there is a need for distributed modelling of the flow exchanges below the river bed and in conjunction with the lateral groundwater inflow.
To contribute to a better understanding of the linkage between the river and the surrounding aquifer, the overall goal of this study is to characterize the mass transport by hydrodynamic dispersion from the vadose zone via the hyporheic zone into the river.The specific objective is to determine a relation for the dispersion parameters so that it can be imple-mented in a more field-realistic model which can be applied under a wide range of spatial and temporal variable saturation and fluxes.

Theory
Transport of material conveyed by the water flow in porous media, like in the riverbed and the surrounding aquifer, can be described by the convective-dispersion equation (CDE).The CDE in one dimension is expressed by the following partial differential equation: where D is the dispersion coefficient (L 2 T −1 ); V is the pore water velocity (L T −1 ); C is the concentration of solute (M L −3 ); t is the time (T) and z is the axial distance (L).Equation ( 1) can be solved analytically (e.g.Lindstrom et al., 1976) for simple geometries and numerically for more complex cases in two and three dimensions.In addition this equation can be expanded for mobile-immobile water in the soil (e.g.Šimùnek et al., 2003) and by sink/source terms representing local in situ degradation or generation.The analytical solution of CDE allows a more straightforward and parsimonious estimation for the underlying physical mechanisms for column-like situations.Alternatively, the parameters can be determined by inverse modelling using more comprehensive numerical models such as HYDRUS 2-D/3-D ( Šimùnek et al., 2006), which allows a more complex geometry and soil layering.At this stage it was preferred to consider only homogeneous sand and using conservative solute.
The hydrodynamic dispersion coefficient is often calculated as a combination of mechanical dispersion and molecular diffusion by: ( where λ represents the dispersivity (L); n is an empirical coefficient ranging between 1 and 2; D e is the molecular diffusion (L 2 T −1 ).At usual flow conditions the mechanical dispersion is much higher than the diffusion, therefore the molecular diffusion is often disregarded and for n equals 1 a linear relation between D and V is obtained (Bear, 1972).The dispersion coefficient is primarily influenced by pore water velocity and dispersivity, which is a function of medium characteristics and water content (Padilla et al., 1999;Nützmann et al., 2002;Toride et al., 2003;Costa and Prunty, 2006).The study of Maraqa et al. (1997) reported that the dispersivity of a soil under unsaturated conditions is higher than when the soil is fully or nearly saturated.
A recent review of dispersivity given by Vanderborght and Vereecken (2007) concluded that for the short travel distance (0 to 30 cm), a clear increase in dispersivity with increasing Hydrol.Earth Syst.Sci., 13, 217-228, 2009 www.hydrol-earth-syst-sci.net/13/217/2009/ flow rate was present, however, this increase was not apparent for long travel distance (>30 cm).Moreover they found that the dispersivity increased when the lateral scale of the experiment increased.Their study also discussed the impact of texture and structure, but little information was given on the influence of soil water content.The unsaturated soil water content is often characterized by the soil-water retention curve.One commonly used parameterization is the van Genuchten (1980) curve: where θ(h) is the soil water retention (L 3 L −3 ); θ r and θ s represent the residual and saturated water content (L 3 L −3 ) respectively; α is the inverse of the air-entry value (L −1 ); n* is a pore size distribution index (>1), both values are considered as empirical coefficients affecting the shape of the hydraulic functions; h is the pressure head (L).The hydraulic conductivity in relation to the soil water retention is given by: where K(h) is the unsaturated hydraulic conductivity (L T −1 ); K s is the saturated hydraulic conductivity of soil; S e is the effective water content (L 3 L −3 ).The numerical model HYDRUS 2-D/3-D uses Eqs. ( 3) and (4) to specify the soil hydraulic properties.Correct interpretation of initial and boundary conditions is required for the analysis of tracer experiments.The prescribed concentration or a Dirichlet boundary condition is adopted by measuring a time-dependent input concentration inside the column, provided the flow is fully developed.The dimensionless column Peclet number, a ratio between solute convection and hydrodynamic dispersion (Bear, 1972), for a given column segment with length L, is defined as: At larger column Peclet numbers (>5) the flow and transport is well developed and the choice of analytical solutions linked to boundary conditions is less critical (van Genuchten and Parker, 1984).As a result at sufficiently high column Peclet numbers the electrical conductivity (EC) measured by Time Domain Reflectometry (TDR) can be used as a prescribed concentration at the upper boundary, and it allows the elimination of uncertainty of the nature of the inlet condition (Avila, 2005).Thus the initial and boundary condition can be set as: Upper boundary condition: End boundary condition: whereC i is the initial concentration; C 0 is the given concentration applied to the system; bothC i and C 0 are assumed as constant.Mojid et al. (2004) developed an efficient transfer-function method based on the Wakao and Kaguei (1982) solution.A short description of the transfer-function method is given in the Appendix.The impulse response in the time-domain becomes: where t is total variable time (T); τ is travel time of the tracer, which is determined by dividing the length of travel by V ; R f is the retardation factor (dimensionless); N is massdispersion number (dimensionless), which is the reciprocal of the column Peclet number.The time-dependent estimated normalized response concentration of solute (C [r.est ) can be predicted by convoluting the input with Eq. ( 9).As shown by Mojid et al. (2006) the analysis method used in this study is not very sensitive to the tail of the pulse and the response.The detection of the end-point is therefore not as critical as that of the start point.

Experiment setup
Soil water content and bulk electrical conductivity was monitored simultaneously by Time Domain Reflectometry (TDR) (Wyseure et al., 1997).Three-rod stainless steel probes, with length of 10 cm, 0.2 cm in diameter, spaced 1 cm apart and attached to a 200 cm coaxial cable were used.The soil volume measured by three-rod sensors is roughly a cylinder limited by the outer rods.Six TDR probes were connected to a Tektronix 1502B metallic cable tester via a Campbell Scientific multiplexer for consecutive scanning.During the experiments the EC (S m −1 ) and the soil water content (cm 3 cm −3 ) were continually measured by using the WinTDR-software, Version 6.1, developed by the Soil Physics Group at Utah State University (Jones et al., 2002 andOr et al., 2004).This software was also used to calibrate the sensors for simultaneous EC and water content measurement.The laboratory model aims at enforcing one dimensional flow-lines within the continuum of unsaturated zone, aquifer and river in accordance to a 1-D column approach.The experimental set-up is shown in Fig. 1.This J-shaped column model was assembled by two vertical columns in transparent Perspex (Polymethyl-methacrylate, PMMA) with two 90 • elbow PVC tubes.The inner diameter of the column was 20 cm, and the height of the left and right column was 100 cm and 50 cm respectively.The two 90 • elbow PVC pipes were supported in a frame.
The J-shaped model was filled with dune sand with bulk density of 1.55 g cm −3 .Clean washed dune sand was preferred for this experiment as higher fluxes can be used and more fluxes can be tested within a reasonable timeframe.Additionally, homogeneous sand allows the investigation for physical processes and relationships without the confounding effect of layering.Texture analysis by sieving gave an average of 100.0%, 97.0%, 51.7%, 7.4% and 1.4% pass-rates through the 2, 0.5, 0.25, 0.1, 0.05 mm sieves respectively.
Section 1, on the left, represents the vadose zone, which remained unsaturated.Section 2 represents the riverbed in conjunction with the adjacent aquifer, which contained the dynamic interface between the unsaturated section 1 and the constantly saturated part.A piezometer was inserted into the top of the saturated section and was connected to a flexible tube.The water levels in the piezometer were compared to the levels in Sect.3, which measures head loss along the saturated zone in Sect. 2. The water level in section 3 at the right hand side was kept constant using an overspill that was connected with a flexible tube.
Three TDR probes were inserted in the Sect. 1 and three in Sect. 2. The TDR probes were numbered from 1 to 6, starting from the top of the unsaturated section and were separated by distances of 20 cm, 20 cm, 20 cm, 30 cm and 40 cm respectively.Each 10 min, a cycle of consecutive measurements for the TDR-probes at all locations and data storage was performed.

Pulse-response experiments
Before the pulse-response experiments, all probes were calibrated for exact length and impedance by WinTDR procedure.Information on the six probes after calibration is shown in Table 1.Firstly, a steady-state water inflow was maintained by a peristaltic pump.The steady-state flux condition was also checked by observing the constant water content by TDR.The coefficient of variation (C V ) for the water content was between 0.08% and 1%.Ordinary tapwater was used as "tracer-free" but had a small background EC.On top of Sect. 1 a paper filter was placed in order to spread water uniformly over the sand.The salt tracer pulse was applied by changing the water source to the pump from tapwater to the potassium chloride (KCl) solution, which was equivalent to a surface application rate of 1.5×10 −3 g cm −2 .The pulse duration was 30 minutes while maintaining the same pumping speed before, during and after the pulse application to ensure a constant pore water velocity.Four different fluxes (i.e.0.5 cm hr −1 , 1 cm hr −1 , 1.5 cm hr −1 and 2 cm hr −1 ) were applied.The transport of the solute per segment, 5 in total, was characterized by monitoring the change in EC at the inlet and outlet of each segment.

Data analysis
The background EC was firstly subtracted from the measured EC responses in order to obtain the increase in EC due to the tracer.The start and the end of the response was determined by a simple and automatic algorithm which had also to avoid a fake start.The slopes of adjacent EC data over the time step were calculated.The start of rising slope was identified by exceeding a specified minimum value.To avoid a false start the slope was taken sufficiently high and the start was set 3 time intervals before exceeding the minimum slope.The EClevel just before the start was also set as the background EC.In this way the algorithm was robust and avoided a false start caused by fluctuation in background EC.The end of a response was set either when the background EC was reached or after a specified maximum duration.Whichever came first was taken as the end, and in most cases the end of response was determined because the background EC had been reached.The EC values after subtraction of the background were summed over the duration.By dividing the EC values by this sum, normalized relative EC values were obtained with sum equal to 1.The normalization also avoided problems due to variation in water content.Imposing the same sum to all response ensures a conservation of the tracer.
The Eq. ( 9) as described by Mojid et al. (2004) was fitted to the pulse-response normalized EC-data.For every segment the signal at the upstream inlet was taken as the input (C in ) while the signal at the downstream outlet was taken as the response (C r ).Pore water velocity and dispersion coefficient were determined for each segment.The method was implemented in the R software, which is a flexible open source and analysis free software under General Public License (GPL).Although R is meant as a language and environment for statistical computing and graphics, it is flexible and allows writing tailor-fitted analysis so that it can be extended through desired packages for specific purposes (Dalgaard, 2004).The automatic calibration was executed by parameter search algorithms available in R. (R-code analysis according to the Mojid et al. ( 2004) is available at simple request).

Testing by HYDRUS 2-D/3-D
HYDRUS 2-D/3-D provides a numerical multi-dimensional solution to the transport equations under variable saturated conditions.It solves the Richard's equation for water flow and solves the CDE for solute and heat transport.The simulation can be displayed graphically and animation can be shown.The HYDRUS program uses Marquardt-Levenberg optimization algorithm for the inverse estimation of soil hydraulic and solute transport parameters.The simulation results for nodes of the mesh can be stored for comparison to the measurements.
A 2-D vertical plane was created to represent the geometry of the J-shaped model in Fig. 1.The initial and boundary conditions were set according to the experimental set-up and pulse-response experiment procedures.Each TDR probe was compared to a node in the finite element mesh corresponding to the middle of the probe.For probe No. 5 and 6 in the bend we also checked extra nodes between the start and the end of the probe.
The simulation started with the establishment of steady state water flux condition with a low background concentra-  tion.Then during the pulse-time of 30 min the tracer concentration was supplied, after which the water inflow with the same low background concentration as before the pulse was used.The simulated concentrations were normalized after subtraction of the background level in the same way as the EC-measurements.
Table 2. Solute-transport parameters for coarse dune sand with four water flux conditions and with concentration in 1.5 g L −1 of the applied pulse of potassium chloride (KCl) for 30 min input.Pore water velocity (V) and dispersion coefficient (D) are determined by the transfer-function method fitted on individual segments (Mojid et al., 2004).
Inverse modelling was performed by using the measured relative EC data from the pulse-response experiment to estimate the soil longitudinal dispersivity.Subsequently the result by inverse modelling was compared with the simulation by using the transport parameters relation identified by fitting the longitudinal dispersivity on segments.Pore water velocity was calculated by HYDRUS afterwhich the longitudinal dispersivity at each sub-region was specified according to Eq. ( 2) with dispersion coefficient determined from the obtained equation.Equation (2) was used while neglecting the molecular diffusion and setting the empirical coefficient n=1.The value of the transversal dispersivity in this simulation was set as 1/10 of the longitudinal dispersivity.Setting transversal dispersivity at 1/20 of the longitudinal dispersivity and zero transversal dispersivity were also tested for checking the parameter sensitivity.

Measured and estimated breakthrough data
Table 2 summarizes the transport parameters obtained from analysing the laboratory experiments by the transfer-function on the segments between probes.
The majority of the measured and estimated response curves (C r and C r.est ) at each segment under different flux conditions indicated a good fit.Coefficients of determination R 2 are very high for almost all segment calculations (majority shows from 0.97 to 0.99).Except to the experiment with water flux of 1 cm hr −1 from probe No. 3 to 4 the response had a shorter duration than the input, which is inconsistent with the dispersion process.Figure 2 illustrates an example of a good and one of a less good fit.Flux condition 0.5 cm hr -1 1.0 cm hr -1 1.5 cm hr -1 2.0 cm hr -1

Solute transport parameters by fitting the transfer function on segments
As shown in Table 2, the average water content between adjacent probes for each segment had different soil water content.In addition water content varied as expected with the flux imposed.The elevation of the piezometric surface as measured by the piezometer increased in function of the flux and was located between the probes of No. 3 and 4. The section delineated between the probes No. 1 and 3 was always unsaturated.All segments downstream of probe No. 4 were always saturated.As shown in Fig. 3 at the unsaturated vertical section (probe No. 1 to 3) the average water contents were less than 30%, and the water contents in the saturated section (probe No. 3 to 6) varied between 37 and 44%, in accordance to the porosity range of the coarse sand (35% to 40%).
As expected the measured pore velocities were much higher in the unsaturated section, as the water content is lower and hence the water filled pore space is less.As shown in Fig. 4 the pore-water velocity decreased dramatically in the saturated section, and the measured dispersion coefficients varied accordingly.
As defined by Eq. ( 2), dispersivity was calculated as the dispersion coefficient divided by the pore water velocity, hereby neglecting the molecular diffusion effect and taking empirical coefficient n=1.Most calculated dispersivities in our experiments were smaller than 1 cm, except for the extreme high value measured at the beginning of the experiment.The values of dispersivity are shown as a function of the volumetric water content in Fig. 5 Flux condition 0.5 cm hr -1 1.0 cm hr -1 1.5 cm hr -1 2.0 cm hr -1 1.0 cm hr -1 1.5 cm hr -1 2.0 cm hr -1 No. 4 to 5 and No. 5 to 6 the dispersivity did not respond to a change in pore water velocity.The latter situation was also reported by Toride et al. (2003); their study observed the occurrence of considerable tailing of breakthrough curves for unsaturated flow.Similarly in our experiment larger dispersivities were observed in the unsaturated soil section as compared to the saturated soil section.For lower fluxes condition the tailing effect of breakthrough curve was not distinct.For higher fluxes the tailing was much more pronounced, especially in the unsaturated section (i.e. at probe No. 1 and 2).
The measurements for the unsaturated section of probe No. 2 to 3 did not perform similarly to the probe No. 1 to 2. This is likely to be caused by a different packing of the sand in that segment.In the segment of probe No. 3 to 4, partly unsaturated and saturated, it was noticed as in segment of probe No. 1 to 2 that dispersivity increased as pore water velocity decreased.
In order to simulate over the entire range of fluxes and water contents as present in the continuum of unsaturated zone, aquifer and river, a more general relationship was identified.After exploring several possible relations the best result was obtained by plotting on a log-log scale the solute dispersion coefficient against the ratio of pore water velocity over water content as shown in Fig. 7.
This relation was similar to the empirical power law shown by Padilla et al. (1999).They compared their data to earlier findings (e.g.those by De Smedt and Wierenga, 1978) and also proposed a log-log relation between the dispersion coefficient and the of pore water velocity over water content.The following equation fitted the data in our study: V/θ v (cm hr -1 )

-1
Figure 7. Log-log relationships between solute dispersion coefficient (D) and ratio of pore water velocity to water content (V/θv); two outliers (*) are not included into the relation Fig. 7. Log-log relationships between solute dispersion coefficient (D) and ratio of pore water velocity to water content (V /θ v ); two outliers (*) are not included into the relation.
The coefficients in Eq. ( 10) are soil specific.During the fitting we excluded 2 outliers.The two outlier data are in the unsaturated section of segment between probe No. 1 and 2 and segment between probe No. 3 and 4 of flux=2 cm hr −1 .The measured dispersivities are relatively low when compared with the data of the same segment in other fluxes.The study of Padilla et al. (1999) used silica sand as material and they found 1.99 as the power coefficient for their data, which is very close to our value of 2.02.A relation between the dispersion parameter, water flux and the soil water content is required for simulating the transport of solutes and other substances in the environment which is imbedded by saturated and unsaturated zones.

Inverse modelling by HYDRUS 2-D/3-D
To the numerical simulation by using HYDRUS 2-D/3-D, initial results showed an anomaly due to different packing of sand in the segment between probe No. 2 and 3, in which the transport was much slower than in the adjacent segments.This anomaly can be improved by changing the parameter α in the soil water retention function for segment between probe No. 2 and 3.The shape of the soil water retention curve according to van Genuchten (1980) and the hydraulic conductivity was modified in that segment, which made it possible to simulate this phenomenon.This result indicated that the shape of the soil water retention curve has a very important effect on the solute transport.The use of "default characteristics" and pedotransfer functions in HYDRUS should therefore be used with caution when comparing to real soils or porous media.
The inverse modelling was therefore performed by optimizing two parameters: λ 1 , the longitudinal dispersivity for Hydrol.Earth Syst.Sci., 13,[217][218][219][220][221][222][223][224][225][226][227][228]2009 www.hydrol-earth-syst-sci.net/13/217/2009/ the region before and after probe No. 2 and 3 with α=0.145, and λ 2 , the longitudinal dispersivity for the segment between probe No. 2 to 3, in which α was set as 0.03.Results of parameters optimization in different fluxes are summarized in Table 3.Initial runs by HYDRUS 2-D/3-D suffered instabilities leading to negative concentrations and oscillating tails of the responses.Also the simulation of the last probes showed initially a slow increase and a fast recession followed by numerical oscillation.The discretization was therefore reduced to 0.5 cm for the vertical section and 0.3 cm for the bend section as finite element size.However, this lead to excessive computational times.For inverse modelling, which requires iteration, this resulted in several days on a PC for one flux only.
The decrease of longitudinal dispersivity for the region before and after probe No. 2 and 3 (λ 1 ) was found as flux increased from 0.5 cm hr −1 to 1 cm hr −1 , subsequently an increase was found as flux increased from 1 cm hr −1 to 2 cm hr −1 .The variation of dispersivity (λ 2 ) at the segment between probe No. 2 to 3 followed a different pattern.Most studies assume that the dispersivity is a constant or an intrinsic property of soils.The inverse modelling showed that longitudinal dispersivity is a variable in function of water flux.This is especially important for porous media under variably saturated conditions.
HYDRUS needed dispersivities as input parameters.This was implemented by using the power relationship of the dispersion coefficient in function of the pore water velocity and water content (Eq. 10 in this study).Different longi-  tudinal dispersivities were therefore specified for each segment between adjacent probes using the pore water velocity and water content as calculated by HYDRUS.The timing of the resulting simulated responses corresponded well with the measured ones.In the bend different velocities were simulated: slower at the inner side and faster at the outer side.One hypothesis was that the atypical behaviour was due to transversal dispersivity between layers with different convective velocities.A sensitivity analysis using different values for relative transversal dispersivity, including zero  transversal dispersivity did not show any visible difference between the simulations.

Comparison of inverse modelling and fitting by segment
The comparison of the measured relative EC with HYDRUS 2-D/3-D simulation by inverse modelling and by using the equation by fitting the transfer function on individual segment is given in Fig. 8a, b, c and d.Three over four simulations by using the equation for dispersion coefficient based on segments presented better correlation than the simulation by fitting dispersivities.For the lowest flux the inverse modelling delivered better results.The fact that inverse modelling delivered a range of dispersivities in function of the flux illustrates that the classical Eq. ( 2) should be used with caution in circumstances with variably saturated conditions.The fitting of dispersion coefficients on individual segments was a lot more straightforward and allows a generalization after exploratory study of the possible relations.A graphical comparison of the responses in the saturated part showed that the inverse modelling overestimated the dispersion, while our Eq.( 10) was leading to a better fit.

Conclusions
The characterization of the lateral inflow from the aquifer to the river in this study contributes to a better insight of riverbed-aquifer flux interactions.In order to determine the hydrodynamic dispersion transport parameters in this variably saturated environment, a laboratory J-shaped column model was designed.By analysing consecutive column segments with the application of the transfer function method proposed by Mojid et al. (2004), the relationship between dispersivities and other physical measurements were determined and discussed.The dispersion parameters were generalized as following: 1.The dispersivity was flux dependent and increased at lower flux in unsaturated section (i.e. in vadose zone); in contrast, dispersivity was fairly constant and independent of flux variation in saturated section (i.e. in riverbed in conjunction with the adjacent aquifer).
2. The longitudinal dispersion coefficient can be best related to the ratio of pore-water velocity over soil water content.This relation can be performed over the range of saturated and unsaturated conditions, and it appears similar to the earlier findings.Further testing of this relation with different soil materials is recommended.
The testing by using HYDRUS 2-D/3-D supports the use of Eq. (10) to determine dispersion coefficients.The result of simulation shows that:

Figure 1 .
Figure 1.Experimental set-up of J-shape model

Figure 2
Figure 2(a).Example of a good fit on the observed pulse-response of segment 4 from probe No. 4 to 5 in flux of 0.5 cm hr -1 Fig. 2a.Example of a good fit on the observed pulse-response of segment 4 from probe No. 4 to 5 in flux of 0.5 cm hr −1 .

Figure 2
Figure 2(b).Example of a less good fit on the observed pulse-response of segment 3 from probe No. 3 to 4 in flux of 1.0 cm hr -1 Fig. 2b.Example of a less good fit on the observed pulse-response of segment 3 from probe No. 3 to 4 in flux of 1.0 cm hr −1 .

Figure 3 .Fig. 3 .
Figure 3. Distribution of average water content (θv) along the model for different fluxes Fig. 3. Distribution of average water content (θ v ) along the model for different fluxes.

Figure 4 .Fig. 4 .
Figure 4. Distribution of pore-water velocity (V) along the laboratory model for different fluxes Fig. 4. Distribution of pore-water velocity (V ) along the laboratory model for different fluxes.

Figure 8 Fig. 8a .
Figure 8(a).Comparison of the measured relative EC with HYDRUS 2D/3D simulated by inverse modelling fitting dispersivity (r 2 = 0.86) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.72) in flux of 0.5 cm hr -1 Fig. 8a.Comparison of the measured relative EC with HYDRUS 2-D/3-D simulated by inverse modelling fitting dispersivity (r 2 = 0.86) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.72) in flux of 0.5 cm hr −1 .

Figure 8 Fig. 8b .
Figure 8(b).Comparison of the measured relative EC with HYDRUS 2D/3D simulated by inverse modelling fitting dispersivity (r 2 = 0.68) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.72) in flux of 1.0 cm hr -1 Fig. 8b.Comparison of the measured relative EC with HYDRUS 2-D/3-D simulated by inverse modelling fitting dispersivity (r 2 = 0.68) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.72) in flux of 1.0 cm hr −1 .

Figure 8 Fig. 8c .
Figure 8(c).Comparison of the measured relative EC with HYDRUS 2D/3D simulated by inverse modelling fitting dispersivity (r 2 = 0.66) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.72) in flux of 1.5 cm hr -1 Fig. 8c.Comparison of the measured relative EC with HYDRUS 2-D/3-D simulated by inverse modelling fitting dispersivity (r 2 = 0.66) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.72) in flux of 1.5 cm hr −1 .

Figure 8 Fig. 8d .
Figure 8(d).Comparison of the measured relative EC with HYDRUS 2D/3D simulated by inverse modelling fitting dispersivity (r 2 = 0.44) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.50) in flux of 2.0 cm hr -1 Fig. 8d.Comparison of the measured relative EC with HYDRUS 2-D/3-D simulated by inverse modelling fitting dispersivity (r 2 = 0.44) and simulated by using equation for dispersion coefficient based on segments (r 2 = 0.50) in flux of 2.0 cm hr −1 .

Table 1 .
Calibrated probe length L p and calibrated characteristic probe impedance Z 0 .