An Integrative Process-Based Model for Biomass and Yield Estimation of Hardneck Garlic (Allium sativum)

We introduce an integrative process-based crop model for garlic (Allium sativum). Building on our previous model that simulated key phenological, morphological, and physiological features of a garlic plant, the new garlic model provides comprehensive and integrative estimations of biomass accumulation and yield formation under diverse environmental conditions. This model also showcases an application of Cropbox to develop a comprehensive crop model. Cropbox is a crop modeling framework featuring declarative modeling language and a unified simulation interface for building and improving crop models. Using Cropbox, we first evaluated the model performance against three datasets with an emphasis on biomass and yield measured under different environmental conditions and growing seasons. We then applied the model to simulate optimal planting dates under future climate conditions for assessing climate adaptation strategies between two contrasting locations in South Korea: the current growing region (Gosan, Jeju) and an unfavorable cold winter region (Chuncheon, Gangwon). The model simulated the growth and development of a southern-type cultivar (Namdo, ND) reasonably well. Under Representative Concentration Pathway (RCP) scenarios, an overall delay in optimal planting date from a week to a month, and a slight increase in potential yield were expected in Gosan. Expansion of growing region to northern area including Chuncheon was expected due to mild winter temperatures in the future and may allow ND cultivar production in more regions. The predicted optimal planting date in the new region was similar to the current growing region that favors early fall planting. Our new integrative garlic model provides mechanistic, process-based crop responses to environmental cues and can be useful for assessing climate impacts and identifying crop specific climate adaptation strategies for the future.

1. Introduction

Garlic (Allium sativum) is a historically important horticultural crop in many countries with global production reaching 30.7 million tons in 2019 after a 40 % increase in production in the last decade (FAOSTAT, 2020). The physiology of garlic has been extensively studied with an emphasis on characteristics as a bulbous crop where a clove of the bulb is planted for growth and a newly grown bulb is harvested for storage and the next round of planting (Takagi, 1989; Kamenetsky, 2007). Some knowledge has been transferred to building crop models specifically targeted for simulating garlic growth and estimating yield at harvest. An early attempt for building a whole-plant garlic model was based on radiation-use efficiency (RUE) to obtain the total amount of carbon assimilates (Rizzalli et al., 2002). CropSyst, which is also a crop model based on RUE, was parameterized for garlic and used for simulating crop rotation between garlic and wheat (Giménez et al., 2016). Other studies focused on a certain aspect of the growth and development of garlic. The crop coefficient (K_c) needed for calculating evapotranspiration with the Penman-Monteith equation was specifically determined for garlic (Villalobos et al., 2004). Photosynthesis and transpiration responses to various environmental conditions were obtained for building a leaf-level gas-exchange model for garlic (Kim et al., 2013). Photosynthetic responses under elevated CO₂ and nitrogen fertilization were further investigated for building a robust model for future climate conditions (Nackley et al., 2016). With a coupled gas-exchange model parameterized for garlic, a process-based model for simulating leaf development and growth of hardneck garlic was developed (Hsiao et al., 2019). The phenology of leaf initiation and appearance was individually tracked by taking account of multiple cues including thermal time accumulation, bulb storage effect, and photoperiod. Individual leaf elongation was translated and aggregated into leaf area expansion at canopy level which was then divided into two layers of sunlit and shaded leaves for accounting for assimilated carbon based on the coupled gas-exchange model. However, carbon partitioning into plant organs such as the bulb was not validated and thus not used for yield estimation. Carbon partitioning is a crucial step in yield estimation modeling for horticultural crops in the sense that the final yield is a result of biomass partitioned into a certain organ, such as the bulb, to be harvested (Marcelis et al., 1998). The early garlic model used a set of multiple partitioning coefficients dynamically varying with developmental stages of the plant (Rizzalli et al., 2002).

Yield estimation for garlic cultivars grown in each region has been crucial to manage growing practices and maintain a stable supply in the market (Põldma et al., 2005; Abdalla et al., 2011; Lee et al., 2011; Portela et al., 2012). Regression models based on weather data or satellite images were often used for yield estimation at a large scale with a minimum set of historical input data (Choi and Baek, 2016; Gómez et al., 2021). Although worldwide garlic production has been steadily increased in recent years, it is unclear if yield will be stable in the current growing regions under future climate conditions. For estimating yield under previously unobserved conditions, it is critical to developing a physiological understanding of how garlic would respond to these environmental cues (Nackley et al., 2016). Regional climate changes may lead to changes in current farming practices and shifts in growing regions to maintain or maximize the yield. The crop model is an important tool for supporting such decisions by enabling simulations of plant growth under diverse climate adaptation strategies (Rosenzweig et al., 2014; Holzkämper et al., 2015; Corbeels et al., 2018).

In this study, our primary objective was to build an integrative process-based garlic model suitable for estimating harvestable biomass as scape and bulb yield under diverse environments including future climate conditions. A new model was developed based on an existing process-based garlic model with an original emphasis on phenology and extended with a focus on biomass accumulation and partitioning (Hsiao et al., 2019). The model was also improved to better reflect physiological responses to temperature by taking into account storage conditions of seed garlic bulbs and cold stress response in terms of leaf-level growth and canopy-level mortality. We used Cropbox modeling framework for reimplementation to take advantage of its declarative modeling language and unified interface for coordinating a large batch of simulations with minimum configurations (Yun et al., 2020). For model testing, a new parameter set was specifically calibrated for Namdo (ND) cultivar (Allium sativum “Namdo”) and validated with multiple datasets with biomass measurements. For demonstrating the yield estimation capability of the model, a climate adaptation strategy was assessed by model simulations for the same cultivar grown in two locations of South Korea. ND is a cultivar originally adapted to the warm climate in the southern region of Korea (Kim et al., 2009). However, the boundary of the growth region between northern and southern types of garlic has moved northward in the past decades due to the warming climate (Heo et al., 2006). Optimal planting dates for achieving maximum yield were discovered through model simulation under current and future climate conditions in the two locations where one is an already established region for growing southern-type garlic and the other has the potential to become a new establishment in the future.

