Long-Term Optimal Management of Rapeseed Cultivation Simulated with the CROPGRO-Canola Model

Abstract

Rapeseed (Brassica napus L.) is an important oilseed crop grown worldwide with a planting area of 6.57 million ha in China, which accounts for about 20% of the world’s total rapeseed planting area. However, in recent years, the planting area in China has decreased by approximately 12.2% due to the low yield and economic benefits. Thus, to ensure oil security, it is necessary to develop high-efficiency cultivation for rapeseed production. Crop growth models are powerful tools to analyze and optimize the yield composition of crops under certain environmental and management conditions. In this study, the CROPGRO-Canola model was first calibrated and evaluated using the rapeseed planting data of four growing seasons in Wuhan with nine nitrogen fertilizer levels (from 120 to 360 kg ha⁻¹) and five planting densities (from 15 to 75 plants m⁻²). The results indicated that the CROPGRO-Canola model simulated rapeseed growth well under different nitrogen rates and planting densities in China, with a simulation error of 0–3 days for the anthesis and maturity dates and a normalized root mean square error lower than 7.48% for the yield. Furthermore, we optimized the management of rapeseed by calculating the marginal net return under 10 nitrogen rates (from 0 to 360 kg ha⁻¹ at an increasing rate of 40 kg ha⁻¹) and 6 planting densities (from 15 to 90 plant m⁻² at an increasing rate of 15 plant m⁻²) from 1989 to 2019. The results indicated that the long-term optimal nitrogen rate was 120–160 kg N ha⁻¹, and the optimal planting density was 45–75 plants m⁻² under normal fertilizer prices. The optimal nitrogen rate decreased with increasing fertilizer price within a reasonable range. In conclusion, long-term rapeseed management can be optimized based on rapeseed and nitrogen cost using long-term weather records and local soil information.

1. Introduction

Rapeseed (Brassica napus L.) is the third largest oilseed crop in the world and the second largest oilseed crop in China [1]. It is mainly used in the mass production of animal feed and vegetable oil [2]. In 2019, the planting area of rapeseed in China was approximately 6.57 million ha, which accounts for nearly 20% of the world’s total rapeseed planting area [1]. However, the rapeseed planting area in China has decreased by about 11.2% in the past five years relative to that in 2010 (about 7.4 million ha) [3]. Moreover, China’s edible vegetable oil self-sufficiency decreased from 54% in 2010 to 31% in 2021 [3,4]. The domestic supply of edible oil has decreased due to the COVID-19 pandemic, which restricted rural labor and edible oil shipments and affected the rapeseed processing enterprises [5]. Recently, China has increased the number of large-scale rapeseed growers and promoted rapeseed production [5]. However, the yield of on-farm rapeseed is typically one-third less than that under regional demonstration projects, which indicates that farmers need more information about optimum management.

The Yangtze River Basin (90°13′~122°19′ E, 24°27′~35°54′ N) is the main production area of winter rapeseed in China, due to flat terrain and sufficient light and water supply. Moreover, winter rapeseed production has great potential in this region, as most fields are typically fallow in the winter [6,7]. Statistical data in this region show that rapeseed yield is approximately 2000 kg ha⁻¹, which is only 65.4% of that in the United States and 58.1% of that in European countries [1]. Demonstration plots under optimum management in this region yield over 3000 kg ha⁻¹, which is higher than on-farm yields [8]. Nitrogen and planting density management play important roles in regulating the development of rapeseed at both individual and population levels [9,10]. Rapeseed responds positively to optimum nitrogen fertilizer and plant population. For example, high nitrogen application can increase the essential components (such as cellulose and lignin) of the stem to improve the resistance of the plants to lodging [11]. However, excessive nitrogen may increase the lodging risk due to decreases in root neck diameter and root dry biomass [12,13]. In addition, a higher planting density can also cause decreases in root neck diameter and root dry biomass [12]. Therefore, optimization of nitrogen and planting population management can help to reduce the application of nitrogen fertilizer as well as improve the yield and overall economic benefits.

