The study area spans four important agricultural areas in Morocco, which have a marked contrast regarding ecological, pedological and climate condition. These regions are: Doukkala, Gharb-Loukous, Tadla and Moulouya (Figure 1).

Doukkala is a region located in the west-central Morocco, covering an area of 1,328,500 ha. Agricultural area is 428,000 ha, divided into rainfed agriculture (75%) and irrigated agriculture (25%) (Doumali and Ibno Namr, 2018). This region is characterized by a semi-arid Mediterranean climate, with rainy winter and dry summer. Rainfall distribution is irregular throughout the months and years, with an annual average of about 322 mm, mostly occurring from October to February. The soils of the region are constituted by six types: poorly developed soils (16.5%), vertisols (30%), calcimagnetic soils (9.4%), isohumic soils (36%), soils with iron sesquioxides (6.7%), and hydromorphic soils (1.4%) (Doumali and Ibno Namr, 2018).

Gharb-Loukous region is constituted by the plain of Gharb and the watershed of Loukous. Gharb -Loukous region is in the northwest Morocco. It is a large plain with very flat relief and an altitude bellow 50 m (Bounouira et al., 2018). Agricultural area is 616,000 ha, divided into rainfed agriculture (40%) and irrigated agriculture (60%). This region is impacted by the oceanic influence, and according to Emberger classification, it is a subhumid bioclimatic stage with temperate winter with a gradual transition to the semi-arid level going east (FAO, 2015). Gharb plain is constituted by the following soil types: poorly developed soils (33%), vertisols (57%), hydromorphic soils (8.5%), and calcimagnetic soils (1.5%).

Moulouya region is located in the lower of the Moulouya watershed, which is the largest watershed in North Africa (Tovar-Sánchez et al., 2016). Climate in the basin is typically Mediterranean to semi-arid, characterized by large variability in precipitation with yearly accumulated values ranging between 230 and 380 mm.

Tadla region is in the central Morocco, with an area of 360,000 ha constituted by the Tadla plain. The average altitude is 400 m. The climate is arid to semi-arid. The dry season extends from April to October, while the wet season occurs during the period from October to March. The summer is very hot because of the scorching winds from the south-west called Chergui which raises the temperatures above 40°C. Tadla plain is constituted by the following soil types: isohumic (83%), calcimagnetic soils (11%), and hydromorphic soils (6%).

Disturbed soil, sampled from topsoil (0–20 cm). We selected to sample the topsoil because this layer presents the layer most affected by the soil management systems principally the alteration in soil organic matter, compared deeper layers. The samples were collected from four agricultural regions in Morocco and georeferenced using a handheld GPS (model Garmin GPSMAP 64s), with a resolution of 3.65 m. These samples were collected in different periods from 2019 to 2022. Total of 331 samples were collected from four regions in Morocco: Doukkala (78 samples), Gharb-Loukous (116 samples), Moulouya (73 samples), and Tadla (74 samples). Location of studied regions is shown in Figure 1. In the laboratory, samples were air-dried, crushed, mixed, and sieved through 2-mm mesh. Soil samples were classified following the USDA soil texture classification, and the textural triangle was plotted using R package soil texture ‘soiltexture’ (Moeys et al., 2018).

Determination of soil organic carbon was performed using sulfochromic oxidation of carbon in a mixture of potassium dichromate (K₂Cr₂O₇) and sulphuric acid (H₂SO₄) at of 135°C, according to Walkley and Black (1934). Particle size distribution was conducted using a hydrometer method (Bouyoucos, 1962), using a sodium hexametaphosphate with a concentration of 5% as dispersing solution, and recording the hydrometer reading at 40 s, 2 h, 4 h, 6h, and 24h. Soil water content at field capacity (FC) and permanent wilting point (PWP) were measured following the methodology proposed by Klute (1986). Disturbed samples using approximately 40 g of soil put in soil rings (0.053 m of diameter and 0.01 m of height) without packeting. Soil samples were saturated by capillarity with water up to two-thirds of the ring height, and submitted after saturation, to the matric potentials of -33 kPa (FC) and -1,500 kPa (PWP) (Klute, 1986; Vomocil, 1986). The use of disturbed samples to determine FC and PWP was used in different previous studies (Shaykewich, 1970; Jong et al., 1983). Available water capacity for plants (AWC) was calculated using the Equation (1). Water content at saturation, FC and PWP is expressed in gravimetric moisture content as mass of water by mass of dry soil, presented by w ₀ at saturation, w ₃₃₀  at FC, and w ₁₅₀₀₀  at PWP_. A W C = w 330 – w 15000 (1)

