Functional Evaluation of Digital Soil Hydraulic Property Maps through Comparison of Simulated and Remotely Sensed Maize Canopy Cover

Abstract

Soil maps can usefully serve in data scarce regions, for example for yield (gap) assessments using a crop simulation model. The soil property estimates’ contribution to inaccuracy and uncertainty can be functionally evaluated by comparing model results using the estimates as input against independent observations. We conducted a functional evaluation of digital maps of soil hydraulic properties of the Zambezi River Basin using a crop growth model AquaCrop. AquaCrop was run, alimented with local meteorological data, and with soil hydraulic properties derived from the digital maps of digital soil mapping (DSM) techniques, as opposed to estimations from the widely used Saxton and Rawls pedotransfer functions. The two simulated time series of canopy cover (CC) (AquaCrop-CC-DSM and AquaCrop-CC-Saxton), which were compared against canopy cover data derived from the remotely sensed Leaf Area Index (LAI) from the MODIS archive (MODIS-CC). A pairwise comparison of the time series resulted in a root mean squared error (RMSE) of 0.07 and a co-efficient of determination (R²) of 0.93 for AquaCrop-CC-DSM versus MODIS-CC, and an RMSE of 0.08 and R² of 0.88 for AquaCrop-CC-Saxton versus MODIS-CC. In dry years, the AquaCrop-CC-DSM deviated less from the MODIS-CC than the AquaCrop-CC-Saxton (p < 0.001), although this difference was not significant in wet years. The functional evaluation showed that soil hydraulic property estimates based on digital soil mapping outperformed those based on Saxton and Rawls when used for simulating crop growth in dry years in the Zambezi River Basin. This study also shows the value of conducting a functional evaluation of estimated (static) soil hydraulic properties in terms of dynamic model output.

1. Introduction

Critical to better soil management is quantitative information detailing the soil resource, its processes, and its variation across landscapes under the broad umbrella of environmental monitoring [1,2]. In this context, digital maps of soil hydraulic properties have been developed using digital soil mapping (DSM) techniques for the Zambezi River Basin (ZRB) [3]. The maps were generated using the machine learning Random Forest model, using covariates related to the SCORPAN concept, whereby soil (S), climate (C), organisms (O), relief (R), parent material (P), and age (A), supplemented with geographic position (N), stand for Jenny’s soil forming factors [2]. The maps were generated for water content at pF0.0 (saturation), pF2.0 (field capacity), pF4.2 (wilting point), and for saturated hydraulic conductivity (Ksat), for three soil depth layers of 30–40 cm, 60–70 cm, and 100–110 cm, and at a spatial resolution of 90 m for the whole ZRB; however, before these maps are used in crop growth or hydrological modelling studies, it is important to check for their performance. In the absence of such confidence, there are uncertainties about the true soil hydraulic properties and processes of the considered soils [4]; hence, it is vital to determine to what extent uncertainties are propagated through the DSM and how these affect subsequent application models relying on these data as input variables.

According to [4,5], in modeling exercises, uncertainty of the model output is a combination of three main sources including model structure uncertainty, model parameter uncertainty and model input uncertainty. These uncertainties in the DSM models can be investigated through functional evaluation; for instance, by checking for the DSM model performance in a crop model such as the FAO’s water driven AquaCrop. Obviously, DSM techniques for estimating soil hydraulic properties at a specific place are not error-free. Depending on the intended model application, the significance of the resulting uncertainty will vary. When the estimation of hydraulic properties is not the final aim, there is scope for evaluating the usability of DSMs by evaluating the accuracy of model output (e.g., from crop growth or hydrological models), to which the DSM-estimates are inputted [6]. This is what is meant by functional evaluation [7,8,9]. References [9,10] recommended using a functional criteria analysis, which analyzes the functional behavior of estimated parameters in various applications. By analyzing how model parameter uncertainties propagate through the stochastic model, a functional analysis combined with a stochastic model of soil hydraulic properties can be used to effectively analyze the relative effectiveness of alternative parameterization methods [11,12]. Some authors have discussed the importance of functional analysis for the uses of pedotransfer functions (PTFs) and have reported the importance of functionally evaluating PTFs for real field-scale applications [9,10,13,14]. Studies have been conducted to evaluate the functionality of PTFs for environmental simulations, such as water flow simulations [14,15], water balance studies [11,12,16,17], irrigation water estimations [18,19] and crop yields [19,20,21], using different models.