In the post-epidemic era, mechanized production, unmanned farms and smart farms are inseparable from efficient management decisions. Long-term optimal management of nitrogen rate and planting density is critical for maximizing the marginal net return (MNR) [14]. However, the optimization must be performed at different spatial scales, considering multiple factors such as fertilizers and seed cost and management cost. Therefore, precise methodologies or tools are needed to combine these factors to make optimal pre-season decisions. The CERES-Rice model has been used to simulate long-term optimal nitrogen management for rice in Northeast China [15]. The CSM-Barley model has been employed to reduce nitrate leaching while improving grain yield and quality in malting barley [16]. Currently, the models used to simulate rapeseed growth include APSIM-Canola [17], CROPGRO-Canola [18] and AquaCrop-Canola [19]. Among them, the APSIM-Canola model has mainly been used to simulate the effect of climate change and soil on crop growth, but this model is not appropriate for simulating the flowering period of spring rapeseed under different sowing dates [20]. The AquaCrop-Canola model is a water-driven model mainly focused on the influence of water on crop biomass and yield [21]. The CROPGRO-Canola model [18] simulates crop growth and development and is open source. This model has been evaluated under irrigation, rainfed and nitrogen stress conditions, as well as calibrated and evaluated to simulate the yield and biomass of rapeseed in eastern Canada (Jing et al., 2016). This model is part of the Decision Support System for Agrotechnology Transfer software [22,23,24]. However, the CROPGRO-Canola model has not been evaluated in China.

In this study, the CROPGRO-Canola model was used to determine the long-term optimal management of rapeseed in the Yangtze River Basin. The specific objectives of this study were to: (1) calibrate the CROPGRO-Canola model for two seasons of N rate and density experiments conducted in Wuhan, China, (2) evaluate the model at the same site for two additional seasons of data, and (3) use the model to estimate the long-term economic optimum N rate and density for this area for different N prices. The results of this work will provide recommendations to farmers on optimum N and population to maximize long-term marginal net return (MNR), and it can serve as a case study to optimize rapeseed management in advance according to rapeseed futures, climate and soil conditions in the future.

2. Materials and Methods

2.1. Experimental Data

Rapeseed experiments were conducted in 2013, 2014, 2015 and 2017 in Wuhan, China, with the variety of Huayouza62, which is a half-winter Brassica rapeseed hybrid bred by Huazhong Agricultural University using Polima cytoplasmic male sterile line “2063A” and the restorer line “05-P71-2”. Huayouza62 has been approved and numbered as National Certified Oil 2011021 and is one of the most common varieties planted in the Yangtze River Basin. Different planting densities and nitrogen rates were applied each year using a split plot design with three replications, with the nitrogen rate being set as the main plot (

Table 1. Datasets for model calibration (I–II) and evaluation (III–IV) of the model.

Dataset	Growing Season	N Rate (kg ha−1)	Planting Density (Plants m−2)	Source of Data
I	2012/10–2013/05	180, 270	15, 30, 45, 60, 75	Yang, 2014 [26]
II	2014/09–2015/04	120, 240, 360	15, 30, 45	Sun, 2016 [11]
III	2013/09–2014/04	120, 240, 360	15, 30, 45	Sun, 2016 [11]
IV	2017/09–2018/05	159	15, 30, 45	Yuan, 2020 [27]

). The anthesis day was recorded as the day when 50% of the plants have open flowers on any node, and the maturity day was recorded as the day when 95% of plants reach a yellow color. Samples were collected from each plot at the overwintering, flowering and the mature stage. After removal of the roots, individual plants were separately bagged, heated at 105 °C and baked at 80 °C to constant weight, and the dry matter weight was determined. The plots were harvested separately, and the actual yield was measured.

The CROPGRO-Canola model requires properties and daily weather data, including daily maximum and minimum temperature and total precipitation, which were collected from a meteorological station near the test site (Figure 1). Daily solar radiation was calculated based on sunshine hours using the method of Angstrom et al. [25]. The soil data used in the study were derived from the Chinese soil dataset in the global soil database HWSD (Harmonized World Soil Database). The soil properties for each soil layer are shown in

Table 2. Soil properties of the experimental site.

Depth (cm)	Clay (%)	Silt (%)	Lower Limit (v v−1)	Drained Upper Limit (v v−1)	Sat. Hydraulic Conduct (cm h−1)	pH	Organic Carbon (%)	Bulk Density (g cm−3)
5	21	50	0.153	0.34	0.447	7.8	1.12	1.22
15	21	50	0.153	0.34	0.447	7.8	1.12	1.22
30	21	50	0.153	0.34	0.447	7.8	1.12	1.22
60	21	45	0.144	0.314	0.414	7.9	0.82	1.31
80	21	45	0.144	0.314	0.414	7.9	0.82	1.31
100	21	45	0.144	0.314	0.414	7.9	0.82	1.31
120	21	45	0.144	0.314	0.414	7.9	0.82	1.31
150	21	45	0.144	0.314	0.414	7.9	0.82	1.31
180	21	45	0.144	0.314	0.414	7.9	0.82	1.31
200	21	45	0.144	0.314	0.414	7.9	0.82	1.31

