1 Introduction

Soil degradation, low inputs of fertilizers, and suboptimal management of pests and diseases escalate the effect of water shortage—the most vital factor for agriculture in Morocco—on cereal production in Moroccan rainfed areas, threatening rural livelihoods and the country’s food security (Balaghi et al. 2013; IPCC 2014). Within the second leg of the national agricultural development plan in Morocco (Plan Maroc Vert 2008–2018) (MAPMDREF 2020a), a major research and development effort was made to improve wheat production in Moroccan rainfed areas (Singh et al. 2009). To continue this effort, appropriate tools are needed to provide farmers with appropriate pre-season and in-season advice for an efficient management of limiting resources, and narrow the gap to water limited potential yield in rainfed areas of Morocco (Shroyer et al. 1990; Jlibene 2009; Attia et al. 2021).

Crop simulation models are useful tools that provide insight in order to understand the interactions between crop genotype, crop management systems, and soil components and their dynamics (mainly water, carbon and nutrients) as affected by climatic conditions (Hoogenboom et al. 2004; Asseng et al. 2013; Yan et al. 2020). They can support recommendations for appropriate crop management practices (Mohanty et al. 2012; Ahmed et al. 2016). The Agricultural Production Systems Simulator (APSIM) (Keating et al. 2003) is among the most widely used models according to the literature. APSIM-wheat, the wheat version of the APSIM model, has allowed an increased understanding of wheat cropping system functioning and improved management practices for wheat crops, such as nitrogen management (De Silva et al. 2021; Berghuijs et al. 2021), phosphorus management (Wang et al. 2014), irrigation management (Chen et al. 2010a; Balwinder-Singh et al. 2011), adjusting sowing dates (Zhang et al. 2012; Hussain et al. 2018), and choosing cultivars and breeding (He et al. 2017; Kawakita et al. 2020). This model has been widely used to estimate potential yield and actual yield gap, as well as to assess and recommend optimal wheat management practices under limited soil and climate conditions, and to estimate wheat yield in various countries with a Mediterranean climate: Australia (Keating et al. 2003), Tunisia (Bahri et al. 2019), Spain, and Italy (Berghuijs et al. 2021).