The field dataset was compared with the regional soil information iSDAsoil dataset, which was from a soil map published by Hengl et al. (2021), with 30 m of resolution. To obtain soil information regarding clay, sand, silt, and total soil organic carbon, digital soil maps were downloaded from the site https://www.isda-africa.com/isdasoil/. For the calculation of the soil organic matter (SOM) from total soil organic carbon, we used the factor of 1.724 by multiplying the total soil organic carbon to obtain SOM (Khatoon et al., 2017). To avoid the effect of spatial autocorrelation, we used average nearest neighbor analyses to remove spatially correlated data (Fortin et al., 2001). After this selection, 282 samples were kept from 331, to compare between both sources of data (field and iSDAsoil); for the modeling all samples (331) were included.

PTFs were evaluated based on gravimetric soil moisture. We selected PTFs used and tested in different studies from literature to predict soil moisture at field capacity (w ₃₃₀ ) and permanent wilting point (w ₁₅₀₀₀ ) (Pidgeon, 1972; Lal, 1978; Arruda et al., 1987; Dijkerman, 1988). The selection of PTFs was based on the selection of PTFs that used soil texture, and soil organic matter as input variables, and gravimetric soil moisture as output. PTFs assessed in this study are summarized in Table 1.

For the evaluation of the PTFs, we followed two procedures. The first procedure is the calibration procedure by building models using the equation developed by other authors, selected in the Table 1 (Contreras and Bonilla, 2018). The second procedure is the validation procedure, by changing the parameters of the PTF equation (Contreras and Bonilla, 2018). Soil database was separated in two subsets, 50% of samples were used for calibration, and 50% of samples were used for validation. Random numbers were generated to split and to select points from each subset as recommended by Minasny and McBratney (2002), using R software (R Core Team, 2017).

In general, the PTFs were based on simple linear regression, to look for the relationship between the estimated variable and an available soil information property (soil texture or soil organic matter). We tested different models of simple linear regressions using different soil properties using the version 4.1.3 of R software (R Core Team, 2017). In addition, we tested new PTFs by using multiple linear regression following a stepwise procedure. For this, we included all soil variables in the model, then we selected models that presented high coefficient of determination (R²) and significant p-value. Furthermore, we evaluated models based on machine learning and non-linear algorithms such as regression tree, Cubist algorithm and random forest.

Regression tree was described by Breiman et al. (1984), as an algorithm using a binary recursive partitioning to split the data into partitions or branches, and then continues splitting each partition into smaller groups as the method moves up each branch.

Random forest is an algorithm constructed using an ensemble of decision trees. The tree is created randomly from different samples of rows and different sample of features is selected for splitting. Each tree performed its own individual prediction, then the prediction of the model is averaged using the prediction of each tree (Breiman, 2001).

Cubist algorithm is a rule-based model based on the predictors used in previous splits, and the prediction is made at the terminal node. Linear models are used at each step of the tree, and the tree is reduced to a set of rules (Kuhn and Johnson, 2013).

These algorithms were performed using the specific R packages and adjusted parameters following the methodology used by Cueff et al. (2021). Regression tree was carried out using “rpart” package (Therneau et al., 2022); the minimum number of observations that must exist in a node was 50 (minsplit = 50). Cubist algorithm was performed using “cubist” package (Kuhn et al., 2022), using the following parameters: number of rules by node (rules = 5). The predictions rules are adjusted to value of extrapolation of 5 (extrapolation = 5), number of the models used (committees = 1). Random forest was performed using “randomForest” package (Liaw and Wiener, 2018), using the following parameters: number of trees in the model (ntrees = 1,000), number of variables in each node (node size = 10), and number of variables randomly sampled in each division (mtry = 2).

The evaluation of the PTFs performances was carried out using different metrics: coefficient of determination (R²) (Eq.2), mean absolute error (MAE) (Eq.3), and root mean square error (RMSE) (Eq.4), BIAS (Eq.5), and the concordance correlation coefficient (CC) (Eq.6). All statistical analysis and modeling equations were performed using the version 4.1.3 of R (R Core Team, 2017). R 2 = ∑ i = 1 n Y ^ i − Y ¯ 2 ∑ i = 1 n Y i − Y ¯ 2 (2) M A E = ∑ i = 1 n Y ^ i − Y ¯ n (3) R M S E = 1 n ∑ i = 1 n Y ^ i − Y ¯ 2 (4) B I A S = ∑ i = 1 n Y ^ i − Y i n (5) where Ŷ _i is predicted value for observation i, Y _i is observed value for i, Ȳ is the average for soil variable selected, n is the number of observations. C C = 2   σ 12 µ 1 − µ 2 2 + σ 1 2 + σ 2 2 (6) where µ 1 and µ 2 are the mean for observed and predicted variables, respectively. σ₁ and σ₂ are the standard deviations for observed and predicted variable, and σ 12 is their covariance.