.

2.2. Model Description and Calibration

The CROPGRO model is a generic crop growth model with a daily time step that computes canopy photosynthesis at hourly time steps using leaf-level photosynthesis parameters and hedgerow light interception calculations. The model provides documentation that defines the parameters of species, ecotype and cultivar traits, making it possible to simulate a new cultivar in a new environment by adjusting the cultivar and ecotype parameters. In this study, a set of appropriate parameters for winter rapeseed growth simulation in Wuhan were obtained by calibrating the default parameters.

The default parameters were defined by the model and updated based on previous research. The anthesis day, maturity day, yield and above-ground biomass were the main outputs for calibration and evaluation. Datasets in 2012 and 2014 (Table 1) were randomly selected to be used for calibration. First, the Genotype Coefficient Calculator and GLUE coefficient estimator within DSSAT were used for preliminary calibration [28,29,30]. In this step, the parameters were adjusted through thousands of iterations to reduce the error of the simulated values. Sensitivity analysis was then performed to determine the main parameters that affect the output variables, and a “trial and error” method was used to adjust these parameters to minimize the error [7]. For example, the anthesis day was the most sensitive to EM-FL (time between plant emergence and flower appearance), OPTBI (minimum daily temperature above which there is no effect on slowing normal development toward flowering), RWDTH (relative width of this ecotype in comparison to the standard width per node (YVSWH) defined in the species file), SLAVR (specific leaf area of the cultivar under standard growth conditions) and SLOBI (slope of relationship reducing progress toward flowering if TMIN for the day is less than OPTBI). Thus, when simulated, anthesis day was calibrated by adjusting combinations of these parameters. Finally, all parameters were checked to determine if they were within the normal expected range. The calibrated parameter values are shown in

Table 3. Calibrated parameters in CROPGRO-Canola for rapeseed with default values and calibrated values for cultivar Huayouza62.

Parameter	Definition	Default	Calibrated
CSDL	Critical short day length below which reproductive development progresses with no daylength effect (h)	24	15.59
PPSEN	Slope of the relative response of development to photoperiod with time (1 h−1)	−0	−0.06
EM-FL	Time between plant emergence and flower appearance (R1) (PD)	29	38.51
FL-SH	Time between first flower and first pod (R3) (PD)	15	10.4
FL-SD	Time between first flower and first seed (R5) (PD)	31	15.9
SD-PM	Time between first seed (R5) and physiological maturity (R7) (PD)	25	26.55
FL-LF	Time between first flower (R1) and end of leaf expansion (PD)	3	0.81
LFMAX	Maximum leaf photosynthesis rate at 30 C, 350 vpm CO2 and high light (mg CO2 m−2 s−1)	1	0.814
SLAVR	Specific leaf area of cultivar under standard growth conditions (cm2 g−1)	250	329.2
SIZLF	Maximum size of full leaf (three leaflets) (cm−2)	100	52.79
XFRT	Maximum fraction of daily growth that is partitioned to seed-shell	1	1
WTPSD	Maximum weight per seed (g)	0	0.003
SFDUR	Seed filling duration for pod cohort at standard growth conditions (PD)	20	15.02
SDPDV	Average seed per pod under standard growing conditions (no pod−1)	22	14
PODUR	Time required for cultivar to reach final pod load under optimal conditions (PD)	10	5.435
SDPRO	Fraction protein in seeds (g g−1)	0.2	0.23
SDLIP	Fraction oil in seeds (g g−1)	0.5	0.48

PD, photothermal days.

.

2.3. Model Evaluation

The model was first run for 2013 and 2017 experiments with the datasets in Table 1 to evaluate its performance. Differences between the simulated and measured anthesis day, maturity day, yield and above-ground biomass were quantified with the following statistical indicators.

RMSE (root mean square error) and normalized RMSE (nRMSE) between simulated ($P_i$) and measured ($O_i$) values were calculated using Equations (1) and (2), where $n$ is the number of observations, and $\overline{O}$ is the average of the measured value $O_i$.

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(O_i-P_i\right)}^2}$$

