Soil Dynamics and Crop Yield Modeling Using the MONICA Crop Simulation Model and Time Series Forecasting Methods

Abstract

Crop simulation models are an important tool for assessing agroecosystem performance and the impact of agrotechnologies on soil cover condition. However, the high uncertainty and labor intensiveness of long-term weather forecasting limits the applicability of such models. A possible solution may be to use time series forecasting models (SARIMAX and Prophet) and artificial neural-network-based technologies (Neural Prophet). This work compares the applicability of these methods for modeling soil condition dynamics and agroecosystem performance using the MONICA simulation model for Voronic Chernozems in the Kursk region of Russia. The goal is to determine which weather indicators are most important for the yield forecast and to choose the most appropriate methods for forecasting weather scenarios for agricultural modeling. Crop rotation of soybean and sugar beet was simulated, with agricultural techniques and fertilizer usage considered as factors. We demonstrated the high sensitivity of aboveground biomass production and soil moisture dynamics to daily temperature fluctuations and precipitation during the vegetation period. The dynamics of the leaf area index and nitrate content showed less sensitivity to the daily fluctuations of temperature and precipitation. Among the proposed forecasting methods, both SARIMAX and the Neural Prophet algorithm demonstrated the ability to forecast weather to model the dynamics of crop and soil conditions with the highest degree of approximation to actual observations. For the dynamic of the crop yield of soybean, the SARIMAX model exhibited the most favorable coefficient of determination, $R^2$, while for sugar beet, the Neural Prophet model achieved superior $R^2$ levels of 0.99 and 0.98, respectively.

1. Introduction

Global population growth is driving the intensification of agricultural production and soil use. The increase in the consumption of agricultural products is determined by many factors, such as the availability of resources, water, and fertilizers, and environmental pollution. Global climate change brings even more uncertainty to this process [1,2]. The efficient use of resources and forecasting the impact of various stress factors on agricultural production are important to ensure food security and reduce the negative impact on the environment [3]. One of the approaches to planning economic activity and forecasting the impact on agroecosystems is the use of dynamic simulation models of agroecosystem performance [4,5,6].

Dynamic crop simulation models are used to predict crop yields and estimate potential yields in different regions [7] and are also implemented in the decision support systems for farmers [8]. Modeling requires the data on soil cover, environmental conditions, weather observations, and crop management practices including sowing and harvesting dates, and the types and amounts of fertilizers used. The most commonly used crop simulation models such as DSSAT [9], APSIM [10], STICS [11], AgroTool [12], and MONICA [13] include soil, climate, and crop properties modules.

The description of soil processes included in most simulation models allows them to be used to predict soil conditions and the dynamics of soil components [14]. In agroecosystem modeling tasks, it is of interest to predict the dynamics of nutrient elements and evaluate their loss as a result of washing out from the soil profile [15,16]. Most models analyze soil profiles by separate horizons, between which substance and energy exchange takes place [17]. The models mostly take into account heat exchange and moisture transfer, as well as C and N cycles [14].

Daily weather observations are necessary to obtain reliable yield forecasts. Many factors limit the possibility of obtaining a long-term weather forecast with daily resolution; for this reason, various methods of weather scenario generation are used in the models. Using a forecast of global climate models as input for crop process-based models is problematic for two reasons: the low spatial resolutions of such models and the notable fluctuations and deviations in crucial processes for crop modeling at the daily time step, such as precipitation, droughts, and others [18,19]. Some of the most commonly used approaches to tackle this limitations are stochastic weather generators [20] and statistical time series forecasting methods [21].

The recent progress in the development of artificial intelligence technologies has made it possible to apply deep neural networks for time series prediction, including weather data [22,23]. The convolutional neural network was used to downscale the climate model scenario from the CMIP6 [24] project, and it can be helpful for agricultural modeling tasks [25]. Machine learning methods are also used for seasonal forecasting, and it is worth noting that the authors have attempted to use interpreted modeling approaches [26,27]. The hybrid use of data-driven and physical-based weather and climate forecast models has become popular in the last decade [28,29].

While stochastic weather generators are quite often used in agroecosystem modeling, machine learning and neural networks have hardly ever been used for weather forecasting in the simulation models. In this paper, we analyze the possibility of using the crop simulation models together with the advanced time series forecasting methods. We also examined how sensitive crop yield models are to the quality of weather forecasts and estimated which weather components impact crop yield more. We assumed that the quality of the weather forecast affects the quality of the estimates of the crop dynamics in different ways. A location in the Kursk Oblast, Russia, was chosen to perform numerical modeling. Modeling was carried out using the MONICA simulation model for sugar beet (sugar beet, Beta vulgaris) and soybean (soyabean, Glycine max) in 2016 and 2017.