However, in order to conduct functional evaluation through a model such as AquaCrop, a functional criterion is needed [18,21,22]. This criterion can come in the form of any reference or measured data such as crop yield, crop water requirements [22] or the crop Canopy Cover (CC). The AquaCrop crop model simulates CC to describe evolving phenological crop development. Through its expansion, aging, conductance, and senescence, CC determines the amount of water potentially transpired, which in turn determines the amount of biomass produced and the final crop yield. If water stress occurs, the simulated CC will be less than the potential canopy cover for no stress conditions [23]; however, this type of approach comes with its own challenges, in that most of the reference data or functional criterion (such as CC) needed to perform functional evaluation is rarely available, especially in places such as the ZRB.

Therefore, remote sensing can be used for field surveys providing reference data such as CC. The collection of spatially explicit crop data for such large areas could be very costly, tedious, and time-consuming, and therefore, the combined use of satellite imagery and ground data collection to monitor crop growth conditions, and subsequently predict yields, is not only convenient, but also economical [24,25]. Among the many available remote sensing products, the Moderate resolution Imaging Spectrometer (MODIS) products are freely available datasets that provide considerable advantages over other remote sensing products in crop monitoring studies [26,27,28,29,30,31]. MODIS provides for many vegetation indices including the normalized difference vegetation index (NDVI), which is an indicator of crop vigor as well as the Leaf area index (LAI), and is a useful biophysical variable defined as the sum of all one-sided green leaf areas divided by the corresponding area on the ground covered by the plant and is generally used for monitoring crop development and estimating crop yields [32,33,34]. The possibility of using remote sensing data instead of field data makes most crop growth models more applicable for use in data-scarce areas such as the ZRB, whereby LAI derived from MODIS data can be assimilated in crop growth modelling studies as a steering factor (reference or observed data and also used as a functional criterion in this study), in which a crop growth model such as AquaCrop is re-parameterized and calibrated to reduce the difference between the model value and the satellite-derived value [29,35,36]. For this study, the Leaf Area Index (LAI, m²m⁻²) derived from MODIS data MCD15A3H Version 6 [37] was used as reference data (functional criterion) for the FAO’s AquaCrop crop growth model after conversion to crop canopy cover (CC, m²m⁻²), for the purposes of performing functional evaluation of the soil hydraulic property maps that were developed using DSM techniques and the machine learning random forest model in another study. The aim of this study, however, was to investigate to what extent the different procedures for mapping the soil hydraulic properties get reflected in a crop-growth model, and to verify if these differences can be noticed by RS data—in this case MODIS. The assumption is that the crop growth model that best reflects the crop growth in the field will also have the highest correspondence with the NDVI value recorded by MODIS by considering the following specific objectives:

• Evaluate and compare the time series of the AquaCrop simulated maize canopy cover (CC) with the time series of the canopy cover derived from the MODIS-satellite LAI-product (MODIS-CC), whereby AquaCrop is alimented by the available digital soil hydraulic property maps (AquaCrop-CC-DSM).
• Investigate whether the AquaCrop-CC-DSM is closer to the MODIS-CC than the CC time series generated by the AquaCrop alimented with the soil hydraulic properties that were estimated by the widely used PTFs of Reference [38].
• Examine whether the performance of the AquaCrop-CC-DSM and the AquaCrop-CC-Saxton depend upon the reference soil group (RSG) and/or upon the rainfall abundance.

2. Materials and Methods

2.1. Study Area and Land Units

The study area coincides with the part of the Zambezi River Basin falling within the Republic of Zambia and is located between 8–18° S latitude and 22–33° E longitude (Figure 1). The spatial resolutions of available soil hydraulic properties (SHPs) data are not sufficient for application in the ZRB, and only a few initiatives such as the SoilGrids information system [39] provide global predictions of SHPs at a resolution of 250 m by 250 m for water only content at pF0.0 and the available water capacity (AWC). Reference [40] also generated a worldwide collection of soil hydraulic properties data and sub-grid variability of soil water retention and hydraulic conductivity curves at 1 km by 1 km resolutions, whereas [41] developed a dataset of the AWC for Africa at a resolution of a 1 km by 1 km scale. Highly weathered soils (Ferralsols) are predominant in the northern section of the ZRB on the higher parts of the plains, and they plateau predominantly in areas constituted of Precambrian basement rocks. The creation of the rift valleys and adjacent mountains has a significant impact on the soils of the eastern section of the basin. These soils are less weathered (e.g., Luvisols, Lixisols) or weakly developed (Cambisols), and are often rocky, stony, or shallow (Leptosols, Regosols). In the driest part, particularly in northern Zimbabwe, soils occur that are affected by the accumulation of soluble salts (Solonchaks, Solonetz). Relatively fertile soils occur on the recent volcanic parent materials (Andosols) in the northern part of ZRB in Tanzania, where older volcanic or calcareous rocks have weathered the fertile black soils of the cracking clay type which have formed (Vertisols).