Running crop models to simulate specific situations requires site- and crop-specific data on climate, soil, crop management, and cultivar characteristics to parameterize, calibrate, and evaluate the model under local conditions prior to actual use (Ahmed et al. 2016; Boote et al. 2016). Parametrization, calibration, and evaluation require large amounts of data that can be highly costly and time-consuming to acquire through experimentation (Makowski et al. 2006; Hsiao et al. 2009; Varella et al. 2010), especially for simulations that cover many different environments. Increasing awareness about the impact of soil health and fertility on plant growth has led researchers to complexify models to integrate more below-ground mechanisms and interactions, such as phosphorus or potassium availability to plants (Delve et al. 2009; Wang et al. 2014; Scanlan et al. 2015). This trend unavoidably results in a larger number of soil and plant parameters. Values for soil parameters can be difficult to acquire, which adds more uncertainty to model outputs (Sinclair and Seligman 1996). Ideally, parameter values for soil or plants should be measured locally, but most frequently soil parameters are estimated from the literature, inferred from proximal situations, and/or calibrated to minimize the error between observed and simulated model outputs on a limited number of trials. Consequently, parameter estimation is responsible for a non-negligible part of model uncertainty (Monod et al. 2006). In many modeling studies conducted over large geographic areas, existing meteorological or soil databases were considered for a direct acquisition of model input parameters (Rigolot et al. 2017; Ojeda et al. 2018; Elli et al. 2020; Liu et al. 2020). The databases’ accuracies, spatially or temporally (Seidel et al. 2018), may impact the simulation results depending on model sensitivity to the corresponding parameters. Little has been done so far to quantify the share of uncertainty on soil input sourced from spatial database within overall model uncertainty. The OCP group (Cherifian Office of Phosphates) in collaboration with the agronomic research institutions consortium and the Moroccan Ministry of Agriculture and Maritime Fisheries has developed the Moroccan soil fertility map, “Fertimap” (http://www.fertimap.ma/map.html). The Fertimap database was constructed from more than 32 000 soil samples collected from different regions of Morocco during the last decade and analyzed for different soil fertility parameters (Fig. 1). Interpolation with a kriging algorithm has been used to establish the full raster map from individual observation points. The Fertimap database could be used in research modelling as a resource for determining soil input parameters. However, some of the soil parameters recorded in Fertimap vary significantly over time, from one season to another. It is therefore necessary to perform for a first assessment of the impact of using the Fertimap database on model accuracy and to estimate model sensitivity to the provided soil parameters in the Fertimap database.

Fig. 1
figure 1

Fertimap advisor (http://www.fertimap.ma/map.html).

In the present study, we propose to assess the impact of sourcing soil parameters from spatialized soil database on model outputs and recommendations that can be driven from simulations. In this perspective, we have assessed the relevance of using the Fertimap database to provide inputs for APSIM-wheat when simulating wheat management scenarios in Moroccan rainfed areas. We proceeded through a three step-approach, estimating (i) the capacity of the APSIM-wheat model to accurately reproduce observed phenological stages, growth state variables, and final yield, but also soil water content (a state variable that is a determinant for plant growth in rainfed areas); and (ii) the effect of using Fertimap data instead of local- and time-specific soil analysis data on APSIM-wheat model accuracy.

2 Materials and methods

2.1 Study area

Data were collected throughout the crop cycle on 126 farmers’ wheat fields during three cropping seasons (2018/2019, 2019/2020, and 2020/2021) across the entire area for rainfed wheat production in Morocco (Fig. 2). The complete database was obtained by combining fields monitored within the framework of the “SoilPhorLife” project during 2019/2020 and 2020/2021 cropping seasons, and farmers’ fields monitored by OCP’s Al-Moutmir program during 2018/2019, 2019/2020, and 2020/2021 cropping seasons.

Fig. 2
figure 2

Study area map: locations of 126 farmers’ fields that were monitored during three crop seasons (i.e., 2018/2019, 2019/2020 and 2020/2021) and used to calibrate and evaluate the models. The colors of filled symbols represent the wheat cultivar. ● ▲ stand for fields whose data were used for model calibration, ■ stands for fields whose data were used for model evaluation. Fields were spread across the four main ISCRAL climate zones for wheat production (Gommes et al. 2009).

Farmers’ fields were evenly and randomly selected across four climate areas (Fig. 2), based on the Strategy for the Conservation and Restoration of Agricultural Land (ISCRAL) zoning developed by the FAO: favorable rainfed areas (i.e., precipitation > 400 mm), intermediate rainfed areas (i.e., 300–400 mm), unfavorable rainfed areas (i.e., precipitation < 300 mm), and mountain rainfed areas. Hence, the final field sample covered a diversity of climatic conditions, soil types, and most importantly different crop management strategies (i.e., cultivars, sowing dates and densities, fertilization, etc.).

Only five different winter wheat cultivars, considered the most widely used in Moroccan rainfed areas, were planted by farmers on these 126 fields: BANDERA, RADIA, FAIZA, ARREHANE, and ACHTAR. Farmers’ choices for wheat varieties were based, essentially, on local recommendations for the given climate area and the availability of seeds on the market.

2.2 Crop and soil monitoring in the field

Phenological stages, yield, above-ground biomass (AGB), leaf area index (LAI), and soil water content (SWC) were monitored on the farmers’ fields. Sampling intensity and observation protocols varied from one field to another. The monitored fields could be classified into three categories: (i) fields with a high measurement intensity, where crop (phenological stages, LAI, and AGB) and soil (SWC) measurements were taken three to four times in-season at critical wheat growth stages (i.e., emergence, tillering, elongation, anthesis, and maturity), in addition to final yield assessment; (ii) fields with medium measurement intensity where only final yield and phenological stages were recorded; and (iii) fields with a low level of measurement where only yield was measured (i.e., end-season measurements).

AGB was assessed through destructive sampling following the sampling design of Bazot et al. (2008): at each measurement date, five samples were collected randomly in each field. For row seeded fields, samples were one linear meter long (0.5 m over two consecutive rows). For fields with broadcast sowing, 0.5 m×0.5 m plots were delimited and sampled. Each sample was then dried at 70°C for 24 h in a well-ventilated oven (Bell and Fischer 1994; Mehrabi and Sepaskhah 2020) and weighed. Wheat yield was measured directly in the field at harvest (i.e., between May and June depending on the climatic zone). Five plots of 1 m×1 m were sampled randomly in each field, according to Bell and Fischer’s (1994) methodology. Phenological stages were recorded using Zadoks scores (Zadoks et al. 1974) and averaged over 10 to 20 individual plants selected across each field using a zig-zag design. LAI was measured using ceptometer (AccuPAR LP80, Decagon Devices Inc.), at 10 stations along a zigzag design in each field at each field visit. SWC was measured in the 0-10 cm top soil layer at every field visit using a capacitive mobile sensor (WET sensor device, Delta-T Devices Ltd). The sensor was calibrated before each crop season using generalized calibrations provided for the most common soil types (Delta-T 2018). In each field and at each visit, 10 to 12 SWC measurements were conducted at different points using a zig-zag design.

2.3 Model description

The APSIM v. 7.10 crop model was used in this study. APSIM-wheat simulates wheat growth and development on a daily time-step on an area basis by responding to the soil components and wheat management practices, and is driven by daily weather data (Asseng et al. 2013; Zhao et al. 2014; Zheng et al. 2015). Crop development was simulated by APSIM-wheat with 11 principal development stages calculated from the accumulation of thermal time, vernalization, and photoperiod sensitivities, and affected by water-nutrient stress (Chen et al. 2010b; Zhao et al. 2020). For the APSIM-wheat model, daily biomass accumulation was calculated based on radiation interception and radiation use efficiency (RUE) which is limited by soil water and nitrogen stress factors. The detailed development history of APSIM-wheat was reported in Holzworth et al. (2014).

2.4 Model parameterization: data sources and calibration procedure

2.4.1 Weather data

Daily weather data for each field and each of the three cropping seasons were obtained from different data sources. Daily precipitation and minimum and maximum temperatures were obtained from ground-based weather stations, either airport stations (available on the website www.tutiempo.net) or weather stations owned by the Regional Offices for Agricultural Development (ORMVA). For each field, weather data were inferred from the closest available ground-based station. Daily global radiation was extracted from NASA’s Prediction Of Worldwide Energy Resources (NASA’s POWER project) (https://power.larc.nasa.gov/data-access-viewer/) by selecting the corresponding grid cell for each field (grid resolution: 0.5°×0.5°, (Zhang et al. 2008; Sparks 2018)).

2.4.2 Soil parameters

Soil parameters for APSIM-wheat calibration and evaluation runs were derived from a combination of on-site measurements, the open access database, and references from the literature.

Parameters related to the water holding capacity were mainly derived from several spatialized soil databases. Soil water thresholds at lower limit of extraction by crops (LL), drained upper limit (DUL), and saturation (SAT)) were extracted from the “Global High-Resolution Soil Profile Database for Crop Modeling Applications” (https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/1PEEY0), which contains soil parameter values for the APSIM model maintained at a resolution of 0.083° (i.e., 10 km) (Han et al. 2019; Liu et al. 2020). Parameters related to soil evaporation dynamics were derived from the “CSIRO protocol for the development of APSoil parameter values for use in APSIM” (Dalgliesh et al. 2016) based on the values for southern Australia, due to the similarity in soil type between southern Australia and Morocco. The USDA-Soil curve number technique was used to parameterize soil run-off and the “CSIRO protocol for the development of APSoil parameter values for use in APSIM” (Dalgliesh et al. 2016) was used to estimate soil albedo and diffusivity parameters, depending on laboratory-determined soil textural classes and crop type (i.e., wheat). However, bulk density (BD) was measured in situ: soil samples were collected each famers’ field, at two depths (soil layer depths sampled depended on soil development in each field).

Soil parameters related to chemical and mineral characteristics of soil were all measured in situ in each field. Soil samples were collected from two layers (layer depths depending on soil depth in each field) and were analyzed before planting for organic matter content, pH, labile P, exchangeable K, electric conductance (EC), soil texture, NH4+, and NO3 using standard procedures.

Finally, parameters related to organic matter dynamics in soil (i.e., the fraction of labile (fbiom) and inert (finert) organic matter) were estimated using the soil type observed in the field and following the CSIRO protocol for the development of ApSoil (Dalgliesh et al. 2016). Furthermore, as APSIM-wheat integrates the simulation of phosphorus dynamics in the soil, phosphorus-related parameters (labile P, P fertilizers (banded and rock P), C/P ratio) were determined in each field.

2.4.3 Crop management parameters

The crop management inputs used in the APSIM-wheat simulations were collected on site based on the real crop practices and strategies recommended by experts and followed by the farmers. Crop management parameters included previous crop, tillage date, sowing factors (type, date, depth, row spacing, densities), cultivar, fertilization (application dates, rates and forms), weed control (weed species, density, and control application dates), and harvesting date.

2.4.4 Crop parameters

Crop parameters were estimated using field observations. Specific parameters related to phenology, biomass accumulation, and yield elaboration had to be calibrated for each of the five cultivars involved in the study. Previous studies on the global sensitivity analysis of APSIM-wheat key outputs (phenology, biomass, yield, and leaf area index) (Zhao et al. 2014; Casadebaig et al. 2014) were utilized to distinguish between the influential and non-influential crop parameters, and to identify the parameters that required calibration (Table 1). Calibration was conducted using wheat data from two cropping seasons (i.e., 2018/2019 and 2019/2020). Data collected across farmers’ fields in the four different agroecological areas and over two climatic years offered a contrasted dataset suitable for a robust calibration procedure. The comparison of climatic data for 2018/2019 and 2019/2020 cropping seasons indicated that, over the four agroecological areas, the 2019/2020 season was characterized by low precipitation and high temperatures compared to the 2018/2019 season. Cereal production recorded on farmers’ fields during 2019/2020 season significantly dropped by 39% compared to the previous crop season of 2018/2019 (MAPMDREF 2020b).

Table 1 Cultivar parameters fitted for APSIM-wheat

The main principle of Boote’s systematic approach (Boote 1999; Li et al. 2018) was followed to calibrate APSIM-wheat crop parameters for the five cultivars considered in this study, using three successive steps as follows: i) calibration of phenological and growth parameters, ii) calibration of leaf area parameters, and iii) calibration of grain yield parameters. Firstly, all crop coefficients were set using default values for the Australian winter wheat cultivar REVENUE since it is widely cultivated in Australian Mediterranean weather conditions, which are similar to the studied Moroccan rainfed areas. Secondly, the phenological parameters were adjusted for each of the five Moroccan cultivars using field observations and local weather data (see 2.4.a) to estimate the duration of phenological stages for each cultivar. Thirdly, a manual trial-and-error method was used to adjust other cultivar parameters (Table 1) for the five winter wheat cultivars, in order to minimize the error (estimated by the Root Mean Square Error, RMSE) between simulated and observed phenological stages, AGB, LAI, and grain yield parameters measured (or observed) during 2018/2019 and 2019/2020 cropping seasons. Fourthly, the consistency of the obtained parameters with the cultivar characteristics (precocity, maximum grain weight) as reported in the literature (Jlibene 2009) was checked and corrected if needed to avoid over-fitting due to the calibration process. Finally, P concentration thresholds (i.e., minimum and maximum P concentrations in plant dry matter) of plant organs at different stages were required as parameters in APSIM-wheat to activate the phosphorus module. The values of P concentration limits were derived from the database of Wang et al. (2014), which was created using experimental measurements (from Elliott et al. 1997; Rose et al. 2007; Bolland and Brennan 2008). For other parameters related to wheat crop (i.e., non-influential cultivar parameters, wheat-water parameters, wheat-nitrogen parameters, wheat-radiation parameters, etc.) the default values, reported in Zheng et al. (2015) and used automatically in the standard wheat-XML file, were conserved and used in this study (see Supplementary Materials).