The purpose of this work was to investigate whether it is possible to use time series forecasting methods for the modeling of soil condition dynamics and agroecosystem performance using the MONICA simulation model.

The main contributions of this paper are:

• We demonstrated the performance of time series models to address weather scenario generation for crop modeling;
• We estimated the efficiency of classical and neural network time series methods for crop and soil modeling;
• We implemented a crop yield forecast based on time series models and a crop simulation model.

2. Materials and Methods

2.1. Research Object

An agricultural field located in the Kursk Oblast, Sovetsky District, Kshensky Settlement, was chosen as a pilot research object. The terrain is a rolling plain crossed by ravines with smooth slopes and shapeless depressions. Climatic conditions correspond to the moderately cold climate with the average annual precipitation of 350–570 mm and the average annual ambient temperature of 5.4 °C. The average temperatures are −7.2 °C in January and +20.1 °C in June. The vegetation period lasts for 188 days. The area of the field plot used in modeling was 48 ha (Figure 1).

Soil was classified as Voronic Chernozem according to the World Reference Base for Soil Resources [30]. An uneven four-field crop rotation was used: barley (2012)—winter wheat (2013)—sugar beet (2014)—winter wheat (2015)—soybean (2016)—sugar beet (2017). Soil samples were collected from the depths of 10, 30, and 80 cm at 51 locations within the field boundaries (N 51 54${}^{\prime }$8${}^{\prime }$, E 37 47${}^{\prime }$34${}^{\prime }$) in October 2016. The key soil characteristics used in modeling are summarized in

Table 1. Soil characteristics used as input data in the modeling.

Horizon Depth, cm	Humus, %	Sand Fraction, %	Clay Fraction, %	pH	C:N	Bulk Density, kg/m3
0–30	5.1	37	9	6.2	12.4	1126
30–50	3.3	45	8.8	6.4	10.5	1056
50–70	2.3	49	8.6	6.6	8.8	1010
70–200	1.6	5.3	9.4	6.7	7.5	1113

.

2.2. Weather Data

Prediction of Worldwide Energy Resource (POWER) data, which is a publicly available NASA project, was used as input data for generating weather scenarios and training machine learning models. The POWER system is a common source of meteorological data in agricultural modeling, including crop yield modeling [31,32], crop disease spread modeling [33], field irrigation optimization [34], and other tasks [35,36]. The data are created by the assimilation of satellite observations into a climate model (Goddard Earth Observing System climate model) and cover the entire land surface [37]. NASA POWER data exhibit a high level of concordance between the obtained estimates and observations at individual weather stations. Other advantages of this data source are the high spatial resolution (1° × 1° degrees) and high speed of updating with the lag of only 5–7 days. The observations starting from 1983 are available at daily resolution. The following observations available in POWER were used as input parameters in modeling: minimum temperature (T2M_MIN, °C), maximum temperature (T2M_MAX, °C), average temperature (T2M, °C), precipitation (PRECTOTCORR, mm), incoming solar radiation (ALLSKY_SFC_SW_DWN, W/m²), wind speed (WS2M, m/s), and relative humidity (RH2M, %). The NASA POWER data were obtained for the Kshensky Settlement, Kursk Oblast (N 51 54′8′, E 37 47′34′), for the years 2010–2017 in CSV format.

2.3. Time Series Forecasting Methods

Here we describe the time series methods used to forecast weather scenarios. Conventional methods and methods based on neural networks were used.

As a baseline we used grouping by mean approach. Grouping by mean value is the simplest method of forecasting a time series with consists in averaging weather observations over several years. The data were grouped by date and averaged over 10 years. As a result, the averaged data were fed to the input of the simulation model.

SARIMAX (Seasonal Autoregressive Integrated Moving Average with exogenous regressors) is among the most common methods for analyzing and forecasting time series [38,39]. SARIMAX combines autoregression (AR), the moving average (MA), and integration (I) for time series modeling. The AR model is used to describe the dependence of the current value on past values, the MA model is used to describe the dependence of the current value on the random errors of the past values, and the I model is used to eliminate non-stationarity of the series:

$$AR\left(p\right):y_t={\varphi }_1{y}_{t-1}+{\varphi }_2{y}_{t-2}+\cdots +{\varphi }_p{y}_{t-p}+{\epsilon }_t$$

(1)

