Soil sampling

Two agricultural fields, labeled 38C and 47, were monitored. These two fields are in the region of La Pocatière, Québec, Canada (47°21′21′′N, 70°01′45′′W and 47°20′37′′N, 70°00′52′′W), and have a fine texture, noted G1 (38C: Kamouraska clay loam) and a medium texture, noted G2 (47: Saint André sandy loam). In this study, field 38C (2.1 ha) was seeded to spring wheat in a binary rotation with soybeans under a conventional 0.15 m deep fall tillage and spring harrowing. Yield levels were 2.4–3.2 Mg·ha⁻¹ for wheat and 1.5–2.2 Mg·ha⁻¹ for soybeans. This 38C field received only mineral fertilization of about 75 kg diammonium phosphate (DAP) and a 100 kg mixture of 80% ammonium nitrate and 20% calcium carbonate (CAN) at seeding (i.e., 40.5 kg N·ha⁻¹ and 15 kg P·ha⁻¹) and 48.6–59.4 kg N·ha⁻¹ as CAN at tillering. Field 47 (1.9 ha), on the other hand, was dedicated to research to test cultivars of different crops. It received crops of quinoa, barley, corn, and wheat from 2009 to 2012 under conventional 0.20 m deep fall tillage. During this period, field 47 received 15 m³ of sheep manure in the fall of 2008 and 23 m³ of swine manure in postemergence of corn (2011). As this field is highly concentrated in P and K, mineral fertilization was limited to a single application of nitrogen in postemergence. No liming was done during the duration of the sampling.

For both fields, soil samples were collected once a week for 33 weeks a year, from mid-April to early December during four years from 2009 to 2012. Since it is neither useful nor practical to collect soils while they are covered with snow, a 19 wk period from early December to mid-April was excluded from the weekly sampling. For each of the two fields, three equidistant soil locations, noted a, b, and c, were sampled along the middle of the field. Each location was clearly identified with flags that remained installed during the four years. Each flag identified the center of a five-meter radius circle. Consequently, 792 samples were collected (2 sites × 3 locations (a, b, and c) × 33 weeks × 4 years). For each location a, b, and c, five sub-samples were taken randomly by a 7 mm diameter Pro-Sonde probe (Khiari et al. 2014) at a 20 cm depth within the sampling circle to form a composite sample of 30 to 40 grams. Precautions were taken to avoid trampling the soil during soil sampling by using wooden planks suspended above the sampling area.

Soil sample testing

Soil samples were tested for pH_W and pH_SMP according to the electrometric method of the CEAEQ (2003a). Soil texture was determined by the hydrometer method (Day 1965). The OM content was determined using the loss on ignition method (CEAEQ 2003b). The soil samples were also tested for their contents in P, K, Ca, Mg, B, Cu, Zn, Fe, Mn and Al extractable at the Mehlich-3 common extractive (Mehlich 1984; CEAEQ 2014). These 10 elements were measured by plasma emission spectrophotometry (Varian model 725-ES ICP-OES, Torch Type: Radial; Australia Pty Ltd.). The cation exchange capacity was estimated by the following equation (CRAAQ 2010):

(1)

Soil phosphorus saturation index (SPS) was calculated by (Khiari et al. 2000):

(2)

Most series were complete, but for some properties, less than 1% of the observations were missing. To have complete time series, these missing values were interpolated by the na.interpolate function of the imputeTS package of the R software (Martitz 2021). The temporal variability of soil properties was expressed as the standard deviation ratio (SDR) and compared with the acceptable quality level for reproducibility (AQLR) of soil sample analyses. The AQLR is defined by the CEAEQ (2015) for analytical quality control of the accredited laboratories in Québec. When a soil test was not on the CEAQ (2015) list, we used the AQLR of NAPT (2001). Of the 90-temporal series (2 sites × 3 localities × 15 fertility indicators), those with an SDR < AQLR were not considered problematic, that is, whatever the sampling time, the sample is considered representative. Time series with SDR > AQLR were decomposed using loess STL of the stlplus package from R software (Hafen 2016). Loess (STL) is an algorithm that has been developed to separate a time series into three components: seasonality, trend, and remainder. Of these three components, only seasonality was used to analyze periods of regularity and irregularity in soil testing reflecting the distribution of observations depending on the sampling period. As for the interannual trend, it only serves to detect whether soil tests may change within a four-year sampling cycle, after correcting for seasonality over the period 2009 to 2012. The third remaining component is a noise signal obtained after eliminating the two effects of seasonality and inter-annual trend and represents the random fluctuations. For the series with an SDR > AQLR, we did an autocorrelation analysis using the acf function of the R forecast package (Hyndman et al. 2021). The auto-correlogram shows the similarity between observations with different time lags. For the seasonal components with SDR > AQLR, we defined a lower limit of validity, noted LL_validity (eq. 3) and an upper limit of validity, noted UL_validity (eq. 4). The AQRL was applied to a central value most representative of the soil tests in Québec, that is, the average calculated over the five weeks (W) of spring (S), noted (eq. 3 and Fig. 1).