2.5 Model evaluation

Model performances were evaluated by comparing plant (yield, phenological stages, LAI, AGB) and soil (SWC) variables as measured in the farmers’ fields during the 2020/2021 growing season (Fig. 2), and the corresponding simulations of APSIM-wheat.

Root mean square error (RMSE, Eq. 1), normalized root mean square error (NRMSE, Eq. 2), model efficiency (ME, Eq. 3), and coefficient of determination (R2) were calculated to evaluate APSIM-wheat’s capacity to simulate yield, phenological stages, LAI, AGB, and SWC adequately. The R2 was determined after modeling a linear regression using SPSS statistical software (SPSS Inc.).

$$ \mathrm{RMSE}=\sqrt{\frac{1}{n}{\sum}_{i=1}^n{\left[{x}_{s_i}-{x}_{m_i}\right]}^2} $$
(1)
$$ \mathrm{NRMSE}=\frac{\sqrt{\frac{1}{n}{\sum}_{i=1}^n\left[{x}_{s_i}-{x}_{m_i}\right]2}}{{\overline{x}}_m}\times 100 $$
(2)
$$ \mathrm{ME}=1-\frac{\sum_{i=1}^n{\left[{x}_{s_i}-{x}_{m_i}\right]}^2}{\sum_{i=1}^n{\left[{x}_{m_i}-{\overline{x}}_m\right]}^2} $$
(3)