where $y_t$—stationary time series with zero mean, ${\varphi }_1,\dots ,{\varphi }_p$—constants $\left({\varphi }_p\ne 0\right),$ ${\epsilon }_t$—Gaussian white noise with zero mean and constant variance ${\sigma }_{\epsilon }^2$.

$$MA\left(q\right):y_t={\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2}+\cdots +{\theta }_q{\epsilon }_{t-q},$$

(2)

where $y_t$—stationary time series with zero mean, ${\theta }_1,\dots ,{\theta }_q$—constants $\left({\theta }_q\ne 0\right)$, and ${\epsilon }_t$—Gaussian white noise with zero mean and constant variance ${\sigma }_{\epsilon }^2$.

In addition, SARIMAX includes the seasonality component, which allows us to take into account the repeating cycles in the data. For this purpose, the seasonality parameters are used, which define the periodicity. The key idea of SARIMAX lies in adding exogenous factors to the ARIMA model. To do this, additional regressors are introduced into the model that can affect the target time series. The general formula of SARIMAX is as follows:

$$Y_t={\beta }_0+\sum _{i=1}^p{\varphi }_i{Y}_{t-i}+\sum _{i=1}^q{\theta }_i{ϵ}_{t-i}+\sum _{i=1}^P{\Phi }_i{Y}_{t-i}^{\left(s\right)}+\sum _{i=1}^Q{\Theta }_i{ϵ}_{t-i}^{\left(s\right)}+X_t^T\beta +{ϵ}_t$$

(3)

where $Y_t$—value at the time point t, ${\beta }_0$—constant, p—order of autoregression, ${\varphi }_i$—autoregression coefficient, q—order of moving average, ${\theta }_i$—moving average coefficient, P—order of seasonal autoregression, ${\Phi }_i$—seasonal autoregression coefficient, Q—order of seasonal moving average, ${\Theta }_i$—seasonal moving average coefficient, $Y_{t-i}^{\left(s\right)}$—value at the $t-i$ time point in season s, ${ϵ}_{t-i}^{\left(s\right)}$—error at the time point $t-i$ in season s, $X_t$—vector of exogenous regressors at the time point t, $\beta$—vector of regression coefficients, and ${ϵ}_t$—error at the time point t.

In this way, SARIMAX takes into account seasonality, autocorrelation, as well as the effects of exogenous variables on the target time series. The maximum likelihood method is used to determine the optimal model parameters.

Prophet is a time series forecasting procedure [40]. It is based on the approach that allows the modeling of various time series with a non-linear trend component, seasonality, and holiday effects.

The model can automatically detect changes in trend and seasonality, as well as take into account additional regressors. In addition, the library can provide model interpretation, which makes it possible to analyze the contribution of each component to the forecast. Prophet is built on the basis of the additive time series model:

$$y\left(t\right)=g\left(t\right)+s\left(t\right)+h\left(t\right)+\epsilon \left(t\right)$$

(4)

where $g\left(t\right)$—trend function, $s\left(t\right)$—seasonality function, $h\left(t\right)$—holiday function, or anomalous data, and $\epsilon \left(t\right)$—error function.

In Prophet, an individual model is fitted for each component:

$$g\left(t\right)=X\left(t\right)\beta +\gamma \left(t\right),s\left(t\right)=\sum _{i=1}^N{a}_{i}cos\left(\frac{2\pi it}{P}\right)+b_{i}sin\left(\frac{2\pi it}{P}\right),h\left(t\right)=\sum _{j=1}^M{k}_{j}I\left(t\in \left[s_j,e_j\right]\right),$$

(5)

where $X\left(t\right)$—regressor matrix, $\beta$—parameter vector, $\gamma$—smoothed seasonality function, N—number of seasonality components, P—season length, $a_i$, $b_i$—seasonality coefficients, M—number of holidays, $k_j$—holiday effect strength, and $s_j$, $e_j$—start and end of the holiday period, respectively.

Neural Prophet is a time series forecasting method based on neural networks [41]. It is an extension of the Prophet method and employs the same idea of time series decomposition into trend, seasonality, and holidays. However, unlike Prophet, Neural Prophet does not rely on the linear models to model trend and seasonality but utilizes neural networks instead.

The model is based on the LSTM (long short-term memory) architecture [42], which allows modeling dependencies in the data based on past values. In the Neural Prophet model, trend is modeled using neural networks, which can manage non-linear dependencies between the time points. Seasonality and anomalies are also modeled based on neural networks, which allows us to capture more complex patterns in the data. The model can also include exogenous variables to account for additional factors affecting the time series.