(1)

$$\mathrm{nRMSE}=\frac{RMSE}{\overline{O}}\times 100$$

(2)

Model simulation efficiency (EF), mean error (ME) and its relative value rME between the simulated and measured values were calculated using Equations (3)–(5).

$$EF=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-P_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-\overline{O}\right)}^2}$$

(3)

$$ME=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left(P_i-O_i\right)$$

(4)

$$rME=\frac{ME}{\overline{O}}\times 100\%$$

(5)

An optimal model gives a good fit to experimental data when RMSE and nRMSE values are small (for winter rapeseed, the nRMSE of the growth period, yield and above-ground biomass were lower than 10%, 15% and 30%, respectively), EF is close to 1, and ME and rME are close to 0 [31].

2.4. Long-Term Optimal Management Strategies

After calibration and evaluation, the model was used to estimate the long-term optimal nitrogen rate and planting density of winter rapeseed in this area. Meteorological data of 1989–2019 from Wuhan were used to simulate the response of yield to 10 nitrogen rates (from 0 to 360 kg ha⁻¹ at an increasing rate of 40 kg ha⁻¹) and six planting densities (from 15 to 90 plant m⁻² at an increasing rate of 15 plant m⁻²). The marginal net return (MNR) was computed for each scenario to evaluate the economic benefits of each treatment:

$${\mathrm{MNR}}_{ij}=\Delta Yield\ast Price-\Delta Nrate\ast Ncost-\Delta Pop\ast Pcost$$

(6)

In the equation, $MN{R}_{ij}$ is the marginal net return in year $i$ and experiment $j$; $\Delta Yield$ is the difference between the simulated yield and yield at the nitrogen rate of 0 kg ha⁻¹ and the planting density of 15 plants m⁻²;$Price$ is the price for rapeseed, which is usually 4.2 CNY kg⁻¹; $\Delta Nrate$ is the nitrogen rate minus 0 (kg ha⁻¹); $Ncost$ is the cost of N (CNY kg⁻¹); $\Delta Pop$ is the planting density minus 15 (plants m⁻²); $Pcost$ is the cost of seed, which is usually 40 CNY kg⁻¹. $Ncost$ is calculated based on the typical price of NPK = 2 CNY kg⁻¹ and urea = 2 CNY kg⁻¹. Based on the historical fluctuations of price, there were six possible combinations of fertilizers, which included three NPK fertilizer prices (2, 4 and 6 CNY kg⁻¹) and two urea prices (2 and 4 CNY kg⁻¹).

$$MNR=\frac{1}{n}{\displaystyle \sum }MN{R}_{ij}$$

(7)

Equation (7) shows the formula used to compute the long-term average MNR of the same experiment from 1989–2019, where $n$ is the number of years. The optimal nitrogen rate and planting density are the levels that maximize the long-term MNR.

3. Results

3.1. Model Calibration and Evaluation

The model cultivar coefficients gave a good fit between simulated and measured anthesis day, maturity day, yield and above-ground biomass for the calibration years (datasets I and II) and evaluation years (datasets III and IV) (Figure 2a,b and Figure 3a,b).

In the calibration years, the RMSE of the anthesis day was 2.90 d and the model gave an nRMSE of 1.85%, an rME of 1.34% and an EF of 0.80, while the maturity day had an RMSE of 2.52 d, an nRMSE of 1.16%, an rME of 1.13% and an EF of 0.29 (Figure 2a). Additionally, to evaluate the model performance for the yield, the simulation results for the yield were presented with the determination coefficient of 1:1 regression plot (Figure 3a). Combined with other evaluation criteria, the nRMSE of yield in the calibration years was 6.25%, with an RMSE of 150.07 kg ha⁻¹, an rME of −0.83% and an EF of 0.94.

The calibrated model was then evaluated using datasets III and IV. In the evaluation years, the RMSE of the anthesis day was 0.55 d, and the model generated an nRMSE of 0.36%, an rME of 0.20% and an EF of 0.91, while the maturity day had an RMSE of 0.88 d, an nRMSE of 0.40%, an rME of 0.21% and an EF of −4.91. Combined with other evaluation criteria, the nRMSE of yield in evaluation years was 6.58%, with an RMSE of 154.76 kg ha⁻¹, an rME of −1.72% and an EF of 0.93. A graph of the results is shown in Figure 2b and Figure 3b. These results suggested that the growth period and yield of crops can be accurately predicted.

