Spatial Variability of Soil Phosphorus Indices under Two Contrasting Grassland Fields in Eastern Canada

: Phosphorus (P) is an essential nutrient for grassland production systems. However, continuous applications of P fertilizers result in soil P accumulations, increasing the risk of P losses in runoff and erosion. This study aims to investigate the ﬁeld-scale variability of soil-test P (STP) in two contrasting grassland ﬁelds using descriptive statistics and geostatistics for accurate recommendations on soil sampling strategy and sustainable approaches to P management. A young grassland (YG; 2 years) and an old grassland (OG; 10 years under permanent pasture) were classiﬁed as humo-ferric podzol and received organic fertilizers. Soil samples were collected in 16-m by 16-m triangular grids at two depths (0–5 and 5–20 cm). They were analyzed for available P and other soil elements extracted using the Mehlich-3 method (M3). The agri-environmental P saturation index (P/Al) M3 was calculated. Phosphorus accumulation was observed in OG (0–5 cm) as a result of long-term manure applications. Repeated applications of organic fertilizers can impact the long-term buildup of soil P, thus decreasing soil P va-riability and spatial dependence in permanent grasslands. A soil sampling strategy focusing on the 0–5 cm layer should be retained in permanent grasslands for sustainable P recommendations in Eastern Canada.


Introduction
Phosphorus (P) is an essential nutrient for most crops, including associated plants of grassland communities. When it is applied in excess of crop requirements, P can accumulate at the soil surface, leading to P runoff via soil erosion and to P contamination of the environment, including eutrophication of surface waters such as lakes [1][2][3]. To prevent this P contamination in the environment, the degree of P saturation of soils (%) was calculated from the ratio of phosphorus extracted with ammonium oxalate (Pox) on the maximum retention capacity of the P (CRP) [4]. The critical P saturation value of 25% has been defined as an international benchmark intended to prevent the risk of P loss and water contamination increases [5].
Grasslands are defined as terrestrial ecosystems dominated by herbaceous and shrub vegetation [6,7] and are distributed according to landscape and site characteristics (soil, climate, water availability) and farm structures. Permanent grasslands include grazed pastures and grasslands harvested for hay and silage, or a combination of the two, and are located primarily in Western Europe, including the United Kingdom and Ireland [8]. Canadian grassland fields covered an area of 1.96 million ha in 2016 [9]. In the Eastern Canadian context, grassland cultivation represents an important economic opportunity for local farmers to increase their revenue. For instance, 50% of agricultural land was under intensive permanent grassland cultivation in the province of Quebec, producing 6.3 In the OG field, calf and pig manures were applied at the rate of 30 Mg ha −1 and 20 Mg ha −1 , respectively, giving approximately 67 kg P ha −1 year −1 over 10 years based on reference values of P published in CRAAQ [26]. At the beginning of May and after the first cut (mid-June), organic manures were surface applied in the OG field using a low-ramps spreading system.
An intensive soil sampling using a grid of 16 m by 16 m, including three cross-shaped designs with several inter-points at about 4 m apart, was performed in spring 2013, providing 151 and 149 georeferenced sampling points for the YG and OG fields, respectively ( Figure 1). A composite soil sample of four cores was taken within a radius of 1 m around each sampling point at the two soil depths (0-5 cm and 5-20 cm) using a 0.05 m-diameter Dutch auger (Eijkelkamp company, Giesbeek, Netherlands). Thus, a total of 302 and 298 soil samples were collected from the YG and OG fields, respectively.