2.4. Forecast Accuracy Metrics

There are a number of quality indices that are mostly used to evaluate the quality of time series prediction achieved by different models [43]. The following quality indices were used in this work: Mean Squared Error (MSE), Relative Root Mean Squared Error (RRMSE), and R-squared ($R^2$). The quality indices were calculated as follows:

$$MSE=\frac{1}{n}\sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2$$

(6)

$$\mathrm{Relative}RMSE=\frac{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}}{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}},$$

(7)

$$R^2=1-\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{{\sum }_{i=1}^n{(y_i-\overline{y})}^2},$$

(8)

where $y_i$—actual value of the target variable in the observation i, ${\widehat{y}}_i$—predicted value of the target variable in the observation i, $\overline{y}$—mean value of the target variable, and n—total number of observations.

Dynamic time warping (DTW) was also used to estimate forecasting accuracy. DTW measures the distance between two time series, accounting for possible distortion shifts and different data scales [44]. This metric provides a more accurate measure of similarity between the time series than the simple Euclidean distance. The approach is based on finding the optimal correspondence between the points in the two time series so that the resulting distance between them is minimal. DTW between the two time series is estimated as follows:

$$DTW(X,Y)=\sqrt{\sum _{i=1}^n\sum _{j=1}^m{\left(X_i-Y_j\right)}^2},$$

(9)

where X and Y—two time series with the lengths of n and m, respectively.

The DTW ($X,Y$) procedure computes the Euclidean distance between each point in the time series X and Y and at each time step, and then finds the optimal correspondence between the points using dynamic programming.

2.5. Crop Model MONICA

MONICA [13] was chosen as the agroecosystem model. MONICA is a one-dimensional dynamic process-oriented simulation model of plant production, which describes the transport and biogeochemical turn-over of carbon, nitrogen, and water in agroecosystems. Using a dynamic iterative algorithm with a daily time step, MONICA simulates the most important processes in the soil and plant, taking into account direct and feedback relations. MONICA is the successor of the previously developed HERMES and AGROSIM models [45,46]. Schematic representation of the MONICA model is presented in Figure 2.

MONICA significantly extended the functionality of its predecessors by adding a full carbon cycle in the agroecosystem, which is a prerequisite for the study of carbon dynamics under climate change. Towards this end, a simple algorithm implemented in HERMES for estimating nitrogen mineralization from soil organic matter was replaced by a comprehensive approach used in the Danish DAISY model [47], which also takes into account the dynamics of soil microbial biomass. MONICA takes into account the impact of the atmospheric CO₂ concentration on photosynthesis and transpiration in crops, as well as the heat stress and possible frost kill of winter crops. The algorithms for calculating plant growth and development implemented in MONICA have been validated for both temperate and tropical latitudes. These experiments identified parameter sets for most typical crops produced in these regions taking into account varietal diversity [48,49]. A sensitivity analysis of MONICA was carried out concerning crop parameters [50] and soil properties [51].

2.6. Numerical Experiments

NASA POWER weather data for the period from 2010 to 2017 were used to build the predictive models of weather scenarios. Two types of experiments were performed to assess forecast accuracy: in the first case, the data were split into a training sample including the years 2010–2015 and a forecast was obtained for 2016 (Figure 3). In the second case, the data were split into a training sample including the years 2011–2016, and a prediction was made for 2017.

Then, the SARIMAX, Prophet, and Neural Prophet forecasting models were trained using the available sample, with the last year of weather observations being used for model validation and hyperparameter tuning. The predictions made for 2016 were used to model soybean growth and development, while the results obtained for 2017 were used to model sugar beet growth and development. An outline of the experiments and used software for modeling is illustrated in Figure 4.

To simulate plant growth dynamics and soil condition, the MONICA model was used in this work, for which the input data were weather forecasts and the soil parameters obtained from field surveys (Table 1). Numerical experiments were conducted for two crops: soybean in 2016 and sugar beet in 2017. For both crops, ploughing to a depth of 25–27 cm was performed prior to sowing. The description of the main agrotechnological events for crops is presented in

Table 2. Crop rotation and fertilization.

Crop	Tillage (cm)	Fertilization (Date–kgN/ha)	Sowing Date	Harvesting Date
Soyabean	25–27 cm	15.04–50	15 April 2016	9 October 2016
Sugar beet	25–27 cm	23.05–90; 15.06–70	2 May 2017	1 November 2017

.