(3)

(4)

Fig. 1.

Example of theoretical variation of an indicator bounded by the validity interval. is the average of first five weeks (W) of spring. Acceptable quality level for reproducibility (AQLR) of soil sample analyses.

These two limits allow the inclusion or exclusion of sampling periods. The average is taken over five weeks since all the fertilization grids and the agronomic and environmental critical thresholds developed in Québec have been designed considering a spring sampling from the third week of April to the end of May.

Fertility indicator values within the inclusion periods are considered stable and representative of spring sampling, those outside the inclusion periods are not representative of spring sampling and should be excluded from the sampling window. Thus, the exclusion period defined in Fig. 1 is a continuous range of weeks not recommended for soil sampling to limit the effect of summer on soil test values. Before using the results from the seasonality of soil tests, it is first necessary to make sure that these effects are present and can be considered. For one or two successive exceedance values, the phenomenon is considered sporadic and does not require the sampling windows to be modified. When seasonal or calendar effects cannot be identified in a time series, the time series is assumed to be deseasonalized (no seasonal effect) even if its SDR exceeds the AQLR. Whether sporadic or nonexistent (SDR < AQLR), the series is free of a sampling calendar effect and no further analysis is required since samples can be taken any time from spring to fall.

Descriptive statistics of indicators

The summary results of SDR (Table 1) show that the overall seasonal variability of the 15 fertility indicators is large and ranges from 1.58% for pH_SMP to 54.97% for Cu_M3. As observed by Díaz-Ravifia et al. (1993), seasonal variations are more important for some. The five indicators pH_W, pH_SMP, P_M3, SPS, and Cu_M3 showed much higher SDRs for the G1 soil than for the G2 soil. Clay appears to contribute to the high random and temporal variability of acidity, phosphorus, and some micronutrients, mainly Cu_M3. Turpault et al. (2008) explained this by the large specific surface areas of clay minerals to react with other soil elements to cause substantial changes over time. Also, the seasonal variability of OM and major nutrients (P, K, Ca, Mg) with SDRs of 9% to 24% is less than that of micronutrients with SDRs of 14% to 55%. Moreover, the two systems NAPT (2001) and CEAEQ (2015) allow more variation on micronutrients (10%–20%) than the other soil tests (5%–15%).

Table 1.

Standard deviation ratio (SDR) of 15 fertility indicators for two soil types (G1 is fine texture; G2 is medium texture) and three locations (a, b, and c), compared with the acceptable quality level for reproducibility (AQLR) of two soil testing laboratory control programs for weekly sampling (from mid-April to early December) during four consecutive years.

Table 1 also shows that SDR values vary very little for a given field for all 15 fertility indicators. Even for Cu_M3, the large SDR values are very close to one other with values of 55%, 55%, and 48%. Therefore, for these two fields, spatial heterogeneity does not seem to interfere with the temporal heterogeneity. For the two soil types, G1 and G2, only indicators pH_SMP, Ca_M3, Mg_M3, Al_M3, CEC_predicted, and B_M3 are within the limits of the accepted variation criteria of CEAEQ (2015). The other indicators, pH_W, P_M3, SPS, Mn_M3, Cu_M3, Zn_M3, Fe_M3, OM, and K_M3, vary more and are above the CEAEQ (2015) limit for G1 and G2. For the soil acidity status, the SDR range of pH_W is 2.4%–4.2%, a variation of 0.2–0.3 pH units, significantly lower than 0.8 units of the temporal variation reported by Colins et al. (1970). For phosphorus status, the SDR range is 13%–27% for P_M3 and 12%–33% for SPS, which upper limits are comparable to the 23% temporal SDR value for bioavailable phosphorus reported by Nyborg et al. (1992). The indicator Cu_M3 has the highest SDR in the range of 23%–55%, with an upper limit is comparable to the temporal SDR of 43.8% for bioavailable Cu reported by Cameron et al. (1994). To define the optimal sampling range, only the latter indicators are considered. Consequently, nine of the 15 fertility indicators were analyzed with time series.