2. Materials and Methods

2.1. Garlic Model

The garlic growth model was extended from a process-based model for leaf development and growth in hardneck garlic (Hsiao et al., 2019). The original model was capable of simulating leaf area expansion at an individual leaf level and estimating carbon assimilated in a canopy calculated by coupled gas exchange, but assessing biomass allocated into a particular organ, i.e., bulb, was not a primary subject of the model at the time. For realistic yield prediction under future scenarios, the model should be able to be reliable under diverse environmental conditions (IPCC, 2014).

Model changes made for the experiments reported in this paper include improved biomass allocation, a dynamically adjusted phyllochron, and cold stress response (Figure 1). For the sake of technical convenience, model code originally written in C++ was reimplemented with Cropbox modeling framework using Julia programming language (Cropbox.jl, 2021). The Cropbox framework allowed a streamlined model development from the model description in a concise declarative form, iterative parameter adjustments within a notebook environment, and batch simulation of large sets of parameters and production of figures included in this paper. The model source code is also available in the public repository (Garlic.jl, 2021).

FIGURE 1

Figure 1. Diagram of key processes updated in the model presented in this paper. Leaf appearance is driven by an accumulation of leaf tip appearance rate (LTAR) which is modified according to dynamic phyllochron (^*). Cold stress (^†) is implemented in two folds; individual leaf-level growth driven by leaf elongation rate (LER) and plant population-level response accounted by plant density (PD). Carbon assimilated through coupled gas exchange forms a carbon pool which determines the amount of carbon partitioned into individual organs. The bulb carbon constitutes a final yield of garlic. For a complete structure of the model, refer to Figure 1 of Hsiao et al. (2019).

2.1.1. Biomass Allocation

The total amount of carbon assimilation was calculated by a C₃ photosynthesis model coupled with a stomatal conductance model and energy balance model as described in the previous paper (Hsiao et al., 2019). Note that the overall model structure including the gas-exchange module was reorganized for taking advantage of domain-specific language provided by Cropbox modeling framework (Yun et al., 2020). The assimilated carbon accumulates and forms a carbon pool (g) ready to be distributed to plant organs. A potential allocation rate of available carbon (g d⁻¹) is driven by the carbon supply rate from the pool (g d⁻¹) excluding maintenance respiration (g d⁻¹). After taking account of carbohydrate synthesis efficiency (Y_g), an actual allocation rate (g d⁻¹) is determined and split into a set of allocation rates for structural organs according to a partitioning table (Figure 2).

FIGURE 2

Figure 2. Partitioning coefficients for plant organs dynamically adjusted according to developmental stage. V, vegetative stage; R1, reproductive stage between scape initiation and appearance; R2, reproductive stage between scape appearance and removal; R3, reproductive stage after scape removal. R3 is not used if the scape is not removed during simulation.

The partitioning table is a 2-dimensional array where each row represents a developmental stage and column represents a destination. Developmental stages span from seed, vegetative, bulb growth before scape appearance, bulb growth after scape appearance, bulb growth after scape removal, and death. Partitioning destinations include root, leaf, sheath, scape, and bulb. For each time step of the simulation, an actual allocation rate weighted by a partitioning coefficient found in the table was used for biomass accumulation of each destination organ.

We determined the initial partitioning coefficients teleonomically for each phenological stage based on our understanding of the crop and literature reviews (Takagi, 1989; Meredith, 2008). We then inspected the model behavior visually and heuristically to adjust the partitioning coefficients to mimic the observations.

2.1.2. Dynamic Phyllochron

An interval between leaf appearance, phyllochron, is not necessarily static throughout plant growth, but can dynamically change depending on the growth condition. However, the maximum leaf tip appearance rate (LTAR_max; d⁻¹) in the original garlic model was determined by thermal time based on storage duration between harvest and planting as well as storage temperature for which seed garlic has been kept during this period. Once initialized, LTAR_max would have remained the same until the end of simulation. This assumption was often held when storage duration was close to an average duration where curve fitting was originally done for. However, sometimes the leaf tip appearance rate may have stayed too high at the end of the growing season when storage duration was longer than usual. The opposite would happen when the storage duration was too short. In other words, scape appearance in the reproductive stage, which is driven by the same mechanism relying on phyllochron that assumes a scape appears after three phyllochrons since the onset of the reproductive stage, could become too sensitive to the initial value of LTAR_max. To alleviate this issue, LTAR_max was adjusted to converge toward half of the asymptote of the maximum leaf tip appearance rate parameter (LTAR_max,a; d⁻¹). The initial value LTAR_max,s is the same as maximal rate of leaf tip appearance rate modified by storage duration SD defined in the original model Hsiao et al. (2019). SD_m is the storage duration that results in half of maximal leaf tip appearance asymptote (LTAR_max,a) and α controls the steepness of the sigmoidal storage function. The same value of parameters was adopted from the original model (Table 1). After each leaf appearance, LTAR_max would linearly increase or decrease depending on the initial rate and the slope toward the asymptote. Convergence is done when leaf rank k reaches the generic leaf number N_g which was set to 10 by default.

$\begin{array}{ll}{r}_0={\text{LTAR}}_{\text{max,s}}=\frac{{\text{LTAR}}_{\text{max,a}}}{1+e^{-\alpha (\text{SD}-{\text{SD}}_m)}}& (1)\end{array}$

$\begin{array}{ll}{r}_1=\frac{{\text{LTAR}}_{\text{max,a}}}{2}& (2)\end{array}$

$\begin{array}{ll}{\text{LTAR}}_{\text{max}}=r_0+(r_1-r_0)\frac{k}{N_g}& (3)\end{array}$