The following MONICA model variables were used as target parameters for the studied crops: yield (kg/ha), aboveground biomass (kg/ha), ecosystem carbon exchange (NEE, kgC/ha), and leaf area index (LAI, –). Organic matter carbon dynamics (kgC/ha), nitrate leaching rate (kgN/ha), ${\mathrm{NO}}_3$ dynamics (kgN/ha), and soil moisture were considered to compare soil condition dynamics.

To assess the quality of the obtained time series of plant growth dynamics and soil condition, the abovementioned quality indices were used, and the forecasts obtained using MONICA based on reanalysis data and experimental data were compared.

3. Results and Discussion

3.1. Forecasting Weather Scenarios

The temperature modeling results are shown in Figure 5. As can be seen on the plot, none of the proposed models can fully describe the weather data. The Prophet and Neural Prophet models smooth the data, resulting in no sharp changes in temperature compared to the real data. If we look at the observations from July to August 2017, it can be seen that the temperature values range from 13 to 26 °C, while the values predicted by the models range from 17 to 20 °C. Most models performed well at capturing the general trend in the data but were rather poor in reproducing the stochasticity of weather events. On the other hand, the grouping by mean and SARIMAX models represented stochasticity; however, the predicted temperature values differed significantly from the real data.

In the case of precipitation modeling, it should be noted that the observations were essentially stochastic and their prediction was nearly impossible. Prophet and Neural Prophet gave smoother results as in the case with temperature but missed the extreme events such as heavy rainfall in June. The grouping by mean model and SARIMAX retained stochasticity, but the peaks predicted by them did not exceed 10 mm/day, while in the real observations, the peaks reached 30 mm/day.

The weather conditions forecast accuracy indices are summarized in

Table 3. Weather variables forecast metrics.

	Average Temperature			Precipitation			Solar Radiation			Wind
	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW
Baseline	0.42	0.82	47.73	1.05	–0.13	64.77	0.52	0.72	50.51	1.00	–0.01	16.67
SARIMAX	0.53	0.71	49.84	1.04	–0.10	68.02	0.57	0.66	46.72	1.06	–0.14	15.88
Prophet	0.42	0.81	61.20	0.99	–0.11	80.34	0.49	0.75	66.57	0.92	0.13	21.94
Neural Prophet	0.41	0.82	61.50	0.99	–0.05	80.06	0.50	0.74	66.95	0.91	0.10	22.64

Bold indicates best model performance.

. The lowest relative error was achieved by the Prophet class models. However, even in their case, the RRMSE for precipitation and wind speed predictions was up to 99%. The coefficient of determination for precipitation and wind speed forecasts was negative for all models. A negative value can be obtained when the model fits the data worse than a simple constant model, or when the predicted dependent variable values deviate strongly from the actual values in the sample. The best values for the other variables were also obtained by the Prophet class models. When using the dynamic time warping algorithm to predict average temperature and precipitation, grouping by mean showed the best results, while in the case of solar radiation and wind speed, the best results were achieved with SARIMAX.

Based on the obtained results, it can be concluded that the Prophet and Neural Prophet models are best suited for predicting temperature and solar radiation, as they have the lowest relative error. However, when forecasting precipitation and wind speed, the grouping by mean and SARIMAX models are more suitable as they retain the stochasticity of the data. In general, it is recommended to use Neural Prophet to achieve the best results in time series forecasting, as it demonstrates the highest quality of prediction combined with the good performance for all the variables considered.

3.2. Crop Growth and Development Modeling

The simulations resulted in the growth and development scenarios for soybean and sugar beet for the different forecasts obtained with the machine-learning-based models. The results for each crop were compared considering the dynamic changes in the main crop growth indicators at one-day intervals. Additionally, quality indices were calculated based on the correspondence between the dynamics of the biophysical crop growth indices obtained using different forecasting methods and based on the real weather observations. The growth curves and dynamics of the key biophysical parameters for soybean are provided in Figure 6. The yield dynamic differs considerably between the forecasts obtained with different models. As can be seen in Figure 6A, Neural Prophet and the grouping by mean model both forecast earlier crop maturation than is predicted based on the observed weather conditions. Prophet underpredicts the total yield by about two times. SARIMAX gives the dynamics that are the closest to the real observations for the whole vegetative period. The aboveground biomass dynamics are similar to the yield dynamics. In this case, SARIMAX also produces the dynamics that are the closest to those observed, and also allows values to be obtained that are the closest to the observed estimate of the aboveground biomass at the end of the season.