Seasonal variability of indicators

Figure 2 presents the time series of pH_W, for soil G1 and location a, and the additive components of seasonal, trend, and remainder. This seasonal component is representative of the three locations of field G1 (Table 2). During the four years, there is clear seasonal variation with four relatively high pH_W peaks, all corresponding to the same period in early spring (Fig. 2b). In mid-summer, between July and August, there is a pH_W trough explained by increased biological activity. These seasonal variations in pH_W are related to temperature and humidity, which are periodic in nature. The results of this analysis show that the average pHw does not show a clear increasing or decreasing inter-annual trend during the four years of study (Fig. 2c).

Fig. 2.

Additive decomposition of pH_W time series for soil G1 and replicate a. [Colour online.]

Table 2.

Averages of fertility indicators, their types (Continuous, Sporadic, or Inexistent), their intervals of exceedance of validity limits, and their autocorrelation periods in the fine textured (G1) and medium textured (G2) soils.

A total of 54 time series were analyzed, that is, the nine indicators selected from the previous section × two soil types (G1 and G2) × three locations (a, b, and c). These nine indicators were first assigned to one of three utility groups: (i) acidity and liming management indicator (pH); (ii) soil nutrient indicators other than phosphorus (K_M3, MO, Mn_M3, Zn_M3, Cu_M3, and Fe_M3); and (iii) agro-environmental phosphorus management indicators P_M3 and SPS.

For each utility group, only one indicator is presented to illustrate its seasonal variation. The three representative variables are pH_W, K_M3, SPS and will be discussed in the following sections. For the other six indicators (P_M3, OM, Mn_M3, Zn_M3, Cu_M3, and Fe_M3), the resulting observations are summarized in Table 2. Analyses of the weekly data show a gradual change in the indicator values, but with occasional sharp fluctuations, as shown in Fig. 3 for pH_W, in Fig. 4 for K_M3, and in Fig. 5 for SPS. The seasonal component is the baseline on which a validity interval for each indicator is applied, as explained in Fig. 1. The trend and remainder series are not required for this study.

Fig. 3.

Time series showing a seasonal pattern of pH_W variation (upper part) over a 33 wk cycle per year, followed by the autocorrelation function (lower part) of pH_W (correlation versus lag). [Colour online.]

Fig. 4.

Time series showing a seasonal pattern of potassium Mehlich-3 (K_M3) variation (mg·kg⁻¹) over a 33-week cycle per year. [Colour online.]

Fig. 5.

Time series showing a seasonal pattern of soil phosphorus saturation (SPS) variation (a; upper figure) over a 33 wk cycle per year, followed by the autocorrelation function (b; lower figure) of SPS (correlation versus lag). [Colour online.]

Variability of acidity diagnostic and liming management indicators