Most of the data preparation and pre-processing for the AquaCrop model inputs [22] were done using a combination of R-software, version 3.5.0 [42] and ArcGIS-software, version 10.7. The AquaCrop model simulations were conducted for 12 locations (land units) in Zambia (Figure 1). These were 12 pure pixels (500 m × 500 m) not under irrigation and were assumed to be a monoculture of maize during the rainy season with fallow land between the maize growing seasons, and they represented mapping units of six different dominant Reference Soil Groups. The land units were generated using ArcGIS-software by combining a climate zone map from the Global Yield Gap and Water Productivity Atlas (GYGA) [43], a land cover map of Africa obtained from http://2016africalandcover20m.esrin.esa.int/ (accessed on 20 November 2021) for extracting agricultural areas, a slope map derived from a Digital Elevation Model (DEM) obtained from http://srtm.csi.cgiar.org (accessed on 20 November 2021), and a soil map based on [44] obtained from the International Soil Reference and Information Centre (ISRIC) SoilGrids web interface at https://soilgrids.org/ (accessed on 20 November 2021).

2.2. AquaCrop Model Inputs and Outputs

AquaCrop is a physical crop water productivity model that simulates crop development and production under various environmental conditions and management practices. The model has been calibrated for more than 30 different crops, including widely cultivated crops such as maize, wheat, and sugarcane, as well as underutilized crops such as bambara groundnut and tef. As AquaCrop keeps an optimal balance between accuracy, robustness, and simplicity, it only requires a limited number of easily obtainable input parameters [23], which makes the model applicable even in data-scarce conditions. AquaCrop has been applied multiple times to assess irrigation, crop, and field management in African cropping systems in the ZRB. The AquaCrop model requires crop input data and observations for parameterization, soil properties, field management practice data, and a time series of weather data, that define the environment in which a crop develops [23].

2.2.1. Climate Data

We used weather data from the Global Yield Gap and Water Productivity Atlas (GYGA) [43] for 6 climate stations in Zambia (Figure 1). The data for the period of 1998–2012 consisted of daily records of rainfall, minimum and maximum temperature, mean relative humidity, mean wind speed, and solar radiation. The time series of daily reference evapotranspiration (ETo), which AquaCrop also requires, was calculated from the minimum and maximum temperature, mean relative humidity, mean wind speed, and solar radiation, according to the FAO Penman Monteith equation [45], by applying the ETo calculator software [23].

2.2.2. Crop Data

We selected maize (Zea mays L.) to simulate crop canopy cover (CC) since it is the most common crop grown by small-holders and large commercial farmers throughout the ZRB during the rainy season. We downloaded MODIS (MCD15A3H Version 6) LAI time series data [37] for the period 2002–2012 by means of the R-MODIStsp package. This data comes as a rasterstack time series dataset that is already pre-processed (i.e., georeferenced and corrected for atmospheric interference). The MCD15A3H Version 6 MODIS Level 4, Combined Fraction of Photosynthetically Active Radiation (FPAR), and Leaf Area Index (LAI) product, is a 4-day composite data collection with a pixel size of 500 m. Within the 4-day period, the compositing algorithm selects the best pixel from all acquisitions of both MODIS sensors on NASA’s Terra and Aqua satellites https://lpdaac.usgs.gov/products/mcd15a3hv006/ (accessed on 20 November 2021) [37], which both have a temporal resolution of 8 days (

Table 1. Description of the standard MODIS LAI/FPAR products.

