In my quest to understand the sources of error from inversion and ways to consider mitigating them I stumbled on a question for which I couldn’t find an answer. Hydrus uses a finite element method (FEM) to simulate water movement in the sample. This is based on Richards’ equation and is used to predict tension from van Genuchten parameters. After a simulation runs, Hydrus compares the predicted tensions to the measured ones that are given at the start. It evaluates the error in the estimation, changes the van Genuchten parameters, and repeats until the error minimized. It is very well documented in the literature that the most critical parameter for a good inversion is theta r. from a practical point of view, this is probably the least important parameter. It has little physical meaning and it is really possible to measure. On the other hand, theta s is one of the most important parameters and can be easily measured. The question is why is the Hydrus model sensitive to theta r and what could be done to reverse this.

Last semester I took a numerical methods class. The last half of the class was on partial differential equations (PDEs) and finite difference method (FDM) solutions. I wasn’t a huge fan of the class because at the time the information seemed intangible and I worked really hard to scrape by. Now that the class is over, I’m glad I took it because I have a much better understanding of modeling. One of the things I couldn’t grasp about Hydrus is that they used a FEM instead of a FDM. Typically, FEM is more complex to implement unless the geometry is complex. The simulations are all 1d so it is not possible to have overly complex geometry.

I spent part of Sunday and all of today implementing a FDM simulation of Richards’ equation. It took several tries because there are some complexities with Richards’ equation that were not in any of the problems I worked on in class. Richards’ equation has two variables, water content and tension, that change in space and time. There are also two functions of tension in the equation. I wanted to use a forward time scheme because this is the easiest to implement. The reason I wanted to use forward time is because I could start at saturation where theta s would be the starting point for the simulation. Because there is a transport component to this phenomena I also used a forward space discretization. After several hours of coding, pondering, and finding errors I arrived at a solution. I have a working code that simulates water movement in a soil sample with water evaporating from the surface. I also now understand why the developers of Hydrus didn’t use this method. The solution has some stability issues. In order to get a stable solution ratio of the resolution of time nodes to the resolution of space nodes needs to be around 1 million. Currently, I have only 10 space nodes but 100,000 time nodes for a 1 day simulation. The code takes about 3.5 minutes to complete the simulation. If I had to iterate the van Genuchten parameters 10 times and extended the time to 10 days this would take 350 minutes or 5.5-6 hours to run 1 sample and the resolution on the space nodes is not adequate. To increase the space nodes to 100 (the same as Hydrus) the simulation would take about 2.5 days to run. This is not practical for my needs.

There are some other backward time schemes that I could implement but that defeats the purpose. I did some searching online for other solutions. Finite volume methods (FVM) look very promising. There have been a few papers published about applying FVM to unsaturated flow problems in porous media. This could be a fun project in the future but as of now I don’t have the required skills to quickly implement the code and I also have a thesis to write (that is not on numerical methods).

Below is a copy of the code. I hope this saves someone some time.

fdm_richards