Modeling biomass and yield production based on nitrogen accumulation in soybean grown in upland fields converted from paddy fields in Japan

ABSTRACT Crop models can help in identifying constraints to crop production and enhancing crop yield. Because one of the major constraints is the availability of nitrogen for seed production, evaluation of nitrogen accumulation is important for soybean crop models. In the present study, we developed a soybean crop model to evaluate biomass production and nitrogen accumulation of Japanese soybean cultivars grown in upland fields converted from paddy fields. The model simplified the effects of leaf nitrogen deficit on biomass production by reducing only the leaf area while maintaining leaf nitrogen concentration, and we introduced a zero-order reaction kinetics model to estimate soil nitrogen supply. The proposed model was observed to simulate the accumulation of aboveground biomass and nitrogen content in both nodulating and non-nodulating soybean cultivars until the mid-seed filling period. The model accounts for 94% of the observed variations in aboveground biomass and nitrogen content and for 69% of the observed variation in seed biomass. The normalized root mean square errors of aboveground biomass, nitrogen content, and seed biomass are 19.3%, 19.6%, and 20.2%, respectively. Because the model does not include the effects of soil water status on pod formation and nitrogen fixation, seed biomass was overestimated in some cases. However, our model quantified the effects of changes in the soil nitrogen supply and biological nitrogen fixation on soybean production and will be useful for identifying and eliminating production constraints in Japanese soybeans.