Official Name	Platform	Raster Type	Spatial Resolution	Temporal Granularity
MOD15A2H	Terra	Tile	500 m	8 Days
MYD15A2H	Aqua	Tile	500 m	8 Days
MCD15A2H	Terra + Aqua Combined	Tile	500 m	8 Days
MCD15A3H	Terra + Aqua Combined	Tile	500 m	4 Days

).

The MODIS LAI rasterstack time series data was overlaid with the land units layer to extract the stack of pixels with the LAI-value corresponding to the center of the 12 land units. As the AquaCrop model produces CC rather than LAI as one of its outputs, the LAI gathered from the MODIS data had to be converted to CC in order to compare the AquaCrop CC outputs. In general, the relation between LAI and CC is based on Beer’s law (commonly used to describe the relationship between the proportion of light penetrating a plant canopy and the leaf area index), and according to [46], it can be expressed as in Equation (1).

$$\mathrm{CC}(\%)=1.005\times {\left[1-{\exp }^{\left(-0.6\mathrm{LAI}\right)}\right]}^{1.2}$$

(1)

with CC as the green crop Canopy Cover percentage and LAI the MODIS Leaf Area Index (m²m⁻²). Before comparing the MODIS CC with the AquaCrop simulated CC, the MODIS CC time series data were filtered using a moving average with parameters (p = 1, n = 3) of the Savitzky-Golay Filter [47,48,49]. The filtering was done because the atmospheric effects on surface reflectance generally cause negatively biased noise within the LAI values of MODIS [48]. The AquaCrop model generates CC outputs using a daily time scale; however, the converted MODIS CC was available on a 4-day time scale. To ensure a fair comparison between the AquaCrop CC and the MODIS CC, 4-day values for the AquaCrop time series were extracted from the daily AquaCrop time series using R by selecting a value of CC every 4 days.

In AquaCrop, crops are characterized by a set of parameters that can easily be calibrated by the modeler to meet specific local growing conditions [22]. We selected default crop parameters for the length of a growing period to reach maturity (1517 growing degree days) and crop planting density (49,231 plants/m²) because this was a comparison study. For the planting date, we used the rainfall criterion used for Southern Africa [50] which was: first day of the first 10-day period since the 1st of July with a cumulative precipitation of at least 25 mm. Prior to the planting date, the actual AquaCrop simulations began on 1 January 1998. At this time, because most portions of the ZRB are already in the peak rainfall period, the initial soil water content was considered to be at field capacity throughout the soil profile.

2.2.3. Soils and Soil Hydraulic Properties Data

The AquaCrop simulations for the 12 land units were conducted on soil profiles representative for the Reference Soil Groups (RSG) according to the 3rd edition of the World Reference Base for soil resources [51]. In the 12 land units, agricultural land occurs on the following six RSGs: Acrisols, Arenosols, Cambisols, Ferralsols, Luvisols, and Podzols. The soil textural classes of each RSG are presented in

Table 2. Soils and soil hydraulic properties data for AquaCrop simulations as generated through DSM (based on ANN derived PTFs and Random Forest regression with Residual Kriging), and the Saxton and Rawls PTF.

Land Unit, RSG and Climate Station	Soil Profile Depth Layer (cm)	FAO Textural Class	Soil Hydraulic Properties Data Used in AquaCrop Simulations
			Digital Soil Mapping (DSM)					Saxton & Rawls PTFs
			pF0.0 vol %	pF2.0 vol %	pF4.2 vol %	AWC vol %	Ksat mm/h	pF0.0 vol %	pF2.0 vol %	pF4.2 vol %	AWC vol %	Ksat mm/h