Ecosystem carbon exchange shows the lowest values in the period from June to August, which is associated with plant biomass growth. Prophet and Neural Prophet give smoother NEE dynamics compared to the actual observations. SARIMAX produces the dynamics that are closer to those observed than all other forecasts, including at the end of the vegetative period when all other models overestimate NEE.

Leaf area index is an important characteristic used to assess crop growth and development during a season. This index is a frequently used variable in the tasks of remote sensing data assimilation into simulation models with the purpose of improving their predictive properties. Therefore, the ability of the model to predict its values during the season is important in assessing its quality. As can be seen on the plot, SARIMAX allows the estimation of LAI dynamics during the season with high accuracy, while the other models have the LAI peak shifted towards the start of the season.

The total soybean yield was 3758 kg/ha in the case of the real weather observations. The closest yield estimate, 3184 kg/ha, was obtained when weather observations were averaged over 10 years, while the worst forecast, 1964 kg/ha, was produced by Prophet. All forecast variants underestimated the final yield. One of the possible reasons for this phenomenon may be the poor quality of precipitation forecasting by the models, which ultimately leads to the underestimation of precipitation amounts over the season. On the other hand, in the grouping by mean model, the total precipitation was 275 mm, which is almost two times lower than the observed value (504 mm), while the predicted total yield is only lower by 5%.

The forecast accuracy indices are provided in

Table 4. Soybean dynamics forecast quality indices.

	Crop Yield			Biomass			LAI			NEE
	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW
Baseline	0.44	0.80	403.71	0.19	0.96	3187.21	0.32	0.89	7.68	0.51	0.73	326.01
SARIMAX	0.12	0.98	969.97	0.04	0.99	1682.33	0.04	0.99	0.58	0.35	0.88	184.21
Prophet	0.49	0.75	4194.16	0.27	0.92	11,062.58	0.41	0.83	4.47	0.9	0.19	280.57
Neural Prophet	0.74	0.40	503.84	0.27	0.92	3435.52	0.43	0.81	8.87	0.64	0.58	280.47

Bold indicates best model performance.

. In the case of the yield prediction, the grouping by mean variant obtained the DTW score of 403.71. However, the SARIMAX model shows the $R^2$ of 0.98, which is higher than 0.80 for the averaging variant. SARIMAX also obtained the highest quality scores in the case of the aboveground biomass dynamics, LAI, and NEE. It is also worth noting that all models were worse at predicting NEE dynamics than all the other parameters.

The growth curves and dynamics of the key biophysical parameters for sugar beet are provided in Figure 7. In the case of sugar beet, unlike soybean, the models overestimate yields compared to the observations. The grouping by mean variant gives earlier crop maturation and a significant overestimation of yield. The Prophet, Neural Prophet, and SARIMAX models also overestimate yields. In contrast to yield, the aboveground biomass dynamics are better estimated by the different models. It is worth noting that unlike yield, the aboveground biomass values are underestimated by SARIMAX. It is also worth mentioning that, unlike soybean, sugar beet yields show less correlation with the aboveground biomass.

The dynamics of ecosystem carbon exchange in the case of real-life observations exhibit a rather high amplitude over the vegetation season. Prophet and Neural Prophet give smoother NEE dynamics compared to the real observations, which was the same for soybean. In turn, the SARIMAX model gives an amplitude that is closer to the observed NEE dynamics.

In the case of leaf area index, the following pattern was observed: all variants produced a forecast that was biased towards the start of the season. In the second half of the season, the Prophet and Neural Prophet models predict the LAI dynamics that replicate the observed dynamics quite well. The SARIMAX model, in its turn, underestimates the LAI values in the second half of the season by about 0.5 units. The Neural Prophet model gives the closest value to the observed estimate of yield with the total value of 24,187 kg/ha, and the worst forecast is produced by the weather averaging approach with the value of 29,395 kg/ha, which is 32% higher than the observed yield. In contrast to soybean, for sugar beet, all forecasts overestimated the total yield. We should also note a rather large variation in the totals of the precipitation values obtained by different models, starting from 291 mm in the observations to 556 mm for the Prophet model.

The forecast accuracy indices for sugar beet are provided in

Table 5. Sugar beet dynamics forecast quality indices.

	Crop Yield			Biomass			LAI			NEE
	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW
Baseline	0.50	0.72	44,301	0.34	0.88	1786	0.35	0.87	0.59	0.18	0.57	284
SARIMAX	0.27	0.92	14,986	0.24	0.94	9664	0.27	0.92	4.40	0.05	0.71	241
Prophet	0.29	0.91	20,650	0.25	0.93	4876	0.25	0.93	1.07	0.53	0.72	342
Neural Prophet	0.22	0.95	11,019	0.22	0.95	4313	0.22	0.4	1.02	0.50	0.75	343