When considering the NAPT (2001) criteria of AQLR, both indicators (pH_W and pH_SMP) would have unacceptable variations for both soil types G1 and G2 (Table 1). The NAPT (2001) system is significantly more stringent for the diagnosis of acidity, with pH tolerances of 0.8%–1.2%, compared with the CEAEQ (2015) system, which has three times greater tolerances: 3.0%–4.0% of pH unit (Table 1). The NAPT (2001) system provides greater accuracy for soil quality control. Therefore, there would be less risk of misinterpretation of active acidity (pH_W), and especially of lime requirements, based on pH_SMP. However, our results are discussed in relation to the criteria of CEAEQ (2015) because it is mandatory in Québec. Based on CEAEQ (2015) criteria, the pH_SMP does not raise any problem of temporal variability for either G1 or G2 soil. However, Lemire et al. (2005) raise the problem of accuracy in the pH_SMP measurement method, whose accepted tolerance in Québec is already high, that is, ± 0.2 pH_SMP units, which would result in an error of estimated lime requirements of ±2.5 t·ha⁻¹. In its official website, the MAPAQ (2021) mentions that this error is unpredictable and difficult to consider when making lime recommendations. On the other hand, for pH_W, the CEAEQ (2015) system detects variations that exceed the AQLR limit (Table 1). As shown in Fig. 3 and Table 2, pH_W decreases below the AQLR limit in the period from weeks 14 and 23, i.e., from July to October. Other studies, such as those by Collins et al. (1970) and Hoskinson et al. (1999), also obtained the same trends of significant decreases in soil pH_W in the summer to early fall. The range of nonacceptable pH_W variation is therefore an average of 10 consecutive weeks, starting in July for G1 soil (Table 2). It is therefore a type of continuous exceedance of the AQLR limit. This drop in pH_W values below the AQLR limit is likely due to climatic conditions during the summer and early fall. This 10-week period is the least suitable for sampling because low pH_W values will lead to poor diagnosis and poor decision-making regarding liming actions. On the other hand, fall and early spring are the best times to sample. Their variations are both similar, nonsignificant and within the acceptable fluctuation according to criteria of CEAEQ (2015). This is consistent with the interpretation grids of soil acidity obtained in Québec from liming tests (Tran and Van Lierop 1982) where sampling was carried out during this period of the year. Since pH_W is the most widely used acidity indicator, it must necessarily be associated with a period for sampling, that is, before July, or in the fall, starting in October. Figure 3 also shows the ACF autocorrelation function reflecting the degree of similarity in pH_W between consecutive weekly samples. Exceeding the critical limits, that is, the two dotted horizontal lines in Fig. 3, derived from the ACF function, shows autocorrelations significantly different from zero (p < 0.05) only for lags of one and two consecutive weeks. In other words, regardless of the sampling period, the statistical similarity of pH_W over time is only maintained for a maximum time lag of two weeks. This low autocorrelation is due to the logarithmic pH scale, which decreases the amplitude of the variation in active acidity and makes repeated measurements over time less autocorrelated. For the other two locations (b and c), the pattern is the same as that in Fig. 3 (not shown, but the main observations are summarized in Table 2) and the autocorrelation is also ≤2. Despite the low reliability of pH_W and pH_SMP, they are still the only investigative and diagnostic tools for the acidity of agricultural soils in Québec, Canada, and in several U.S. states. Since the 1980s, Follett and Follett (1980) have suggested adding other indicators such as soil texture, clay type, organic matter content, exchangeable aluminum, and cation exchange capacity to better diagnose acidic soils and assess their lime requirement more accurately.

Variation in diagnostic indicators of soil nutrients

Five soil fertility indicators (Ca_M3, Mg_M3, Al_M3, CEC_predicted, B_M3) were not affected by seasonality since their SDRs were ≤ AQLRs (Table 1). However, six indicators (K_M3, OM, Mn_M3, Zn_M3, Cu_M3, Fe_M3) had a SDR ≥ AQLR for G1 and (or) G2 soils. A representative case is shown in Fig. 4 for K_M3 in soil G2 at location a. Result summary is shown in Table 2 for all six indicators, for both soils G1 and G2 and for all three locations (a, b, and c). Figure 4 shows a range of variation of K_M3 that is not acceptable or below the AQLR limit because of three values, one isolated at week 12 and two at week 20 and 21. This variation of 100 mg K·kg⁻¹ between week 12 (summer) and weeks 20 and 21 (early autumn) is very important. The causes and mechanisms leading to this situation remain unclear due to the complex environment of the crops (Hoskinson et al. 1999). For this group of indicators, only Cu_M3 showed a significant and continuous exceedance in the fine textured soil, G1, compared with the average of the five weeks of spring (Table 2). However, this Cu_M3, which has the greatest seasonal variation with SDRs of 47%–55% (Table 1), also shows that it is the least concentrated element, at 1.4–1.5 mg Cu·kg⁻¹. These extremely low values make estimates of the concentration of Cu_M3 usually less reliable than the measurement of the other indicators. Therefore, none of the five indicators, pre-examined by descriptive statistics (Ca_M3, Mg_M3, Al_M3, CEC, and B_M3), and none of the six indicators examined by time-series decomposition (K_M3, OM, Mn_M3, Zn_M3, Cu_M3, and Fe_M3) should be subject to a representative sampling calendar effect.