where \( {x}_{m_i} \) are the measured values, \( {x}_{s_i} \) are the simulated values, xm is the mean of the observed values, and n is the number of observations.

2.6 Assessment of the impact of Fertimap database uncertainty on APSIM-wheat model output

Fertimap is a spatial soil fertility database of Moroccan agricultural areas developed from 32 000 soil samples collected from different regions of Morocco during the last decade and analyzed for different soil fertility parameters. Fertimap advisor gathers four soil parameters for each grid cell, in addition to soil texture: organic matter content, labile phosphorus content, potassium content, and pH (averaged over the whole soil profile). As K availability and its effect on plants are not represented in APSIM-wheat’s model, the influence of using Fertimap data to determine soil parameters could be assessed on the following model parameters: organic carbon, labile P, and pH.

First, input values for this three parameters when using one or another source of information were compared. Organic carbon, labile P, and pH values from local- and time-specific soil analyses (averaged over soil layers) conducted in 100 farmers’ fields were compared to their counterparts extracted from the Fertimap database. RMSE between the two values obtained from these two sources were calculated over the 100 farmers’ fields for each of the three parameters.

Second, a simple local sensitivity analysis was conducted to assess the effect of varying these three soil parameters (i.e., organic carbon, pH, and phosphorus) on model output for two farmers’ fields in contrasting climate and management conditions. The sensitivity analysis was focused on one single model output: crop yield. Indeed, crop yield was the variable of interest in this study and the three soil parameters were expected to influence simulated crop growth variables and eventually yield as reported in Attia et al., (2021). Local sensitivity analysis, which consists in varying the value of a few parameters while keeping the others constant, has several defects including the inability to detect interactions between parameters (Saltelli et al. 2008). To overcome this drawback, we repeated the sensitivity analysis for two contrasted situations: one field located in the favorable rainfed area with fertilizer application, and one field in an unfavorable rainfed area without any fertilizer application. We intended to explore the effect of weather and wheat management conditions in the studied farmers’ fields on the results of this sensitivity analysis with this approach.