Bold indicates best model performance.

. The Neural Prophet model revealed the highest quality index values in the case of the yield forecast, in particular, $R^2$ = 0.95. The Neural Prophet model also obtained the best $R^2$ values for all other biophysical indicators compared to other models.

3.3. Soil Dynamics Modeling

The soil conditions in the plough horizon obtained in the soybean simulations are presented in Figure 8. As can be seen on the plot, most models gave poor quality weather forecasts when reproducing soil moisture. Prophet and Neural Prophet underpredicted soil moisture compared to the real weather observations. Using SARIMAX and the grouping by mean model allowed soil moisture values to be obtained that were closer to the actual observation, although a certain underestimation can also be seen. In the summer season, a sharp increase in soil humidity can be seen in the case of the real-life observations with peaks exceeding 0.35, while the predicted peak values in the corresponding period are not higher than about 0.3.

The nitrate leaching dynamics reveal a pattern similar to soil humidity. None of the weather forecasting methods obtained the same high leaching values as in the real observations. All the models and grouping by mean produce similar dynamics and underestimate the maximum leaching values, from 1.0–1.2 to 0.3–0.4 kgN/ha. It is also worth noting that the Prophet and Neural Prophet models produce a peak of nitrogen leaching at the start of the season, while two peaks can be distinguished in the real weather observations, one of them occurring in early May and the other in early June. This may be due to the presence of high precipitation, as can be seen in the Figure A1. SARIMAX predicts the first peak in early May, similar to the real observations, but the second peak in early June is not predicted.

The amount of precipitation and its distribution over the season play a significant role in forming the dynamics of soil moisture and nitrate leaching. As mentioned above, precipitation is one of the parameters that are challenging to predict because of its stochasticity, which in turn impacts the possibility of modeling those soil conditions for which precipitation is the main affecting factor. For this reason, smoothed precipitation predictions underestimate the moisture levels and may lead to the underprediction of nitrate leaching.

The nitrate content dynamics do not show such pronounced fluctuations during the season as moisture and nitrate leaching. For both model forecasts and real weather observations, the nitrate content level appears to be very low throughout the season and increases only after crop harvesting in early September. The closest to the observed dynamics was achieved with the SARIMAX model. As can be seen on the plot, grouping by mean and the Neural Prophet model overestimated the nitrate content at the end of the season, while the Prophet model gave nitrate content levels that were lower than the observed ones. The organic matter content showed a dynamic that is similar to the nitrate content. The plot shows an increase in the organic matter content after harvesting, which is associated with the residual plant parts left in the soil horizon. Similar to the case with the nitrates dynamic, SARIMAX allowed us to obtain a forecast closely replicating the real-life observations.

Table 6. Soil parameters dynamics forecast quality indices for soybean.

	Soil Organic Carbon			Soil Moisture			Nitrate Leaching			Nitrate Content
	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW
Baseline	0.18	0.96	0.003	0.69	0.46	1.590	0.92	0.07	2.27	0.82	0.32	0.180
SARIMAX	0.05	0.99	0.001	0.55	0.68	0.510	0.86	0.19	2.24	0.25	0.93	0.030
Prophet	0.28	0.92	0.004	0.71	0.31	1.750	0.95	0.00	2.35	0.73	0.41	0.170
Neural Prophet	0.13	0.98	0.003	0.67	0.46	1.710	0.92	0.08	2.12	0.76	0.40	0.180

Bold indicates best model performance.

compares the quality indices for the soil condition forecasting models in the case of soybean. As the table shows, the RRMSE and $R^2$ for soil humidity and nitrate leaching are worse than those for the organic matter and nitrate content dynamics, which can be also seen in Figure 8. The same tendency was observed for all the weather scenario modeling methods in this study. When the models were compared to each other, the SARIMAX model achieved the best quality metrics values for all soil parameters.

The soil conditions in the plough horizon obtained in the sugar beet simulations are presented in Figure 9. For sugar beet soil moisture, significant differences between the observations and model forecasts were also obtained.

