*2.4. Calculating Net Soil N Mineralization*

Mn was calculated from the end of winter to the beginning of autumn from the mineral N mass balance of a maize crop not fertilized with N, as follows:

$$\text{Mn} = \text{Nf} - \text{Ni} + \text{N}\_{\text{uptake}} + \text{N}\_{\text{leached}} \tag{1}$$

with Ni and Nf corresponding to the soil mineral N content in the 0–90 cm soil profile in March and October, respectively; Nuptake corresponding to N uptake by the plant (kg N ha−1); and Nleached corresponding to nitrate leaching that may occur in spring, after measurement of Ni (kg N ha<sup>−</sup>1).

Ni, Nf, and Nuptake were measured in triplicate. Nleached was estimated using the STICS model [41], which was parameterized with the soil properties of each field, and initialized at the measurement date of Ni. Equation (1) is a simplified approach for estimating the mineral N mass balance, but it is valid in situations without N fertilization. Gaseous N losses can be assumed to be very low and compensated by atmospheric deposition and symbiotic fixation of N.

Since N mineralization depends strongly on weather conditions, we controlled for the influence of weather on mineralization by converting true time into normalized time (nday), using functions integrated into STICS [41], in order to model the effects of temperature and soil water content on N mineralization. The effect of soil temperature on mineralization is described by a logistic function, which is roughly exponential from 0–25 ◦C. This function is similar to an Arrhenius function, with an activation energy of 78 kJ mol−<sup>1</sup> K−<sup>1</sup> from 0–35 ◦C, and also equivalent to a Van 't Hoff function, with a Q10 coefficient of 3.15 from 0–35 ◦C from 0–25 ◦C. The effect of soil water content on mineralization is described by a linear function. Mineralization in temperate soils peaks when soil water content equals field capacity and stops when the ratio of soil water content to field capacity is less than 0.3 [38]. The Mn estimated from the N mass balance for each year and field was then divided by the normalized time for each year, which gives a daily "normalized" rate of mineralization Vn (kg N ha−<sup>1</sup> nday<sup>−</sup>1). Mean Vn (Vnmean) was calculated for each field by averaging Vn for all three years of measurements.

## *2.5. Data Screening*

Data were screened to exclude the fields in which agronomic or measurement problems had occurred. They were represented in a decision tree with three nodes (Figure 3): (i) the presence of weeds, which can compete with maize and bias estimates of Mn; (ii) non-homogeneous maize cover at harvest, which can bias measurements of Nf and N uptake and, thus, estimates of Mn; and (iii) excessively high variability in Mn among the three subplots. We considered, as outliers, fields whose coefficient of variation of Mn was greater than or equal to 25% (i.e., the 95th percentile) in 2012, 2013, or 2014. To retain the most consistent data and, thus, decrease uncertainty in estimates of Mn, we excluded these fields from further analysis and modeling.

**Figure 3.** Decision tree used to filter the data from each field. Conditions had to be true for 2012, 2013, and 2014 to pass to the next step. CV: coefficient of variation, Mn: soil net nitrogen mineralization.

#### *2.6. Soil and Plant Analysis*

Initial soil mineral N content (Ni) was measured at the end of winter (March), and final soil mineral N content (Nf) was measured at the beginning of autumn (October), before resumption of nitrate leaching. A composite sample of 10 soil cores was created for each subplot, for three layers (0–30, 30–60 and 60–90 cm). Soil mineral N was extracted in a 1M KCl solution using a soil/KCl ratio of 1:2, and the NH4 <sup>+</sup> and NO3 − contents of the soil extracts were then determined by continuous flow colorimetry by the methods developed by [42] and [43], respectively.

The soil of the upper layer was sampled in March 2013 to estimate soil properties. Samples consisted of 10 cores, which were pooled. Total carbon (C) and N were determined by the Dumas dry-combustion method. Cation exchange capacity was established using the Metson method [44], and pH was obtained in water [45]. Soil texture was based on measuring the particle size of five fractions: clay (<2 μm), fine silt (2–20 μm), coarse silt (20–50 μm), fine sand (50–200 μm), and coarse sand (200–2000 μm) [46]. POM was determined using a simple fractionation method, by wet sieving under water with a 50 mm sieve [16]. POM was recovered on the sieve, dried, weighed, and finely ground before analyzing the C and N contents by the Dumas dry-combustion method using a Thermo Finnigan Flash EA 1112 Series analyzer.