TABLE 1

Table 1. Parameters for dynamic phyllochron.

2.1.3. Cold Stress

Two types of cold stress response commonly found with garlic plants grown in the field during the winter season were added to the model: cold injury and cold damage. Cold injury represents impeded leaf growth under below normal air temperatures. The actual leaf elongation rate (LER; cm d⁻¹) is scaled down from potential leaf elongation rate (LER_p; d⁻¹) by cold injury effect [E(t)] at time t. The potential cold injury effect (E) is based on an apparent cold injury effect at the moment (C) attenuated by the period of cold the plant has experienced (D). The number of cold days (D) increases when the air temperature (T) is below a critical temperature for cold injury (T_c,i; °C) which assumed to be 0 °C. Once the temperature rises above the threshold T_c,i, the cumulative days D decreases until it resets to zero. The actual cold injury effect at the time [E(t)] is preserved to be less than or equal to the degree of effect experienced previously [E(t − 1)] while the same episode is in effect (D > 0). It assumes that the effect of cold injury is prolonged for a certain period of time while recovering under normal temperature even after the cold episode is over. The longer a plant is exposed to the cold temperature, the more difficult it is to recover from the injury.

The potential cold injury effect (E) was derived from a growth chamber experiment (Supplementary Figure S1). ND cultivar was planted and grown in pots at room temperature until the third leaf emerged. Plants were then placed in a growth chamber subject to a constant temperature of 0°C, –5°C, –10°C, and –15°C, respectively. The leaf area of the third leaf blade from surviving plants was measured multiple times for 8 days. Each treatment had three replicates. For each sampling point, normalized leaf area was calculated from the final leaf area measured under room temperature before the onset of treatment. The size of the potential cold injury effect (E) was assumed to be the same as the normalized leaf area at a given time (Supplementary Figure S2). A range of apparent cold injury effects (C) was then obtained from each line of fitted E (Supplementary Figure S3). With the obtained values, C is now formulated as a logarithmic function with two parameters (a, b) (Table 2).