1 Acrisol Kabwe	30 60 100	Sandy loam Sandy Clay loam Sandy Clay loam	31.7 29.6 31.1	17.3 17.1 17.4	6.7 6.1 7.0	10.6 11.0 10.4	204.6 201.8 185.8	38.8 39.5 40.4	20.0 26.2 29.3	12.0 16.8 19.2	8.0 9.4 10.1	15.7 5.8 3.4
2 Ferralsol Kabwe	30 60 100	Sandy Clay loam Sandy Clay loam Clay loam	36.6 36.5 37.3	22.6 21.5 21.8	9.1 9.7 10.3	13.5 12.0 11.5	100.8 66.4 69.3	40.0 41.3 42.9	26.2 31.1 34.6	15.7 19.8 22.7	10.5 11.3 11.9	6.6 2.7 1.45
3 Ferralsol Mpika	30 60 100	Sandy Clay loam Sandy Clay Sandy Clay	37.7 36.5 38.3	23.1 22.2 21.3	9.8 9.6 9.4	13.3 12.6 11.9	120.5 98.5 101.9	40.3 41.3 42.4	28.1 31.7 34.6	18.1 21.1 23.4	10.0 10.6 11.2	4.5 2.2 1.2
4 Luvisol Mumbwa	30 60 100	Sandy Clay loam Sandy Clay loam Sandy Clay	37.0 36.6 37.4	25.6 24.9 23.5	11.6 10.2 10.7	14.0 14.7 12.8	69.8 71.8 73.8	39.7 41.0 42.1	25.1 30.0 33.0	15.1 19.3 21.6	10.0 10.7 11.4	7.8 3.3 1.9
5 Acrisol Mumbwa	30 60 100	Sandy Clay loam Sandy Clay loam Clay loam	38.2 38.8 39.5	26.5 26.6 23.8	11.3 11.3 11.5	15.2 15.3 12.3	68.3 50.3 48.0	40.0 41.6 43.0	26.7 31.3 34.6	16.3 19.9 22.8	10.4 11.4 11.8	5.9 2.8 1.5
6 Acrisol Choma	30 60 100	Sandy loam Sandy Clay loam Sandy Clay loam	30.5 29.3 30.6	18.3 17.9 17.6	6.8 6.7 7.2	11.5 11.2 10.4	163.1 190.4 197.1	38.7 38.9 39.5	19.0 23.6 26.2	11.4 15.0 16.8	7.6 8.6 9.4	18.4 8.6 5.7
7 Acrisol Choma	30 60 100	Sandy loam Sandy Clay loam Sandy Clay loam	28.5 28.3 28.7	21.1 21.0 20.8	7.2 7.0 7.3	13.9 14.0 13.5	208.4 232.5 242.2	38.7 38.8 39.0	18.6 22.5 24.1	11.4 14.4 15.6	7.2 8.1 8.5	19.3 10.3 7.9
8 Arenosol Livingstone	30 60 100	Sandy loam Sandy loam Sandy loam	33.3 32.1 32.0	21.2 21.4 21.0	9.3 8.8 7.9	11.9 12.6 13.1	187.8 180.0 225.8	38.9 39.1 39.5	22.8 24.5 26.1	14.4 15.6 16.8	8.4 8.9 9.3	10.1 7.6 5.9
9 Arenosol Livingstone	30 60 100	Sandy loam Sandy loam Sandy loam	30.3 31.5 31.4	20.5 20.8 20.3	7.8 8.0 8.0	12.7 12.8 12.3	182.5 198.8 233.5	38.9 39.1 39.4	22.5 24.9 26.5	14.4 16.2 17.4	8.1 8.7 9.1	10.6 6.9 5.2
10 Podzol Mongu	30 60 100	Sandy loam Sandy loam Sandy loam	26.3 26.0 26.1	17.5 16.6 16.0	5.2 5.2 5.4	12.3 11.4 10.6	336.7 351.4 342.2	39.2 38.5 38.5	12.6 17.1 18.1	6.9 10.6 11.3	5.7 6.5 6.8	44.7 22.8 19.9
11 Arenosol Mongu	30 60 100	Sandy loam Sandy loam Sandy loam	25.7 25.7 25.8	15.8 14.7 14.0	4.6 4.4 4.6	11.2 10.3 9.4	386.1 392.5 383.2	39.4 38.7 38.5	12.1 13.4 14.8	6.4 7.5 8.7	5.7 5.9 6.1	49.2 38.8 31.6
12 Cambisol Mongu	30 60 100	Sandy Clay loam Sandy Clay loam Sandy Clay	36.9 36.1 36.8	21.8 21.4 21.5	8.3 9.0 9.8	13.5 12.4 11.7	84.7 78.1 75.2	39.7 40.8 42.2	26.9 30.8 34.4	17.5 20.4 23.4	9.4 10.4 11.0	5.2 2.5 1.2

.