Statistical and Geostatistical Analyses
Descriptive statistics were carried out using SAS Software version 9.4 [29] to investigate Both soils were classified as humo-ferric podzols [24] and attributed to the Beaurivage soil series [24]. Soil surface texture for the YG and OG fields varied from sandy loam to loamy sand. They were classified to the same G3 soil texture group (i.e., texture with clay ≤ 20%) utilized in the province of Quebec recommendations [11]. Both fields were moderately well-drained. The annual winter mean temperature is −9.4 • C, and the mean growing season temperature is 16.6 • C [25]. Mean annual precipitation is 1120 mm, of which 320 mm is snow [24].
The main grass species consisted of Ladino white clover (Trifolium repens) cultivated in association with timothy (Phleum pratense L.) in the YG field and tall fescue (Festuca arundinacea Schreb.) in the OG field. In 2004-2005, 2008-2009 and 2012-2013, cattle manure was applied and incorporated twice a year in the first 5 cm using a low-ramps spreading system in the YG field at the rate of 18,700 L ha −1 , equivalent approximately to 20 kg P ha −1 year −1 . Cattle manure was applied in spring and later after the first cut (mid-June).
In the OG field, calf and pig manures were applied at the rate of 30 Mg ha −1 and 20 Mg ha −1 , respectively, giving approximately 67 kg P ha −1 year −1 over 10 years based on reference values of P published in CRAAQ [26]. At the beginning of May and after the first Agronomy 2021, 11, 24 4 of 18 cut (mid-June), organic manures were surface applied in the OG field using a low-ramps spreading system. An intensive soil sampling using a grid of 16 m by 16 m, including three crossshaped designs with several inter-points at about 4 m apart, was performed in spring 2013, providing 151 and 149 georeferenced sampling points for the YG and OG fields, respectively ( Figure 1). A composite soil sample of four cores was taken within a radius of 1 m around each sampling point at the two soil depths (0-5 cm and 5-20 cm) using a 0.05 m-diameter Dutch auger (Eijkelkamp company, Giesbeek, Netherlands). Thus, a total of 302 and 298 soil samples were collected from the YG and OG fields, respectively.
Pearson's correlation coefficients (r) between soil P indices (P M3 , (P/Al) M3 ) and other soil chemical properties (Al M3 , Ca M3 , Fe M3 , TC, and pH water ) were also calculated for each soil layer using SAS. Data were analyzed using the CORR and univariate procedures. Multiple regression equations were generated for the soil P indices for each soil depth.
Geostatistical analyses were performed for all soil parameters. Data were achieved for each of the soil layers using the GS+ version 9 software [31] geostatistical computations and model validations. The spatial structure of the different soil chemical properties was evaluated via isotropic and anisotropic semivariograms. Experimental semivariograms, the main component of kriging, are an effective tool for assessing spatial variability [32]. Semivariogram parameters for each theoretical model (spherical, exponential and Gaussian) were generated. The corresponding nugget (C 0 ), partial sill (C), sill (C 0 + C), and range values of the best-fitting theoretical model were calculated. The partial sill ratio (C/(C 0 + C)), expressed as the percent of total semivariance, was used to define the spatial dependency or structure of the soil variables for the two fields. Semivariograms with a partial sill ratio of < 25%, 25% to 40%, 40% to 60%, 60% to 75%, or > 75% were considered to have a low, low-moderate, moderate, strong-moderate, or strong spatial dependency, respectively [33]. After the suitable theoretical model for each data set and the corresponding semivariogram parameters were selected, spatial variability maps were generated using 1 m × 1 m block kriging interpolation. Kriging map reliability was evaluated using cross-validation analysis (R 2 CV ) [34].

Descriptive Statistics of Soil Phosphorus Indices and Other Chemical Properties
For the YG field, the mean values of soil P indices and other soil chemical properties were similar in both soil layers (0-5 cm and 5-20 cm), with the exception of Ca M3 (Table 1). Overall, the CVs ranged from 5% to 81%, with the highest obtained for the (P/Al) M3 , and were classified as follows: low for soil pH water and Al M3 , moderate for Fe M3 , high for Ca M3 and TC, and very high for the two soil P indices. For the OG field, soil P indices and other soil chemical properties were higher in the 0-5 cm layer than in the 5-20 cm layer, with the exception of Al M3 , for which the mean concentrations were similar (Table 1). On average, the CVs ranged from 6% to 64% in both soil layers, with the highest value obtained for the (P/Al) M3 , and were classified as follows: low for soil pH water and Al M3 , moderate for Fe M3 and TC, and high for Ca M3 . The CVs of the soil P indices varied from high to very high in both soil layers. The intensity of variability of the soil P indices was higher in the YG compared to OG for both soil layers.