$\begin{array}{ll}d=\{\begin{array}{cc}1,& T<T_{c,i}\\ -1,& T\ge T_{c,i}\end{array}& (4)\end{array}$

$\begin{array}{ll}D={\displaystyle \sum }max(d,0)& (5)\end{array}$

$\begin{array}{ll}C=max\left(0,min\left(log\left[a·\left(T-T_{c,i}\right)+b\right],1\right)\right)& (6)\end{array}$

$\begin{array}{ll}E=1-\frac{C}{e^{\frac{1}{D}}}& (7)\end{array}$

$\begin{array}{ll}E(t)=\{\begin{array}{cc}1,& D=0\\ min(E,E(t-1)),& D>0\end{array}& (8)\end{array}$

$\begin{array}{ll}\text{LER}=E(t)·{\text{LER}}_p& (9)\end{array}$

Cold damage represents the death of the plant under more extreme temperatures, resulting in a reduced plant density (PD). The current plant density (PD) starts from an initial planting density (PD₀) and then may decrease over time as the air temperature (T) drops below a critical temperature for unrecoverable cold damage (T_c,d). Unlike cold injury where leaf growth only slows down under low temperature and eventually recovers once the temperature rises back up, cold damage leads to permanent wilting from which plants no longer recover. Mortality due to cold damage (M) is a logistic curve representing a relative portion of plants that survived at the end of the cold damage treatments. Cold damage shape factor (s) determines how quickly the plants will die off under such extreme conditions. Due to the lack of experimental datasets available for parameter fitting, we instead looked up episodes of cold damage and corresponding temperatures reported in the newspaper during the past years to determine parameter values. The model assumed mortality (M) started building up when the temperature went below around −10 °C and quickly reached peak damage at around −20 °C (Supplementary Figure S4).

$\begin{array}{ll}M=\frac{e^{-s(T-T_{c,d})}}{1+e^{-s(T-T_{c,d})}}& (10)\end{array}$

$\begin{array}{ll}S=1-M& (11)\end{array}$

$\begin{array}{ll}\text{PD}=S·{\text{PD}}_0& (12)\end{array}$

TABLE 2

Table 2. Parameters for cold stress response.

2.2. Parameter Estimation

A parameter set for ND cultivar was mostly based on existing parameter sets calibrated for Korean Mountain (KM) and Shantung Purple (SP). Only a few parameters were modified to reflect how ND grows in general compared to the other cultivars (Table 3). Minimum length of longest leaf (LM_min) was increased to 100 cm which was the largest value found in our datasets. The potential maximum elongation rate (LER_max) was accordingly adjusted to 5.56 cm d⁻¹ assuming full leaf expansion takes 18 days under optimal conditions. Initial leaves at harvest (ILN) was set to 6 as found in a dissected seed garlic clove. Stay green (SG) was set to 1.5 d which was a value calibrated for SP cultivar. Storage temperature (ST) was assumed constant at 8 °C during the entire storage period. Seed bulb was assumed harvested on June 30th and storage duration (SD) was calculated accordingly. Note that it was our intention to select parameters that allow overestimation of biomass and leaf area for tracking potential growth while calibrating parameters to keep the timing of phenology as close as possible to the observation.

TABLE 3

Table 3. Parameters for leaf development of Namdo (ND) cultivar.

We validated our new parameter set for ND cultivar using three datasets. The first dataset (D1) was collected from an experiment plot located at the Research Institute of Climate Change and Agriculture (RICCA), Jeju, South Korea. ND cultivar was planted on October 8th, 2014, and harvested on June 19th, 2015. Growth and development measurements, such as leaf count, leaf area, biomass for each part, were recorded from the mid-vegetative stage until harvest. This dataset was used for evaluating the overall response of the ND parameter set we obtained above. The second dataset (D2) was collected from a temperature gradient greenhouse (TGG) located at the same site in Jeju. TGG is a glass house equipped with heaters on the one side of the wall that keeps a small, but constant gradient in the enclosed planting zones. Five planting zones were set up with 2 °C to 3 °C of temperature differences kept from end to end. ND cultivar was planted on October 7th, 2014, and harvested on May 17th, 2015. Similar measurements to the first dataset were recorded from the early vegetative stage until harvest. This dataset was used for evaluating the response of the model to modest temperature changes in a similar way that the model would be subject to future climate conditions. Scape was not removed in both datasets. The last third dataset (D3) was obtained from a farm field located at Jeongsil (JS) neighborhood in Jeju, South Korea. ND cultivar was planted on September 9th, 2009, and harvested on June 18th, 2010. Types of recorded variables were the same as other datasets. Scape was assumed to be removed shortly after its appearance out of the whorl. Only the visible portion of the scape was cut and measured for biomass. Hourly time-series of weather data for each dataset were obtained from nearby Jeju station (184) operated by Korean Meteorological Administration.

2.3. Future Climate Projections

For assessing garlic yield changes in future climate conditions, we specifically chose two locations in South Korea (Figure 3). The first site was Gosan, Jeju located on a Southern island with a humid subtropical climate off the Korean peninsula, which is in Zone 9b in terms of USDA Hardiness Zones. Gosan is where ND cultivar is currently grown on a large scale for commercial purposes where stakeholders are interested to see how the cultivar would be performing and if any adjustments in growing practice such as if planting date shifts would be required in the long term. The second site was Chuncheon, Gangwon which is located in the middle of the Korean peninsula with a humid continental climate with cold winter. Chuncheon is in Zone 6b according to USDA Hardiness Zones and the current climate is not favorable for growing a Southern cultivar, such as ND, which has been historically more adapted to warm climates. However, it is not clear whether future climate conditions could allow the growing of new crops in a region previously not suited for production. We assumed two RCP (Representative Concentration Pathway) scenarios, RCP4.5 and RCP8.5, for future climate projections (IPCC, 2014). Daily weather dataset for the two locations under RCP4.5 and RCP8.5 scenarios were obtained from the AgClimate data portal (AgClimate, 2019) operated by Rural Development Administration in South Korea. Elevated CO₂ concentrations for RCP scenarios were used accordingly for driving the integrated coupled gas-exchange model. For comparison with current climate conditions of each location, we also obtained 30-years of normal weather data from 1980 to 2010 provided by the same data portal. The daily weather dataset was then scaled down to hourly time-series by applying an interpolation method specifically designed for different types of weather variables (Moon et al., 2019). Elevated CO₂ levels for each year under future scenarios were obtained by linear interpolation between CO₂ abundance projection reported in the Climate System Scenario Tables (IPCC, 2013).

FIGURE 3

Figure 3. Geographic locations of Gosan and Jeju in South Korea. Gosan is located on a Southern island with a humid subtropical climate off the Korean peninsula. Chuncheon is located in the middle of the peninsula with a humid continental climate with cold winter. A color overlay represents the average minimum temperature (°C) in January during past 30 years from 1981 to 2010 provided by the National Institute of Agricultural Sciences, South Korea. According to USDA Hardiness Zones, Gosan is Zone 9b and Chuncheon is Zone 6b.

For each location, multiple runs of the simulation were executed with combinations of treatments. In the case of future climate projections, 8 time windows from the 2020s to 2090s with 10 years interval multiplied by 10 repetitions from different random seeds produced 80 weather datasets for each RCP scenario. Historical weather was organized in a single time window referred to 1980s with 10 random repetitions, leading to 10 weather datasets. The planting date was adjusted from 240th to 350th DOY (day of the year) with 10 days intervals, requiring 12 runs for each weather condition. In turn, the current normal scenario involved 120 runs of simulation while each RCP scenario required 960 runs. Immediate scape removal after appearance was assumed. The harvest date was set fixed to be May 15th considering practical harvest dates of ND cultivar grown in Jeju. The yield was based on fresh biomass of bulb at the harvest date. In the model simulation, fresh biomass was calculated with 85 % moisture content in bulb. The optimal plating date was then determined by finding out a planting date that gives maximum yield for a given year. For the sake of yield comparison between planting dates, years were grouped into three time periods, namely 1980–2010s, 2020–2050s, and 2060–2090s. A data point in the 1980–2010 period was composed of 10 samples originated from a set of stochastic weather data pooled for the entire period whereas a data point in the future climate scenario was composed of 50 samples where the result of 10 stochastic weather datasets was combined with 5 periods of 10 years interval.

3. Results

3.1. Model Evaluation

3.1.1. Dataset 1 (RICCA Field)

Our parameter for ND cultivar was first evaluated with measurements from field grown garlic in the dataset D1. Green leaf area from the model exhibited consistent overestimation with an average RMSE (root mean square error) of 275 cm² while keeping a similar trajectory of leaf expansion and senescence pattern to the observation (Figure 4). The maximum green leaf area was estimated to be 1.193 cm² on April 28th while the maximum value in the dataset was 1.150 cm² observed on April 30th. The peak was reached in simulation shortly after the scape appeared on April 23rd during early vegetative growth.

FIGURE 4

Figure 4. Comparison of model simulated green leaf area and observation from dataset D1 which was from field measurements and primarily used for evaluating new parameters.

Simulated biomass allocation to each part of plant structure, namely bulb, leaf, stalk, and scape, was close to the observation with an average RMSE of 5.0, 1.6, 2.5, and 1.7 g, respectively (Figure 5). The maximum biomass of bulb from the simulation was 22.6 g reached on June 14th while the observed maximum was 26.0 g measured on June 4th which was the final date recorded in the dataset. The maximum living leaf biomass from the simulation was 8.4 g reached on April 27th and the maximum observed value was the same 9.9 g on April 30th. The maximum biomass of stalk, which consists of leaf sheath and scape, from the simulation was 15.3 g close to the observed maximum of 15.9 g. However, the timing of the peak was shifted by almost 3 weeks later on May 27th compared to the observed May 7th. Since leaf sheath withered away at the end of the growth stage, the only remaining part of the living stalk would be a scape which was the reason why a convergence between two values occurred in mid-June. The simulated biomass of the scape reached its maximum value of 9.5 g on June 20th which was almost identical to the observed maximum which, however, occurred earlier on June 4th. Note that the two slopes of biomass increase for simulation and observation were almost parallel to each other.

FIGURE 5

Figure 5. Biomass allocation to leaf, bulb, stalk, and visible portion of scape simulated over time compared with measurements from dataset D1 which was from field measurements and primarily used for evaluating new parameters.

Leaf development phenology simulated by the model was in agreement with fresh leaf count recorded in the dataset with an average RMSE of 1.6 leaves (Figure 6). According to the model, leaf initiation was completed on March 22nd with the total number of 16 leaves and the 16th leaf appeared on April 2nd. The final leaf was completely senesced and dropped after about 2.5 months of growth on June 15th.

FIGURE 6

Figure 6. Simulated leaf development phenology compared with fresh leaf count recorded in dataset D1 which was from field measurements and primarily used for evaluating new parameters.

3.1.2. Dataset 2 (RICCA TGG)

Plants grown in higher temperature zone 1 had a larger green leaf area of 1.521 cm² and reached its peak earlier on April 20th when compared to plants grown in lower temperature zone 5 whose maximum green leaf area was estimated to be a smaller 1.303 cm² reached on May 9th (Figure 7). The simulation result was in agreement with observation recorded in dataset D2 in which zone 1 had a larger leaf area of 1.176 ± 0.239 cm² on April 17th, in contrast to zone 5 having a smaller leaf area of 1.124 ± 0.090 cm² on the same date. After 3 weeks on May 9th, the difference was flipped over where plants in zone 1 had senesced faster, resulting in a smaller green leaf area of 512 ± 55 cm² while zone 5 was still greener with 701 ± 70 cm². The mean temperature recorded in dataset D2 for zone 1 was 15.5 °C and zone 5 was 12.6 °C.

FIGURE 7

Figure 7. Comparison of model simulated green leaf area and observation for two temperature zones (Zone 1: high, Zone 5: low) recorded in dataset D2 from a temperature-controlled glass house.

3.1.3. Dataset 3 (JS Field)

Simulated biomass allocation compared with another independent dataset D3 came with an average RMSE of 7.4, 3.6, and 4.2 g for bulb, leaf, and stalk, respectively (Figure 8). The maximum bulb biomass from the simulation was 37.4 g reached on June 22nd while the observed maximum was 28.7 g measured on May 26th. The maximum living leaf biomass from the simulation was 10.0 g reached on May 1st and the maximum observed value was 9.1 g on April 28th. The maximum biomass of stalk from the simulation was 13.5 g on May 1st when scape just appeared and was subsequently removed. The maximum stalk biomass recorded in the dataset was 12.3 g observed on April 28th. Once the scape was removed and no longer a part of stalk composition, stalk biomass gradually decreased as the mature sheath stopped growing and started senescence. According to the dataset, scape was removed during the reproductive stage, but an exact date of the removal was not recorded. Thus, we assumed that scape removal took place as soon as a tip of the scape became visible out of the whorl which occurred on May 1st by model estimation. The scape biomass at the time of removal was 5.7 g.

FIGURE 8

Figure 8. Biomass allocation to leaf, bulb, stalk, and visible portion of scape simulated over time compared with measurements from dataset D3 which was from field measurements with scape removal.

3.2. Yield Projection

3.2.1. Gosan

In the current climate condition, the fresh yield on ND cultivar in Gosan, Jeju was estimated to be maximum at an average of 6.8 kg m⁻² with an SD of 0.8 kg m⁻² when planted in late August (240 DOY) which was closely followed by early September planting (260 DOY) with an average of 6.7 ± 0.7 kg m⁻² (Figure 9). A similar level of high yield was maintained until late September then yield gradually decreased with later planting dates. The difference between maximum yield from early planting date and minimum yield from later planting date was 3.7 kg m⁻². In the near future from the 2020 to 2050s, a similar pattern was observed that high yield was achieved with early planting in September. An overall yield was increased to 7.5 ± 0.8 kg m⁻² under RCP4.5 and 7.6 ± 0.8 kg m⁻² under the RCP8.5 scenario when compared to the current climate. There was no clear difference between RCP4.5 and RCP8.5 scenarios in terms of estimated yield for a given planting date. In the distant future from the 2060s to the 2090s, the RCP4.5 scenario still maintained a similar pattern where a high yield was achieved in early planting dates, but then the range of potential high yield was expanded to later planting dates. For example, planting dates from late August (240 DOY) to mid-October (290 DOY) all resulted in an average yield that closely ranged from 7.6 to 7.7 kg m⁻². The difference between the maximum and minimum yield estimated in the range of planting dates was reduced to 3.0 kg m⁻². Under the RCP8.5 scenario, the yield curve was more flattened that the difference was only 1.8 kg m⁻² between all planting dates. The maximum yield was 8.3 ± 0.6 kg m⁻² in early-October (280 DOY).

FIGURE 9

Figure 9. Yield estimation within a range of planting dates under current and future climate projection in Gosan, Jeju which is a region that grows ND cultivar at a commercial scale. The current climate indicates simulation results with 30-years of normal weather data from 1980–2010. Shades represent a range of ±1 standard deviation from the mean fresh yield.

When optimal plating dates were assessed from current to future climate conditions by 10-year intervals, a clear trend of shifting toward later planting date was found in the future projection although its strength varied depending on scenarios (Figure 11). Both scenarios began with optimal planting date in early to mid-September then showed a strong divergence at the end of the century that the optimal planting date surfaced from late September to early October under the RCP4.5 scenario, while mid-November became a possibility under the RCP8.5 scenario.

3.2.2. Chuncheon

In the current climate condition, the yield of ND cultivar in Chuncheon, Gangwon was estimated close to nil with a maximum yield of 0.06 ± 0.05 kg m⁻² due to impeded leaf growth and increased mortality from cold winter temperature during early vegetative growth (Figure 10). However, under future projections, yield became more viable thanks to a warmer climate. In the near future from the 2020 to 2050s, the maximum estimated yield was 1.5 ± 0.7 kg m⁻² under RCP4.5 and 1.3 ± 0.9 kg m⁻² under RCP8.5 scenario when both planted in early September (250 DOY). In the distant future from the 2060 to 2090s, more yield was achievable with 2.5 ± 0.8 kg m⁻² under RCP4.5 and 4.3 ± 1.1 kg m⁻² under RCP8.5 scenario with early to mid-September plating dates, 250 DOY and 240 DOY, respectively. The difference between the maximum and minimum yield estimated in the range of planting dates was 0.6 kg m⁻² in the early period and 1.0 kg m⁻² in the late period under RCP4.5 scenario. RCP8.5 scenario initially had a similar range of estimated yield over multiple planting dates with 0.5 kg m⁻² in the early period, but then showed a much higher variance of 1.7 kg m⁻² in the late period.

FIGURE 10

Figure 10. Yield estimation within a range of planting dates under current and future climate projection in Chuncheon, Gangwon where the current climate is not favorable for growing ND cultivar. The current climate indicates simulation results with 30-years of normal weather data from 1980–2010. Shades represent a range of ±1 SD from mean fresh yield.

In terms of optimal planting dates for each period, planting in around early September consistently turned out to work best regardless of climate scenarios and estimated yield gradually declined with later planting dates in all scenarios (Figure 11).

FIGURE 11

Figure 11. Optimal planting dates estimated for each year under current and future climate projections in Gosan and Chuncheon. The year 2010 indicates simulation results with 30-years normal weather data from 1980–2010 representing current climate. For RCP scenarios, each year indicates a 10-year average, e.g., 2090 refers to 2090–2100. Error bar and shades represent a range of ±1 SD from mean optimal planting date.

4. Discussion

4.1. Model Evaluation

4.1.1. Parameter Calibration

Our parameter set for ND cultivar was mostly derived from existing parameters for Korean Mountain (KM) and Shantung Purple (SP) with small changes (Hsiao et al., 2019). The final parameter set used in our simulation turned out to be similar to one calibrated for the SP cultivar. It would not be unexpected given that ND was originated from a Chinese cultivar and hence there is a chance SP might be also related to this original cultivar due to geographical proximity (Kim et al., 2009).

4.1.2. Overestimation

Simulation with the calibrated parameter set often resulted in overestimation in the sense that model output tends to touch upper boundaries of observed data points (Figures 4, 5, 7, and 8). This overestimation is intended and expected because our model estimates for biomass and yield potential for a given environmental condition with assuming no limiting factors such as water and nutrient while the field conditions are likely less than optimal. Some types of overestimation also presumably came from the difference between physical measuring methods and modeled algorithms. For example, the green leaf area recorded in the dataset was measured by taking only non-withered leaves at the time of observation which is prone to lose leaves with some green parts intact. On the other hand, the green leaf area from the model was calculated for individual leaves at a fractional scale by tracking the current green portion of the leaf for each time step. Such difference could have been one of the reasons leading to overestimation in leaf area simulation (Figures 4, 7).

4.1.3. Temperature Regimes

Evaluation of model response under slightly separate temperature regimes as conducted in temperature gradient house provided an insight into how elevated air temperature under RCP scenarios would affect plant growth (Figure 7). In short, the plant will grow larger and faster and die earlier in warmer conditions. An average of 3 °C difference in temperature led to more than 2 weeks of the shift in growth peak and 15 % change in total green leaf area. Subsequent changes in senescence timing would imply a need for finding out optimal harvest dates which we assumed as a rather constant management decision in our simulation.

4.1.4. Storage Condition

Storage duration (SD) affects maximum leaf tip appearance rate (LTAR_max) via the dynamic phyllochron model (Equations 1, 2, 3). Storage temperature (ST) along with storage duration then decides number of leaves initiated inside a seed bulb at the time of planting. Despite the importance on phenology, experimental records rarely exist for when seed bulbs were harvested and in which condition they were stored until the date of planting. With controlled experiments on storage conditions, we could have a better understanding of how leaves are initiated during the storage period and whether more sophisticated approaches like a dynamic plastochron would be necessary.

4.1.5. Other Environmental Cues

While our assessment for yield projection under future condition was primarily driven by model response to temperature regimes, there is still room for considering other environmental cues. For instance, we did not have a separate vernalization process in the model but assumed winter chilling was always at an adequate level for triggering continuous development in the following warm spring. While it worked reasonably well in most cases, some processes like mortality from cold stress could capitalize on this additional cue to make dormant plants less vulnerable to cold damage, for instance. Soil water availability would be another important cue in terms of assessing irrigation requirements under future climate conditions given that garlic plots are often irrigated to prevent water deficiency during the reproductive stage.

4.2. Yield Projection

4.2.1. Gosan

Jeju, including Gosan, is a region where southern-type garlic cultivars like ND have been grown at a commercial scale thanks to the warm climate. Northern regions in the Korean peninsula often have colder winters which prevents growing southern-type cultivars. Growers instead choose northern-type cultivars for cold hardiness. As southern cultivars have the advantage of higher yield and shorter growing seasons, whether they could be adapted to northern regions has been an important question to many growers and stakeholders. According to our simulation results using a process-based model, ND cultivar would continue growing well with a slightly higher yield (Figure 9). The optimal planting date in early September estimated for the current climate condition was indeed close to the current practice held in Jeju. This practice is, however, expected to be delayed in the near future if growers keep trying to maximize yield under both climate scenarios (Figure 11). The amount of planting date delay depends on the scenarios and the model estimated that the RCP4.5 scenario could see a shift less than a month while the RCP8.5 scenario could potentially result in 2 months of late planting. Note that although optimal planting dates were seemingly changed drastically, yield estimation curve vs. planting date was actually getting more flattened and the difference between the maximum and minimum yield became smaller in warmer climate conditions. Such changes imply that growers would have more options to choose planting date better suits their own needs. For example, October planting in the future would still promise a yield level close to the potential maximum while taking a relatively shorter growing season to help reduce the overall cost of labor and resources. High yield does not necessarily lead to high profitability when a longer growing season increases overall water demand entailed by higher irrigation cost for compensation (Lobell, 2014). Looking further in the distant future, it was more clear to expect a shift in optimal planting to later dates although yield difference to earlier planting dates was not significant.

4.2.2. Chuncheon

Chuncheon is located in the middle of Korean peninsula and its continental climate frequently experiences below freezing temperatures during the winter season. Hence, only northern-type cultivars, adapted to cold climates have been grown in this region. Our simulation results confirmed expected yield for the current climate was almost non-existent as most plants subject to mortality (M) could not survive after episodes of extreme cold. By contrast, future climate conditions were projected to become more favorable in both scenarios so that expecting some tangible yields would be at least feasible (Figure 10). Generally, early planting dates were favored in terms of optimal yield as similar in the current growing region like Jeju (Figure 11). Although the level of estimated yield was still lower when compared to the yield currently obtainable in other established regions, growers may be able to take advantage of this new opportunity for expanding crop portfolio during the winter season.

4.3. Climate Adaptation

It is important to tease apart an effect of adaptation introduced by planting date shift in the future from an impact of planting date shift in the current condition (Lobell, 2014; Challinor et al., 2018). According to the result from Gosan, planting date shift in the established growing region has little effect on adaptation as evidenced by a similar or smaller range between maximum and minimum yield when comparing the current and future climate scenarios (Figure 9). Higher estimated yield was a result of increased productivity under more favorable conditions primarily due to the higher temperature and elevated CO₂, but not from potential phenological changes that occurred by planting date shift. In non-established regions like Chuncheon, however, an effect of adaptation increased virtually from zero to a substantial amount when transitioning to new climate conditions as evidenced by a rising slope of yield estimation curve (Figure 10). In other words, the importance of planting date as a climate adaptation strategy depends on how crops are adapted to the local condition.

Underlying uncertainties of crop models often hinder model driven assessment of climate change adaptation (Rosenzweig et al., 2014; Holzkämper et al., 2015; Corbeels et al., 2018). We adopted a process-based model with coupled gas-exchange to minimize uncertainties in physiological responses to temperature and CO₂, but still there remain many processes that can use improvements. For instance, the water stress response was solely dependent on leaf water potential disconnected with soil water status albeit a less practical implication when garlic plants are generally irrigated to meet high water demand during the reproductive stage. Nitrogen response is another important factor that the model currently does not take into account and instead assumes non-limiting due to extensive use of fertilizers in practice, but will be critical for assessing an exact cost of production in finding optimal yield.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

KY and S-HK contributed to conception and design of the study, performed the modeling and analysis, and wrote sections of the manuscript. KM and MS collected and organized the dataset. KY wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

The information, data, or work presented herein was funded in part by the Cooperative Research Program for Agricultural Science and Technology Development, Rural Development Administration, Republic of Korea under grant no. PJ015124012022.

Author Disclaimer

The views and opinions of the authors expressed herein do not necessarily state or reflect those of the funding agency.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.783810/full#supplementary-material

References

Abdalla, M. M. A., Aboul-Nasr, M. H., and Metwaly, A. K. (2011). Growth and yield of fifteen garlic ecotypes. Assiut J. Agric. Sci. 36, 52–65.

AgClimate (2019). Available online at: https://agclimate.epinet.kr

Challinor, A. J., Müller, C., Asseng, S., Deva, C., Nicklin, K. J., Wallach, D., et al. (2018). Improving the use of crop models for risk assessment and climate change adaptation. Agric. Syst. 159, 296–306. doi: 10.1016/j.agsy.2017.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, S., and Baek, J. (2016). Garlic yields estimation using climate data. J. Korean Data Inf. Sci. Soc. 27, 969–977. doi: 10.7465/jkdi.2016.27.4.969

CrossRef Full Text | Google Scholar

Corbeels, M., Berre, D., Rusinamhodzi, L., and Lopez-Ridaura, S. (2018). Can we use crop modelling for identifying climate change adaptation options? Agric. Forest Meteorol. 256, 46–52. doi: 10.1016/j.agrformet.2018.02.026

CrossRef Full Text | Google Scholar

Cropbox.jl (2021). Available online at: https://github.com/cropbox/Cropbox.jl

FAOSTAT (2020). Available online at: https://www.fao.org/faostat

Garlic.jl (2021). Available online at: https://github.com/cropbox/Garlic.jl.

Giménez, C., Stöckle, C., Suárez-Rey, E., and Gallardo, M. (2016). Crop yields and N losses tradeoffs in a garlic-wheat rotation in southern spain. Eur. J. Agron. 73, 160–169. doi: 10.1016/j.eja.2015.11.016

CrossRef Full Text | Google Scholar

Gómez, D., Salvador, P., Sanz, J., and Casanova, J. L. (2021). Regional estimation of garlic yield using crop, satellite and climate data in mexico. Comput. Electron. Agric. 181, 105943. doi: 10.1016/j.compag.2020.105943

CrossRef Full Text | Google Scholar

Heo, I., Kwon, W.-T., Chun, Y., and Lee, S. (2006). The impact of temperature rising on the distribution of plant – in case of bamboos and garlic. J. Environ. Impact Assess. 15, 67–78.

Google Scholar

Holzkämper, A., Calanca, P., Honti, M., and Fuhrer, J. (2015). Projecting climate change impacts on grain maize based on three different crop model approaches. Agric. Forest Meteorol. 214, 219–230. doi: 10.1016/j.agrformet.2015.08.263

CrossRef Full Text | Google Scholar

Hsiao, J., Yun, K., Moon, K. H., and Kim, S.-H. (2019). A process-based model for leaf development and growth in hardneck garlic (Allium sativum). BAnn. Bot. 124, 1143–1160. doi: 10.1093/aob/mcz060

PubMed Abstract | CrossRef Full Text | Google Scholar

IPCC (2013). Climate change 2013: The physical science basis. contribution of working group i to the fifth assessment report of the intergovernmental panel on climate change. Technical report.

Google Scholar

IPCC (2014). Climate change 2014: Synthesis report. contribution of working groups i, ii and iii to the fifth assessment report of the intergovernmental panel on climate change. Technical report.

Google Scholar

Kamenetsky, R. (2007). BGarlic: Botany and Horticulture, Chapter 2, John Wiley & Sons, Ltd.

Google Scholar

Kim, C.-H., Seong, K.-C., Lee, J.-S., Kang, K.-H., Um, Y.-C., and Suh, H.-D. (2009). Production of seed garlic by sawing bulbils of southern type garlic in Jeju island. Protect. Hortic. Plant Factory 18, 74–80.

Google Scholar

Kim, S.-H., Jeong, J. H., and Nackley, L. L. (2013). Photosynthetic and transpiration responses to light, CO₂, temperature, and leaf senescence in garlic: analysis and modeling. J. Am. Soc. Hortic. Sci. 138, 149–156. doi: 10.21273/JASHS.138.2.149

CrossRef Full Text | Google Scholar

Lee, J.-W., Park, G.-A., Joh, H.-K., Lee, K.-H., Na, S.-I., Park, J.-H., et al. (2011). Analysis of relationship between vegetation indices and crop yield using KOMPSAT (KoreaMulti-Purpose SATellite)-2 imagery and field investigation data. J. Korean Soc. Agric. Eng. 53, 75–82. doi: 10.5389/KSAE.2011.53.3.075

CrossRef Full Text | Google Scholar

Lobell, D. B. (2014). Climate change adaptation in crop production: beware of illusions. Glob. Food Security 3, 72–76. doi: 10.1016/j.gfs.2014.05.002

CrossRef Full Text | Google Scholar

Marcelis, L., Heuvelink, E., and Goudriaan, J. (1998). Modelling biomass production and yield of horticultural crops: a review. Sci. Hortic. 74, 83–111. doi: 10.1016/S0304-4238(98)00083-1

CrossRef Full Text | Google Scholar

Meredith, T. (2008). The Complete Book of Garlic. Portland, OR: Timber Press.

Google Scholar

Moon, K., Seo, H., Shin, M., Song, E., and Oh, S. (2019). “Generation of hourly weather data using daily data,” in Proceedings of The Korean Society of Agricultural and Forest Meteorology Conference (Jeju), 195–196.

Google Scholar

Nackley, L. L., Jeong, J. H., Oki, L. R., and Kim, S.-H. (2016). Photosynthetic acclimation, biomass allocation, and water use efficiency of garlic in response to carbon dioxide enrichment and nitrogen fertilization. J. Am. Soc. Hortic. Sci. 141, 373–380. doi: 10.21273/JASHS.141.4.373

CrossRef Full Text | Google Scholar

Põldma, P., Merivee, A., Pae, A., and Justus, K. (2005). Influence of planting time on the development, yield and quality of garlic (Allium sativum L.) in Estonia. Acta Hortic. 688, 333–338. doi: 10.17660/ActaHortic.2005.688.49

CrossRef Full Text | Google Scholar

Portela, J. A., Sidoti, B., Reybet, G., Bellaccomo, C., and Astorquizaga, R. (2012). “Yield stability of ten garlic (Allium sativum) clonal cultivars in northern Patagonia, Argentina,” in VI International Symposium on Edible Alliaceae, Vol. 969 (Fukuoka), 107–112.

Google Scholar

Rizzalli, R., Villalobos, F., and Orgaz, F. (2002). Radiation interception, radiation-use efficiency and dry matter partitioning in garlic (Allium sativum L.). BEur. J. Agron. 18, 33–43. doi: 10.1016/S1161-0301(02)00094-1

CrossRef Full Text | Google Scholar

Rosenzweig, C., Elliott, J., Deryng, D., Ruane, A. C., Müller, C., Arneth, A., et al. (2014). Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison. Proc. Natl. Acad. Sci. U.S.A. 111, 3268–3273. doi: 10.1073/pnas.1222463110

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, H. (1989). Onions and Allied Crops, Volume 3, Chapter Garlic Allium sativum L., Boca Raton, FL: CRC Press.

Google Scholar

Villalobos, F., Testi, L., Rizzalli, R., and Orgaz, F. (2004). Evapotranspiration and crop coefficients of irrigated garlic (Allium sativum L.) in a semi-arid climate. Agric. Water Manage. 64, 233–249. doi: 10.1016/S0378-3774(03)00198-7

CrossRef Full Text | Google Scholar

Yun, K., Timlin, D., and Kim, S.-H. (2020). Coupled gas-exchange model for C₄ leaves comparing stomatal conductance models. Plants 9, 1358. doi: 10.3390/plants9101358

PubMed Abstract | CrossRef Full Text | Google Scholar