"too many steps in qromb" when using PFAS module

[siminhydrus]

Dear support team,
I’ve developed a 3D PFAS transport model using HYDRUS 5.04.0230. The model represents about 30 meters of unsaturated zone and 6 meters of saturated zone at the bottom. I’ve defined five different soil materials based on observed lithology at the site.
I was able to successfully run the model for flow using atmospheric BC at the top and variable pressure head at the bottom. Then, I added solute transport and sorption to the solid phase using linear sorption with Kd values specific to each material. For solute transport, I used measured sorbed concentrations in the soil as the initial condition by multiplying the measured values by the average dry bulk density. I set the initial conditions as "in total concentrations." These initial sorbed values act as the source of PFAS, and since there’s no more application of PFAS at the site, I applied zero concentration flux at the top and bottom boundaries. My goal is to calculate the solute mass flux at the bottom boundary over time.
The solute transport model with linear sorption ran without any issues. However, when I activated the PFAS module and defined linear Kh_AWI values for each material (I have one solute and considered equilibrium sorption to the air-water interface), the model refuses to run. Two observations are:
1) Right after reading the solute transport information, I can see the message "too many steps in qromb," and the model doesn’t start running (I waited for about 5 minutes).
2) If I press Enter, the same error line appears twice, followed by "error when writing to an output file." (screenshot attached)
This seems to be related to calculating the air-water interface area for each material, but I’m not sure what’s causing it or how to fix it. I haven’t changed anything else in the model other than defining Kh values after activating the PFAS module.
Are there any additional parameters or options I need to adjust for the PFAS module to run? Any thoughts on what might be causing this issue? Your help is greatly appreciated.
Best,
Simin

[Jirka]

The subroutine qromb (from Numerical Recepies) is used in HYDRUS for the integration of a function over a closed interval using Romberg Integration (and it pauses when it does not converge). It is used only in two places:
1. to calculate the air-water interfacial area
2. to calculate soil matric flux potential
and thus, clearly, the first one applies to you. The air-water interfacial area is calculated by integrating the retention curse (see the manual) from the soil hydraulic parameters (alpha, n, Qr, and Qs). I have tested the routine for all 12 textural classes of Carsel and Parrish (1988), and I have never encountered any problems. I would, thus, guess that you must be using some strange (unrealistic) parameters.

For some clay soils, and for a very dry conditions (not likely to occur), huge numbers of the air-water interfacial area may be obtained. The algorithm simply pauses when it does not converge, and once you hit any keyboard key, it continues. As a result, this should not affect the rest of the simulation. On the other hand, such huge AWI values will undoubtedly cause trouble for transient conditions (with destruction or creation of AWI).

I looked more closely on the code and adjusted the parameters in the Qromb routine. I have recompiled the computational module and copied it to DropBox. You can download it at this link: https://www.dropbox.com/scl/fo/nmeekfsh ... 0awk8&dl=0. You can test it by replacing the file in this folder of the HYDRUS installation (HYDRUS 5.x 64-bit\Bin), or you can run it externally of GUI (see FAQ 1 at PC-PROGRESS - FAQ). This exe will be in the upcoming version 5.05 (also for other dimensions), and will be deleted once this version is released.

J.

[siminhydrus]

Dr. Simunek, thanks for the updated .exe file. It ran successfully.
I have another question regarding the same model (a different topic, but since I've already explained the project here, I decided to post under the same thread). I am using a constant head boundary condition at the bottom of the domain. After running the model for four years, the constant boundary flux mostly shows negative values (meaning water is added to the system to maintain the constant head at the bottom, I assume), except after two major precipitation events.

My questions are:
1) Is it still possible to estimate groundwater recharge with the current model setup (i.e., a constant head boundary at the bottom)? Or are recharge estimates only possible when using free drainage at the bottom?
2) Checking the constant boundary solute flux, there is always a positive solute flux, meaning mass is being removed from the bottom boundary. I am a bit confused about why solute mass is leaving the system while the water flux at the bottom is negative (indicating water is being added to the system). Let me know if I am interpreting this incorrectly. I should also mention that I have sorbed PFAS in all model layers, including the saturated portion of the soil—I’m not sure if that is relevant here.

I am attaching screenshots of both water and solute fluxes for your reference. Thanks a lot in advance!

[Jirka]

Hydrus calculates bottom water fluxes (usually referred to as recharge) for any specified boundary conditions. For a free drainage BC, the bottom flux is always downward (i.e., negative). For a pressure head BC, it can be both downward (i.e., negative) or upward (i.e., positive), depending on the conditions in the profile (i.e., whether there is upward or downward pressure head gradient).

Since I do not know what your bottom solute BC is, I cannot comment on bottom solute fluxes.

J.

[siminhydrus]

Thanks, Jirka, for the clarification. I now realize I was completely off about the positive and negative fluxes. A negative bottom flux indicates downward flow, meaning water is leaving the system. Based on your explanation and the constant boundary flux plot, can we say that the average recharge rate is approximately 0.003 m³/day?

Apologies for not mentioning the solute boundary condition earlier. I am using a third-type (flux) boundary condition with zero concentration at the bottom of the domain. As you can see, there is a positive solute flux from the bottom boundary. Can we interpret this as the mass of PFAS leaving the bottom boundary (or, more precisely, reaching the groundwater) due to desorption occurring in the vadose zone from the initially sorbed PFAS concentration?

Thank you,
Simin

[Jirka]

Average recharge: Visually it seems so. However, it would be better to take cumulative bottom flux and divide it by time to get the average recharge rate.

Solute BC. Solute fluxes as positive for inflow and negative for outflow. I do not know why you don’t find this information in the manual or help, rather than wasting my time! I do not think what you write is correct. If GW concentration is zero, then the upward bottom solute flux (positive) should be zero as well, while downward bottom solute flux should reflect solute concentrations in the domain.

J.

[siminhydrus]

I did check the manual, and I just checked it again. This is a 3D model. According to the manual, 'Boundary fluxes are positive when water is removed from the system' and 'Values of the solute flux and the cumulative solute flux are positive when the solute is removed from the system.' This is exactly how I explained it when asking my question. thanks for the explanation anyway.

Simin