Fitting the parameters of the RIA parameterization of the soil water retention curve

RIAfitter is a  Fortran2008 code that fits the five parameters of the RIA parameterization of de Rooij et al. (2021) and de Rooij (2022) of the soil water retention curves on observed values of the volumetric water content for an arbitrary number of matric potentials.

Version 2.0, September 2024

Changes with respect to version 1.2.2

- The user can opt between

1) creating a map with regular grid point of the objective function in the parameter space, and having a limited number of parameter fits (hard-coded as 3). The first complexes of two of those runs are filled with points based on that map.

The parameter correlation matrix is calculated for each run, based on a random sample of the last iterations.

2) having a large number of fits (hard coded as 50), and using those to calculate both the covariance and correlation matrix, as well as histograms for each parameter.

The number of fits for each option is set in Module ParameterInformation (RunsWithMap and RunsWithCovariance). To select option 1, set UseMap to T. For option 2, set UseMap to F.

- The user can provide a complete set of fixed parameters and fill the rest of the input files with dummy values. The code then generates a table of the soil water retention curve. In that case, it is assumed that the fitted value of hd is provided. This means that this value of hd is already adjusted in such a way that the target value, at which the water content reaches zero, is equal to (1 + c) * hd.

- The first convergence criterion was relaxed. In earlier versions, the parameter range had to satisfy both the absolute and the relative tolerances. In this version, it has to satisfy at least one of them. This is more realistic for parameters whose permitted ranges span more than an order of magnitude.

- The user can now specify for which parameters the log-transformed value needs to be fitted. It is recommended to select this option for alpha and hd. For the selected parameters, the map of the objective function has an equidistant grid on the log-transformed axes representing these parameters. The calculation of reflection and contraction points also uses the log-transformed axes, and so does the choice of the coordinate values of a random point.

- Parameter statistics, correlations, and covariances are computed using the logarithms of the selected parameters.  

N.B.1 The user-specified parameter range is for the untransformed parameter. The code log-transforms the upper and lower limit.

N.B.2 The user-specified tolerances for the parameters are applied to the log-transformed values directly, so they must be valid for the log-transformed values.

- The error standard deviations for the observations are scaled. The scaling factor that is applied to each error SD is the same, regardless of the variable (water content or matric potential) of which the error SD is scaled. It is hoped that this will create an objective function that is less flat further away from its minimum.

- Routine hjChecker was added to avoid that combinations of parameters can be generated that that would result in hj (the junction point) lying between hae (the air-entry value) and zero.

- Several other minor routines added.

- Compliance with the Fortran-2008 standard is better enforced.

- The calculation of the RMSE is now aligned with the mainstream method, and no longer based on the Chi-squared statistic.

- The implicit equation for the correction factor c imposes limits on the valid combinations of parameter values. These limits are now fully implemented.

- Bugs fixed in BisectionForC, routine rewritten to increase efficiency.

- Correlation coefficients and covariance matrices are based on sample (co)variances instead of population (co)variances.

- If all parameters are fixed, the corresponding table of the retention curve is generated. Still, both input files have to be available, and be filled with numerous dummy values and a few dummy observations.

- Some comments were improved for clarity.

- Sections or lines of code specific to RIAfitter (which would have to be replaced or rewritten when used for other parametric models or entirely different applications) is labeled 'specific to RIAfitter' in the comments.

- Pure functions and subroutines are declared as such.

Version 1.2.2.1, May 2023

Changes with respect to version 1.0:
- Small error in the partitioning of the liquid water content for plotting data fixed.
- Clarity of the comments improved, references to equations in de Rooij (2022) added.
- Calls to WCRIA with dummy variables now use multiple dummies instead of using a single dummy multiple times because not all compilers can handle repeated invocations of the same dummy in a single subroutine call.
- FUNCTION hjRIA now slightly decreases the calculated value of hj so that it is not quite at its theoretical limit.

Changes with respect to version 1.1:
- Estimates of the saturated hydraulic conductivity calculated according to Eqs. (1) and (15) of Timlin et al. (Soil Sci. Soc. Am. J. 1999).
- During every iteration, c and hd are iteratively determined from the value of alpha, n, and the target value of hd to ensure that (1 + c)*hd closely approximates the target hd. Therefore, c and the corrected value of hd are no longer computed after the fitting process is completed but are instead updated during the parameter optimization process by subroutine WCRIA. See the User Manual for details.

Changes with respect to version 1.2 and 1.2.1:
- Bug removed in subroutine BiSectionForC
- Reduced the value of the upper limit for the first iteration of the Bisection algorithm and moved it from subroutine BiSectionForC to MODULE SCEBisParameters.
- Akaike's information criterion calculated and written to output.

Change with repect to version 1.2.2: a bug in the calculation of Akaike's Information Criterion was fixed.

The version of the equation that is fitted is that of de Rooij (2022), in which the matric potential at oven-dryness (hd) and alpha are among the fitting parameters, and the matric potential at the junction point (hj) is a derived parameter.

Each data point needs to have an estimate of the standard deviation of the error of the matric potential, and of the water content. These are used to calculate the weight assigned to each data point.

The code minimizes the root mean square error (RMSE) of the fitted vs. the observed values using shuffled complex evolution. It also  generates a regular grid of map points with the associated RMSE that covers the parameter space. The code evaluates ten convergence criteria for each fitting parameter. Five of these apply to all parameters simultaneously. See Subroutine ConvergenceCheck for details.

To guard against local minima, three runs are performed. If so desired, two of these runs have starting points based on the map of the parameter space. Based on the RMSE, the best fit is determined. If at least two runs return the same RMSE, the run with the largest total number of times a parameter met a convergence criterion is selected. If at least two runs score the the same on both criteria, the run with the smallest rank number is picked (Run 1 over run 2, run 2 over run 3).

The correlation matrix of the fitting parameters is calculated based on a random sample of a prescribed number of fits prior to and including the final fit.

A table with points on the soil water retention curve (and its components) is written to output.

Together with the code, a document that desribes the code and its use is made available.

N.B. The theory underlying version 1.2.2 and up is given in:

de Rooij, G.H. 2022. Technical note: A sigmoidal soil water retention curve without asymptote that is robust when dry–range data are unreliable. Hydrology and Earth System Sciences, 26, 5849-5858, doi 10.5194/hess-26-5849-2022.