Dear Hydrus team and Forum Members,
We simulated a simple conservative tracer experiment in a homogeneous soil column (L=100cm), with constant hydraulic head boundary conditions inducing a saturated column with flow from top to bottom.
We are simulating a constant injection of a conservative tracer on the top. We are interested in the breakthrough curves along the soil profile. We would like to compare the Hydrus output to the OgataBanks solution, expecting they should be identical for the entire column length.
We played with the solute boundary conditions and have two cases:
 Case 1, where we apply a constant concentration boundary condition in Hydrus1D
 Case 2, where we apply a constant flux boundary condition
The breakthrough curves are indeed identical but only at the distance from 0cm to approximately 92cm from the source of injection (case 1). Starting from 92cm till the column outlet, or the final 8cm or so of the profile, the Hydrus solutions and OgataBanks solutions are increasingly deviating (case 1). We have tried several different sytem settings and conditions, and still got the same kind of differences between the numerical and analytical solutions.
But when the upper BC for the solute transport is concentration flux BC (case 2), the solutions are deviating along the entire column length but gradually converging to the OgataBanks solutions at the column outlet. So the Hydrus breakthrough curve at the outlet looks more similar to the analytical solutions than in the constant concentration BC case. We all know that it’s the constant concentration BC in OgataBanks.
We suspect what might have caused these differences in Hydrus solution is in the way how Hydrus implements the boundary conditions which we are not sure of. We would appreciate any comments of clarification and suggestions for minimizing the difference between the analytical and the Hydrus solution.
A summary of our model is below, and I also attached the Hydrus codes in the zip file.
L = 100 cm
i = 0.01
K_sat = 8.64 cm/day
Porosity = 0.25
Dispersivity = 4.423 cm
Diffusion coefficient = 0
Simulation time = 600 days
Water flow upper and lower BCs: Constant pressure head
Solute transport upper BC: Constant concentration of 1 mmol/cm3 (case 1); Concentration flux of 1 mmol/cm3 (case 2)
Solute transport lower BC: Zero concentration gradient
Thank you.
Hue
Boundary condition implementation in Hydrus1D
Boundary condition implementation in Hydrus1D
 Attachments

 cases12.zip
 (696.87 KiB) Downloaded 16 times
Re: Boundary condition implementation in Hydrus1D
The solute transport part of HYDRUS1D has been extensively and successfully tested against the analytical solutions from CFitm, CFitIm, and CXTFIT. Note that these analytical solutions are valid for different surface (Dirichlet or Cauchy) BCs, as well as for different bottom boundary conditions (finite and infinite columns). Obviously, the solution must to be different for Dirichlet and Cauchy BCs (as more solute enters the column for the Dirichlet BC), and for different bottom BCs. You show only one analytical solution, which you call the OgataBanks solution (which I’m not familiar with), without actually specifying for which BCs it was derived. You clearly cannot compare two different numerical solutions for different surface BCs with a single analytical solution. J.
Re: Boundary condition implementation in Hydrus1D
Thank you, Jirka, for the quick comment. We will be looking at the 3 analytical solutions you mentioned.
With our cases, I should have been clearer.
As an analytical solution to the advectiondispersion equation (ADE), the OgataBanks solution (1961) we used has the following form:
C(x,t) = (Co/2) * [erfc[(xvt)/2sqrt(Dt)] + exp(vx/D) * erfc[(x+vt)/2sqrt(Dt)]]
Where erfc is the complementary error function, and:
C = solute concentration [M/L^3]
x = distance of interest from the source of solute injection, along the flow path [L]
t = time since the release of the solute [T]
Co = initial solute concentration [M/L^3]
v = average linear groundwater velocity (= q / porosity) [L/T]
D = longitudinal dispersion coefficient [L^2/T]
The OgataBanks solution was derived from the following initial and boundary conditions:
Initial condition: C(x,0) = C_i
Upper BC: Dirichlettype (or concentrationtype BC); C(0,t) = Co
Lower BC: Semiinfinite
So our Case 1 is comparable to the OgataBanks solution. Comparing OgataBanks to case 2 was just our curiosity.
Do you have additional comments? We would greatly appreciate it.
Ref: Ogata, A., and Banks, R. B. 1961. A solution of the differential equation of longitudinal dispersion in porous media. U.S. Geological Survey Professional Paper 411A, A1A9.
With our cases, I should have been clearer.
As an analytical solution to the advectiondispersion equation (ADE), the OgataBanks solution (1961) we used has the following form:
C(x,t) = (Co/2) * [erfc[(xvt)/2sqrt(Dt)] + exp(vx/D) * erfc[(x+vt)/2sqrt(Dt)]]
Where erfc is the complementary error function, and:
C = solute concentration [M/L^3]
x = distance of interest from the source of solute injection, along the flow path [L]
t = time since the release of the solute [T]
Co = initial solute concentration [M/L^3]
v = average linear groundwater velocity (= q / porosity) [L/T]
D = longitudinal dispersion coefficient [L^2/T]
The OgataBanks solution was derived from the following initial and boundary conditions:
Initial condition: C(x,0) = C_i
Upper BC: Dirichlettype (or concentrationtype BC); C(0,t) = Co
Lower BC: Semiinfinite
So our Case 1 is comparable to the OgataBanks solution. Comparing OgataBanks to case 2 was just our curiosity.
Do you have additional comments? We would greatly appreciate it.
Ref: Ogata, A., and Banks, R. B. 1961. A solution of the differential equation of longitudinal dispersion in porous media. U.S. Geological Survey Professional Paper 411A, A1A9.
Re: Boundary condition implementation in Hydrus1D
I think that the answer is obvious. The numerical model does not consider semiinfinite profile (as the analytical solution does); the soil profile is always finite and the boundary conditions are specified at a given locations and this BC affect the concentration profile in its near vicinity. If you want to compare the two solutions, you need to apply the BC (e.g., a zero gradient) at a location far from locations for which you compare the results. J.
Re: Boundary condition implementation in Hydrus1D
Thanks Jirka. That really solved the puzzle.