Soil texture distribution is presented in Figures 2, 3 following the USDA classification. Clay is the dominant class (38.07%), followed by Clay Loam (20.24%) and Sandy Clay Loam (9.97%), representing the different soil textures which could be present in arid and semi-arid regions (Figure 3). Figure 4 shows the results of hydraulic properties in different sampled regions. Soil hydraulic properties presented the same trend comparing the four regions, increasing from left to right side of the textural triangle (Figure 3). However, the lower values of w ₀ and w ₃₃₀ were in Doukkala region, 0.55 g g⁻¹ and 0.15 g g⁻¹ respectively (Figure 4). Though, the highest values were in Gharb-Loukous region with value of 0.76 gg⁻¹ and 0.30 gg⁻¹, respectively. Moulouya region presented the highest value of AWC of 0.11 gg⁻¹. In other words, the results show that both the w ₃₃₀ and w ₁₅₀₀₀ increase with decreasing of soil particle size.

In general, field data have a high variability when comparing the coefficient of variation (CV), and the classification of soil variability established by Oku et al. (2010), which classifies soil variability in the following classes: low variability (CV: 0%–15%), medium variability (CV: 16%–35%), high variability (CV>36%).

For the Doukkala region, iSDAsoil dataset presented low variability, with CV less than 15% in all soil properties studied. In Table 2, it is observed that the CV is higher than 35% for clay, silt, and SOM in the collected samples of the Doukkala region. Gharb-Loukous region presented the same trend regarding the low variability in soil information used from iSDAsoil dataset. However, the field data presented high variability with values of CV of 85.16%, 38.16%, and 40.01% in the sand, silt, SOM, respectively. In the Moulouya region, moderate variability was observed in the clay fraction in field data and iSDAsoil dataset. Regarding the field data, this region presented lower soil variability when compared to other regions. In the Tadla region, low variability in all soil properties in iSDAsoil data, contrariwise moderate to high variability was observed in the field data. In general, in all studied regions, iSDAsoil dataset tended to underestimate clay and silt content, and overestimate sand content and SOM (Table 2).

Regarding the performance of iSDAsoil data among studied regions, Gharb-Loukous region presented higher value of RMSE and MAE in all soil variables, principally for clay with RMSE 27.32%. According to Minai et al. (2021), clay prediction in digital soil mapping is based on the position on landscape using the information derived from digital elevation model. High clay content was predicted in flat depositional area and low clay contend in the alluvial area characterized by slopy and erosional landscape. We noted high clay content in the studied soil compared with Mediterranean soils, which could be explained by the position in landscape of the sampling points, as described in previous sections the studied regions constituted by plains and flat landscape. Field dataset with similar results were obtained by Eljebri et al. (2019) in the Doukkala region.

ISDAsoil data presents an ambitious project of mapping African soil with 30 m of resolution. This project was modelled using a machine learning framework with a compilation of 150,000 soil samples and different covariate layers representing information related to land cover, land use, land surface temperature, digital terrain parameters and bioclimatic variables.

Different criteria of evaluation of the performance of iSDAsoil data (predicted) and field data (observed) are presented in Table 3, comparing studied variables in the four regions. Silt showed lower values of RMSE and MAE. RMSE ranged from 11.22% to 27.32%, 9.34%–20.92%, 16.44%–40.31%, and 1.37% to 1.86, for clay, silt, sand, and SOM, respectively. Gharb-Loukous region presented highest values of RMSE and MAE in all soil variables, and the lower values are presented in the Moulouya region and Doukkala region. Regarding the SOM, Tadla region presented low values of RMSE (1.37%) and MAE (1.17%).

Variables used in the present study had a good accuracy and performance (Hengl et al., 2021). Using this iSDAsoil dataset allowed to predict hydraulic properties with different accuracies. However, iSDAsoil data tended to overestimate soil moisture content using PTFs, that could be explained by the overestimation of the soil properties used as input variables in the models. In addition, PTFs developed and selected presented the same trend regarding the performance in both data (iSDAsoil and field data).