The fumigation-extraction method was used to estimate soil microbial biomass (SMB) [47], using 40 g of an oven-dried equivalent of soil, shaken in 200 mL of 0.025M K2SO4 for 45 min. Non-fumigated soils were also extracted in the same way. Oxidizable C in the fumigated and non-fumigated K2SO4 was determined using a V-WS SHIMADZU TOC analyzer. A kEC conversion factor of 0.38 was applied to convert the flush of oxidizable C into SMB. EON was determined using the method of [48]: 4 g of soil were steam distilled in 40 mL of phosphate borate extractant (buffered at pH 11.2) for 8 min, and ammonium in the distillate was back titrated using 0.0025 M H2SO4. The amount of organic N hydrolyzed, which corresponded to EON, was obtained by subtracting the total NH4-N extracted from native NH4-N. Mean soil properties, POM-N, SMB, and EON are summarized in Table 1.

**Table 1.** Measured mean (±1 standard deviation) soil variables for the 137 fields of the network and the 67 fields selected. Physico-chemical properties of the 0–30 cm soil layer include texture, C and N contents, pH, Metson cation exchange capacity (CEC), soil microbial biomass (SMB), extractable organic N (EON), particulate organic matter (POM-N), and soil organic nitrogen (SON) stocks.


Bulk density was measured once, in triplicate, in 2011, for each experimental field and each of the three soil layers (0–30, 30–60 and 60–90 cm), using an 8 cm diameter root auger, which cored undisturbed samples of a known volume. The bulk density of the fine-earth fraction calculated from the dry mass and the core volume was used to convert the mineral N content of the samples to kg N ha−<sup>1</sup> and to calculate the stocks of SON and POM-N (t N ha<sup>−</sup>1), as well as those of EON and SMB (kg N ha−1).

Aboveground biomass and N content of the maize crop were quantified at harvest, when maize plants were harvested and weighed in all subplots. Total N uptake of maize was calculated by multiplying N in the aboveground biomass by 1.15 to estimate its belowground N at harvest [49].

### *2.7. Calculating an Indicator of the Cropping System*

An indicator of the cropping system (I\_Sys) was calculated to integrate the diversity of field management (i.e., crop rotation and organic waste application) among fields, considering a period of 15 years before the year of interest. I\_Sys was calculated by summing an indicator of the effect of the N returned to the soil in crop residues and an indicator of the effect of repeated applications of organic waste on soil mineralization. See [40] for details on calculation of the indicator. To assess the influence of I\_Sys on the N mass balance and soil N mineralization, we classified its values into three levels using k-means clustering: low (≤63 kg N ha−1), moderate (63–98 kg N ha−1), and high (>98 kg N ha<sup>−</sup>1). These three classes contained 32%, 53%, and 15% of the fields, respectively.

#### *2.8. Statistical Analysis*

All statistical analyses were performed with R, v 4.1.0 [50]. Pearson correlations were used to assess relations between variables. Analysis of variance (ANOVA) was used to compare means when data were normally distributed (according to the Shapiro–Wilk test) and had homogenous variances (Levene's test); if not, the Kruskal–Wallis test was used.

Generalized additive models (GAMs) were used to predict Vn [51]. This innovative method consists of selecting soil properties for a given GAM, which can best explain the variation in N dynamics [51,52]. GAMs can distinguish the relative effects of covariates by allowing for nonlinear relations between them and the variable studied (i.e., Vnmean). Variables were chosen to obtain models with the smallest mean square error of prediction (MSEP) [53]. This criterion was calculated by applying a "leave-one-out" strategy, which represented an internal validation of the model and avoided over-fitting the model to the data. The agreement between predictions and observed Vnmean was evaluated statistically by calculating the coefficient of determination (R2) and the root mean squared error (RMSE) [53]. Additionally, the ratio of performance to inter-quartile distance (RPIQ) was calculated as the ratio of the inter-quartile range to the square root of the MSEP [51,54]. RPIQ represents the degree to which the dispersion of the response variable exceeds the model's prediction error. We also calculated the indicator dMSEP, which represents the proportional increase in prediction error (MSEP) when a given variable is removed from the model, as a metric of the relative importance of the variable.

#### **3. Results**