The calibration and evaluation results demonstrated that the model could duplicate historical test results, indicating that the model is appropriate for simulating the long-term optimal management of Huayouza62. Moreover, the model produced an RMSE of 3184.52 kg ha⁻¹, an nRMSE of 29.62%, an ME of −1717.69 kg ha⁻¹, an rME of −15.98% and an EF of −0.12 for above-ground biomass in the evaluation years. The dynamic simulation of the above-ground biomass in this study resulted in a lower nRMSE than that reported by Deligios et al. [18]. In addition, these results are consistent with those of other evaluations with the CROPGRO-Canola model in colder environments [32], indicating that the model can simulate N yield response under a warm climate as well. Model evaluation criteria (

Table 4. Model statistical criteria.

Attribute	Statistical Criteria
	N	RMSE	nRMSE	ME	rME	EF
	Calibration years
Anthesis day (D)	19.00	2.90	1.85	2.11	1.34	0.80
Maturity day (D)	19.00	2.52	1.16	2.47	1.13	0.29
Seed yield (kg ha−1)	19.00	150.07	6.25	−20.05	−0.83	0.94
Above-ground biomass (kg ha−1)	19.00	3051.07	30.53	−293.42	−2.94	−0.47
	Evaluation years
Anthesis day (D)	13.00	0.55	0.36	0.31	0.20	0.91
Maturity day (D)	13.00	0.88	0.40	0.46	0.21	−4.91
Seed yield (kg ha−1)	13.00	154.76	6.58	−40.38	−1.72	0.93
Above-ground biomass (kg ha−1)	13.00	3184.52	29.62	−1717.69	−15.98	−0.12
	Different nitrogen applications
Anthesis day (D)	8.00	2.00	1.29	1.00	0.65	0.84
Maturity day (D)	8.00	2.09	0.96	1.63	0.75	0.20
Seed yield (kg ha−1)	8.00	138.23	5.39	−64.75	−2.53	0.95
Above-ground biomass (kg ha−1)	8.00	4582.25	36.24	−3454.00	−27.31	−1.25
	Different planting densities
Anthesis day (D)	15.00	2.37	1.50	1.60	1.02	0.76
Maturity day (D)	15.00	1.95	0.89	1.67	0.76	0.26
Seed yield (kg ha−1)	15.00	182.10	7.48	22.47	0.92	0.61
Above-ground biomass (kg ha−1)	15.00	3138.82	30.54	−1330.13	−12.94	−0.50

) showed a good agreement between the simulated and observed anthesis day, maturity day, yield and above-ground biomass.

3.2. Impact of Rapeseed Management on Model Prediction

In addition to reasonable simulation of the yield, the model could also accurately simulate responses to nitrogen rate and planting density (Table 4).

Figure 4 shows good agreements between the simulated and observed anthesis day, maturity day and yield for different nitrogen rates. Specifically, under different nitrogen rates, the RMSE of anthesis day was 2.00 d, and the model generated an nRMSE of 1.29%, an ME of 1.00 d, an rME of 0.65% and an EF of 0.84 (Figure 4a). The maturity day had an RMSE of 2.09 d, an nRMSE of 0.96%, an ME of 1.63 d, an rME of 0.75% and an EF of 0.20 (Figure 4b). The yield had an RMSE of 138.23 kg ha⁻¹, an nRMSE of 5.39%, an ME of −64.75 kg ha⁻¹, an rME of −2.53% and an EF of 0.95 (Figure 4c). These results indicated that the model simulated the response of growth period and yield to nitrogen rate well.

The model also predicted the response to planting density well (Figure 5a,b). Specifically, under different planting densities, the RMSE for anthesis day was 2.37 d, and the model generated an nRMSE of 1.50%, an ME of 1.60 d, an rME of 1.02% and an EF of 0.76 (Figure 5a). The maturity day had an RMSE of 1.95 d, an nRMSE of 0.89%, an ME of 1.67 d, an rME of 0.76% and an EF of 0.26 (Figure 5b). In addition, the yield had an RMSE of 182.10 kg ha⁻¹, an nRMSE of 7.48%, an ME of 22.47 kg ha⁻¹, an rME of 0.92% and an EF of 0.61 (Figure 5c).