For each of the two fields, ranges of uncertainties for OC, pH, and labile P parameters were determined, using the Fertimap database, as the limits (minimum and maximum) of possible parameter values in the climatic zone represented by each if the two fields (i.e., favorable areas for field (a) and unfavorable areas for field (b)). Ranges of parameter variations were therefore defined as follows: for field (a) OC [0.5 to 3.5%], labile P [5 to 50 ppm], and pH [5.5 to 8.5], and for field (b): OC [0.2 to 1.6%], labile P [5 to 50 ppm], and pH [6.5 to 8.8]. After checking the quasi-linearity of the variation of simulated yield when varying each input parameter (see Supplementary Materials), and to better detect the parameters interaction effect on APSIM-wheat simulated yield, a complete randomized experimental plan was built by varying the three parameters within their respective possible range at three equidistant values. All other input parameters were fixed according to the field characteristics. Simulations were run for the 27 possible combinations of values for the three parameters for each field (soil organic carbon, labile P, and pH) and simulated yield results were analyzed using a one-way analysis of variance (ANOVA) to assess the effect of each parameter on APSIM-wheat model output variability.

Finally, the distributions and medians of simulated yield and AGB when using Fertimap data were compared to the simulated results when using time and local-specific soil analysis data to parameterize the model for 100 farmers’ fields monitored in this study. A Mann-Whitney U non-parametric variance test (Mann and Whitney 1947) was performed to evaluate the difference between model outputs when using one or another source of data. All the statistical analyses were executed using the SPSS statistical software (SPSS Inc.).

3 Results and discussion

3.1 Model calibration: effects of rainfed climate conditions

A good agreement was shown between observed and simulated values used to calibrate the APSIM-wheat model and to estimate the cultivars’ genetic coefficients (Table 2). The calibration of APSIM-wheat provided an accurate prediction of phenological stages for the five cultivars with an overall RMSE of 4.4 stages (out of 100 stages of Zadoks growth scale), and NRMSEs ranging between 5 and 16% depending on the cultivar. Similarly, the overall error on yield prediction for the five cultivars was within 0.4 t.ha–1, with variation of the quality of prediction between cultivars (NRMSEs between 13% and 23%). Measured AGB values for the five cultivars were comparable to APSIM-wheat simulated values with an overall RMSE of 0.8 t.ha–1. In general, calibration results indicated a reasonable estimation of phenological development, biomass, and yield-related genetic coefficients for the five studied wheat cultivars.

Table 2 APSIM-wheat model calibration results: goodness-of-fit statistics; using root-mean-square error (RMSE) and normalized root-mean-square error (NRMSE) indices, after a comparison between APSIM-wheat simulated and observed Zadok’s phenological stages, yield (t.ha-1), above-ground biomass (AGB, t.ha-1) and leaf area index (LAI, m2.m-2).

The LAI predictions were of lesser quality, despite the model calibration processes. NRMSEs on LAI exceeded 35% for all cultivars except for BANDERA (NRMSE = 19%). Low quality in LAI prediction is a common feature of modelling studies (DeJonge et al. 2012; Falconnier et al. 2019; Saddique et al. 2019; Attia et al. 2021). A first possible explanation for this poor prediction of LAI lies to low accuracy in LAI measurement. LAI was not measured directly from destructive sampling but assessed from light interception measured in the field with a ceptometer, which could result in overestimation of LAI as light and dry leaves also reduce light transmission through the canopy as measured by the ceptometer. However, variations in the quality of LAI prediction between cultivars, fields, and year within the study dataset indicate that poor LAI simulations could rather from limitations in APSIM-wheat algorithms for Moroccan conditions. In APSIM-wheat, daily LAI increase is calculated as a function of the number of nodes; afterwards, several stress factors (water stress, N stress, and P stress) apply to potential LAI increase, as well as a carbon stress factor that applies when the amount of leaf dry matter produced on the preceding day is not enough to support the LAI increase of the day. The daily number of nodes depends on thermal time accumulation and on the achievement of the phenological stages, which were properly predicted according to model calibration results. Daily leaf dry matter depends on daily biomass production and partitioning equations that monitor the allocation of daily biomass to the different organs (Zheng et al. 2015; Ahmed et al. 2016). While daily biomass seems to be properly calculated according to calibration results, it is more likely that model errors originate in the simulation of dry matter partitioning or of the stresses reducing leaf area expansion or inducing early senescence. Better quality of simulations was obtained for BANDERA cultivar that was mainly used in favorable areas where water stress is less frequent. This finding supports a possible inaccuracy in the calculation of water stress effects on the LAI dynamic in the APSIM-wheat model algorithm (Mohanty et al. 2012). Alternatively, poor LAI simulation could result from an insufficient description of the tillering and tiller senescence processes. The APSIM-wheat algorithm describes the wheat plant as uniclum. Tillering and tiller senescence are encompassed in the coefficients allowing the calculation of daily LAI from the number of nodes on the main stem (Wang et al. 2003). Moreover, as in reality, tillering highly depends on plant density, which is also not accounted for in the model algorithm (Saddique et al. 2019). However, in Morocco, plant density can be very heterogeneous and low, especially in rainfed conditions, compared to plant density in mechanized wheat systems around the world. This could explain the limited capacity of the model to finely reproduce LAI dynamics in Moroccan wheat fields.

