Hi Board,
I am modeling drainage out of a seepage face at the bottom of a hillslope that has been irrigated to just above field capacity. There is a water table at the bottom of the hillslope and a wedge of water that creeps up the hillslope along the bedrock.
I am specifying initial conditions. At first, I set the entire hillslope to a given water content (just below theta.s), but that didn't yield enough water out the seepage face. Then I used pressure head to specify initial conditions, and defined the elevation of the water table (pressure head = 0) using empirical data, with increasing pressure toward the ground. This yielded much more water, which is closer to empirical results, but still not enough. Now I would like to add more water (ie, the saturated wedge). However, I don't have spatially fine resolution data to know the exact dimensions/shape of the wedge, and the wedge is an irregular shape (unlike the water table, which is just a horizontal line). Does anybody know how to add the saturated wedge to the initial conditions?
Another option is to create an initial conditions file by irrigating my model at a high rate until saturation, which would give a more realistic shape of the water table with saturated wedge. However, this would introduce some artifacts (such as not accounting for macropore flow in my model), I think?
Thanks,
Ray
how to model saturated wedge as initial condition?

 Posts: 2
 Joined: Wed Apr 20, 2016 1:22 pm
 Location: USA
Re: how to model saturated wedge as initial condition?
Hi Ray, all:
I am attempting to use a similar set up to you. I am modeling a 2D hillslope where the downslope edge of the domain is a combination of seepage face in the upper meter and gradient boundary below (down to bedrock at 10 m depth). I have observation well data that tells me the location of the groundwater table at the beginning of my simulation at the edge of my simulation domain. I am trying to initialize a saturated wedge by setting the measured groundwater level to h=0 cm in the simulation with linearly increasing pressure head downwards from there, as you tried.
However, I am running into numerical stability issues. When I initialize any part of the domain with h>=0 cm, the model will crash within its first time step. If I set the bottom of the domain to be h<0 cm (even 1 cm), it will run, but my results are not accurate as I know there is groundwater within my simulation. I have tried many things to get Hydrus to run with a portion of the domain initially saturated but no luck so far (increased pressure head tolerance to 2 cm; decreased minimum time step to 10^12 days; increased spatial discretization).
Any thoughts on why Hydrus crashes immediately when initialized with a saturated wedge?
Mark
I am attempting to use a similar set up to you. I am modeling a 2D hillslope where the downslope edge of the domain is a combination of seepage face in the upper meter and gradient boundary below (down to bedrock at 10 m depth). I have observation well data that tells me the location of the groundwater table at the beginning of my simulation at the edge of my simulation domain. I am trying to initialize a saturated wedge by setting the measured groundwater level to h=0 cm in the simulation with linearly increasing pressure head downwards from there, as you tried.
However, I am running into numerical stability issues. When I initialize any part of the domain with h>=0 cm, the model will crash within its first time step. If I set the bottom of the domain to be h<0 cm (even 1 cm), it will run, but my results are not accurate as I know there is groundwater within my simulation. I have tried many things to get Hydrus to run with a portion of the domain initially saturated but no luck so far (increased pressure head tolerance to 2 cm; decreased minimum time step to 10^12 days; increased spatial discretization).
Any thoughts on why Hydrus crashes immediately when initialized with a saturated wedge?
Mark
Re: how to model saturated wedge as initial condition?
The initial condition in the entire domain must be physically realistic, otherwise the program may not be able to resolve the first time step. Not only below the water table but also above. The code clearly can handle saturated bottom conditions, see a couple of tutorials at: , https://www.pcprogress.com/en/Default. ... tutorials , e.g., " 3.03  Furrow Infiltration with a Solute Pulse" or " 3.04  Flow and Transport in a Transect to a Stream" for 2D or "5.02  A Simple Reservoir" for 3D.
You may also try, for example, to get some initial pressure head profile with simplified boundary conditions (e.g., no fluxes at boundaries) and then import this pressure head profile as initial conditions into the project with more complex BCs.
J.
You may also try, for example, to get some initial pressure head profile with simplified boundary conditions (e.g., no fluxes at boundaries) and then import this pressure head profile as initial conditions into the project with more complex BCs.
J.
Re: how to model saturated wedge as initial condition?
Thanks Jirka,
I have been using the hydrostatic equilibrium from lowest point initial condition and am still running into the same numerical issues when having h>=0 cm at the bottom and the vertical boundary on the downslope edge of the domain set as gradient or seepage face. When using no flux or constant flux at the boundary the simulation runs ok, but these are obviously not appropriate for my project. Also when setting the lowest nodal point h<0 it runs fine.
I have also noticed that it runs when setting the gradient boundary value to a positive number. However, water is instead being forced into the domain since my domain slopes down from left to right. The calculated slope angle and gradient value should thus be negative.
The final thing I have tried is to set up the model with no flux conditions as you suggested and use those final conditions as the initial conditions for my problem. The final conditions reached in the no flux simulation eventually reach hydrostatic equilibrium and I am having the same issues described above.
Mark
I have been using the hydrostatic equilibrium from lowest point initial condition and am still running into the same numerical issues when having h>=0 cm at the bottom and the vertical boundary on the downslope edge of the domain set as gradient or seepage face. When using no flux or constant flux at the boundary the simulation runs ok, but these are obviously not appropriate for my project. Also when setting the lowest nodal point h<0 it runs fine.
I have also noticed that it runs when setting the gradient boundary value to a positive number. However, water is instead being forced into the domain since my domain slopes down from left to right. The calculated slope angle and gradient value should thus be negative.
The final thing I have tried is to set up the model with no flux conditions as you suggested and use those final conditions as the initial conditions for my problem. The final conditions reached in the no flux simulation eventually reach hydrostatic equilibrium and I am having the same issues described above.
Mark