All the above results suggest that the model can simulate the response of different indices to nitrogen rate and planting density well, particularly the yield. Overall, the model simulated most of the variations in yield under different nitrogen rates (nRMSE = 5.39%) and planting densities (nRMSE = 7.48%) (Figure 4c and Figure 5c). In contrast to the simulation of yield, the simulation of above-ground biomass was rather insufficient, with an nRMSE of about 30% (Table 4). However, the results regarding the yield and growth period were sufficient to support our optimization of rapeseed management.

3.3. Long-Term Optimal Management

Table 5. Yield, marginal net return (MNR) of winter rapeseed and soil N under different nitrogen rates. The values are mean of 30 seasons.

Treatments	Yield (kg ha−1)	MNR (CNY ha−1)	Soil N (kg ha−1)
N0	1453.90	0	153.13
N40	1691.90	645.8	151.87
N80	1857.20	986.26	153.58
N120	1970.97	1110.28	156.22
N160	2076.33	1199.02	159.36
N200	2094.93	923.34	163.20
N240	2166.70	870.36	167.73
N280	2254.43	885.64	172.75
N320	2333.87	865.46	176.77
N360	2393.73	763.1	183.39

The MNR is for rapeseed price of 4.2 CNY kg⁻¹ and N price of 2 CNY kg⁻¹. N₀: 0 kg N ha⁻¹; N₄₀: 40 kg N ha⁻¹; N₈₀: 80 kg N ha⁻¹; N₁₂₀: 120 kg N ha⁻¹; N₁₆₀: 160 kg N ha⁻¹; N₂₀₀: 200 kg N ha⁻¹; N₂₄₀: 240 kg N ha⁻¹; N₂₈₀: 280 kg N ha⁻¹; N₃₂₀: 320 kg N ha⁻¹; N₃₆₀: 360 kg N ha⁻¹.

shows the yield, soil N and MNR of winter rapeseed Huayouza62 under ten nitrogen rates averaged over 30 seasons. Figure 6 shows a sensitivity analysis of yield and MNR to different N prices. Generally, an increase in nitrogen rate led to increases in the yield of winter rapeseed and soil N. However, an increase in yield did not necessarily increase the MNR. Hence, MNR should be fully considered to maximize the economic and ecological benefits.

Changes in the price of N fertilizer were also considered. The lines with different colors in Figure 6 represent the changes in MNR under different nitrogen fertilizer prices. As the nitrogen rate increased, the MNR at normal N fertilizer price (NPK = 2 CNY kg⁻¹; urea = 2 CNY kg⁻¹) first rose and reached the maximum at 120–160 kg N ha⁻¹. However, the MNR was strongly influenced by the price of N fertilizer. With increasing price of the N fertilizer, the difference between the minimum MNR (the lowest point of the curve) and maximum MNR (the highest point of the curve) increased. For instance, when the price of NPK compound fertilizer was < 2 CNY kg⁻¹ and that of urea was < 4 CNY kg⁻¹, the MNR tended to reach the maximum at the nitrogen rate of 120–160 kg N ha⁻¹. However, when the price of N fertilizer increased, the optimal nitrogen rate was lower than 120 kg N ha⁻¹, and MNR even became negative with increasing nitrogen rate.

The model was also run for ten nitrogen rates, six planting densities and six N fertilizer prices over thirty seasons of historical weather to determine the long-term optimal nitrogen rate for this area. In Figure 7, the range in between the red dotted lines represents the optimal nitrogen rate. At the NPK fertilizer price of 2 CNY kg⁻¹ and the urea price of 2 CNY kg⁻¹, the MNR reached the maximum at the nitrogen rates of 120–200 kg N ha⁻¹ and plant densities of 75–90 plants m⁻², and then decreased slightly with further increases in nitrogen rate; at the NPK fertilizer price of 2 CNY kg⁻¹ and the urea price of 4 CNY kg⁻¹, the optimal range of nitrogen rate was 80–160 kg N ha⁻¹. Moreover, a further increase in the price of NPK fertilizer and urea would lead to a decrease in the optimal nitrogen rate. Overall, the optimal planting density of rapeseed in this area is 45–75 plants m⁻² with a reasonable nitrogen rate, because rapeseed is prone to lodging at high planting densities, and dense planting may not necessarily bring greater economic benefits. Yield reduction due to lodging is not considered in the model.

4. Discussion

4.1. Simulation Performance of the CROPGRO-Canola Model