Variability of agro-environmental diagnostic indicators

The two indicators, P_M3 and SPS, are used to prescribe phosphorus rates (CRAAQ 2010) and to prevent pollution and eutrophication of surrounding surface waters (MDDELCC 2017). The average contents of P_M3 available phosphorus for soils G1 and G2 are 122 and 617 mg P·kg⁻¹, respectively, and those are 110 and 571 mg P·kg⁻¹, respectively (Table 2). The two latter averages indicate soils with high phosphate fertility, even for P-demanding crops, since they exceed the agronomic critical thresholds of 90 mg P_M3·kg⁻¹ (CRAAQ 2010; Khiari et al. 2000). The SPS averages of the G1 and G2 soils are, respectively, 5.7% and 27.9%, and the values are 5.0 and 26.6%. These values show that the fine-textured G1 soil is clearly below the environmental critical thresholds of 7.6% (MDDELCC 2017) and 8% (CRAAQ 2010). However, soil G2 is supersaturated in phosphorus and greatly exceeds the environmental critical threshold of 11% (CRAAQ 2010). As an example, the seasonal component of the SPS time series is shown in Fig. 5 for soil G1 and location a. Results are shown in Table 2 for both indicators, the two soils and the three locations. Figure 5 shows quasi-periodic oscillations characterized by alternating periods of SPS stability within the lower LL_validity and upper UL_validity limits (eqs. 3 and 4) and periods of increases in SPS above UL_validity (above the blue line). In Fig. 5, SPS is characterized by a long period of 16 weeks, between week 4 and 19, where the SPS values exceed UL_validity. These high SPS values cause an environmental risk of soil P saturation (CRAAQ 2010) going from a lower medium risk class (4.0 ≤ SPS < 6.5) in spring, to a higher medium risk class (6.5 ≤ SPS < 8.0) in summer for 30% of the samples, and up to a high-risk class (8.0 ≤ SPS < 14) for 14% of the samples. When compared with the critical value of 7.6%, above which the regulation (MDDELCC 2017) considers the risk of diffusion of P to surface water to be high, this results in an 18% exceedance. To select sampling windows for collecting representative data, period of high risk must be avoided. The twelve graphs of the agro-environmental indicators: two indicators (SPS and P_M3) × two soils (G1 and G2) × three locations (a, b, and c) showed patterns of seasonality similar to those of Fig. 5. For field G1, the six seasonal time series resulted in nonrepresentative sampling windows between weeks 4 and 19 (Table 2). In Québec (CRAAQ 2010) sampling should be done in the spring (weeks 0 to 5). Outside this interval, SPS is likely to vary considerably with an upward trend up to week 19 and then a return to the spring values during the fall (Fig. 5). During summer, soils warm up, which stimulate the biological activity and make phosphorous more available (Habibiandekordi et al. 2015). For the G1 field, spring SPS values of 5.4 (Table 2) are not problematic and can increase by 2.5 % by week 11 (Fig. 5) and exceeds the environmental threshold of 7.6 %, established by the regulation (MDDELCC 2017) for soils with a clay content >30%. For SPS values above the 7.6% threshold, the phosphate fertilization strategy must aim at reducing soil saturation below the threshold. A sample taken during the wrong time period may under-estimate phosphate fertilization amounts. An over-estimation of SPS would falsely reduce the areas available for manure application, and force farmers to rent land to dispose of their manure, resulting in extra time and costs. For the G2 medium textured field, SPS of almost 28% above the environmental threshold of 13% (MDDELCC 2017). In Table 2, when comparing the intervals of exceedance for P_M3 and SPS, we observe an early seasonal variation starting in week 4 or 5 for G1 soil. But for G2 soil, it is shifted between weeks 8 and 14. Since the G2 soil is much more saturated with phosphorus than the G1 soil, this could explain such a delay for the seasonal effect of P_M3 and SSP. Fig. 5b shows significant autocorrelations for five consecutive weeks. Therefore, regardless of the sampling period, the statistical similarity of SPS across time is only guaranteed for a maximum of 5 wk. From one sampling cycle to another, it is therefore more consistent and representative to always sample at the same time within 5 weeks of each other Locations b and c yielded a pattern similar to that of Fig. 5b (not shown), and almost the same 3–5 wk intervals where the SPS measures are significantly autocorrelated (Table 2).