3.2 Model evaluation

Observed phenological stages, AGB, yield, LAI, and SWC from 27 farmers’ fields during the 2020/2021 cropping season were compared to corresponding variables simulated using the calibrated crop parameters (Fig. 3).

Fig. 3
figure 3

APSIM-wheat validation results for the five wheat cultivars used in the Moroccan rainfed areas. Comparison between simulated and observed a Zadok’s phenological stages, b yield (t.ha-1), c above-ground biomass (AGB, t.ha-1), d leaf area index (LAI, m2.m-2), and e soil water content (SWC, m3.m-3). Red line represents (x=y) equation. R2, root-main-square error (RMSE), normalized root-main-square error (NRMSE) and model efficiency (ME) are the statistical indices used to evaluate APSIM-wheat model accuracy.

When compared to observed data independent from those used for calibration, APSIM-wheat accurately predicted phenological stages for each of the five cultivars (NRMSEs lower than 15% and ME exceeding 0.95). The APSIM-wheat model also reasonably reproduced observed yield for the five cultivars with overall NRMSEs equal to 14%. Also, the model acceptably simulated AGB with an overall RMSE equal to 1.2 t.ha–1 and ME equal to 0.91. AGB was better simulated by the APSIM-wheat model for BANDERA, FAIZA, and RADIA cultivars, that were used in favorable and intermediate rainfed areas. The agreement between observed and simulated AGB was better during vegetative stages (before Z70) than during grain filling (Z70 to Z90). These two findings were consistent with potential possible model error originating in the calculation of water stress factors and their impact on daily leaf expansion and dry matter production in the APSIM-wheat model, as suggested by calibration results.

APSIM-wheat’s predictive quality for yield, AGB, and phenology as evaluated in this study was consistent with the literature on model evaluations for wheat production in rainfed areas. Calculated RMSEs and yield were in the same range or smaller than in most other studies using the APSIM-wheat model (Makowski et al. 2006; Hussain et al. 2018; Zhao et al. 2020). Ahmed et al. (2016) reported higher model quality scores in simulated yield and AGB and a similar performance in predicting phenological stages during the calibration and evaluation of APSIM-wheat under favorable rainfed conditions in Pakistan. However, in the present study, the calibration and evaluation process was carried out in 126 farmers’ fields with various soil, management practices, and climatic conditions which offset slightly higher errors. This diversity of conditions in the calibration and evaluation datasets warrants a robust model evaluation for wheat simulations in Moroccan rainfed areas. In addition, the quality of model prediction for phenology, AGB, and yield that was found in this study matched that of other crop models (WOFOST and CropSyst) that have been used previously in Morocco (Bregaglio et al. 2015). This study therefore enlarges the range of models suitable for simulating wheat in Morocco and provides an opportunity for simulating more complex scenarios such as P fertilization scenarios that are only possible with APSIM-wheat.

As in the calibration step, APSIM-wheat was found to have a limited capacity to reproduce LAI observations recorded in the evaluation dataset. APSIM-wheat model error on LAI prediction was acceptable only for the BANDERA cultivar (NRMSE = 26%) and to a lesser extent FAIZA (NRMSE = 36%) that were cultivated in favorable rainfed areas. In contrast, for other cultivars that were cultivated in intermediate and unfavorable areas, LAI was underestimated by the model, which seemed to overestimate the stress effect intensity. Nevertheless, acceptable to high coefficients of determination (R2 > .80) were observed between the simulated and observed LAI values for ACHTAR cultivar in unfavorable rainfed areas. Similarly, Ahmed et al. (2016) and Seyoum et al. (2018) found a high agreement (NRMSE < 9%) between observed LAI and LAI as simulated by APSIM-wheat only in areas with cumulative rainfall exceeding 600 mm in the crop season in Pakistan and Ethiopia, respectively. Additional improvements in the models’ response to drought could be achieved by improving the simulation of soil water state and/or adapting the calculation of the water stress coefficient in the model.

APSIM-wheat showed poor accuracy in simulating soil water content in the top soil layer (SWC) with RMSE = 0.08 m3.m–3 and tended to overestimate SWC. The model even seemed to have a decreased capacity to simulate SWC in drier environments: in favorable rainfed areas (> 400 mm) NRMSE on SWC ranged between 22% and 24% while NRMSEs exceeded 40% in most cases in intermediate (300–400 mm) and in unfavorable (< 300 mm) areas. These findings were consistent with results reported in Mohanty et al. (2012) which showed that the prediction of SWC by APSIM-wheat was relatively accurate (RMSEs in the 0–15 cm soil layer were below 0.04 m3.m–3) in non-water stress conditions. In contrast, a weak prediction of SWC in the soil surface layer by APSIM-wheat was reported in Australia, in dry conditions with only 160 mm in-season rainfall (Zeleke 2020). In Zeleke’s work, APSIM-wheat was also found to overestimate SWC in the surface layer (i.e., 0–15 cm) in unirrigated and irrigated plots.