In general, the simulation error of the model for the growth stages was less than 3 days, which is consistent with the measured growth period in the field experiments, indicating an accurate simulation of the growth stages by the model. For the yield, the simulated nRMSE of the model was below 8%, which is lower than 15% (acceptable standard). The model also simulates the response to nitrogen rate and planting density well. The CROPGRO-Canola model may be highly appropriate to simulate optimum management of rapeseed. However, the nRMSE for the simulation of the model on the above-ground biomass was nearly 30%, which was not as good as the simulation of yield. Although the simulated above-ground biomass under different nitrogen rates and planting densities also showed a similar trend to the measured data (Figure 8), the accuracy of this model for the above-ground biomass needs to be further improved.

4.2. Reduction in Uncertainty in Calibration

Uncertainty has always been a major problem in the application of crop growth models. In the calibration of the CROPGRO-Canola model, it is necessary to determine the optimization targets and strategies, which requires a better understanding of the mechanism of the model. We conducted a sensitivity analysis by employing the Extended Fourier Amplitude Test (EFAST), and found the key parameters for growth period, yield and above-ground biomass [7]. We then used the Genotype Coefficient Calculator and GLUE coefficient estimator within DSSAT to iteratively adjust the parameters to make the output variables as close to the measured data as possible [28,29,30]. Finally, the main parameters were slightly adjusted to improve the calibration. In this process, if the simulated anthesis day was later than the measured anthesis day, it was necessary to adjust the parameters related to the development, such as EM-FL and PPSEN. Generally, the parameters are adjusted in the order of development, distribution of carbon and nitrogen (the biomass of different parts) and yield components, which can reduce the uncertainty of the simulation process to some extent.

Simulating growth and development of winter rapeseed could be further enhanced by measurements such as LAI (leaf area index), plant height and biomass allocation [18]. However, these measurements were not conducted in this experiment. To improve qualitative simulation of crop growth, crop models have been combined with remote sensing, which can be a surrogate for LAI measurements [21,33,34]. Thus, incorporating remote sensing with model calibrations could likely improve the calibration process.

4.3. Practical Value of Long-Term Optimization Results

In previous studies, the optimal management was usually judged based on the yield. For instance, most recent studies have indicated that the optimal nitrogen rate in this area is 120–200 kg N ha⁻¹ [11,26], which is consistent with our results at conventional fertilizer prices. However, when the price of the fertilizer increases, this optimal nitrogen rate may not be adopted by the farmers, since their prior concern is still the economic benefits, which may be affected by many factors such as investment in seeds, fertilizer, irrigation, labor and machinery cost for field operations [35]. Because the irrigation cost is generally low for the wet climate and low temperature in the growing season of rapeseed [36], the mechanization of rapeseed production has not been widely promoted, and the rapeseed farming labor is highly mobile and difficult to quantify. Thus, we only considered the most important inputs in the cultivation of rapeseed in future studies. However, more factors affecting MNR may be incorporated into Equation (6) later to make the optimization more reliable.

The results of this study would be affected by other possible external factors. Specifically, fertilizer prices were defined considering normal fluctuations in this study, but many other factors affect fertilizer prices, such as pandemic and regional conflicts [37]. Our results indicated the relationship between fertilizer inputs and fertilizer prices in Figure 7. Thus, it could be served as a basis for reducing fertilizer supplies when fertilizer prices rise above the upper end of the range in the study. The effect of climate change and extreme climate were not considered in this study; however, crop yield is highly sensitive to this [38]. Thus, more possible external factors could be considered and make the results of long-term optimal management have greater reference value to practical application.

5. Conclusions

In this study, the long-term optimal management of winter rapeseed cultivation was studied. The results considered long-term historical climate, yield, fertilizer (i.e., price, amount and types), rapeseed population, and seed prices. The major results can be concluded as follows: (1) the CROPGRO-Canola model was robustly calibrated and evaluated in Wuhan, China; (2) fertilizer input increased rapeseed yield within a certain range, but decreased when it exceeded the range; (3) long-term optimal nitrogen rate was 120–160 kg N ha⁻¹, and the optimal planting density was 45–75 plants m⁻² under normal fertilizer prices.

Although the main factors were considered, some possible external factors (i.e., pandemic and climate change) may influence the results. Some uncertainty of the model results exists due to incomplete datasets, calibration method and potential errors in the model growth and development theory.

In conclusion, this study optimized the CROPGRO-Canola model and provided researchers and decision-makers with references to take adaptive management to ensure oil security in future.