We used the International Soil Reference and Information Centre (ISRIC) or the SoilGrids web interface to download basic soil datasets such as soil granulometry fractions (sand, clay, and silt), soil organic carbon content, and bulk density, with a spatial resolution of 250 m at depth layers of 30, 60, and 100 cm to compare two methods of estimating soil hydraulic properties data [39] at (https://soilgrids.org/, accessed on 20 November 2021). Each of the downloaded datasets was then overlaid with the 12 land units, and using zonal statistics with the mean function, the extracted values were used as input for the Saxton and Rawls PTFs to estimate soil hydraulic properties (Table 2) for the first method. In the second method, the soil hydraulic properties data at depth layers of 30–40 cm, 60–70 cm, and 100–110 cm were obtained from the digital soil maps of hydraulic properties derived based on Digital Soil Mapping (DSM) techniques. Each map layer (Figure 2) was also overlaid with the 12 land units, and using zonal statistics with the mean function, we extracted soil hydraulic properties values (Table 2) that were used as input in the AquaCrop model.

2.2.4. Simulation Outputs

Two daily CC time series covering the 1 January 2002 to 31 December 2012 period were generated for each of the 12 land units by means of AquaCrop. One was obtained from AquaCrop operated with soil hydraulic properties data (Figure 2 and Table 2) derived from DSM, hence referred to as AquaCrop-DSM-CC, whereas the other one, referred to as AquaCrop-Saxton-CC, was generated from AquaCrop operated with soil hydraulic properties data (Table 2), which were derived from the [38] PTFs. The two daily time series were eventually converted to a 4-day time step to match the temporal resolution of the MODIS CC time series.

2.2.5. Statistical Analysis

The MODIS CC time series were used as the reference so that the two comparisons considered were (a) AquaCrop-CC-DSM versus MODIS-CC and (b) AquaCrop-CC-Saxton versus MODIS-CC. The AquaCrop CC-simulations were evaluated against the MODIS-CC-time series using the root mean squared error (RMSE) and the co-efficient of determination (R²). To determine whether the deviation between AquaCrop-DSM-CC and MODIS-CC on the one hand, and AquaCrop-Saxton-CC and MODIS-CC on the other hand are statistically significant, an Analysis of Covariance (ANOCOVA) was conducted. The ANOCOVA model in Equation (2) was put forward. The quantitative dependent variable was the deviation of either the DSM-CC or Saxton CC time series from the MODIS CC time series, all based on a 4-day time step, whereas the independent variables were the Reference Soil Group (RSG), the rainfall (also on a 4-day time step) and interactions between the RSG and rainfall. The rainfall variable was primarily included in the equation because the maize crop simulations in the AquaCrop model were rainfed rather than irrigation based.

$$\mathrm{Deviation} \mathrm{of} \left({\mathrm{CC}}_{\mathrm{DSM} \mathrm{or} \mathrm{Saxton}} \mathrm{from} {\mathrm{CC}}_{\mathrm{MODIS}}\right)=\mathrm{RSG}+\mathrm{Rainfall}+\mathrm{RSG}\times \mathrm{Rainfall}$$

(2)

Furthermore, in order to identify dry, normal, and wet years from the six meteorological stations, the RAINBOW-a software package for hydrometeorological frequency analysis and for testing the homogeneity of historical data sets was used [52].

3. Results

3.1. Maize Canopy Cover (CC)

During the 2002–2012 growing seasons, there was a good agreement of the MODIS CC time series with both the DSM and Saxton CC time series for all the 12 land units (Figure 3). The deviations of either the DSM or Saxton CC time series from the MODIS CC time series were mainly observed in the growing seasons with lower rainfall events or dry spell periods. For the land units 1 (Acrisol) and 2 (Ferralsol), the deviations in the dry spell periods were not as pronounced as those for the rest of the land units.

3.2. Statistical Analysis

For the 12 land units, a mean RMSE = 0.07 and a mean R² = 0.93 was observed for the MODIS CC time series versus the DSM CC time series, whereas a mean RMSE = 0.08 and mean R² = 0.88 was obtained for the MODIS CC time series versus the Saxton CC time series (Figure 4). Although, the AquaCrop model performance was slightly better with the soil hydraulic properties data obtained from DSM maps than with those from the Saxton and Rawls PTFs. Figure 4 seems to suggest that there was no marked difference between the DSM and the Saxton methods, but rather, that there are two different population effects, one leading to a straight line and the other leading to an S-shaped curve.

Table 3. Intercepts, coefficients, and p values for the ANOCOVA tests based on Equation (2).

	Digital Soil Mapping (DSM)		Saxton & Rawls PTFs
Variable	Coefficients	p Values	Coefficients	p Values
Intercept	0.312	0.754	−0.879	0.379
Rainfall	−6.288	3.330 × 10−10 ***	−6.990	2.910 × 10−12 ***
Arenosol	4.459	8.320 × 10−06 ***	1.138	0.255
Cambisol	−0.349	0.727	−3.407	6.580 × 10−4 ***
Ferralsol	−5.044	4.620 × 10−07 ***	−3.726	1.960 × 10−4 ***
Acrisol	−1.131	0.258	−2.050	0.040 *
Luvisol	4.883	1.060 × 10−06 ***	4.036	5.480 × 10−05 ***
Podzol	10.266	<2.000 × 10−16 ***	6.797	1.120 × 10−11 ***
Arenosol * Rainfall	6.794	1.150 × 10−11 ***	4.777	1.800 × 10−06 ***
Cambisol * Rainfall	2.717	0.007 **	−0.340	0.734
Ferralsol * Rainfall Acrisol * Rainfall Luvisol * Rainfall	−3.429 0.650 5.941	6.080 × 10−4 *** 0.516 2.910 × 10−09 ***	−2.336 −0.726 5.188	0.019 * 0.467 2.16 × 10−07 ***
Podzol * Rainfall	9.499	<2.000 × 10−16 ***	7.443	1.060 × 10−13 ***

Significant codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

for the ANOCOVA tests indicates that the RSG, the rainfall, and interactions between the RSG and rainfall variables each have a significant impact on the deviations of either the DSM CC or the Saxton CC from the MODIS CC data points. Most of the p-values of the RSG, rainfall, and interactions between the RSG and rainfall variables were all p < 0.001 for the deviations of either the DSM CC or the Saxton CC from the MODIS CC. Furthermore, using the intercepts and coefficients in Table 3, on average, for a every 3 mm of rainfall received on a RSG, the deviations of the DSM CC from the MODIS CC were about −4.68%, whereas those for the deviations of the Saxton CC time series from the MODIS CC time series were about −13.31%. Moreover, the deviations of the Saxton CC time series from the MODIS CC were more pronounced and statistically significant (p-values < 0.001) in the dry years (2002–2003, 2005–2006, and 2011–2012 of Figure 3), implying that overall, the DSM outperformed the Saxton and Rawls PTF method when it came to estimating the soil hydraulic properties data.

4. Discussion

Time series of canopy cover is an important intermediate output of water-driven crop growth models such as the AquaCrop model, and therefore, an accurate estimation of this variable is essential for the model to produce good/acceptable estimates of crop evapotranspiration, biomass, and yields [53]. In this study, two methods of estimating soil hydraulic properties data were functionally evaluated by using the SHPs in AquaCrop and analyzing the simulated CC. AquaCrop model simulations were conducted for 1998 to 2012 with soil hydraulic properties data derived either from DSM or from the Saxton and Rawls PTFs; hence, we obtained from the AquaCrop model two simulation outputs: a maize crop CC time series based on the DSM derived soil hydraulic properties data, and a maize crop CC time series based on the Saxton and Rawls PTFs soil hydraulic properties data. The AquaCrop model performance was evaluated using the RMSE and R² obtained by comparing each of both obtained CC time series with a MODIS CC time series for 12 considered land units.

Whereas the statistical indices indicate a relatively good performance of the AquaCrop model, the results seem to suggest that there was no marked difference between the DSM and the Saxton methods, but rather that there are two different groups of years, one leading to a straight line, and the other leading to an S-shaped curve (Figure 4); however, the ANOCOVA tests showed that the deviations of the Saxton CC time series from the MODIS CC were more pronounced than the DSM-CC from the MODIS CC and statistically significant with (p-values < 0.001) in the dry years. This was most likely (Table 2) due to the crops’ water stress, which was more severe in dry years than in the wet years. In the wet years, the difference between the DSM and the Saxton technique for calculating soil hydraulic properties is likely to be negligible because no water stress occurs in a wet year, so the difference in SHPs is of lesser importance for CC (water content does not drop below thresholds for canopy expansion and early leaf senescence, below which growth reduction occurs). Furthermore, when compared to the DSM method, the significant deviations for the Saxton and Rawls PTFs were much larger in dry years because Saxton produces lower values of the available water capacity and saturated hydraulic conductivity, but slightly larger values of the water content at pF0.0, pF2.0, and pF4.2 (Table 2).

We may conceive a few possibilities for why the S-shaped curves in Figure 4 appear: (1) the rainfall variable’s cyclic change between dry and wet seasons; (2) a discrepancy between the AquaCrop CC and the MODIS CC at the start of the growth season; (3) a chosen land unit might not be a pure agricultural pixel; (4) the crop on a particular land unit may not have been maize, and/or the NDVI may be spiked by, e.g., weeds; and lastly (5) in the months between May and sometime in November, MODIS CC may show higher values in the beginning of the season, because for AquaCrop simulations, there were no crops on the field; hence, there were zero values for CC, but MODIS generates some CC values in this period, which probably came from other vegetation such as weeds, and other crops that farmers plant before the start of the rain season, sometime in the month of November.

The findings of this study reveal that the DSM maps produce estimates of soil hydraulic properties, which, when fed into the AquaCrop crop model, produce maize CC estimations that are closer to the MODIS CC, compared with the AquaCrop CC, simulated with the SHPs that were estimated by the Saxton and Rawls PTF. The DSM estimates of the SHPs differ from the PTF estimates at three levels: (i) the deterministic estimation model type (Random Forest versus Multiple Linear Regression), (ii) the predictor variables (DEM-derived, climatic, SoilGrids-derived etc., versus SoilGrids only), (iii) the accounting for the spatial autocorrelation of the deterministic model residuals when they are present, versus ignoring this spatial autocorrelation.

Apart from using the AquaCrop simulated CC as the criterion for the functional evaluation of the DSM versus Saxton in crop growth modelling, other variables can be considered, with the condition that the reference data are available, e.g., biomass and yield of the considered crop [20,21]. Such information was not available for the studied land units, nor can it easily be derived from available remotely sensed imagery. Other possible variables for functional evaluation relate to the terms of the soil water balance, which in a water driven crop growth model, have a direct relationship with model outputs such as CC, biomass, and yield. Again, reference data must be available for such variables to assess the usability of the SHP estimates. One such variable is the water content of the topsoil which is simulated by AquaCrop for each time step and can be derived through microwave remote sensing [54]. It can be argued that the water content of the topsoil is more directly influenced by the SHPs (especially Ksat) than the canopy cover, and hence, would be a more appropriate variable for functional evaluation. Our reasoning was that the canopy cover (and biomass and yield) is the more pertinent variable of interest in the context of crop growth modelling, whereas for this context, topsoil water content is a rather intermediate indicator.

5. Conclusions

In modelling, researchers often limit themselves to comparing their modelled data with observed data, but this gives little insight into the relevance of the modelled data for the land manager. By performing a functional evaluation, we aimed to gain insight into the applicability of the crop modelled results; however, in order to conduct functional evaluation through a soil water balance-based model such as AquaCrop, a functional criterion is needed. This criterion can come in the form of model output deviation from any reference or measured data such as crop yield, the temporal profile of crop water requirements, or crop canopy cover. As CC is closely related to biomass production and yield, and because of the availability of reference data through the MODIS-LAI time series, we used the CC time series of the maize crop as the functional criterion. The location specific CC time series were derived from the time series of remotely sensed LAI, as available in the MODIS archive. AquaCrop alimented with local meteorological data and with soil hydraulic properties data derived from digital maps, and the widely used standard Saxton and Rawls PTFs, were run to generate two time series of AquaCrop simulated CC for each considered location so that pairwise comparison of time series could be made. These comparisons resulted in a RMSE of 0.07 and an R² of 0.93 for AquaCrop-CC-DSM versus MODIS-CC, and a RMSE of 0.08 and R² of 0.88 for AquaCrop-CC-Saxton versus the MODIS-CC time series. An ANOCOVA test showed that the deviations of the Saxton based CC time series from the MODIS CC were more pronounced and statistically significant (p-values < 0.001) in the dry years than the DSM based CC time series. The findings of this study reveal that the developed DSM maps produce estimates of soil hydraulic properties data, which, when fed into the AquaCrop crop model, produce maize CC estimations that are closer to the MODIS CC; therefore, this study recommends the use of soil hydraulic properties data derived from our DSM maps in crop growth models in favor of other available data products; however, instead of using the CC as the functional criterion in future functional evaluation studies, other variables such as the topsoil moisture content should be used.