Low-quality SWC estimation in the top soil layer combined with a limited capacity to estimate LAI in the driest region was both consistent with possible errors originating in the model’s water module. Overestimation of the top soil water content could result from an underestimation of soil evaporation or infiltration in the heavy cracking soil that is common in Morocco and known as “Tirs” (Moussadek et al. 2017). Conversely, some authors have previously suggested that low-quality simulation of the soil water content could result from inaccurate calculation of the canopy cover and water stress factors, hence of the transpiration flux (Paredes et al. 2014; Montoya et al. 2016; Mehrabi and Sepaskhah 2020), which was also supported by the medium to low quality of LAI prediction found in the present work. Deviations between observed and measured SWC could also result from errors on the soil parameters that we inferred from the “Global High-Resolution Soil Profile Database” (Han et al. 2019), in particular DUL and LL to which the model has a strong sensitivity. Uncertainties on the climatic input in particular on rainfall that may differ between the field location and the weather station could also participate in the error on SWC.

3.3 Impact of Fertimap database uncertainties and parameters on APSIM-wheat output

As reported in previous studies, the uncertainty in a crop model comes from three sources: model structure, model parameters, and input data (White et al. 2011; Challinor et al. 2013; Seidel et al. 2018). The use of existing large databases affected by temporal and spatial variability (e.g., Fertimap) as model input parameters constitutes a source of uncertainty (Xiong et al. 2020). In this study, we conducted a first assessment of the impact of integrating the digitalized Moroccan soil fertility database (Fertimap) as a source of APSIM-wheat input parameters on the model output (i.e., AGB and yield).

The values of soil pH, OC, and labile P inputs extracted from Fertimap were first compared to those obtained from site specific measurements. The comparison between the two data sources of the distributions of values obtained for these parameters for the 126 framers field in the study showed a small difference between the medians for the three parameters (OC medians = 1.2% and 1.1%, labile P medians = 28 ppm and 21 ppm, and pH medians = 7.9 and 7.8, for Fertimap and soil specific analysis respectively) (Fig. 4). RMSEs between inputs from the two sources were equal to 0.63% for OC, 10 ppm for labile P, and 0.67 for pH. Therefore, uncertainties on model parameters imputable to data source were small compared to the variation observed between fields and remained within the ranges of variation used for the local sensitivity analysis presented afterwards. Differences between parameter values from two sources could originate from several causes. First, spatialized data available in Fertimap were obtained from interpolation of localized measurements; differences could therefore result from the fact that the location of farmers’field used in this study did not exactly coincide with one of the location where measurements were taken to elaborate Fertimap. Second, the temporality of parameter value assessment differs between the two sources of inputs: soil chemical characteristics evolve under the effect of soil management and fertilization practices applied to a given field. While values recorded in the Fertimap database were measured a few years ahead, parameter values measured in situ were assessed a few days before planting date for each year. Finally, the variation and complexity of soil organo-chimical processes (e.g., P sorption and desorption, P immobilization and mineralization, organic matter decomposition etc.) overlap with human management practices in complicated mechanisms and processes that contribute to in those variations. Despite these differences, parameters contained in Fertimap database appeared to remain rather close proxies of actual soil characteristics when measured at the beginning of the copping seasons, revealing that i) soil chemical properties would vary little over short distances within the wheat production area in Morocco, and ii) soil fertility characteristics vary slowly in time, maybe due to low organic matter and reduced fertilization inputs in Moroccan wheat production area.

Fig. 4
figure 4

Comparison of Fertimap data and specific soil analyses data values (left) and their distributions (right) for three soil parameter values: a organic carbon, b labile P, and c pH. Black line in left graphs represents (x=y) equation.

The effect of the variation of each of the three parameters on wheat grain yield as simulated by APSIM-wheat was then assessed for field type (a) (> 600 mm in-season rainfall and added fertilizers) and field type (b) (< 200 mm in-season, no applied fertilizers), with a local sensitivity analysis for the two fields. For field (a), the ANOVA revealed a significant effect of OC [p < .001] on simulated grain yield with an increase of 274 kg.ha-1 between the use of the minimum (0.2%) and maximum (1.6%) observed OC values, while there was no significant impact of labile P [p = .076], and pH variations had no effect on simulated yield [p = 1]. This finding was consistent with real-world experimental evidence showing that wheat is a crop that has low demand for phosphorus (Olsen and Sommers 1982; Tsadila et al. 2012) and is tolerant to pH variation (Franc de Ferrière 1933; Froese et al. 2015). In the case of field (a), phosphorus demand may be satisfied with applied fertilizers (130 kg.ha–1 of N and 60 kg.ha–1 of P fertilizers added during the crop season). Conversely, the increase of simulated yield when OC was varied most likely only results from an improved N status of the plant as simulated by the model, as the OC parameter is not involved in the calculation of the soil water balance in APSIM. Since water stress levels were moderate in field (a), it is likely that N was the limiting resource in the corresponding simulations; an increase in the OC soil parameter most likely allows an increase in available N along the simulated crop cycle and, consequently, a slight increase of yield.

