Water Infiltration in 1D soil column

[wurunjian]

Hi all,

I did a simple water infiltration in 1D soil column simulation. The
input file is revised based on "demo_simple_flow" example. I also did
the similar simulation in Hydrus-1D, but I found the two results of soil
water content are different. I listed the main parameters below. Can
someone help me explain the difference or point out the mistakes I made
in unit conversion? Thanks in advance!

=== Soil Column Property ===
Ks = 0.0468 cm/min (PERM_Z = 7.96d-13)
thetas_s = 0.39 (POROSITY = 0.39)
theta_r = 0.035 (LIQUID_RESIDUAL_SATURATION = 0.0897d0)
alpha = 0.02 cm-1 (ALPHA = 2.04d-4)
n = 2.26 (M = 0.558d0)
tortuosity = 0.5 (TORTUOSITY = 0.5 )

=== Boundary Condition ===
Top BC: Flux, 12mm/h (Flux 3.333d-6 m/s)
Bottom BC: Free Drainage

=== Initial Condition ===
The whole column, theta_i = 0.1 (PRESSURE = 1.89d4)


Runjian

[Gautam Bisht]

Runjian,

There where following two errors in your PFLOTRAN inputdeck:

- Initial condition pressure was such that theta_i /= 0.1, and

- A dirichlet bottom boundary condition was used, which was set to initial condition. While the correct bottom boundary condition should have been a 'unit_gradient' 

Attached below is the updated PFLOTRAN inputdeck and time series plot. A qualitative comparison of Hydrus-1D and PFLOTRAN results now look reasonable to me.

-Gautam.

[Peter Lichtner]

Runjian: if possible I would like to use your problem as a benchmark for the PFLOTRAN User_Manual if you agree. I would need the data from the Hydrus run.

…Peter

<pflotran_theta_time_series.pdf><pflotran_gbisht.in>

[wurunjian]

Hi Gautam,

Thanks for your help! I have two questions about the changes you made.
(1) I used the equation (A-6) in the user manual to convert initial water content to pressure. But I got the different value (1.89d4) with you. How can you got the pressure (8.4d4)?

(2) The FLOW_CONDITION (unit_gradient) is meaning free drainage? Actually I cannot find the "unit_gradient" type in the user manual.

============

Hi Peter,

You are free to use this example. I attached the result file. Hydrus-1D uses FE method and the water content is calculated on the node (In this example there is total 101 nodes). If you need more information, please ask me.

Thanks,

Runjian

[Gautam Bisht]

Hi Runjian,

(1) I used the equation (A-6) in the user manual to convert initial water content to pressure. But I got the different value (1.89d4) with you. How can you got the pressure (8.4d4)?

An equation is missing in the manual to relate capillary pressure (in section A.2) with liquid pressure (in section A.1). 

P_liquid = P_reference - P_cap; 

where P_reference = 101,325 [Pa]

The value of 1.89d4 you obtained from equation A-6 is capillary pressure, while the primary unknown variable in RICHARDS mode is liquid pressure. Thus, initial pressure should be

P_liquid = 101325 - 1.89d4 = 82425 [Pa] ~ 8.4d4

 

(2) The FLOW_CONDITION (unit_gradient) is meaning free drainage? Actually I cannot find the "unit_gradient" type in the user manual.

Yes, unit_gradient is meant for free drainage and it is one of the unadvertised capabilities in PFLOTRAN. The reason a particular capability is unadvertised is because it may not be thoroughly tested or implemented for only to a limited number of PFLOTRAN modes.

-Gautam

[Peter Lichtner]

I will update the user_manual.

…Peter

[wurunjian]

Hi Gautam,

What's the difference between "unit_gradient" and "zero_gradient" in pressure condition?
How PF deals with "seepage face" condition? Just sets the pressure head = 0 if saturation = 1?

Thanks,

Runjian

[Hammond, Glenn E]

A zero gradient flow condition would be "no flow".  I believe that unit gradient allow drainage if the pressure gradient exceeds hydrostatic.  Someone else can correct me.

The seepage face prevents inflow if the boundary pressure Is below the reference pressure (usually atmospheric).  Outflow is always possible as long as the gradient is in that direction.

Glenn

[Gautam Bisht]

Hi Runjian,

Darcy flux (eqn A-2 in the user manual) is given as

q = - (k * kr)/mu * dphi/||dist||

dphi = (P - W*rho*g*z)_up - (P - W*rho*g*z)_dn,

||dist|| = norm of separation distance.

where 'up' and 'dn' corresponds to upwind and downwind control volumes.

For the unit_gradient BC, the difference in liquid pressure between up and dn cells is not accounted for, so

dphi = - (W*rho*g*z)_up + (W*rho*g*z)_dn (see richards_common.F90)

A zero-gradient BC on a region for RICHARDS mode = Not specifying any BC on a region for RICHARDS mode.

As I understand it, the advantage of ZERO_GRADIENT_BC is for mode other than RICHARDS. e.g. In TH mode, you could specify a ZERO_GRADIENT_BC for pressure, but a dirichlet/neumann for temperature on a region. But, if you one doesn't specify a BC for a region, "no flow" BC is assumed by default for both pressure and temperature.

-Gautam.