Geostatistical Parameters of Soil Phosphorus Indices and Other Chemical Properties
In the YG field, soil P indices and the other soil chemical properties (except for Ca M3 ) were fitted with geostatistical spherical models, and pure nugget models (PNs) were used for Ca M3 in both soil layers. Spatial ranges varied from 28 to 81 m. Spatial ranges of soil chemical properties extended beyond the 16 m × 16 m sampling grid, indicating that the sampling grid used is appropriate. Spatial dependence ratios of soil P indices and other soil chemical properties ranged from 58% to 100% for the two soil layers and were classified according to Whelan and McBratney [33] as follows: Fe M3 ranged from moderate to strong-moderate, pH water and soil P indices were classified as strong-moderate, and TC and Al M3 were classified as strong ( Table 2).
In the OG field, soil P indices and the other soil chemical properties (except for Al M3 and Ca M3 ) were fitted with geostatistical spherical models, and PNs were used for Ca M3 and for Al M3 in the 0-5 cm soil layer ( Table 2). Spatial ranges varied from 27 to 76 m and extended beyond the 16 m × 16 m sampling grid, indicating that the sampling grid used is appropriate. Spatial dependence ratios of soil P indices and the other soil chemical properties ranged from 26% to 69% in both soil layers. According to Whelan and McBratney [33], the spatial dependence of Al M3 was low-moderate in the 5-20 cm layer. Spatial dependence values were moderate for soil pH water and Ca M3 , strong-moderate to moderate for Fe M3 , low-moderate to moderate for soil P indices, and strong-moderate for soil TC in both soil layers ( Table 2). Spatial dependency of soil P was lower under OG compared to YG for both soil layers.
In the YG field, cross-validation coefficients (R 2 CV ) ranged from 0.31 to 0.54 for soil chemical properties and from 0.31 to 0.40 for soil P indices. In the OG field, R 2 CV were generally lower. The R 2 CV values ranged from 0.20 to 0.36 for the soil P indices. High R 2 CV values (R 2 CV > 0.6, [35]) indicate a good fit for kriged map reliability. Consequently, the R 2 CV values obtained for the YG field may indicate relative mean fits for map reliabilities for most soil chemical properties, including soil P indices, while the fits were weaker for map reliabilities of soil chemical properties in the OG field.

Spatial Distribution of Soil Phosphorus Indices and Other Chemical Properties
The P M3 and (P/Al) M3 spatial distribution maps showed visual similarities for each soil layer studied (Figure 2a,b,g,h and Figure 3a,b,g,h) . In the YG field, spatial distribution values of soil P indices, TC and Fe M3 were relatively similar in both soil layers (Figure 2a,b,dg,i,j). Spatial distribution values of Al M3 were lower in the 0-5 cm soil layer compared to the 5-20 cm soil layer (Figure 2c,i). Table 2. Geostatistical parameters of the soil chemical properties for two soil layers (0-5 cm) and (5-20 cm) from the young and old grassland fields. However, in the OG field, spatial distribution values of soil PM3, (P/Al)M3, TC, and FeM3 were higher in the 0-5 cm soil layer than in the 5-20 cm soil layer ( Figure  3a,b,d,e,g,h,j,k). For instance, spatial distribution values of these soil chemical properties varied from 46 to 232 mg kg −1 , 2.8% to 13%, from 20 to 55 g kg −1 and from 135 to 324 mg kg −1 against 15-139 mg kg −1 , 1.1% to 8%, 20-40 g kg −1 , 135-297 mg kg −1 in 0-5 cm and 5-20 cm OG soil layers, respectively. In the OG field, no visual similarities were noted between AlM3 0-5 cm and 5-20 cm or between CaM3 0-5 cm and 5-20 cm. This was not surprising, as the PN models were obtained with the semivariogram analysis.    However, in the OG field, spatial distribution values of soil PM3, (P/Al)M3, TC, and FeM3 were higher in the 0-5 cm soil layer than in the 5-20 cm soil layer ( Figure  3a,b,d,e,g,h,j,k). For instance, spatial distribution values of these soil chemical properties varied from 46 to 232 mg kg −1 , 2.8% to 13%, from 20 to 55 g kg −1 and from 135 to 324 mg kg −1 against 15-139 mg kg −1 , 1.1% to 8%, 20-40 g kg −1 , 135-297 mg kg −1 in 0-5 cm and 5-20 cm OG soil layers, respectively. In the OG field, no visual similarities were noted between AlM3 0-5 cm and 5-20 cm or between CaM3 0-5 cm and 5-20 cm. This was not surprising, as the PN models were obtained with the semivariogram analysis.   However, in the OG field, spatial distribution values of soil P M3 , (P/Al) M3 , TC, and Fe M3 were higher in the 0-5 cm soil layer than in the 5-20 cm soil layer ( Figure  3a,b,d,e,g,h,j,k). For instance, spatial distribution values of these soil chemical properties varied from 46 to 232 mg kg −1 , 2.8% to 13%, from 20 to 55 g kg −1 and from 135 to 324 mg kg −1 against 15-139 mg kg −1 , 1.1% to 8%, 20-40 g kg −1 , 135-297 mg kg −1 in 0-5 cm and 5-20 cm OG soil layers, respectively. In the OG field, no visual similarities were noted between Al M3 0-5 cm and 5-20 cm or between Ca M3 0-5 cm and 5-20 cm. This was not surprising, as the PN models were obtained with the semivariogram analysis.
In general, spatial distribution maps also revealed visual associations between soil P indices and other soil elements, indicating the influence of these different soil chemical properties on the soil-available P. For instance, in the YG field, the Al M3 and TC concentrations from the soil layers were inversely related to the soil P indices (Figure 2a-d,g-j).
In the OG field, the Al M3 concentrations in the OG soil layer (5-20 cm) were inversely related to the soil P indices in OG (5-20 cm) (Figure 3g-i). The high Ca M3 concentrations correspond to soil P accumulation regions in the OG soil layer (5-20 cm), confirming P M3 -Ca M3 accumulation (Figure 3g,h,l). These results showed that Al M3 and TC content might influence soil P retention capacity in both YG soil layers, whereas Al M3, Fe M3 and Ca M3 may affect soil P in both OG soil layers.