However, unlike the case with soybean, the sugar beet models tended to overestimate soil moisture levels. The Prophet and Neural Prophet models allowed smoothed moisture dynamics with values of about 0.3 to be obtained, while the real observations showed more pronounced fluctuations and ranged between 0.1 and 0.3 on average. The SARIMAX model and the grouping by mean variant showed dynamics close to the observed one. Nitrate leaching shows a significant increase at the start of the season and then a gradual decrease, which is associated with the use of nitrogen fertilizers, which are taken into account in the modeling. The grouping by mean variant produced a dynamic that was closest to the real-life observations. Prophet and Neural Prophet showed some anomalous increases at the end of the season, which may be mentioned as shortcomings of these forecasts. The nitrate content dynamics showed a sharp peak in the first half of the season after nitrogen fertilizers were applied. The SARIMAX and the grouping by mean models produced forecasts that were quite close to the real content, while the Prophet and Neural Prophet models underestimated the content by about two times. The dynamics of the organic matter content are similar to the soybean case and show a growth after crop harvesting resulting from the plant residues remaining in the plough horizon.

The best quality soil condition forecasts for sugar beet were obtained using weather averaging over several years (

Table 7. Soil parameters dynamics forecast quality indices for sugar beet.

	Soil Organic Carbon			Soil Moisture			Nitrate Leaching			Nitrate Content
	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW	RRMSE	R2	DTW
Baseline	0.17	0.97	0.0007	0.54	0.7	0.52	0.79	0.34	0.26	0.29	0.91	0.06
SARIMAX	0.24	0.94	0.0016	0.71	0.5	0.75	1.29	–0.69	0.29	0.35	0.87	0.08
Prophet	0.16	0.97	0.0007	0.74	0.4	1.82	1.8	–25.25	4.31	0.81	0.25	0.25
Neural Prophet	0.13	0.98	0.0006	0.61	0.57	1.38	2.02	–9.56	1.65	0.86	0.16	0.28

Bold indicates best model performance.

). Neural Prophet obtained the best quality indices for the soil organic matter forecast only, although the differences between this model and the averaging variant were not significant. We should also note the poor quality of the nitrogen leaching forecasts obtained with the models, with only the grouping by mean variant obtaining a positive $R^2$ value equal to 0.34. It can be assumed that the quality of the soil dynamics prediction depends on the quality of the crop dynamics prediction.

As part of our research, we generalized machine learning and statistical analysis methods for working with climate data. We identified which variables can be accurately predicted using these methods and which contain uncertainty that cannot be predicted. In addition, based on the identified results, it is possible to conclude the contribution of individual variables to subsequent modeling using simulation methods.

We hope that our research might be helpful for the development of methods and techniques for the cheap generation of surrogate weather scenarios and risk estimations. However, such research might require an additional scale of computations related to one more mathematical model for risks. Nevertheless, our research shows the advantages and limitations of popular time series forecasting methods for different weather and soil parameters. We see that some parameters can be predicted accurately with our approach, and some require another approach to be found.

One of the main limitations of most seasonal forecast approaches is the complexity of modeling extreme events. This drawback has a powerful effect on the ability to predict precipitation [52]. However, this is important for modeling nitrogen leaching for soybeans. As can be seen, for observed weather, there are several peaks of nitrogen leaching due to extreme precipitation. At the same time, only one peak is observed for the model variants since the models did not predict extreme precipitation. Thus, there may be an underestimation of essential factors, such as nitrogen leaching, which directly depend on the possibility of predicting extreme events. It is important to note that within the framework of the presented study, meteorological parameters were evaluated as independent time series data. In subsequent studies, it makes sense to consider the possibility of conducting a correlation analysis between the abovementioned variables.

Future studies could explore the reasons behind the underestimation of crop yield, explicitly investigating the role of poor precipitation forecasting in model predictions. Additionally, using transformer-based methods that are popular in NLP tasks can improve the quality of the forecast and the use of hybrid approaches based on statistical weather generators and machine learning models [53,54].

4. Conclusions

The numerical experiments conducted in the present work demonstrated the possibility of using machine learning methods for time series forecasting to predict weather scenarios, which could be used in the simulation modeling of agrosystems. The SARIMAX and Neural Prophet models allowed soybean and sugar beet yield estimates to be obtained that were rather close to the observed weather conditions. However, the time series forecasting methods used smooth the series of weather observations, which seems less critical in the case of predicting soil nutrient scenarios, but can make modeling soil moisture dynamics quite challenging since dynamics smoothing produced by these models does not allow extreme events (e.g., heavy rainfall) to be taken into account. Currently, the effectiveness of long-term weather forecasts is limited due to the high level of uncertainty. However, as forecasting technology continues to advance and become more accurate, it is expected that agricultural systems will benefit from improved forecasting performance.