Study site

Field sample collections were performed in Dianzi (118.29 °E, 37.06 °N, 13 m Altitude), Boxing County of Shandong Province, China. The region is characterized by a temperate continental monsoon climate. The annual mean rainfall and air temperature were 601 mm and 12.5 ℃, respectively. The precipitation during winter wheat growth season only accounts for 20–30% of the annual precipitation, which is far less than the demand of wheat grows (Si et al., 2020). The mean rainfall during 2017–2020 winter wheat seasons were 253.5, 95.5 and 162.7 mm, respectively. The mean air temperature during 2017–2020 winter wheat seasons were 8.2, 8.4 and 8.8℃, respectively. The basic information of soil properties are showed in Table 1. The mean soil available nitrogen (N), phosphorous (P), potassium (K) and soil organic carbon (SOC) were 59.69, 10.67, 127.57 mg kg− 1, and 1.33%, respectively.

Table 1

Basic information of soil properties in 0-100 cm layer

Soil depth (cm) | Clay (%) | Silt (%) | Sand (%) | Field capacity (g g− 1) | Saturated water content (g g− 1) | Bulk density (g cm− 3) |

0–20 | 5.25 | 78.55 | 16.20 | 0.25 | 0.32 | 1.41 |

20–40 | 6.33 | 77.74 | 15.94 | 0.21 | 0.26 | 1.52 |

40–60 | 7.63 | 76.89 | 15.48 | 0.21 | 0.24 | 1.54 |

60–80 | 8.23 | 76.92 | 14.85 | 0.21 | 0.25 | 1.53 |

80–100 | 8.79 | 78.90 | 12.31 | 0.21 | 0.24 | 1.52 |

The winter wheat variety in 2017–2018 season was Jimai-22, and the variety of Jimai-23 was used in 2018–2020 seasons due to variety alternative. The wheat was sown in early October and harvested in next early June. Irrigation of 90 mm was applied on April 2, 2018, and on March 28 and May 10, 2019, and on March 20 and May 7, 2020 in the three wheat seasons, respectively. Samples were collected in three plots (9 m × 20 m). The base fertilizer rates of N: P2O5: K2O were 120-120-120 kg ha− 1, and the topdressing N fertilizer 120 kg ha− 1 was applied at early jointing stage.

Sample collection and isotopic analyses

Wheat stem and soil samples were collected for stable isotopes measurement at the different growth stages (jointing, heading, flowering, filling, and harvest). In addition, at 3 days after both irrigation and rainfall and before irrigation, the xylem and soil samples were collected with three replicas. The wheat xylems between the soil surface and the first node about 2–3 cm in length were collected. Around the sampled wheat plants, soil cores from 0 to 150 cm depth (at 0 ~ 5, 5 ~ 10, 10 ~ 20, 20 ~ 30, 30 ~ 40, 40 ~ 60, 60 ~ 80, 80 ~ 100, 100 ~ 120 and 120 ~ 150 cm) were extracted using a soil auger in the three plots. All soil and stem samples were sealed with parafilm and kept stored at -15 ℃ before isotopic measurements. Precipitation or irrigation water were collected at 4 ℃ in each rainfall or irrigation event using a 50 ml polyethylene bottle and sealed with parafilm.

Soil water content (SWC, cm3 cm− 3) was measured by the Insentek sensors (Beijing Oriental Ecological Technology Ltd., Co., Beijing, China) at 10 cm increment in 0–60 cm depth, 20 cm increment in 60–120 cm depth and 30 cm increment in 120–150 cm depth. Accuracy and reliability of the sensors for measuring SWC were assessed in the previous study (Qin et al., 2021). Further, we also use the oven-dry method to correct the SWC measured by sensors. Based on the soil water balance, the daily crop evapotranspiration (ETtotal, mm) and the daily changes in soil water storge of each soil layer (ETeach, mm) can be obtained using Insentek sensors.

A cryogenic vacuum distillation system (LI-2000, LICA United Technology Limited, Beijing, China) was used to extract water from stem and soil samples. To avoid evaporation fractionation, all water samples were stored in a refrigerator at 4 ℃ before isotopic analysis. The isotopic measurements for D and 18O content were conducted using a Continuous-Flow Isotope-Ratio Mass Spectrometer (Delta V Advantage, Thermo Fisher Scientific, Bremen, Germany). All water samples were calibrated against the Vienna Standard Mean Ocean Water standard (V-SMOW) and then convert to δD and δ18O, expressed as \({\delta }\left(\text{‰}\right)=\frac{{\text{R}}_{\text{s}\text{a}\text{m}\text{p}\text{l}\text{e}}-{\text{R}}_{\text{s}\text{tan}\text{d}\text{a}\text{r}\text{d}}}{{\text{R}}_{\text{s}\text{tan}\text{d}\text{a}\text{r}\text{d}}}\times 1000=\left(\frac{{\text{R}}_{\text{s}\text{a}\text{m}\text{p}\text{l}\text{e}}}{{\text{R}}_{\text{s}\text{tan}\text{d}\text{a}\text{r}\text{d}}}-1\right)\times 1000\).

Model description