Relationships Between Soil Phosphorus Indices and Other Chemical Properties
In the YG field, Pearson's correlation coefficients (r) with the highest significance were observed between P M3 and Al M3 , followed by TC in both soil layers (Table 3). For (P/Al) M3 , the highest r values were obtained with TC in both soil layers. In the OG field, significant relationships were observed between P M3 and Al M3 , Ca M3 and soil pH water in both soil layers. However, the highest r was observed between (P/Al) M3 and Ca M3 in each soil layer (Table 3). Thus, soil P M3 in the YG field may be influenced significantly by Al oxides and hydroxides (soil Al), soil organic matter (TC) and soil solution (soil pH water ), whereas in the OG field, Al oxides/hydroxides (soil Al), calcium carbonate (soil Ca), soil solution (soil pH water ) and iron oxides/hydroxides (soil Fe) may have a significant influence on soil P accumulation. Table 3. Pearson's correlation coefficients (r) of the soil-available phosphorus (P) indices (P M3 , and (P/Al) M3 in relationship with soil chemical properties for two soil layers (0-5 cm) and (5-20 cm) from the young and old grassland fields. Multiple regression equations were generated for the soil P indices in the YG and OG soil layers to investigate the influence of these soil chemical properties (Table 4). In the YG field, the regression coefficients (R 2 ) were very significant, ranging from 0.789 to 0.893, which means that 79% to 89% of soil-available P variations may be explained by Al M3 , TC and soil pH water in both soil layers. Similarly, in the OG field, the R 2 varied between 0.532 and 0.659, indicating that 53% to 66% of soil-available P accumulations in both soil layers may be affected significantly by Al M3 , soil pH water , Ca M3 and Fe M3 . Consequently, permanent grassland may accumulate soil P while significantly affecting the effect of these soil chemical properties.

Soil Phosphorus Stratification and Its Environmental Implications
In this study, soil P accumulation was observed in the 0-5 cm OG soil layer, mainly due to multiple applications of organic fertilizers (calf slurries, pig manure, and cattle excreta) over a long time period. Indeed, soil P accumulation in the 0-5 cm OG top layer may also be due to the elimination of soil disturbance (no tillage, compared to the conventional tillage) under permanent grasslands.
Many studies [3,19,36,37] have shown that soil P accumulation in the surface layer from intensively cultivated lands is linked to increased applications of P fertilizers, organic inputs, soil land-use conversion and management. For instance, Sun et al. [3] reported that significant increases in soil-available P (P Olsen) were attributed to intensive local applications of swine manure over a 15-year period (1982-1997) in the Rugao region (China). Other studies [19,38] found that soil P indices (P-Bray I, P-ammonium-acetate-EDTA (P AAE ): referred to as the Swiss reference method) varied significantly in different land uses and crop systems, including permanent grassland fields.
This study revealed that repeated application of organic fertilizers raised soil P accumulation in the first five soil centimeters, increasing the potential agri-environmental risk for P pollution in permanent grasslands. The mean (P/Al) M3 indicator raised up to twice (7%) in the 0-5 cm OG soil layer, relative to the 0-5 cm YG soil layer (3%). These results have several environmental implications for the sustainable management of soil P in grassland fields. Indeed, the (P/Al) M3 mean values obtained from the OG site were close to the national critical threshold P saturation (15%) elaborated by the Province of Quebec for coarse-textured soils, meaning that (P/Al) M3 values higher than 15% increased up to water pollution or groundwater contamination in this Eastern Canadian region.
Moreover, these results also revealed potential agri-environmental impacts for the management of permanent grasslands in the Eastern Canadian region. According to the national Agricultural Operations Regulation (REA) adopted by the province of Quebec in 2002, these obtained results were close to the threshold level (13%) for P saturation (P/Al) M3 for agricultural lands, including grassland fields. Soil P stratification or buildup is an important agri-environmental concern because the excess surface P may be lost in runoff [39]. At the same time, however, lower concentrations at depth in the rooting zone may reduce crop yields [40].

Variability of Soil Phosphorus Indices and Other Soil Chemical Properties
The highest CVs were obtained for soil P indices and the lowest CVs for soil pH measurements. These results were similar to other studies [20,[41][42][43][44]. A low CV for soil pH water may be due to the logarithmic scale of pH measurement [45] or may be attributable to management practices, including regular soil liming activities in intensive agriculture [46]. West et al. [47] reported that a high CV (47%) for soil-available P could be explained by chemical forms of measured P that are highly labile.
Many studies [20,[48][49][50] have likewise reported very high CVs for soil-available P (CV > 50%), particularly on Irish soils under intensive grassland farming. McCormick et al. [50] reported a CV of 57% (n = 334) for a 50-ha parcel of land encompassing 15 small grassland fields ranging in size from 0.94 ha to 5.65 ha in Northern Ireland. Fu et al. [49] obtained a CV of 68% (n = 537) for a 52-ha grassland field located in southeastern Ireland.
This study shows that a reduction in the CV of soil P was observed in permanent grasslands. For example, Jia et al. [18] reported a large decrease in the CV of soil P from 151% to 55% after 25 years of intensive farming practices in China. However, some authors [18,19,30,48] found that long-term fertilization practices or intensive land-use changes, including permanent grasslands, may influence CV values of soil P indices. Roger et al. [19] reported extremely high CV values for soil-available P (P AAE ) at permanent grasslands sites (CV = 107%) compared to mountain pasture sites (CV = 87%) and cropland sites (CV = 70%). From their study, the CVs of soil P saturation [(P ox /(Al+ Fe) ox ) × 100] ranged from 35% to 38%. In general, CV values are good indicators of the intensity of soil variability, but not of the nature of this variability (structured or randomized variability), nor do they provide an accurate model of soil variability [42].

Spatial Dependence of Soil Phosphorus Indices and Other Soil Chemical Properties
Spatial dependence (R = (C/C 0 + C)) describes the spatial structure of a soil property [51]. Some authors [20,52,53] have indicated that soil chemical properties with strong spatial dependence (R > 75%) are mainly influenced by soil intrinsic factors (soil formation factors, parental material, soil types or series), those with low spatial dependence (R < 25%) are strongly influenced by extrinsic factors (agronomic/cultural or fertilization practices, intensive use of organic fertilizers or mineral amendments), and those with moderate spatial dependence (25% < R < 75%) result from the interaction between intrinsic and extrinsic soil factors, which codetermine the soil property analyzed. The PNs models are characterized by a spatial structure smaller than the sampling grid because of an inappropriate sampling scale [44].
The spatial dependence of soil P indices in both YG and OG soil layers was moderate (31% to 71%) as a result of the interaction between soil intrinsic and extrinsic factors, such as intensive fertilization or farming practices. Many authors [18,20,33,[54][55][56][57][58] have found moderate spatial dependence for soil P, while others [3,19,42,50] have reported strong spatial dependence for soil P (R > 75%) in their respective studies conducted on different scales. The strong spatial dependence of soil-available P observed in these studies may be attributable to the strong influence of intrinsic factors, such as different soil types or soil series, on the spatial distribution of soil-available P at different agricultural, local, or regional scales.
Spatial dependence values for soil pH water and Fe M3 in the YG and OG soil layers were moderate (41% to 74%), owing to the interaction between soil intrinsic and extrinsic factors. Arrouays et al. [44], Cambouris et al. [42], and Di Virgilio et al. [41] obtained similar results for soil pH water . Spatial dependence values for soil Al M3 and TC ranged from strong to moderate in the YG and OG soil layers, respectively. Spatial dependence values obtained for Al M3 in YG soil layers were lower than those reported by Cambouris et al. [42]. For soil TC, Kokulan et al. [59] and Arrouays et al. [44] reported similar results for YG and OG soil layers (0-5 cm and 5-20 cm). Consequently, Al M3 and TC in both YG soil layers may result from the strong influence of soil intrinsic factors, whereas Al M3 and TC in both OG soil layers result from the interaction between soil intrinsic and extrinsic factors.
This study also revealed that a combination of intensive fertilization practices over a longer period reduced the spatial dependence of soil P in both OG soil layers. This may be explained by the increasing influence of soil extrinsic factors, such as P-rich fertilization, which modified the long-term spatial structure of soil P in the OG field. Other authors [3,18,48] reported a reduction of spatial dependence for soil P on Chinese and Irish soils at different regional and agricultural scales. These studies demonstrated that a weak spatial dependence of P distribution could be explained by various intensive agricultural farming practices (i.e., increasing fertilization, long-term P-rich mixed fertilization, and strong intensification of sheep grazing pasture activities).
From an agronomic perspective, the reduced spatial structure of soil P may affect the use of delineated MZs for sustainable management of soil P in grassland fields. A strong spatial structure can have a positive effect in terms of delineating MZs, while a weak spatial structure makes it difficult to use MZs, owing to the lack of spatial dependence.
The use of MZs is an agricultural approach based on the existence of within-field spatially structured soil and crop variability [42,45]. Subdividing a field into MZs is an effective way of understanding or controlling spatial variability of soil and crop properties in order to optimize soil nutrients and crop management [35]. Delineation of MZs has been tested successfully for various soil chemical properties [17,60] and crop fields [42,45], including old grassland fields [43], as well as different field sizes [35,45]. However, there are few examples of the use of this approach with the aim of optimizing the management of P variability in grassland fields.
The present study on soil P variability was carried out in 2.5-ha fields. Consequently, it would be mitigated to delineate MZs for these fields due to (1) the reduced spatial structure of soil P, which is caused by an increase in extrinsic soil factors, and (2) the size of the fields with respect to the crop studied. Further research on larger fields is needed to better define P spatial structure in grasslands.

Spatial Distribution of Soil Phosphorus Indices and Other Soil Chemical Properties
The R 2 CV values provide an indication of the reliability of spatial distribution maps. Weak R 2 CV values were obtained for soil P indices and other chemical properties in the OG field compared to the YG field. Use of R 2 CV values in grassland fields aims to evaluate the reliability of kriged spatial distribution maps of soil P indices for developing better agri-environmental strategies for sustainable P management in the context of precision agriculture. Thus, R 2 CV should also be considered as an important reliability factor for mapping soil P indicators in grassland fields to prevent P pollution. Some authors [58] reported lower R 2 CV values (0.0-0.20), while others [20,35] found higher R 2 CV values (0.59-0.62) for soil-available P. In this study, the lowest R 2 CV values obtained may be explained by the weak spatial structure of the soil P obtained, as suggested by Kravchenko [61].
Despite their low reliability (R 2 CV < 0.4), spatial distribution maps of soil P indices (Figure 3a,b) confirmed that the soil P accumulation observed in the OG soil layer (0-5 cm) is attributable to the repeated application of organic fertilizers. Furthermore, in this study, high P concentration areas were found in the YG (81-129 mg kg −1 ; 5% to 8%) and OG (139-232 mg kg −1 , 8% to 13%) fields, probably due to uniform applications of P-rich organic fertilizers. The YG field was fertilized with cattle manure under corn and soybean rotation and included 2 years of grassland over a period of 10 years, whereas the OG field was fertilized for 10 years with multiple applications of organic fertilizers (e.g., calf slurries, pig manure, and cattle excreta).
Spatial distribution maps for soil Al M3 in both sites also confirmed the strong influence of high available Al content on soil P adsorption capacity in grassland soils in Eastern Canada. In a previous study, Tran et al. [62] highlighted this possibility in relation to Quebec agricultural soils. Taking into account soil Al M3 concentrations, these authors classified Quebec soils into three groups according to their capacity for soil P adsorption: (1)  revealed that grassland fields in Eastern Canada have an average-high capacity for soil P adsorption. From these results, the spatial distribution of soil Al M3 could be evaluated as a potential indicator for soil P adsorption in grassland fields. This approach would aim to implement sustainable management of P as mineral or organic fertilizer, based on VRA of P fertilizer for preventing pollution. For instance, other studies [51] found that Al M3 could be considered as an important indicator of maximum P adsorption capacity for a VRA of P fertilizer in the agricultural fields from Eastern Canada.
The spatial distribution maps for soil P indices and other soil elements (Al M3 , TC, Fe M3 and Ca M3 ) in YG and OG soil layers developed in this study confirm the observed correlation, particularly between soil P indices and soil Ca M3 and Fe M3 in the 5-20 cm OG soil layer. The southern area of the OG field (5-20 cm) may reflect soil P accumulation linked to Fe oxides or hydroxides and neoformation of Ca phosphates as a result of long-term fertilization (Figure 3k,l).

Relationships Between Soil Phosphorus Indices and Other Soil Chemical Properties
The results showed differences between YG and OG. Soil P may accumulate or increase in OG, and there may be significant changes of different soil parameters (Al M3 , TC, Ca M3 , Fe M3 ) on soil P indices, which are probably due to a decrease in soil pH.
Decreasing soil pH is common in intensively farmed soils due to heavy applications of inorganic N fertilizers [36,63,64]. It is well known that soil pH plays an important role in P solubility and adsorption processes, such as P adsorption onto soil oxides, minerals and soil organic matter [65,66]. For instance, Fu et al. [48] mentioned a significant relationship between soil-available P and soil pH (r = 0.69, n = 396, p < 0.01) in a grassland field used for grazing for 16 years. In this study, we found significant relationships between soil P indices and soil pH (r = 0.43 for P M3 , p < 0.001; r = 0.46 for (P/Al) M3 , p < 0.001) in the 0-5 cm OG soil layer.
Results also show that soil P accumulation in the OG soil layers (0-5 cm and 5-20 cm) may be related to soil-available concentrations of Al, Ca and Fe resulting from a decrease of soil pH. Many studies [3,36,[67][68][69] reported that various agricultural and farming practices affect soil P accumulation. Other studies [70,71] observed relationships between soil P buildup and Fe/Al/-P, attributed to the intensive use of organic and inorganic P fertilizers in cultivated acid soils.
Positive relationships observed between soil P M3 buildup and soil Ca M3 in the OG field may be indicative of soil P accumulation as Ca-bound P or neoformation of Ca phosphates after long-term P fertilization. This possibility has previously been advanced by several authors [36,72]. These earlier studies showed that the P-Ca relationship also depends on the relationship between soil pH and Ca. In the present study, we reported a significant relationship between soil pH and soil Ca M3 in the 0-5 cm (r = 0.65; p < 0.001) and 5-20 cm (r = 0.57; p < 0.001) OG soil layers.

Agronomic and Environmental Implications
Repeated applications of organic fertilizers and other agronomical practices such as no-till soil may increase soil P accumulation in the first 5 cm of soil in agricultural fields [58,73,74], including grasslands [43]. This points to the importance of developing an effective soil sampling optimization strategy to provide better guidance for practical grassland management. Fertility management in grasslands is related to sampling depth. Although some studies [75][76][77] demonstrated that the 0-5 cm layer might provide the majority of N and P in grassland fields, a deeper soil profile may be important to sustaining healthy grasslands.
Soil sampling strategy in grassland fields represents a balancing act between fertility management and minimizing environmental impacts (i.e., P pollution, water quality). The current 0-17.5 cm sampling strategy may be more suitable for yield optimization, while the 0-5 cm sampling strategy would be better for more sustainable P recommendations in grassland fields to reduce P pollution in the environment. According to the results of this study, we recommend the use of a soil sampling strategy focusing on the 0-5 cm layer in permanent grassland fields, given the high risk of soil P stratification associated with P-rich fertilization, instead of the soil sampling strategy (0-17.5 cm) currently used in Eastern Canada, in order to derive more sustainable P recommendations for grasslands. Furthermore, other studies [1,50,76,78] recommended soil sampling strategies in grassland fields, focusing on the 0-5 cm, 0-7.5 cm or 0-10 cm layers in other countries.
In this study, reduction of soil P variability under permanent grasslands was observed, inversely to other previous studies achieved within grassland fields [8,47]. Therefore, sustainable management of soil P in grassland fields should aim to better control withinfield soil P variability. Spatial distribution maps of soil P indices should be used to guide variable-rate P fertilizer applications in order to prevent under-and over-fertilization with P in grassland fields.
Furthermore, new technology advances (reflectance and evaluation of P) and the democratization of prices for VRA will make possible the potential use of this technique in the grasslands in the near future. Indeed, the implementation of local or environmental policies about agri-environmental and sustainable management of P will encourage grassland farmers to turn to these new technologies in the context of PA. Hence, although it is true that the additional economic costs can be high, the VRA approach would be an advantage or benefit for Eastern Canadian grassland farmers in reducing P losses in the environment.
Delineating MZs in grassland fields represents another option that can be used in the future to optimize the management of soil P variability in grassland fields. However, this agronomical practice has some limitations in adopting such an approach. These constraints are (1) a more reduced spatial structure of soil P, caused by a strong influence of extrinsic factors comparatively to the intrinsic factors, and (2) the field size with respect to the crop studied, which may influence the spatial distribution of MZs and their reliability. To overcome these constraints, MZs should be based on stronger soil P spatial structures obtained from reliable semivariograms, which are needed to calculate sampling intervals and develop an accurate site-specific application scheme for soil P in grasslands. In this context, field size is an important soil-related factor that needs to be taken into consideration. Further research on larger fields is needed to better define P spatial structure in grasslands.

Conclusions
This study represents, to our knowledge, the first one on the spatial variability of soil P in grassland fields in Eastern Canada. The main goal of this study was to evaluate the spatial variability of soil-available P in grassland fields managed under intensive farming systems, using descriptive statistics and geostatistical tools for accurate recommendations on soil sampling strategy, and developing potential future and sustainable approaches of P management in grassland fields.
In this study, repeated applications of organic fertilizers under grassland production with minimum soil disturbance increased soil-available P buildup [P M3 , (P/Al) M3 ] at the 0-5 cm OG soil layer. This soil P accumulation increased the potential agri-environmental risk for P pollution in permanent grasslands. Under permanent grasslands, reducing variability and spatial dependence of soil P indices were observed in both soil layers (0-5 cm and 5-20 cm), as a result of the increasing influence of extrinsic soil factors, such as repeated applications of organic fertilizers over a long time. From an agronomic perspective, the reduced spatial structure of soil P may affect the use of delineated MZs for sustainable management of soil P in grassland fields. Despite their low reliability, spatial distribution maps of soil P indices should be used to guide variable-rate P fertilizer applications in order to prevent under-and over-fertilization with P in grassland fields.
In keeping with the results of this study, we recommend a soil sampling strategy focusing on the 0-5 cm layer in permanent grassland fields, given the high risk of soil P stratification caused by P-rich fertilization, as opposed to the soil sampling strategy (0-17.5 cm) currently used in Eastern Canada, in order to generate more sustainable P recommendations for grasslands. Moreover, this primary essay on the spatial variability of soil P was only restricted to two selected grassland fields. Further studies in other and larger grassland fields are required for a better understanding of the spatial variability of P to improve agri-environmental P fertilization for grassland fields.