For field (b), in conditions with lower precipitation (< 200 mm in-season) and no fertilization, labile P [p < .0001] and OC [p = .001] significantly influenced the variability of simulated yield, whereas pH did not [p = 1]. Yet, the sensitivity of simulated yield to OC parameter variations was much smaller than in field (a), as simulated yield increased only by 1.6 kg.ha–1 when OC increased from minimum to maximum observed values. For labile P, yield increased by 96 kg.ha–1 when labile P increased from minimum to maximum observed values. This result consolidates the interpretation given above: when water stress increases and water becomes the main limiting resource, the effect of higher OC in soil resulting in higher N availability is less, as N uptake is limited in the model by water stress. Actually, in APSIM’s Soil-N component, soil water and soil temperature are factors that impact the decomposition of organic matter, nitrification, and urea hydrolysis processes. This interaction between water and N resources on plant N uptake and growth has been largely observed in reality and reported in the literature. Cramer et al. (2009) and Nawaz et al. (2012) reported that early or late drought in the cropping season decreased nutrient uptake, which would be due to both reduced mass flow in the soil and sap flow and vascular conductivity in the plant as evidenced in controlled conditions. The higher sensitivity of model output to variation of labile P parameters in field (b) compared to field (a) may result from the absence of applied fertilizer, and, consequently, the plant relying only on soil P to satisfy its P demand, but also from reduced labile P input from mineralization which is also dependent on soil humidity in the model (Delve et al. 2009). Finally, the low sensitivity of simulated yield to the pH parameter in field (a) as well as in field (b) was consistent with the fact that, in real-world contexts, wheat is tolerant to alkaline and acid soil and that fertilizing practices as entered in the model may output the effect of soil pH on crop yield.

The results outlined by the sensitivity analysis were confirmed by the comparison of model outputs when using one or another data source to infer soil parameters in APSIM-wheat. No significant difference was found between the median of simulated outputs when using Fertimap or local data as soil inputs with APSIM-wheat, neither for yield (median: 0.21 t.ha–1; p = 0.936 ) nor AGB (median: 1.42 t.h–1, p = 0.939). Distributions of model outputs over the 100 farmers’ fields and the three cropping seasons were similar whatever data source was used, for the three parameters (Fig. 5). Whether parameter values were extracted from Fertimap or local measurement, yield and AGB predictions varied within similar ranges over the entire study zone and the 3 years: yield varied within [6509–149 kg.ha–1] and [6509–167 kg.ha–1] ranges with site-specific measurements and Fertimap data respectively, and AGB varied within [8.23–0.09 t.ha–1] and [8.23–0.09 t.ha–1] with site-specific measurements and Fertimap data respectively. Therefore, the source of data used to parameterize OC, pH, and labile P did not seem to impact APISM-wheat model outputs overall in rainfed Moroccan conditions.

Fig. 5
figure 5

Overall comparison of APSIM-wheat simulated output values (left) and their distributions (right); a yield (kg.ha-1) and b above-ground biomass (AGB, t.ha-1), when using two soil data types as model soil input: specific soil analyses and Fertimap database. Black line in left graphs represents (x=y) equation.

4 Conclusion

This study delivered crucial results that open avenue to use APSIM for pre-season and in-season advising of farmers in Morocco: (i) it delivered verified parameters for the five main wheat cultivars grown in the rainfed area of Morocco. (ii) It demonstrated a satisfactory capacity of the APSIM model to predict wheat yield for different climate and fertilization management within the same area. Although APSIM-wheat was found to fail to predict LAI with an NRMSE lower than 26% for the five cultivars, several hypotheses to explain this lack of accuracy have been identified. Indeed, the low predictive quality of the model for LAI was partly explained by errors in the prediction of SWC, especially in the driest conditions, where SWC was systematically overestimated by the model. The present study allowed the identification of possible ways to adapt APSIM-wheat to Moroccan conditions and improve the quality of LAI prediction, such as the adaptation of soil evaporation algorithms for cracking soils. (iii) It validated the Fertimap database as a data source for parameterizing soil in APSIM, especially for labile P, OC, and pH parameters. We showed that despite the possible variation in time and space of the variables it compiles, Fertimap can be a valuable source of data to determine soil parameters for crop models and particularly to simulate P availability in soil. Using verified APSIM-wheat model parameterized with Fertimap for the five cultivars should allow to simulate appropriate management scenarios spatialized across the wheat production area in Morocco and to advise farmers to adapt their planting and fertilization strategies at sowing and in season, depending on climate experienced and/or expected during the on-going season.