Three Bayesian mixing models, including MixSIAR, MixSIR and SIAR, were used to determine plant RWU proportion. There is no isotope fractionation occurred during water transfer from the soil to the roots and through xylem vessels, thus the fractionation values were set to 0 in the Bayesian models (Donovan and Ehleringer, 1994). The average and standard deviation value of the stem and each source isotopic compositions (δD and δ18O) were input into the three Bayesian models. Meanwhile, these Bayesian models output a source water’s most likely contribution to the stem water (Evaristo et al., 2017; Stock et al., 2018; Wang et al., 2019). For SIAR, the model was run with 500,000 iterations (discarding the first 50,000) using Markov chain Monte Carlo (MCMC) built on R software. The MixSIR used sampling importance resampling (SIR) algorithm, and the model was run with 500,000 iterations. The run length of the MCMC was set to ‘long’ and the error option was set to ‘residual only’ in MixSIAR. The more detailed information of three models can be referred to literatures (Moore and Semmens, 2008; Parnell et al., 2013; Stock et al., 2018; Wang et al., 2019).

Data preparation and model performance assessment

In this study, three different soil moisture conditions, i.e. normal, dry and wet condition, were used to evaluate the performance of different models.

Normal condition: the soil moisture ranged between 70% and 90% of field capacity without soil water movement by gravity in the 0–100 cm soil column (i.e. one week after irrigation).

Dry condition: the soil moisture ranged between 40% and 60% of field capacity due to the consumption of crop development.

Wet condition: the soil moisture ranged between 80% and 100% of field capacity and the gravity water flows within the 0–100 cm soil column after irrigation (within a week after irrigation).

Although we can currently analyze the process of water absorption by roots through models and experiments, it is still difficult to obtain the true value of the proportion of water contributed from each source to plant. Therefore, we indirectly evaluate the accuracy of different models by evaluating the matching degree between observed and predicted values of isotopic compositions of stem water. Based on the theory that the stem isotopic compositions (δD and δ18O) can be considered as the mixture of different water sources and thus the calculation of predicted value (Pi) was expressed as follows:

$${\text{P}}_{\text{i}}=\sum _{\text{i}=1}^{\text{n}}{\text{f}}_{\text{i}}{{\delta }}_{\text{i}}$$

4

where \({\text{f}}_{\text{i}}\) represents the proportion of water contributed from each source i (n = 4 in this study), and \({{\delta }}_{\text{i}}\) indicates the measured isotopic composition of the corresponding water source i.

Model performance indicators such as the coefficient of determination (R2), the root mean square error (RMSE), the index of agreement (IA), the Nash-Sutcliffe efficiency coefficient (NS) were adopted to evaluate the model prediction accuracy by comparing the matching degree between observed and predicted values of stem water isotopic compositions.

The equations used for RMSE was:

$$\text{R}\text{M}\text{S}\text{E}=\sqrt{\frac{1}{\text{n}}\sum _{\text{i}=1}^{\text{n}}{\left({\text{P}}_{\text{i}}-{\text{O}}_{\text{i}}\right)}^{2}}$$

5

where n is the number of validation samples. Regarding the value of RMSE, lower values implied the less error and accuracy of the models.

The equations used for R2 was:

$${\text{R}}^{2}={\frac{\sum _{\text{i}=1}^{\text{n}}({\text{P}}_{\text{i}}-{\text{P}}_{\text{a}\text{v}\text{e}})({\text{O}}_{\text{i}}-{\text{O}}_{\text{a}\text{v}\text{e}})}{\sqrt{{\sum }_{\text{i}=1}^{\text{n}}{({\text{P}}_{\text{i}}-{\text{P}}_{\text{a}\text{v}\text{e}})}^{2}}\sqrt{{\sum }_{\text{i}=1}^{\text{n}}{({\text{O}}_{\text{i}}-{\text{O}}_{\text{a}\text{v}\text{e}})}^{2}}}}^{2}$$

6

where the Pave and Oave are the mean of the predicted and observed values, respectively. The closer the R2 to 1, the more perfect agreement between Pi and Oi value.

The equations used for IA and NS was:

$$\text{I}\text{A}=1-\frac{\sum _{\text{i}=1}^{\text{n}}|{\text{P}}_{\text{i}}-{\text{O}}_{\text{i}}|}{\sum _{\text{i}=1}^{\text{n}}{\left(\right|{\text{O}}_{\text{i}}-\text{O}}_{\text{a}\text{v}\text{e}}|+\left|{\text{P}}_{\text{i}}-{\text{O}}_{\text{a}\text{v}\text{e}}\right|)}$$

7

The value of IA, based on a calculation of dispersion, ranged from 0 to 1. An IA value of 1 meant more awesome prediction in the models. An IA value close to 0 indicated a worst agreement between Pi and Oi value.

The equations used for NS was:

$$\text{N}\text{S}=1-\sum _{\text{i}=1}^{\text{n}}{({\text{P}}_{\text{i}}-{\text{O}}_{\text{i}})}^{2}/\sum _{\text{i}=1}^{\text{n}}{({\text{O}}_{\text{i}}-{\text{O}}_{\text{a}\text{v}\text{e}\text{r}\text{a}\text{g}\text{e}})}^{2}$$

8

The value of NS also measured the effectiveness of the model. It was reported the simulation results were regarded as perfect prediction in water-limited ecosystems if the NS value is close to 1 and as less reliable if the NS value is negative (Wang et al., 2019). In addition, we used field data from our previous study at the same site to verify which model predicts best for water source apportionment.

## Statistical analysis

According to the similarities of SWC, soil layer classifications was analyzed using a hierarchical cluster analysis (HCA) method described in Zhao et al., 2018. The statistical analysis was performed in SPSS 25.0 (SPSS Inc., Chicago, IL, USA). One-way analysis of variance (ANOVA) was used to compare the differences in RWU source partitioning by the different models at α = 0.05. All figures were drawn using Origin 2016 (Origin Lab, Northampton, MA, USA).