Introduction
The improvement in the average soybean yield in Japan has stagnated over the past 30 years, and the yield was approximately half of that in two major soybeanproducing countries -the United States and Brazil -in 2017(FAO, 2020. Several reports have identified genetic and environmental factors as the cause of the low yields. Regarding genetic factors, Japanese soybean cultivars have shown lower seed yields than the United States cultivars regardless of location (Kawasaki et al., 2018;Rao et al., 2002). Regarding environmental factors, in Japan, soybean is mostly cultivated in upland fields converted from paddy fields under high temperature and humidity conditions; these factors often limit the production owing to conditions including waterlogging caused by excessive rainfall, poor drainage, and a high water table (Kokubun, 2013). The poor drainage condition restricts nitrogen acquisition by soybeans via the reduction of biological nitrogen fixation and soil nitrogen absorption (Maekawa et al., 2011;T. Takahashi et al., 2006). In addition, continuous soybean production in upland fields converted from paddy fields reduces soil nitrogen fertility, and it has been suggested to cause a long-term reduction in soybean yield in northern Japan (Sumida et al., 2005). Because soybeans require larger amount of nitrogen for seed production owing to the high protein content of the seeds, enhancement of total nitrogen uptake is essential to increase soybean yields (Rotundo et al., 2014;Y. Takahashi et al., 1992). However, the effects of nitrogen limitation on soybean yield have not been quantitatively evaluated for Japanese soybeans cultivated in upland fields converted from paddy fields.
Crop models are useful tools to quantitatively evaluate crop productivity because they numerically integrate multiple processes of crop growth and consider interactions between genotypes, environment, and management at the field scale (Boote et al., 2013;Jeuffroy et al., 2014). Various crop models have been developed over the past few decades with different levels of complexity, depending on the objectives of the models (Jones et al., 2017). However, most soybean crop models have been developed in common upland fields under relatively dry conditions (Boote et al., 2008;Sinclair et al., 2003) and may be unsuitable to evaluate production constraints in upland fields converted from paddy fields. Although some models have been developed for soybean production in Japan, they only simulated biomass production without nitrogen accumulation in plant (Abe & Momoya, 1987;Sameshima, 2000). Therefore, developing a soybean crop model that simulates nitrogen accumulation in the crops cultivated in upland fields converted from paddy fields is essential to quantify the production constraints and improve the soybean yield in Japan.
Soybean acquires nitrogen from biological fixation as well as from soil nitrogen mineralization. The fixation amount is determined by highly complex processes involving rhizobial bacteria in nodules and an association between plant and soil conditions. These complexities alter the description of nitrogen fixation among crop models (Liu et al., 2011). The CROPGRO model estimates the potential nitrogen fixation rate based on nodule biomass including the carbon cost for nodule activities (Boote et al., 2008). By contrast, the Sinclair model estimates the potential nitrogen fixation rate based simply on the linear relationship between nitrogen fixation and aboveground biomass (Sinclair, 1986). This linear relationship was confirmed by experimental results for soybean production in upland fields converted from paddy fields in Japan (Shiraiwa et al., 1994). Considering the difficulty of measuring nodule biomass in the field, the APSIM model uses this linear relationship to estimate the nitrogen fixation for legume crops (Chen et al., 2016).
The descriptions of soil nitrogen mineralization are different among the abovementioned models, although the Sinclair model does not consider nitrogen acquisition from mineralized soil nitrogen. The CROPGRO and APSIM models estimate the soil nitrogen mineralization rate (RNSM) based on the dynamics of water, carbon, and nitrogen in soil systems (Gijsman et al., 2002;Probert et al., 1998). However, these models are highly complicated, and numerous parameters are required to adapt to a particular condition. Moreover, obtaining direct measurements of some parameters associated with organic matter pools is difficult (Benbi & Richter, 2002). By contrast, Sugihara et al. (1986) simply expressed RNSM via a reaction model using the Arrhenius law because soil nitrogen mineralization is a result of microbial activity. Nira and Hamaguchi (2012) develop a zeroorder reaction kinetics model by employing the Arrhenius law and estimated the mineralized nitrogen accumulation in soybean cultivated in upland fields converted from paddy fields.
The objective of the present study was to develop a soybean crop model that simulates crop growth and production in upland fields converted from paddy fields by incorporating nitrogen components of biological fixation and soil nitrogen mineralization. The soybean crop model reported by Sinclair (1986) and Sinclair et al. (2003), which estimated nitrogen fixation based on its linear relationship with aboveground biomass, was used as base model and the zero-order reaction kinetics model was introduced to estimate soil nitrogen mineralization (Nira & Hamaguchi, 2012). For validation of the model, we compared the model results to the experimental data of nodulating and non-nodulating soybeans cultivated in upland fields converted from paddy fields. The developed model can separately estimate nitrogen acquisition from mineralized soil nitrogen and that from biologically fixed nitrogen and can help identify and quantify the constraints to soybean production in Japan.

Field experiments and measurements
Field experiments were conducted in an experimental upland field, which had been converted from a paddy field, of the National Agriculture and Food Research Organization (NARO), Tsukubamirai (36°00′ N, 140°01′ E, gray lowland soil), Ibaraki, Japan, from 2013 to 2015. In 2013, basal chemical fertilizer (2.5, 10.0, and 10.0 g m −2 of N, P 2 O 5 , and K 2 O, respectively) was applied. In 2014 and 2015, 2.7 g m −2 of N and 3.3 g m −2 P 2 O 5 were applied to the soil before sowing, but K 2 O was not applied based on the soil analysis results. Three determinate Japanese soybean cultivars, Ryuho, Enrei, Fukuyutaka, and a non-nodulating mutant derived from Enrei (En1282, Francisco & Akao, Francisco & Akao, 1993) were planted at two sowing dates each year in plots with 0.7-m row width. In addition, the plots with 0.35-m row width were added at the second sowing date in 2015. All plants were planted 0.15 m apart from each other within the row. Ryuho was planted for 2 years -in 2013 and in 2015; Enrei and Fukuyutaka were planted for 3 years from 2013 to 2015; and En1282 was planted for 2 years -in 2013 and 2014. Each cultivar was sown as follows: on June 12 and July 31 in 17 June 2013 and July 17 in 2014, and June 4 and June 30 in 2015. The experiments were conducted as a split-split plot design with three replicates where planting dates were the main plots, cultivars were the subplots, and plant densities were sub-sub plots. In 2013 and 2014, each plot comprised six rows of 6 m long (25.2 m 2 ). In 2015, each 0.7-m row plot comprised six rows of 7 m long (29.4 m 2 ), and each 0.35-m row plot comprised eight rows of 7 m long (19.6 m 2 ). Intertillage and ridging were conducted in the 0.7-m row plots before flowering, and no irrigation was applied. Weather condition in 2013 was warm and sunny compared with the 20-year mean during the soybean growing season ( Table 1). The precipitation in October 2013, in June and October 2014, and in July and September 2015 was higher than the 20-year mean, whereas that in July and November 2013, in September 2014, and in October 2015 was lower. The solar radiation in August and September 2015 was lower than the 20-year mean but that in October was higher.
A total of 8 plants from the central four rows in the 0.7-m row plots and 12 plants from the central six rows in the 0.35-m row plots were destructively sampled at the vegetative, flowering (R2), initial seed filling (R5), mid-seed filling, and full maturity (R8) stages. The growth stages were determined according to Fehr and Caviness (1977), and all plants were manually harvested above the cotyledon node. Four representative plants were selected to measure the dry weight of the stem, leaf, petiole, pod, and seed to calculate the ratio of the dry weight of each organ. The total dry weight of the entire plant was measured for all collected samples, and the dry weight of each organ was calculated by multiplying the dry weight of the entire plant with the ratio of each organ's dry weight. Soybean biomass was presented based on the dry weight. In 2013 and 2014, the nitrogen concentration of each organ was measured using the dry combustion method (Sumigraph, NC-22 F, Sumika Chemical Analysis Service CO, Ltd, Japan) for one replicate. A fraction of canopy radiation interception (FRI) was measured using digital image analysis as previously reported (Nakano et al., 2020). Soil samples were obtained from a depth of 0.05-0.15 m at five randomly selected points before the basal fertilizer was applied. Thereafter, the samples were mixed and analyzed for NO 3 -N and NH 4 -N concentrations.

Model development
The soybean crop model proposed in this study simulates daily soybean growth and yield using soil properties and meteorological data, i.e. air temperature, day length, and solar radiation. The model consisted of six submodels for predicting phenology, leaf area, biomass production, biomass partitioning, nitrogen supply, and nitrogen partitioning ( Figure 1 and Table 2). The phenology and leaf area submodels specific to Japanese soybean cultivars have already been developed (Nakano et al., 2015(Nakano et al., , 2020. However, the previous leaf area submodel has only covered leaf area growth before seed filling. Therefore, we extended the present model to include the adaptations after seed filling taking into consideration the reduction of leaf area growth and promotion of leaf abscission with seed growth.

Extending the leaf area model to include the seed filling period
To extend the leaf area submodel for including the period after seed filling, we first determined the timing when the leaf appearance and leaflet expansion were stopped by the transition from the vegetative to reproductive phase. Leaf area growth at the main stem mostly stops at R5 (Sameshima, 2000;Tenorio et al., 2017). However, although the rate of increase was considerably reduced, it was observed that the primary stem node, where the leaf is attached, continued increasing after R5 (Torrion et al., 2012). Moreover, leaf area growth increases with a higher CO 2 concentration even during the seed filling period (Jin et al., 2017). Therefore, our model assumed that leaf area growth continued after seed filing to account for the increase in leaf and petiole biomass attributed to biomass partitioning for seed production. In addition, we assumed that the timing of the complete shift of the soybean growth stage from a vegetative to reproductive phase was the same as the timing of the complete cessation of a new pod formation. Consequently, leaf area growth was also considered complete until this time. The time course of flower production of determinate soybeans indicates bimodal peaks, and the termination of flowering at higher racemes is observed to coincide approximately with the initiation of seed growth at lower racemes (Saitoh et al., 1998). Therefore, our model calculated the period from flowering to the last pod set for the entire plant as twice the period from R2 to R5 based on the assumption that the termination of the last flowering coincided with the timing of R5. Nitrogen translocation from the leaf and petiole to seed during the seed filling period induces a simultaneous decline in leaf nitrogen concentration and leaf area (Oikawa et al., 2013). However, to simplify the model, we quantified the reduction in leaf nitrogen content by reducing only the leaf area following Soltani et al. (2006) for  Table 2. Table 2. List of the equations of the soybean crop model proposed in the present study.

Phenology submodel
Details of soybean phenology model are provided in Nakano et al. (2015).
Leaf area submodel Details of leaf area growth model are provided in Nakano et al. (2020). Leaf area growth continues until twice the period from flowering (R2) to initial seed filling (R5) after R2, when there is an increase in leaf and petiole biomass.
Biomass production submodel (2) Biomass partitioning submodel  Day of year of initial seed filling (R5) RBSDmx: Potential rate of seed biomass increase (g m −2 d −1 ) BSD: Seed biomass (g m −2 ) RBSPmx: Potential rate of stem and pod biomass increase (g m −2 d −1 ) BSP: Stem and pod biomass (g m −2 ) BPS: Biomass partitioning ratio for stem and pod (= 0.4) RBLPmx: Potential rate of leaf and petiole biomass increase (g m −2 d −1 ) BLP: Leaf and petiole biomass (g m −2 ) BLPnm: Leaf and petiole biomass at maximum nitrogen concentration (g m −2 ) RBVG: Rate of available biomass for vegetative growth (g m −2 d −1 ) Nitrogen supply submodel Nitrogen partitioning submodel ) estimating daily biomass production as described in the subsequent sections.

Modeling biomass production
Daily biomass production was estimated from canopy intercepted radiation and radiation use efficiency (RUE). The canopy intercepted radiation was calculated from solar radiation and FRI, which is dependent on leaf area index (LAI). RUE is affected by environmental factors, such as air temperature, soil water deficit, and atmospheric CO 2 concentration. We only considered the temperature effects on RUE as our field experiments were not conducted under conditions of severe soil water deficit or high CO 2 concentration. Thereafter, RUE was estimated using a three-segment function with cardinal temperatures and the maximum RUE (RUEmx). Our model defined RUEmx as a model parameter, and it was determined by minimizing the differences between the observed and estimated aboveground biomass before R5. In the optimization process for RUEmx, FRI measured using the digital image analysis was used to avoid including estimation errors of LAI in calculating canopy intercepted radiation. Further, biomass production was observed to be affected by the nitrogen content in leaves (Sinclair & Horie, 1989), and the reduction of leaf nitrogen content, in turn, induces the reduction in both the nitrogen concentration in leaves and leaf area. The reduction of the leaf nitrogen concentration then decreases RUE. However, to simplify the model structure, our model quantified the effects of leaf nitrogen reduction based on leaf area. Therefore, we used the LAI at maximum nitrogen concentration (LAInm) for calculating biomass production to compensate for the effects of the reduction of leaf nitrogen concentration (Soltani et al., 2006). The value for LAInm was calculated using the nitrogen content at the maximum leaf and petiole nitrogen concentration (CLNPmx), which was identified as model parameter, and the specific leaf area, which included the petiole, because we considered the leaf and petiole as one group.

Modeling biomass partitioning
Vegetative biomass was divided into two parts in the model: the first was the stem and pod and the second was the leaf and petiole. Subsequently, the biomass distribution in these two parts was calculated using a biomass portioning ratio for stem and pod biomass (BPS), which was defined as a model parameter. We treated the leaf and petiole as one group because they dropped with senescence during the seed filling period, with the translocation of nitrogen to the seed.
The harvest index (HI) -the ratio of seed biomass to total plant biomass -is used to estimate seed growth in crop models because the HI linearly increases with time during the seed filling period (Brisson et al., 2003;Robertson et al., 2002;Sinclair, 1986). This linear increase has also been reported for Japanese soybean cultivars (Sameshima, 2000). However, HI for soybean increases with increasing seed biomass as well as with increasing fallen leaf and petiole biomass, resulting in an overestimation for biomass partitioning to seed (Hay, 1995). Therefore, to eliminate the uncertainty associated with the abscission of leaves and petioles from the HI, we defined a modified HI (HIM) as the ratio of seed biomass to total plant biomass excluding the leaf and petiole biomass. The rate of increase in HIM was determined using the observed data. However, because almost all leaves and petioles were fallen at R8, the maximum value of HIM was set to 0.6 according to the value of HI (Sameshima, 2000). In addition, seed growth required 1.3 times the biomass products as vegetative growth when considering the difference in the carbon costs between seed and vegetative organs (Sinclair, 1986).

Modeling nitrogen supply
Our model assumed three sources of nitrogen acquisition for soybeans: mineralized soil nitrogen, fertilizer nitrogen, and biologically fixed nitrogen. The mineralized soil nitrogen consisted of two sources: the nitrogen contained before sowing (INSM) and that mineralized from organic soil nitrogen during the soybean growth period. The value of INSM was determined from the measured soil inorganic nitrogen concentration and soil mass of the sampled soils before plowing. The total nitrate and ammonium inorganic nitrogen concentrations in the upper soil layers (0.05-0.15 m) were 1.90 mgN 100 g −1 soil in 2013, 1.65 mgN 100 g −1 soil in 2014, and 1.85 mgN 100 g −1 soil in 2015, respectively. We used the average value (1.8 mgN 100 g −1 soil) as a model input parameter throughout the research period. The soil mass of rooting depth was calculated from the root depth and bulk density. We have used the values of maximum rooting depth as 0.2 m and the bulk density as 115 kg m −3 from the field measurements, and the root growth rate was set as 0.03 m d −1 according to Soltani and Sinclair (2012). The fertilizer nitrogen was applied as basal fertilizer: 2.5 gN m −2 for 2013 and 2.7 gN m −2 for 2014 and 2015, and the nitrogen utilization rate for chemical fertilized nitrogen was set as 30% in the model (Kato et al., 1983).
RNSM during the soybean growth period was estimated using a zero-order reaction kinetics model, which is reportedly applicable to Japanese soybeans grown in an upland field converted from a paddy field (Nira & Hamaguchi, 2012). Because our study was conducted at the same experimental farm as that of Nira and Hamaguchi (2012), we used the same model parameters: 0.496 mgN kg −1 d −1 for the mineralization rate constant at 25°C (k 25 ) and 56,900 J mol −1 for the apparent activation energy (Ea). Although soil temperature is required to estimate RNSM using the reaction kinetics model, we used air temperature instead of soil temperature. As one example, the measured mean air temperature and soil temperature at the middle of rooting depth were 22.56° C and 22.51°C, respectively, during the soybean growth period (from August to October) in an upland field converted from a paddy field in 2016. This difference was observed result in an overestimation of ~2.3% of mineralized soil nitrogen. Then, the nitrogen utilization rate for mineralized soil nitrogen (UNSM) was defined as a model parameter.
Considering the carbon cost of biological nitrogen fixation, the biologically fixed nitrogen was assumed to be utilized for soybean growth when the nitrogen requirement in the soybean crop exceeded the nitrogen uptake from the soil (Liu et al., 2011). Moreover, because most of the absorbed nitrogen from the soil is not directly used for seed production (Ohyama, 1983), our model assumed that only the fixed nitrogen was directly used for seed growth. The absorbed soil nitrogen was assumed to be transported to the vegetative organs and was subsequently translocated to seed. The potential rate of biological nitrogen fixation (RNFX) was estimated from the linear relationship between the nitrogen fixation rate and the aboveground vegetative biomass using the nitrogen fixation coefficient (CNFX). The value of CNFX was set as 0.0017 gN g −1 according to Shiraiwa et al. (1994). In addition, considering the competition for the carbon source between seed growth and nodule activity, nitrogen fixation was assumed to cease when daily biomass production was lower than the daily potential seed growth rate (Sinclair et al., 2003).

Modeling nitrogen partitioning
Our model has calculated the daily nitrogen requirements from the daily biomass increase and nitrogen concentration in each organ. The growing organs were assumed to contain the maximum nitrogen concentrations, which were defined as model parameters for the stem and pod (CNSPmx), the leaf and petiole (CNLPmx), and the seed (CNSDmx). Subsequently, the actual nitrogen accumulation was determined by comparing the daily nitrogen requirements of soybean with the daily nitrogen acquisitions from soil and from the biological nitrogen fixation. In addition, the model assumed that the soybean crop utilized the nitrogen contained in the sown seed and recovered the nitrogen from fallen leaves and petioles via canopy shading prior to using the soil nitrogen and biologically fixed nitrogen. The sown seed was assumed to have a dry weight of 0.3 g seed −1 and contain the maximum nitrogen concentration of seed (CNSDmx).
When the seed nitrogen requirement was not fulfilled using the acquired nitrogen content in the model, the nitrogen accumulated in the vegetative organs was translocated to the seed. The fraction of the nitrogen translocated from each vegetative organ -the stem and pod and the leaf and petiole -was determined by the relative amount of translocatable nitrogen in each part (Sinclair et al., 2003). The translocatable nitrogen was estimated from the difference between the present nitrogen content and minimum nitrogen content, which was determined from the minimum nitrogen concentration for the stem and pod (CNSPmi) and the leaf and petiole (CNLPmi).

Model calibration and validation
Similar to a previous study (Nakano et al., 2020), we estimated the model parameters using the data from 2014 and the 0.7-m row plots for 2015 and further validated the model results using the data from 2013 and the 0.35-m row plots for 2015. Model performance was evaluated based on root mean square error (RMSE) and the normalized RMSE (NRMSE, %), which was calculated as the RMSE divided by the mean observed value. The model included the following 10 parameters: RUEmx, BPS, two parameters for estimating RHIM (RHIMa and RHIMb), UNSM, CNSPmx, CNLPmx, CNSDmx, CNSPmin, and CNLPmi. These parameters were directly determined or estimated by minimizing the error between the estimated and observed values. Parameter optimization was performed using Microsoft Excel Solver (Microsoft Excel for Office 365 MSO, Microsoft Corporation, Redmond, WA, US).

Calibration of the model parameters
The aboveground biomass was determined to be closely related to the cumulative radiation intercepted by the canopy, and the value of RUEmx was calculated as 1.0 g MJ −1 by fitting the estimated aboveground biomass before R5 (Figure 2). The stem and pod biomass linearly increased with increasing aboveground biomass before R5, and the value of BPS was determined to be at 0.4 (Figure 3). These parameters were only evaluated before R5 because the leaves and petioles start to fall during the seed filling period, and this affects the parameter estimation. Our model assumed that the values of RUEmx and BPS were defined as constant during the entire growing period. The value of HIM linearly increased with increasing days after R5, although the occurrence of some fluctuations was observed (Figure 4). We   defined the values of RHIMa and RHIMb -the rate of HIM increase and the lag period for HIM increase after R5 -as 0.0142 d −1 and 4 d, respectively. The outlier point in Figures 2 and 3 was collected from the June sowing of Fukuyutaka in 2014 at R5, which showed detached leaves and petioles due to severe lodging and selfshading before the seed filling period.
The stem and pod nitrogen concentration was ~0.03 gN g −1 at the initial growth stage ( Figure 5). Thereafter, the upper limit was decreased and became stable at approximately 0.02 gN g −1 as biomass increased ( Figure 5). The upper limit of the leaf and petiole nitrogen concentration also showed a similar tendency to the stem and pod, and it was stable at approximately 0.04 gN g −1 . By contrast, the seed nitrogen concentration was stable at approximately 0.07 gN g −1 and 0.06 gN g −1 for nodulating and non-nodulating cultivars, respectively. Near maturity, the lower limits of the nitrogen concentration of the stem and pod and the leaf and petiole were found to be below 0.005 gN g −1 and 0.01 gN g −1 , respectively. Based on these results and previous reports (Oikawa et al., 2013;Soltani & Sinclair, 2012), we defined the values of CNSPmx, CNLPmx, CNSPmin, and CNLPmi as 0.02 gN g −1 , 0.04 gN g −1 , 0.005 gN g −1 , and 0.01 gN g −1 , respectively; meanwhile, the value of CNSDmx was defined as 0.07 and 0.06 gN g −1 for nodulating and nonnodulating soybeans. The value of UNSM was calculated as 100.6% when fitting the aboveground nitrogen content for En1282, and we adjusted this to 100% for the model.

Validation of the model
Examples of the seasonal patterns of estimated and observed aboveground biomass and nitrogen content for Enrei and En1282 in 2013 are presented in Figure 6.
The aboveground biomass of En1282 was similar to that of Enrei until the mid-seed filling period for the June sowing and throughout the growing season for the July sowing. By contrast, the aboveground nitrogen content of En1282 was considerably smaller than that of Enrei after R5, although it was similar until R2. The model estimated aboveground biomass and nitrogen content were consistent with the observed values between Enrei and En1282 until the mid-seed filling period. However, both the aboveground biomass and nitrogen content of Enrei were found to be overestimated at maturity (R8). Notably, the rapid decline of the estimated aboveground biomass and nitrogen content before R8 was caused by the abscission of leaves and petioles. Our model simply estimated the dropping of leaves and petioles when the calculated leaf and petiole nitrogen concentration from the biomass and nitrogen contents decreased below CNLPmi. This simplification did not affect the seed biomass estimation because our model calculated the potential increase in seed biomass from HIM, which excluded the leaf and petiole biomass from total plant biomass.
The model accounted for 94% of the variation in aboveground biomass and nitrogen content during the entire growth period and 68% of the variation in seed biomass at R8 (Figure 7). The calculated relative RMSE was 19.3%, 19.6%, and 20.2% for the aboveground biomass, aboveground nitrogen content, and seed biomass estimations, respectively. The graph shows the systematic bias for seed biomass, indicating that the model tended to overestimate at low yield and underestimate at high yield ( Figure 7). This seed biomass overestimation induced the aboveground biomass and nitrogen content overestimations of Enrei at R8 (Figure 6). Figure 5. Relationship between the biomass and nitrogen concentration in the stem and pod, leaf and petiole, and seed for the model calibration data. Solid lines indicate the maximum and minimum nitrogen concentrations in the stem and pod (CNSPmx = 0.02 gN g −1 , CNSPmi = 0.005 gN g −1 ) and leaf and petiole (CNLPmx = 0.04 gN g −1 , CNLPmi = 0.01 gN g −1 ) as well as the maximum nitrogen concentration in the seeds (CNSDmx = 0.07 gN g −1 ).

Discussion
The present study aimed to develop a soybean crop model that estimates biomass production and nitrogen accumulation specific to Japanese soybeans. We have applied three modifications to the base model (Sinclair, 1986;Sinclair et al., 2003) in an effort to simplify and improve the model, as described in the Materials and Methods section. The first modification was to quantify the effects of the leaf nitrogen deficit on biomass production by reducing only the green leaf area. This simplification compensated for the effects of reducing the leaf nitrogen concentration on RUE by overestimating the loss in leaf area. Model analysis indicates that the effects of nitrogen deficit on maize production was not sensitive to the difference in approaches -reduction of both RUE and leaf area or reduction of only leaf area. (Li et al., 2009). Moreover, for soybean crops, the loss of leaf nitrogen during the seed filling period was observed to induce a larger reduction of green leaf area than  leaf nitrogen concentration (Muchow et al., 1993). Therefore, the simplification of the effects of the leaf nitrogen deficit by reducing only the leaf area will be reasonable for estimating the biomass production of soybeans.
The second modification of the model involved the use of a value for HIM that excluded the leaf and petiole biomass from HI to estimate seed growth. The promotion of leaf abscission is highly complicated because it is affected by several internal and external factors, including the translocation of nitrogen from leaves to seed, soil nitrogen deficiency, and canopy shading (Oikawa et al., 2013). This modification avoids the need to precisely estimate the senescence process of detached leaves for considering the increase in HI due to loss of leaf and petiole biomass. A close relationship has been observed between HIM and actual HI, which includes abscised leaves and petioles, for Japanese cultivars because of the close relationship between leaf biomass and stem biomass (Kakiuchi, 2013). In our model, the HIM was only used to estimate the potential seed growth, and the actual seed growth was calculated by accounting for the nitrogen available for seed production.
The third modification to the model was the addition of the zero-order kinetics model to estimate the mineralized soil nitrogen. The calibration procedure and parameter data sets of the zero-order kinetics model have been presented for numerous soils in Japan (Furue & Uwasawa, 2001;Sugihara et al., 1986). This modification enabled our soybean crop model to simulate the growth of nodulating and non-nodulating cultivars ( Figure 6). Using the calibration data of the non-nodulating soybean cultivar, we estimated the value of UNSM as 100%, implying that all the mineralized nitrogen in the soil was used for soybean growth. The values of UNSM are estimated to be in the range 82%-98% for wheat (Nira & Nishimune, 1998) and 68%-97% for corn (Saito & Ishii, 1987), which are comparable to our results; this finding further suggests that the zero-order kinetics model can be applied for estimating the soil nitrogen supply for soybean production. Because our soybean crop model can separately estimate the amounts of nitrogen uptake from soil and from biological nitrogen fixation, the model can be used to assess the contribution of these nitrogen variables to soybean yield. This could be useful in evaluating the effects of the inhibition of nitrogen fixation caused by excess water stress and soil nitrogen deficit on soybean production, which have been emphasized as the causes of the low yield of Japanese soybeans (Sumida et al., 2005;T. Takahashi et al., 2006).
Our model has accurately estimated the aboveground biomass and nitrogen content regardless of cultivars, sowing dates, and plant densities (Figure 7). However, seed biomass was overestimated for Enrei sown in June 2013 and all the nodulating cultivars sown in July 2013. During the June sowing in 2013, the observed seed biomass of Enrei was smaller than that of Ryuho by 92.6 g m −2 , although the phenologies of both cultivars were similar. In 2013, the field had experienced temporary drought from early-to mid-August, with a total precipitation of 4 mm from July 30 to August 19, during which both Enrei and Ryuho sown in June reached the initiation of pod formation. However, we could not identify the factors that induced a larger yield reduction in Enrei than in Ryuho. On the other hand, during the July sowing in 2013, the seed biomass overestimation was only found for nodulating cultivars and not for the non-nodulating mutant En1282. Therefore, this overestimation may be associated with the nitrogen fixation. Waterlogging occurred in mid-October in 2013 with a total precipitation of 183.5 mm from October 14 to October 15, during which the soybeans sown in July reached their mid-seed filling period. Because nitrogen fixation has been determined to be sensitive to the soil water stress of excessive wetness (Maekawa et al., 2011), the effect of soil water stress should be considered in the model for future studies.
By contrast, it is difficult to explain the cause of the seed biomass underestimations for Fukuyutaka, which was sown in July 2015 in the 0.35-m row plot. The underestimation was not apparent for Enrei cultivated in the 0.35-m row plots and the calibration data of Fukuyutaka in the 0.7-m row plot with the same sowing date. Therefore, the model underestimation for Fukuyutaka in the 0.35-m row plot might be caused by a combination of the model adaptability for the cultivar, sowing date, and plant density.
Finally, to the best of our knowledge, our model is the first to estimate soybean growth by balancing biomass production and nitrogen accumulation for Japanese soybeans. Moreover, this model separately estimates the nitrogen acquisition from mineralized soil nitrogen and that from biologically fixed nitrogen by soybeans cultivated in upland fields converted from paddy fields. Although further studies are needed, our model will contribute in identifying the yield-limiting factors in Japanese soybeans and can thus lead to better management.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was partly supported by the Council for Science, Technology, and Innovation (CSTI), Cross-ministerial Strategic Innovation Promotion Program (SIP), 'Technologies for creating next-generation agriculture, forestry and fisheries' (funding agency: Bio-oriented Technology Research Advancement Institution, NARO), and JSPS KAKENHI Grant Number 20H03110.