The comparison between field and iSDAsoil data showed high values of RMSE and MAE in the Gharb-Loukous region. The clay soil in Gharb-Loukous region is a Vertisol formed on marl or on alluvial and colluvial material, as described by Moussadek et al. (2017), with clay content that could reach 80%. The solum formed by marl can reach several meters before finding a calcareous crust or a sandy horizon (Moussadek et al., 2017). However, the soils in Doukkala region are clay sandy soil as confirmed by ours results (Moussadek et al., 2017). Low accuracy to predict SOM could be explained by under-representation of the studied regions, and soil organic carbon dynamic in the soil, affected by the soil management system, crop systems, and climate change (rainfall, temperature, and mineralization) (Hengl et al., 2021).

Table 9 shows results of calibration and validation of new pedotransfer functions using different algorithms. Overall, calibration dataset exhibited slight improvement compared to validation dataset for all hydraulic properties. Nevertheless, both variables w ₃₃₀ and w ₁₅₀₀₀ presented the best performance compared to w ₀ in all algorithms.

The results revealed better performance of random forest compared to other algorithms for all properties (w ₀ , w ₃₃₀ , and w ₁₅₀₀₀ ). Prediction of w ₀ showed low performance compared to w ₃₃₀ and w ₁₅₀₀₀ with high value of RMSE and low value of R². Clay presented the most important variable for the prediction of w ₀ , w ₃₃₀ , and w ₁₅₀₀₀ , with variable importance of 35.42%, 33.27% and 37.40%, respectively, followed by sand, with variable importance of 33.16%, 32.18% and 31.48%, respectively (Figure 6). Variable importance order presented almost the same behavior for all predicted properties. Clay has the most important variable. Finally, SOM was less important variable.

Regarding regression tree algorithm, the results showed R² from 0.646 to 0.789 in the calibration dataset and ranged from 0.736 to 0.852 in the validation dataset: higher value of R² associated with w ₃₃₀ prediction, medium value with w ₁₅₀₀₀ prediction, and lower value with w ₀ prediction. Concerning Cubist algorithm, R² ranged from 0.626 to 0.786 in calibration dataset and from 0.782 to 832 in validation dataset. In Cubist algorithm, w ₃₃₀ , and w ₁₅₀₀₀ presented the same performance and low R² values were observed in w₀ prediction.

Several studies used models based on machine learning to predict water retention properties (Lamorski et al., 2008; Al Majou et al., 2018; Kotlar et al., 2020; Cueff et al., 2021). In general, random forest presented the best performance compared to other models based on machine learning or not. For example, Kotlar et al. (2020) compared random forest with linear regression using 192 samples collected in different land used (native forest, integrated crop livestock systems and grazing pastures. Similar results were obtained by Cueff et al. (2021) regarding the performance of random forest compared to other algorithms. In the present study, random forest showed the best performance compared with regression tree and Cubist algorithm.

Performance variability and low-quality prediction of published PTFs, could be explained by different reasons: 1) lack of harmonization of the protocol used in the determination of soil moisture at field capacity and wilting point, 2) difference in sampling, quantity and representativity of soil used, 3) difference in sample size and number of replicates to perform the analysis (Román Dobarco et al., 2019; Cueff et al., 2021), and 4) difference in selected soil properties (chemical, physical and mineralogical properties) used to develop these PTFs.

PTFs developed is this study can estimate soil hydraulic properties from soil particle size distribution, and soil organic matter, and to understand soil hydraulic pattern in Moroccan conditions. This is particularly the case in the current context of climate change impact and the increasing demand for agricultural production. Sustainable water management is still a challenge to reconcile the increasing demand for irrigation water and the optimization of the use of resources. Therefore, sustainable water management can be performed using crop water productivity models (e.g., AquaCrop), to define the net irrigation and water requirements of the crops. The estimation of soil hydraulic properties using regional PTFs could improve the accuracy and precision of these type of models, increase water use efficiency, avoid waste of water by evaporation or deep percolation. In addition, crop water productivity models request a large databases and soil information as inputs for to give good results. Soil properties estimated by PTFs could be an alternative to fill the information gap, and to run these models with low cost.

The present research allowed to develop PTFs for Moroccan soils, and to constitute the first database of soil hydraulic properties in the region. However, further studies are needed to improve and to upgrade the models developed in the present study to integrate other soil properties, pedological information and soil management systems, using a whole approach from the modeler to field user of soil information. In addition, there is a need to harmonize methodologies, and to develop an easily accessible database of soil information in the region. This database that could be integrated to regional or continental geographic information system, including databases of soils information, natural and climatic variables, which could be easily updated to respond to various requests from users at different levels (famers, agronomists, researchers, and decision makers).