Developing perennial fruit crop models in APSIM Next Generation using grapevine as an example

A new model for grapevines ( Vitis vinifera ) is the first perennial fruit crop model using the Agricultural Production System sIMulator (APSIM) Next Generation framework. Modules for phenology, light interception, carbohydrate allocation, yield formation and berry composition were adapted or added into APSIM Next Generation to represent the nature of fruit-bearing vines. The simulated grapevine phenological cycle starts with the dormancy phase triggered by a critical photoperiod in autumn, and then goes through the subsequent phenophases sequentially and finally returns to dormancy for a new cycle. The canopy microclimate module within APSIM Next Generation was extended to allow for row crop light interception. The carbohydrate arbitrator was enhanced to consider both sink strength and sink priority to reflect carbohydrate reserve as a concurrent competing sink. Weather conditions and source–sink ratio at critical developmental stages were used to determine potential grapevine yield components, e.g. bunch number, berry number and berry fresh weight. The model was calibrated and tested extensively using four detailed data sets. The model captured the variations in the timing of measured budburst, flowering and véraison over 15 seasons across New Zealand for five different varieties. The calculated seasonal dynamics of light interception by the row and alley were consistent with field observations. The model also reproduced the dynamics of dry matter and carbohydrate reserve of different organs, and the wide variation in yield components caused by seasonal weather conditions and pruning regimes. The modelling framework developed in this work can also be used for other perennial fruit crops.


IN TROD UCTION
Perennial crops represent a long-term investment by any landowner.There are fewer opportunities to change location, genotype and plant configuration to adapt to climate compared with annual crops.As such, reliable models for evaluating options at establishment and for ongoing management would be valuable aids for decision-making.However, in contrast with annual crops where a large amount of modelling work has been undertaken on assessing the yield and adaptation options under climate change (Schlenker and Roberts 2009;Challinor et al. 2014), less work has been undertaken on the perennial crops (Santos et al. 2011;Bock et al. 2013).The research on perennial crop models lags behind the annual crop models largely owing to the physiological complexity that regulate yield formation and the difficulties and high cost of data collection.
Grapevine (Vitis vinifera) is one of the most economically important perennial fruit crops worldwide, used both as table grapes and for wine production (OIV 2020).It has played an important cultural role in many parts of the world.Because of its global economic importance and being found in diverse growing climates, the grapevine has emerged as a model perennial fruit crop species.There have been a large number of scientific studies from genomics to production practices carried out on grapevines over the years.The grapevine reaches reproductive maturity in 4-5 years after planting, and may remain economically productive for 20-60 years depending on management.
The key attributes required of a grapevine model include: (i) A reproductive cycle from inflorescence initiation to harvest, lasting about 15-18 months.Grapevines initiate inflorescence primordia in the axial buds of each shoot in the summer of the year prior to which they flower.(ii) The buds undergo a dormancy phase in winter and this breaks in spring.A period of intense vegetative growth follows after budburst.Initial vegetative growth is supported by the internal carbohydrate reserves until the vine achieves a net carbohydrate gain around flowering.(iii) Vegetative and reproductive growth occur concurrently.The reproductive stage occurs during the flowering process of the inflorescences.Some growers perform a top pruning around this time to reduce the carbohydrate competition between vegetative and reproductive growth.(iv) The number of berries that set and continue to grow is determined shortly after flowering.Fruit clusters gradually become the primary sink for photosynthates during development.
Several functional-structural plant (FSP) models with explicit descriptions of grapevine canopy architecture have been developed.These often focus on certain functional processes such as vegetative growth (Louarn et al. 2008b;Schmidt et al. 2019), light interception (Louarn et al. 2008a), canopy photosynthesis and transpiration (Prieto et al. 2012;Zhu et al. 2018), plant water status (Zhu et al. 2018;Albasha et al. 2019) and the effects of whole-plant carbohydrate and water status on berry water and sugar uptake (Zhu et al. 2019).Those models allow us to better understand the functioning of a grapevine plant, and often the important processes highlighted by FSP models can be summarized by functional relationships for including in other models, e.g.farming system models.A farming system can be defined as an arrangement of land, crops, livestock, other capital goods and labour put together for the primary purpose of producing plant and animal products for consumption (Woodward et al. 2008).Farming systems models contain interconnected biophysical and management modules which are able to simulate agricultural systems comprising soil, crop, tree, pasture and livestock, and have the flexibility to integrate non-biological farm resources such as availability of water or farm machinery (Holzworth et al. 2014).Farming systems models are essential tools for investigating the effects of climate and management on crop and livestock growth, productivity and the system's water and nitrogen balances (Holzworth et al. 2018;Keating 2020).Thus, compared to FSP models and crop models which only focus on crop growth and development, farming system models have advantages in assessing the behaviour at farm level, such as large-scale water and nutrient balances with the consideration of spatial variability of farming types, soil and climate conditions.
A number of farming or cropping system grapevine models have been developed (Moriondo et al. 2015), namely VineLogic (Godwin et al. 2002), STICS grapevine model (Brisson et al. 2003(Brisson et al. , 2009)), Cropsyst (Stöckle et al. 2003), Gutierrez et al. (1985), Williams et al. (1985), VIMO (Wermelinger et al. 1991), Bindi et al. (1997), NVINE (Nendel and Kersebaum 2004), SWAP (Ben-Asher et al. 2006), Vitisim (Lakso et al. 2008), WALIS (Celette et al. 2010).The STICS grapevine model has the most functionalities among all the listed models.It can simulate the light interception by the grapevine rows based on geometry description, water balance of the wholevineyard and the grapevine yield (Valdés-Gómez et al. 2009;Fraga et al. 2016).However, it has some limitations as well.For instance, it does not capture seasonal variations of yield, or consider the availability of carbohydrate on leaf area and root expansion (see details in the Discussion section).Furthermore, a crop model is generally built to address certain objectives of a study (Van Ittersum et al. 2003) and therefore would need to be adapted for other applications.Such an approach may be difficult in existing grapevine farming system models.
The Agricultural Production System sIMulator Next Generation (APSIM Next Generation) (Holzworth et al. 2018) platform has advantages over other modelling platforms, allowing scientists to construct and configure a crop model visually, by dragging and dropping components in the user interface (UI) (Brown et al. 2014).It combines a range of tools in the UI to assist model developers to build, run, test, evaluate and auto-document the model (Fig. 1; see Supporting Information-Fig.S1).The APSIM Next Generation plant modelling framework provides a collection of generic software classes, e.g.plant organs and functions that can be assembled and parameterized in different ways to create different plant models (Brown et al. 2018).It means that dynamic composition of lower-level processes and organ classes (e.g.photosynthesis, leaf) into larger constructions, e.g.wheat model, can be achieved by the model developer without additional coding.Graphs and statistics for assessing model performance can also be constructed in the UI.However, no formal perennial fruit crop model has been developed in APSIM Next Generation using the plant modelling framework.
The aim of this study was to: (i) develop a perennial fruit crop model in APSIM Next Generation using grapevine as a model plant; (ii) improve the carbohydrate allocation scheme for perennial crops; (iii) improve the light interception calculation for row-planted trees/ vines; (iv) improve the representation of phenology, leaf cohort and perennial organ set-ups in APSIM Next Generation in general.The APSIM grapevine model developed in this study can be used to explore the long-term dynamics of phenology, canopy development, and yield and carbohydrate reserves under different pruning systems at different locations, and later be used to assess the environmental footprint of the viticultural practices, e.g.water and nutrient footprint, and the suitability of integrated grapevine-livestock farming systems.

Model description
A series of developments were undertaken to adapt the APSIM Next Generation source code to enable simulation of perennial fruit crops.The purpose of the model presented here was to simulate the phenology, canopy development, light interception, carbohydrate assimilate and allocation, carbohydrate storage and remobilization, water uptake, organ growth and yield formation in a row-planted vineyard (Fig. 2; see Supporting Information-Method S1 APSIM Next Generation grapevine model documentation).The basic principle of the model is that phenology is largely driven by temperature and growth is dynamically driven by both environmental factors and differential allocation of carbohydrate among organs.It includes the following novel features: • A new carbohydrate allocation method to distribute carbohydrate to structural and non-structural organ components simultaneously based on sink strength and priority • A yield module that uses the weather conditions in both the previous and current season as well as the carbohydrate status to determine potential bunch number per shoot, berry number per bunch and berry weight • A row and alley strip configuration to represent a typical vineyard set-up with rows of grape vines growing in a strip of bare ground, interspersed with alleys of short grass [see Supporting Information-Fig.S2].Roots from each strip are present beneath their neighbouring strip and compete for water and nutrients • An adapted row light interception method based on that of Goudriaan (1977) to partition light interception between the row canopy and the alley • Flexibility in simulating different cane-or spur-pruned training systems by controlling the cane or cordon number, bud number per vine, dynamics of canopy height, canopy width and depth in the UI (Fig. 2).New training systems can be added in the UI through management scripts.• An adapted phenology method that defines the start of endo-dormancy by critical photoperiod and then goes through each phenology phase progressively.Budburst was used for triggering canopy growth instead of crop emergence as in annual crops (Fig. 1) • A flexible parameter optimization procedure that calls APSIM Next Generation through R bash command and optimizes model parameters through using an R package.

Phenology module.
The phenology module starts with dormancy phase of the buds.This is followed by budding, flowering, fruit setting, berry development and leaf death phases (Fig. 1).The dormancy phase include two stages: endo-dormancy and eco-dormancy.Endo-dormancy is defined as the period when buds are dormant because of physiological conditions, and eco-dormancy is defined when buds remain dormant because of unfavourable environmental conditions (Sarvaš 1974).The budburst stage in the budding phase occurs when eco-dormancy is broken.For consistency in calculating the phenology stage in both the Southern and Northern Hemispheres, we use days after winter solstice (DAWS) as the standard time in the model.We assume the endo-dormancy stage of the bud starts in autumn when photoperiod drops below a critical value (Camargo-Alvarez et al. 2020).The critical value depends on grape variety and was set to 13.1 h in the current model based on the findings of Camargo-Alvarez et al. (2020).Three budburst approaches were tested using a phenology data set: (i) forcing-only approach, which only describes the release of eco-dormancy; (ii) chilling-forcing, which describes the release of endo-and eco-dormancy in sequence; (iii) chilling-forcing-overlap, which describes the release of endo-and ecodormancy in sequence plus the effect of chilling units on forcing accumulation through implementation of an exponential decrease function.The chilling-forcing approach was selected in the current grapevine model because it enabled the best predictions of the timing of budburst.In this approach, endo-dormancy is completed when the chilling requirement is satisfied, and eco-dormancy is completed when sufficient heat units are accumulated to force budburst (Chuine 2000).This approach has two major assumptions: (i) Chilling must be fulfilled before budburst.However, it has been shown that buds with insufficient chilling can still burst, although this often results in a delayed, desynchronized and low percentage of budburst (Dokoozlian 1999).(ii) Temperature is the only driver for budburst while all other factors are not limiting budburst, e.g.carbohydrate reserves and water supply.Bonada et al. (2020) showed that severe water deficit in winter and spring in Australia can delay the budburst.Thus, the budburst functions may need to be adjusted if other factors were limiting budburst based on their local production conditions.A negative sigmoid function was selected for simulating the responses of chilling units (R c , Eq. ( 1)) and a growing degree approach was used to calculate forcing heat unit (R f , Eq. ( 2)).
where T a is the hourly air temperature and T b is the base temperature above which forcing heat unit will be accumulated.R c (T a ) and R f (T a ) are calculated at hourly intervals, and then summed and averaged daily to give the daily chill and heat unit accumulation.In APSIM Next Generation, hourly temperature is interpolated based on the method described in Goudriaan and Van Laar (1994).Values for c, d and T b were fitted for each variety.
Configurations can be easily adapted based on the needs of the users.More screenshots of the UI relating to the basic simulation set-up for a vineyard, list of model simulations and validations, sensibility to water stress and a simulation example of bringing sheep into the vineyard are in Supporting Information-Fig.S3.The flowering, fruit setting and berry development phases were simulated using a physiological development day (D s ) target which is the number of days the phase would take to complete at optimal temperatures.The daily increment of development days (dD) is calculated using the Wang-Engle temperature response function (Eqs (3) and (4)), which returns a value of 1 when temperatures are optimal and a value between 0 and 1 at sub-and supra-optimal temperatures (Wang and Engel 1998).Development days were used for all other temperature-dependent processes, e.g.leaf appearance and expansion.T min ≤ T a ≤ T max (3) where T a is air temperature, T min is the minimum temperature, T opt is the optimum temperature and T max is the maximum temperature.dD is calculated at hourly intervals, using interpolated hourly T a data, and these values are averaged to give dD.Wang and Engel (1998).Parameters were optimized for achieving the best prediction for flowering and véraison simultaneously.
After the summer solstice, the photoperiod decreases.When photoperiod decreases to a critical value, the bud for the following season's growth enters into endo-dormancy.However, the whole-plant phenology cannot be set to endo-dormancy at this time because the current year's shoots and berries are still growing.Although it would be tidy to develop two different phenological clocks, one for the bud and one for the plant in growth, we decided to work around it by resetting the plant phenology phase into endo-dormancy at the winter pruning event, which was defined in the management script with user input time.The chilling accumulation between the critical photoperiod and the pruning event was captured by a function and then added to the total chill accumulation in the endo-dormancy phase.

Structure module.
The structure module keeps track of the structure of the canopy, e.g.bud number per vine, main-stem leaf appearance, branching and canopy height (Fig. 1).Main-stem leaf appearance is modelled with the appearance of a new leaf every time dD accumulation exceeds the current phyllochron (the D s interval between two successive leaves).Leaf appearance occurs between budburst and véraison with a phyllochron of 1.64 for the first five leaves to appear, increasing to 1.82 for the 10th leaf onwards (Greer and Weston 2010).
Branching rate was calculated as the product of three components: potential branching rate, the effect of leaf area per vine and the effect of retained bud number per metre row.Parameters were optimized based on the retained node number experiment (see detailed description of the experiment in Data sets used for model calibration and validation).

Row-planted crop light interception method.
The row-planted crop light interception method developed by Goudriaan (1977) and Pronk et al. (2003), and later widely used for intercrop studies (Zhang et al. 2008;Gou et al. 2017) was implemented in APSIM Next Generation.The light interception method assumes that on a daily basis, the light intensity from all incoming directions is the same (Goudriaan 1977).It can be considered as a uniform overcast sky where only diffuse radiation is presented.This assumption enables easy spatial integration over all directions without the need to consider time of day.Furthermore, it represents the row canopy as a cube and assumes leaves are randomly distributed within the cube (Pronk et al. 2003).Based on those assumptions, the light interception by the row canopy can be expressed as a function of canopy width, canopy depth, the distance between two rows, leaf area of the row crops and light extinction coefficient [see Supporting Information-Fig.S2].This light interception method is used when the field is configured by two rectangular zones.Row orientation is not accounted in the light interception calculation.If the canopy of a training system cannot be assumed to be a cube (e.g.tatura training system), new light interception methods can be added in the APSIM Next Generation source code (https://github.com/APSIMInitiative/-ApsimX/),see details in 4.1 Canopy development and radiation interception.
Several modifications were made in the calculation compared with the method described in Gou et al. (2017) in order to capture the vineyard/orchard canopy features.For instance, for a vertical-shootpositioned (VSP) grapevine, all the leaves were accumulated in the top part of the canopy with a bare trunk at the bottom.The original model assumed that the base of all canopies was at ground level.We removed this constraint by extending equations to allow it to start from any height.
Canopy depth and canopy width were used to calculate the light interception by the grapevine row instead of using canopy height and row width [see Supporting Information-Fig.S2].Agricultural Production System sIMulator Next Generation distinguished canopy height from canopy depth.Canopy height was defined as the distance from the top of the canopy to the ground, while canopy depth was defined as the distance from the top of the canopy to the bottom of the canopy excluding the trunk part where there is no leaf.Canopy width was estimated by the width used in the summer trimming practice and can be viewed as the width where 95 % of leaves were covered.

Carbohydrate allocation method. All organ classes in APSIM
Next Generation use the same biomass subclass to keep track of their live and dead biomass status.The biomass subclass currently deals with dry matter (DM) and nitrogen (N) content of organs.Each of these components of biomass were separated into three types of listed, with arrows indicating the start and end time of the process and a drawing indicating the dynamics of the process.DAWS represents days after winter solstice and the values indicated are the mean DAWS for each phenological event for Sauvignon blanc in Marlborough, New Zealand.DM represents dry matter and NSC represents non-structural carbohydrate.Perennial organs include cane, trunk and structural root.pools represented by the properties (see Supporting Information-Method S1 3.12 Arbitrator and Fig. 2): • StructuralDM and N are essential for the growth of the organ.They remain within the organ once they have been allocated, and pass from Live to Dead pools as the organ senesces.• MetabolicDM and N are essential for crop growth and their concentrations can influence the function of organs, e.g.MetabolicN content is a main factor to determine leaf photosynthetic efficiency.MetabolicDM and MetabolicN may be reallocated (moved to another organ upon senescence of this organ) or translocated (moved to other organs at any time when supplies do not meet their structural and metabolic DM demands).• Non-structural DM and N are non-essential to the function of an organ and they may be reallocated or translocated.
The DM allocation arbitrator steps through a series of calculations to determine the supply of biomass available for allocation each day.Biomass supply includes reallocations from senescing organs, photosynthesis fixation and remobilization of non-structural reserves, and they are used in sequence (see detailed description in Brown et al. (2019)).
The existing allocation methods in APSIM Next Generation allocate DM to non-structural reserves only when all other carbohydrate demands are satisfied.This is not true for perennial crops, where trunk and root reserves are important competing sinks all the time.Therefore, a new allocation method, Q-Priority-then-Relative-Allocation, was developed based upon the existing Priority-then-Relative-Allocation method, where priority used only maintenance demand in APSIM Next Generation (Brown et al. 2019).The Q-Priority-then-Relative-Allocation method is based on the common carbohydrate pool assumption that the ability of developing plant organs to attract carbohydrates is independent of the topological position of the organs.However, the priority for different DM demand varies.The use of the priority parameter was adapted from the mathematical kiwifruit model (Buwalda 1991) and the functional-structural kiwifruit model (Cieslak et al. 2011).In the Q-Priority-then-Relative-Allocation method, DM was allocated in two steps.In the first step, an allocation was calculated based on the relative DM demand and priority parameter (q) of each biomass type (c) for each organ (o, Eq. ( 5)): where Prop[o] c is the proportion of DM allocated to organ o and biomass type c, and q is the corresponding priority parameter.
If there is still biomass unallocated and demand unsatisfied, a second step is carried out where the remaining biomass is allocated based on the relative demand without the priority to ensure every organ has a chance to fulfil their full demand.Thus, the priority factor limits growth only when there is a carbohydrate deficit.
Growth respiration coefficient was set as 0.2 gDM gDM −1 for all organs except berry, which had a much smaller growth respiration (0.02 gDM gDM −1 ) (Dai et al. 2009).Maintenance respiration for different components was set as in Zhu et al. (2019).The exponential response of maintenance respiration to temperature was set according to the Arrhenius equation (Poni et al. 2006).The temperature sensitivity coefficients Q 10 [-] were 1.9 for stem, 1.77 for berry and 1.66 for leaf (Poni et al. 2006).Q 10 was defined as the increase in the respiration rate resulting from a temperature increase of 10 °C.

Leaf module. The grapevine model uses the phytomer-based
Leaf organ class that simulates the expansion and senescence of cohorts of leaves at each position on the primary shoot and estimates the number of leaves in each cohort based on branching rates.Cohort leaves are treated the same as main-stem leaves appearing at the same time.
The maximum area of each cohort is currently set as a function of leaf rank.The growth duration, lag duration, senescence duration and minimum and maximum specific leaf area (SLA [m 2 g −1 ]) were parameterized based on our field experiment (see Section 2.2.4).For correctly modelling leaf senescence and the dynamics of leaf area, a photoperiod acceleration effect was added.The photoperiod acceleration effect also reflected the effects of the time of harvest on leaf senescence, as we observed that leaves would stay green much longer if the fruit were not harvested (see Supporting Information-Method S1 3.5 Leaf module).Other leaf properties were parameterized to ensure no nitrogen on leaf expansion in the current grapevine model, as the nitrogen dynamic parts require further calibration.
The potential total DM demand of each leaf was calculated as the product of the daily leaf area increase (constrained by stress) and mean specific leaf area.The total DM demand was then divided into structural and metabolic demands based on the fraction of structural DM in total leaf weight.The priority factor q for leaf structural and metabolic demand was set to a high value (1.6) as it is the immediate source of photosynthates and has the highest priority early in the season (Buwalda 1991;Lakso et al. 2008).Nonstructural reserves were not considered in the current leaf model.At leaf senescence, metabolic leaf DM was reallocated to other organs (DMReallocationFactor = 1).
The dynamics of seasonal canopy width and canopy depth were defined in the leaf class based on our field observation.A canopy pruning event was simulated each summer, which removed 20 % of leaf and shoot biomass.The pruning time and intensity can be configured easily in the management script for each simulation.For current simulations in the test set, the same pruning dates and intensities were used, as we do not have the exact pruning information for each site and year.

Cane, trunk and structural root module.
Cane, trunk and structural root models were all represented by the Generic Organ class, which has properties of biomass status and daily biomass demand and supply (Fig. 2; see Supporting Information-Method S1 3.7 to 3.9).The structural DM demand of Cane and Trunk (Eq.( 6)) were calculated by their radius (r [m]) derived from their DM and length, daily growth rate of the radius (dr/dt [m day −1 ]), length (l [m]) and wood density (ρ[g m −3 ]) (Cieslak et al. 2011): where i is organ type, q str,i is the priority factor for structural DM demand for a certain organ, determined as 0.7 for all three organs (see Materials and Methods section).Daily growth rate of the radius was calculated based on the trunk circumference increment in 14 years measured on 800 field-grown irrigated Sauvignon blanc vines.
The structural carbohydrate demand of the structural root was calculated by the structural demand of the trunk times the structural root/trunk ratio (set to one in this model).The structural root/trunk ratio may vary between vineyards and training systems and further work is required to capture this in the model.The non-structural carbohydrate (NSC) DM demand was modelled as an active competing sink (Cieslak et al. 2011), and the parameters were kept the same for those three organs.NSC included total soluble carbohydrates and starch.The rate of non-structural DM synthesis depended on organ size and limited by overloading (Eq.( 7)): The equation limited carbohydrate storage owing to non-structural DM overloading and the organ storage capacity, which was proportional to the current organ structural biomass, C max,NSC is the maximum non-structural DM per unit of structural DM, set to 0.26 based on field trunk and root samples of Sauvignon blanc (Greven et al. 2016).q NSC was set to 0.61 throughout the growing season.k syn is the daily nonstructural DM synthesis rate per unit of structural biomass, which is represented by a beta growth function in Eq. ( 8) (Yin et al. 2003) for capturing the fast recovery of carbohydrate reserves after flowering (Greven et al. 2016;Seleznyova et al. 2018).Parameters were optimized based on the measured dynamics of NSC in trunk and root by Greven et al. (2016) (see Materials and Methods section): where k max is the maximum non-structural DM synthesis rate, 0.03 [g g −1 ].For potential carbohydrate retranslocation, 3.7 % of the total carbohydrate storage in each organ per day (DM retranslocation factor) was set available to be retranslocated.Actual retranslocation depends on the daily carbohydrate supply and demand differences.Carbohydrate supplies from organ senescence and photosynthesis will be used first before using carbohydrate storage.
The biomass of retained cane for the following season was reset by the winter pruning events based on input cane diameter in the UI for the cane-pruned vines.The biomass of trunk and structural root keeps growing each year.In addition, retained cane can be parameterized as cordon and keep growing each year as well.

Annual shoot module. Shoot was represented by a Generic
Organ class as well.The structural biomass demand of the shoot was calculated as the daily delta of a beta growth function, which was fitted to individual shoot weight as a function of development days, and then multiplied by the shoot population per square metre and priority factor for shoot demand (Eq.( 9); see Supporting Information-Method S1 3.5).Shoot population per square metre was calculated internally based on shoot number per vine and vine density.
where StructuralDMs(t) [g m −2 ] is the structural DM of all shoots at development day t from budburst.Here shoot refers to the first year shoot stem and does not include leaves.SP is the shoot population per square metre.StructuralDM max [g] is the fitted maximum structural DM per shoot under four-cane pruned vines.t e is the time when StructuralDM max is reached (79.2 D s ), t m is the time when maximum growth rate occurs (43 D s ).Note the shoot DM includes both the primary and secondary shoots.The priority factor for shoot structural DM demand was set to a very low value (4e-3) to capture the reduction of shoot biomass under high shoot density per metre caused by both smaller primary shoots and fewer lateral shoots (Greven et al. 2014).Parameters for shoot non-structural DM demand, e.g.k syn and priority factors were set as half the values for Cane.

Root module.
The root module was parameterized to represent the fine roots of the plant, which are responsible for the uptake of water and nitrogen (see Supporting Information-Method S1 3.11).The root module calculates root growth in terms of rooting depth, biomass accumulation and subsequent root length density in each soil layer.A proportion of the fine roots can be allocated to different zones as well, e.g.inter-row.
Only structural DM demand was considered for the fine root, which was estimated by current root DM and relative root growth rate (gDM gDM −1 day −1 ).The daily loss of roots was calculated using a SenescenceRate function.All senesced material was automatically detached and added to the soil fresh organic matter.Fine roots often represents less than 20 % of the total root biomass (Buwalda 1993).Based on the data of Lakso et al. (2008), we set the initial fine root DM to be 60 g per vine (converted into g m −2 in the model based on population density).The annual fluctuation of the fine root DM was captured by dynamics of the relative growth rate and senescence rate.The dynamics of relative growth rate and senescence rate throughout the growing season were added based on the data presented in Comas et al. (2005).The priority factor for root DM demand was set to the same as for leaf DM demand.The perennial structural roots were separated from the root module because they have different functions and the separation makes the parameterization easier, e.g.root length growth.
Root length growth was calculated using the daily DM partitioned to fine roots and specific root length (SRL [m g −1 ]).Root proliferation in layers was calculated using an approach similar to the generalized equimarginal criterion used in economics.The uptake of water and N per unit root length was used to partition new root material into layers of higher 'return on investment' .

Berry module.
The berry module was represented by a reproductive organ class.Final yield per vine was calculated as the product of berry number per vine and berry fresh weight.Both berry fresh weight and dry weight were calculated dynamically (see Supporting Information-Method S1 3.10).
The berry number per vine was calculated by shoots per vine (as a function of retained node number), bunches per shoot and berries per bunch.Strong seasonal variations could be found in bunch number, berry number and berry fresh weight especially in cool-climate viticulture regions.The seasonal variation was captured by including the effects of weather conditions at critical periods around flowering of the previous and current season (see details in Zhu et al. (2020)).Parameters were readjusted in the current model as the calculation of development days was updated.
The methods for determining bunch number per vine in Zhu et al. (2020) were further generalized into the product of bunch number per shoot and shoot number per vine to account for different numbers of retained nodes per vine.Mean inflorescence number per node counted on 80 nodes per plot during the flowering assessment in the regional phenology and yield monitoring data set (see Section 2.1) were used for model calibration.Mean daily maximum temperature (TmaxIni) and daily radiation during the inflorescence initiation were used for predicting mean bunch number per shoot for the following season.TmaxIni and mean daily maximum (TmaxFlow) and minimum temperature around the flowering period of the current season were used to predict berry number per bunch.TmaxFlow, mean daily radiation intensity and cumulated rainfall around flowering, and cumulated rainfall after flowering and before véraison were chosen to predict the potential berry fresh weight.Furthermore, in the current model, carbohydrate effects represented by total carbohydrate supply and demand were added in the calculation of bunch number, berry number and potential berry fresh weight.The yield estimation method assumes that the potential yield components can be estimated based on the weather and carbohydrate conditions at certain critical periods, e.g.inflorescence initiation, while other periods have no effect.It has been shown that weather and carbohydrate conditions pre-budburst have strong effects on flower number (Eltom et al. 2017), and potentially on bunch number as well.However, our data set is insufficient to quantify such responses.
The dynamic of berry fresh weight was calculated based on the development days from flowering using the beta growth function (Eq.( 8)).Potential berry fresh weight (k max in Eq. ( 8)) was adjusted at véraison based on weather conditions and source/sink ratios in each season.
Total berry potential DM demand was determined based on berry number and the mean DM demand of individual berries described using the beta growth function (Eq.( 8)).Parameters were determined based on Section 2.1.The priority factor for berry DM demand was set to one.
Total soluble solids (TSS in °Brix) after véraison was calculated based on the ratio of berry dry weight to fresh weight through an empirical relationship between water content and TSS derived by Garcia de Cortazar-Atauri et al. (2009).Total titratable acid concentration (TA in g L −1 tartaric acid equivalents) was simulated through a negative exponential equation (Eq.( 10)): where TA max is the TA observed at véraison, 30.3 g L −1 in our data set, TA min is the fitted minimum TA, 7.8 g L −1 .Four 32-berry samples were collected weekly immediately prior to véraison and this continued until harvest, to determine the berry weight and TSS for each site and cultivar combination.TA was measured only for Marlborough Sauvignon blanc.A threshold soluble solids concentration of 8 °Brix was used as an alternate measure of the midpoint (50 %) of véraison (Parker et al. 2014).Yield, TSS, and TA prediction were calibrated only for the four long-term Sauvignon blanc sites.

Data set 2: retained node number trial.
The effects of retained node number per vine on leaf appearance, leaf area and biomass allocation were calibrated with the data from the retained node number trial (Greven et al. 2014;Greven et al. 2015).Three-year-old Sauvignon blanc grapevines (clone UCD1MS on Schwarzmann rootstock, Vitis riparia × V. rupestris) growing on a well-drained silt-loam soil (Wairau series) at the Rowley Crescent vineyard [41.49S; 173.95E] were selected for the trial in July 2006.Vines were planted with 2.8 ×1.8 m row-by-vine spacing.The vines, which had been previously pruned to four canes of 12 nodes (48 nodes), were cane-pruned to five different numbers of nodes per vine: 24, 36, 48, 60 or 72 by retaining either two, three, four, five or six 12-node canes per vine, respectively.The node treatment was repeated on the same vines each winter during pruning time for four consecutive years.
Early-season primary shoot vigour was assessed by taking detailed pre-flowering measurements of shoot length, shoot diameter (at third internode) and leaf number.Furthermore, shoot lengths were continuously monitored by a timelapse outdoor PlantCam (Wingscapes WSCA04) and manually measured from the photograph based on a reference ruler in the images.Leaf area per vine was measured by point quadrat after fruit set and by whole-vine destructive sampling at harvest.A random 30-berry sample was collected weekly from 50 % véraison and this sampling continued until harvest to determine the berry weight and TSS from a four-vine bay (plot).At harvest, all bays were hand harvested, with yield and total bunch number per plot (bay) recorded.Final berry weight was determined from a randomly collected 100-berry sample.Shoot number, shoot and cane fresh weight were assessed after harvest in winter.
Each season, during late winter, samples were taken from the trunk of one grapevine in each research plot for the assessment of overwinter carbohydrates reserves.A small 5-mm trunk corer was used to extract a piece of trunk wood that was quickly returned to the laboratory and freeze dried.Once freeze dried, samples were ground into powders and kept at −20 °C until analyses were performed.

Data set 3: post-harvest defoliation trial.
Starting in October 2008, a trial for studying the effects of post-harvest defoliation on carbohydrate reserves was conducted in the neighbouring rows of the retained node number trial.Two pruning treatments were compared: 48 nodes and 72 nodes (Greven et al. 2016).Eight carbohydrates measurements were taken throughout the season, starting at five-leaf stage, then at bloom, lag phase, véraison, mid-ripening, harvest, leaf fall and dormancy.These measurements were conducted for four years.Samples from both trunk and root were collected.Root samples were taken from a mixture of old and young roots varying from 1 to 5 mm in diameter at a depth of ~150 mm.The results of this trial were grouped with the data from the retained node number trial for simplifying the simulation set-up.

Data set 4: MRL_Rapaura_Irrigation.
The canopy development, dynamics of leaf and shoot biomass were validated by data from a regulated deficit irrigation trial at Central Rapaura Marlborough (noted as MRL_Rapaura_Irrigation) [41.47 S; 173.89 E] from 2001 to 2005 (Greven et al. 2005).The experimental results showed little differences in canopy development, dynamics of leaf and shoot biomass among different irrigation treatments; thus, we grouped the results and only used the mean treatment values.The comparison of the dynamics of soil water content in vine row and in alley between observation and simulation is shown in Supporting Information-Method S1 6.1.5and 6.1.6soil water.
The trial was conducted on four-cane pruned Sauvignon blanc vines (Telek 5C rootstock) grown on Wairau silt loam.In total, 72 vines were monitored for leaf, shoot and berry development.Canopy density was measured using the point quadrat with biweekly or monthly intervals.One or two destructive harvests of all leaves and shoots to establish the total leaf area and shoot biomass on the vines were done in mid-January or February.Fruit composition changes (TSS, TA and pH) The regional phenology and yield data set were used for calibration and validation of the phenology and yield module. b The retained node number and post-harvest defoliation data sets were used for calibration and validation of the yield estimation, leaf area, carbon allocation between different components, and the dynamics of organ biomass and NSC.c The MRL_Rapaura_Irrigation data set was used for calibrating the fraction of radiation interception, the water use and water uptake competition between vine row and alley.It was also used for validation of leaf area, carbon allocation between different components and the dynamics of organ biomass.

Model calibration methods
Model parameters were determined based on data sets one to three where possible, e.g.parameters linked with phenology, leaf appearance rate and the dynamics of biomass demand of shoot and berry (Table 2).Data set 4 was mainly used for validation.When parameters were not measured directly or linked with whole-plant function, e.g.q str,i (priority coefficient) they were optimized through the DEoptim function in R (Ardia et al. 2011).
A link between APSIM Next Generation and R programming language was created through the APSIM Next Generation edit feature.The APSIM Next Generation edit feature enables users to edit the parameters through bash commands.An R script was developed where we can call APSIM Next Generation, pass the parameters and assess the performances of the model under new parameter values.The script was then coupled in the optimization function (DEoptim) for improving the model performances through parameter optimization.The criteria for optimization were set to minimize the sum of squares between observed values and predicted values when one observed variable was of concern.When there was more than one variables of concern, the criteria were set to minimize the sum of squares of normalized differences between observed and simulated values (Eq.( 11)): where i is the sample number (or time point), n the total number of measurements, j is the variable number, m the total number of variables in concern.X sim,j,i is the simulated value for variable j and sample number (time point) i, and X obs,j,i is the observed value for variable j and sample number (time point) i.
Variables that were less influenced by other processes were calibrated first.In general, the model calibration and validation were conducted in sequence of phenology (budburst, flowering and véraison), canopy height and width, leaf appearance rate, leaf area via changing branching rate and leaf extension duration, the biomass of different organs via radiation use efficiency (RUE) and carbohydrate allocation priority (with yield components as input in management script), NSC, bunch number, berry number and berry weight.Several iterations were undertaken, as there are correlations between different parameters.
Goodness-of-fit was assessed by root mean square error (RMSE) of the concerned variables and coefficients (intercept, slope and R 2 ) of linear regressions between observed and simulated values: Values and units Source

DM-Retranslocation-Factor
The maximum percentage of NSC in each organ is available to be retranslocated each day 3.7 % (unitless) Data set 2 and 3 Annual shoot module SP Shoot population per square metre m −2 (Eq.( 9)) Structural-DM max The fitted maximum structural DM per shoot under four-cane pruned vines g (Eq.( 9)) Data set 2 and 3 where i is the sample number, n the total number of measurements.X sim,i is the simulated value and X obs,i is the observed value.The unit of RMSE is the same as that of the investigated variable.

Simulation set-up
An initiation method (treeInitiation) was developed in each simulation in a management script for initializing the vineyard and vine conditions.Users can specify the variety, the starting date of the simulation, vine age, vine distance, bud number per vine, cane (or cordon for spur-pruned vine) number per vine, cane diameter, cane length, trunk diameter, trunk length, ratio of structural root biomass to trunk biomass, fraction of NSC in trunk, structural root and cane.The width of vine row and alley was defined in the soil zone.
Bunch number per node and berry number per bunch were determined/influenced by vine carbohydrate status and weather conditions in the season prior to the season of harvest.They were input for the first year of a simulation.As such, the yield result from the first year should be treated with caution as it does not reflect the model's predictive mechanisms.We suggest starting the simulation one season earlier than the first observation season.
For the moment, the yield module has only been calibrated for irrigated and fertilized Sauvignon blanc grapes.Cane and trunk biomass were estimated based on the input dimension and default wood density in the management script.Non-structural carbohydrate content was calculated accordingly based on the fraction of NSC input for each organ.The age of the vine can be used in the calculation of some agerelated attributes, e.g.canopy height and width.After initiation, the simulation can run continuously for the whole experimental period, e.g. 15 seasons for the phenology monitoring trial.
Soil physical parameters, water, nitrogen and organic matter properties can be specified for each simulation.For the current data sets, the root module was conducting calculations of water and nitrogen uptake, and the nitrogen arbitrator was automatically conducting nitrogen allocation and nitrogen-limited growth for each organ in the background.However, water and nitrogen stress on canopy and fruit growth were minimized by irrigation and fertilizer input in the currently simulation.A detail sensitivity analysis of the performance of the model outputs to days of interval between two irrigation events without the influence of rain (rainfall was removed in the script) were shown in Supporting Information-Method S1 8.1 Irrigation Response.

Phenology
Budburst, flowering and véraison dates for Sauvignon blanc in eight sites across New Zealand were simulated well (Fig. 3), although the temperature profile varied greatly among sites owing to local terrain and distance to the ocean.For the calibration data set, the RMSEs were 4.1 days (based on date) for budburst, 3.1 days for flowering and 6.8 days for véraison for the calibration data set.For the validation data set, the RMSEs were 4.9 for budburst, 4.2 for flowering and 5.3 for véraison.
The model was parameterized for four other varieties: Chardonnay, Merlot, Pinot gris and Pinot noir [see Supporting Information-Fig.S3].The overall RMSEs for those four varieties were 6.5 for budburst, 4.85 for flowering and 6.5 for véraison.

Leaf number and leaf area dynamics
The rates of leaf appearance, leaf expansion and branching had considerable effects on the leaf area dynamics.The performance of the model in those three aspects was evaluated by the dynamics of leaf number per main shoot (Fig. 4A), leaf number per whole shoot [see Supporting Information-Fig.S4], individual leaf size [see Supporting Information-Fig.S5], leaf area of the main shoot and leaf area of the lateral shoot (Fig. 4C and D).The rate of leaf appearance was similar across different retained node number treatments except for the 72-node treatment, which had fewer leaves on the main shoot in the 2009/10 season (Fig. 4A).This was not captured by the model, as only a constant leaf appearance rate was used in the model.The leaf area per vine did not differ greatly between different treatments in either observed or simulated values.However, a big drop in leaf area of the lateral shoot under high retained node number per vine was revealed by the model (Fig. 4D).This was captured through branching rate in the Structure Module.The sudden drop in leaf area around 200 DAWS was caused by the summer pruning events.

Radiation interception
Radiation interception is an important aspect of the model performance, and together with RUE they drive the DM production and carbohydrate accumulation.The newly implemented row-planted crop light interception method captured the seasonal dynamics of light interception in different years in the MRL_Rapaura_Irrigation trial (Fig. 5).The corresponding dynamics of leaf area are shown in Supporting Information-Fig.S4.Measured and simulated fractions of light interception increased rapidly from budburst (~100 DAWS) and reached a maximum fraction of light interception of 53 % around 200 DAWS.The measured fraction of light interception by light sensor in the alley was slightly higher than the simulated values.The deviation between the simulated values and observed values could be caused by the light interception by trunks and canes, which would cause a reduction in the radiation arriving at inter-row but would not increase DM production.The optimized RUE based on the observed leaf, shoot and berry biomass before flowering was 0.8 gDM MJ −1 , and after flowering was 0.9 gDM MJ −1 .

Dynamics of the DM of different organs
Simulated leaf DM per shoot (Fig. 6A, R 2 = 0.87), shoot DM per shoot (Fig. 6B, R 2 = 0.84), and berry DM per vine (Fig. 6C, R 2 = 0.82) were well aligned with the observed values for both retained node number trial (model calibration) and the MRL_Rapaura_Irrigation trial (model validation).The DM per shoot decreased rapidly with increasing number of nodes per vine.The mean shoots DM reduced from 35 g per shoot under 24 nodes treatment to 12 g per shoot under 72 nodes treatment at leaf fall.
The seasonal patterns of DM dynamics of different organs were reproduced by the model (Fig. 7).Initially, about 90 % of the DM were stored in trunk and structural root.After budburst, the stored nonstructural DM was translocated into the newly growing leaf, shoot and berry (Fig. 7A).With photosynthesis and the support of stored DM, leaf and shoot DM increased rapidly.Leaf DM reached its peak at the end of January.Total berry DM increased rapidly after flowering and became the largest biomass pool at harvest (Fig. 7B).
The DM retranslocation before flowering decreased NSC concentration in the trunk and structural root rapidly (Fig. 8).The NSC pool started to replenish after flowering, mainly in competition with berry growth, as leaf and shoot DM demands had greatly decreased (Figs 7  and 8).The rate of NSC replenishment slowed down during véraison but increased quickly when the berry DM demand went down, and leaf DM retranslocation accelerated during leaf senescence.

Yield components
Wide variations in bunch number per vine (35-110), berry number per bunch (40-90) and yield per vine (3.5-13 kg) caused by seasonal fluctuations and retained node number per vine were captured by the model (Fig. 9).For the calibration data set, the RMSEs were 5.6 for bunch number per vine, 7.7 for berry number per bunch and 1.3 kg for yield per vine.For the validation data set, the RMSEs were 6.3 for bunch number per vine, 9.8 for berry number per bunch and 2.0 kg per vine for yield per vine.The model underestimated these attributes at high values, e.g.10-13 kg yield per vine, compared with the observations.
Two seasons (2016/17, and 2017/18) were selected to have a close check of the simulated dynamics of berry fresh weight, TSS and TA [see Supporting Information-Figs S6-S8].The season of 2016/17 had slightly higher cumulated growing degree days (°Cd) (>10 °C) from October to April than the long-term average (1299 °Cd versus 1228 °Cd).However, this season had more rainfall between véraison    perennial fruit crop to be implemented into APSIM and it provides a useful template for the development of models for other perennial fruit crops.

Canopy development and radiation interception
Apart from two structural grapevine models that simulated leaf area at the individual leaf level (Lecoeur et al. 2011;Schmidt et al. 2019), all other process-based grapevine models simulated leaf area at the canopy scale through different methods.For example, Bindi et al. (1997) calculated total leaf area based on the number of leaves that had appeared.In the STICS grapevine model, leaf area growth was driven by temperature, phenological stages and stresses, with an empirical plant densitydependent function representing inter-plant competitions.A few other grapevine models have calculated the increase in the leaf surface based on the biomass allocated to leaves and SLA (Wermelinger et al. 1991;Nendel and Kersebaum 2004).
The current APSIM grapevine model took the 'leaf-to-leaf ' approach, using phyllochrons to control the leaf initiation, leaf appearance and duration of leaf expansion to simulate the dynamics of potential individual leaf area.Actual leaf area depended on biomass allocation, SLA and other abiotic stresses, e.g.water stress.The variation in SLA caused by the ratio of structural and non-structural weight was captured by maximum and minimum SLA.Furthermore, for correctly simulating leaf area on lateral shoots, we controlled the branching rate through an empirical relationship between bud number per metre row and percentage of reduction in branching rate (Fig. 4).Further development could replace this empirical relationship with sounder physiological controls, e.g.source/sink ratio, and add the branching pattern caused by succession of phytomer types (P0-P1-P2) (Louarn et al. 2007).
The dynamics of leaf area was tightly linked with radiation interception, and the Beer-Lambert law was often used for the interception calculation (Bindi et al. 1997).However, grapevines are usually planted in rows with substantial alley space in between.This violates the homogenous assumption of the Beer-Lambert law.Several approaches have been developed based on canopy height and width, and canopy porosity, to account for such heterogeneity (Brisson et al. 2003;Oyarzun et al. 2007), but are not yet available in APSIM.Here, we added the row-planted crop radiation interception methods developed by Goudriaan (1977) and Pronk et al. (2003) into APSIM.We also adapted the method to account for three conditions: (i) canopy width is the same as row width, and leaves are distributed evenly across the plant height vertically (Zhang et al. 2008;Gou et al. 2017).(ii) Canopy width is smaller than the row width and leaves are concentrated in the top part of canopy, e.g.hedge-pruned vines.(iii) Canopy width is larger than the row width and leaves are concentrated in the top part of the canopy, e.g.apple orchards.The current method is also applicable to conditions where gaps within the row are not negligible, by considering the canopy as an array of cuboids as shown in Pronk et al. (2003) on the conifer eastern white cedar (Thuja occidentalis 'Brabant').For conditions where the canopy shape largely violates the cube assumption, other light interception methods for heterogeneous canopy need be implemented in the APSIM source code.de Castro and Fetcher (1998) developed a method to deal with three-dimensional (3D) canopy analytically.They divided the canopy into many cubic cells, and then calculated the probability that a beam would penetrate to any given cell without being intercepted by the foliage in the path, using an exponential extinction function.Another promising approach was to transform the 3D canopy light interception problem into a two-dimensional (2D) problem by using the shadow area and intensity on the ground at any given time (Rojo et al. 2020).

Carbohydrate allocation and carbohydrate reserves
Three approaches have mainly been used for simulating carbohydrate allocation in existing grapevine models: (i) Harvest index or fruit biomass index.For instance, Bindi et al. (1997) assumed the ratio of fruit dry weight and total aboveground biomass increased linearly during the reproductive phase.(ii) Partitioning coefficient.For example, The DM was expressed in g m −2 , which was the default unit in APSIM.One vine occupies 1.62 m 2 in this simulation set-up (1.8 m vine spacing × 0.9 m width of the vine row).The lines of trunk were overwritten by the lines for structural root, as we use a biomass ratio of 1 between structural root and trunk, and the same DM demand parameters for them.
Nendel and Kersebaum ( 2004) used a partitioning coefficient for biomass allocation and the partitioning coefficient changes according to phenological stages.(iii) Relative sink strength (Gutierrez et al. 1985;Wermelinger et al. 1991;Lecoeur et al. 2011).Carbohydrate reserves in perennial structures were rarely considered and only a few models considered the perennial reserves in different organs as a whole (Brisson et al. 2009;Pallas et al. 2011).
We expanded the existing carbohydrate allocation approaches as follows by: (i) simulating both NSC synthesis and hydrolysis explicitly in separated organs, e.g.cane, trunk and structural root, based on the approach presented in L-kiwi (Cieslak et al. 2011).(ii) Carbohydrate reserve (starch synthesis) was simulated as an active competing sink instead of allocating carbohydrate to reserve only when there was surplus carbohydrate.To achieve this, a new allocation method (Q-Priority-then-Relative-Allocation) was added to the PMF arbitrator to account for sink priority and sink strength simultaneously (Buwalda 1991;Cieslak et al. 2011).The new approach effectively captured the seasonal dynamics of NSC concentrations in both trunk and structural root compared with our field observations (Fig. 8).The pattern of the dynamics of NSC concentration depends on the rate of NSC synthesis and hydrolysis, but is strongly influenced by the rate of NSC consumption in spring and the rate of carbohydrate assimilation by leaf as well.The transition point of NSC from decrease to increase indicates when the vine achieves a net carbohydrate gain.

Variations in yield
Until now variations in yield and yield components, namely bunch number per vine, berry number per bunch and berry weight, were not well represented in process-based grapevine models (Moriondo et al. 2015).As discussed above, the DM accumulation of the fruit compartment was often simplified and the variation in berry number per vine was not considered.The STICS grapevine model did calculate the berry number per square metre during the fruit setting period (a thermal time duration).During this period, on each day, the number of set fruit was calculated using the potential number of set fruit per inflorescence and per degree day, the daily development rate (degree day), the number of inflorescences per plant (input or calculated based on plant carbohydrate status at the early stage), the plant density, the trophic stress index and the frost stress index acting on fruit from flowering (Brisson et al. 2009).As thermal time would linearize the temperature effect, the final total berry number will be similar each year if the carbohydrate status is similar.Thus, the STICS grapevine model did not account for the effects of weather conditions around flowering on bunch number per vine, berry number per bunch and potential berry fresh weight.This may not be a problem when applying the model in hot climate conditions where temperature does not limit the inflorescence initiation (Sadras and Soar 2009).However, for cool-climate viticulture, remarkable seasonal variations in yield, up to 100 %, can be found (Zhu et al. 2020).Much of this can be attributed to temperatures at inflorescence initiation and flowering, although inadequate carbohydrate reserves can result in inflorescence abortion and reductions in bunch weight (Eltom et al. 2013;Greven et al. 2014).We simulated the seasonal variations in bunch number per vine, berry number per bunch and berry weight by integrating the statistical yield model developed by Zhu et al. (2020) with the process-based model which added additional carbohydrate effects to the statistical yield model.The approach accounted for 74 % of the variations in bunch number per vine, 46 % in berry number per bunch and 54 % in yield per vine across seven sites with different weather conditions and retained nodes per vine, even with the same parameter set.
When validating the model, we found that some high-yield sites and years were underestimated.This was probably caused by underestimating the shoot number per vine.High-yield vineyards often have high vine capacity with large trunks, and a high proportion of vigorous shoots growing from adventitious buds (uncounted at winter pruning) or growing from secondary and tertiary buds of the counted buds.Including more detailed site information like the soil physical properties, soil nutrients and water information into the model for accounting for variations in vine capacity, shoot number per vine, and leaf area development, are likely to further improve the prediction.For correctly simulating those aspects, further work needs to be done in calibrating the water and nutrient uptake and competition with the inter-row cover crops in the current grapevine model.Furthermore, updating the model with real-time information, e.g.shoot and inflorescence number in spring, flowering and véraison date, would largely improve the model accuracy for the predictions of leaf area, yield, berry TSS and berry TA and could assist in-season management decisions.The yield module was calibrated with the four Marlborough long-term sites, namely Upper Brancott, Western Wairau, Seaview Awatere and Central Rapaura, and validated with the remaining sites (Table 1).Hawke's Bay Bridge pa and Waipara Valley Glasneven were not included in the yield comparison because retained number of canes per vine changes from year to year and even within the same year owing to management decisions on yield control.

CON CLUSIONS A ND FU T UR E PROSPECTS
and their effects on grapevine development.The modelling framework developed in this study can also be able to be applied to other perennial fruit crops as shown by our preliminary modelling work on apple (Malus domestica), peach (Prunus persica) and kiwifruit (Actinidia deliciosa).Once the grapevine model is validated across a broader range of training systems, climatic conditions and varieties, it will become a promising tool for illustrating the impacts of climate change on current viticulture practices and for bridging the gap between real-time field data and decisions that growers often face.

SUPPORTIN G INFOR M ATION
The following additional information is available in the online version of this article-Method S1.Agricultural Production System sIMulator (APSIM) Next Generation grapevine model documentation.

Figure 1 .
Figure 1.Illustration of the APSIM Next Generation UI.This illustrates the basic structure of the grapevine model (left panel), the configuration of the leaf and berry module (middle panel) and the configuration of structure and trunk module (right panel).

Figure 2 .
Figure 2. Illustration of the main processes happening in different modules according to the phenology timeline.The model follows the traditional calculation loop: phenology, potential organ growth, light interception, carbohydrate assimilation and allocation, and actual organ growth in a row-planted vineyard.The basic principle of the model is that phenology is largely driven by temperature and growth is dynamically driven by both environmental factors and differential allocation of carbohydrates among organs.Users first need to initialize the vineyard and vine condition in the tree initiation method (organ box) before running any simulation.Green boxes on the left list the main modules in the model.The main processes in each module are −Tmin) α (Topt−Tmin) α −(Ta−Tmin) 2α (Topt−Tmin) 2α

Figure 3 .
Figure 3. Verification and validation of the simulated budburst (A), flowering (B) and véraison (C) date expressed as DAWS for eight Sauvignon blanc sites across New Zealand.The solid black line is the one-to-one line between predictions and observations, and the red line is the linear regression between the two.HB represents Hawke's Bay, MRL represents Marlborough, WV represents Waipara Valley.The phenology module was calibrated with the four Marlborough long-term sites, namely Upper Brancott, Western Wairau, Seaview Awatere and Central Rapaura, and validated with the remaining sites.

Figure 5 .
Figure 5. Model validation of the fraction of light interception by grapevine canopy in the Central Rapaura Vineyard in 2003-2004 (A) and in 2004-2005 (B).Points were observed values and lines were simulated ones.The decrease in the fraction of light interception around 200 DAWS in both observed data sets was due to summer pruning events.

Figure 4 .
Figure 4. Model verification of vineyard leaf number per main shoot (A), leaf area per vine (B), leaf area on main shoots (C) and lateral shoots (D) in the retained node number trial in the season 2009/10.The sudden drop in leaf area around 200 DAWS was caused by the summer pruning events.

Figure 6 .
Figure 6.Verification and validation of leaf DM per grapevine shoot (A), shoot DM per shoot (B), berry DM per vine (C).Shoot DM is the stem part and does not include leaf DM.The dynamics of DM were calibrated based on the data from the retained node number trial (MRL_Rowley), and validated with the data from the MRL_Rapaura_Irrigation trial.The solid black line is the oneto-one line between predictions and observations, and the red line is the linear regression between the two.

Figure 7 .
Figure 7.The dynamics of the dry weight of each grapevine organ (A) and proportion of each organ in the total plant DM in the 48-retained node number treatment (four-cane) (B).The DM was expressed in g m −2 , which was the default unit in APSIM.One vine occupies 1.62 m 2 in this simulation set-up (1.8 m vine spacing × 0.9 m width of the vine row).The lines of trunk were overwritten by the lines for structural root, as we use a biomass ratio of 1 between structural root and trunk, and the same DM demand parameters for them.

Figure 8 .
Figure 8.The dynamics of the NSC concentration in grapevine trunk (A) and root (B) with different retained node number per vine.Lines are predictions and points are observed data.The trial started in 2006 when the vines were 3 years old.Detail NSC measurement started in 2009 when the vines were 6 years old.Root samples were taken from a mixture of old and younger roots varying from 1 to 5 mm in diameter at a depth of ~150 mm.

A
grapevine model has been developed and is the first perennial fruit crop model in APSIM Next Generation.It has all the structures and functions to simulate the phenology, radiation interception, leaf growth, DM production, carbohydrate and nitrogen allocation and storage, and yield formation of a perennial plant.The model can be easily adapted in the UI to reflect different process formulations (e.g.yield formation), and to simulate different cultivars or VSP training systems without changing the source code.This extends the usability of the model under different conditions.The model performed well accounting for the effects of weather conditions on phenology across eight sites in New Zealand and the yield components in the Marlborough region.It can be used to explore the long-term suitability of different crop loads, based on the stability of the yield between years and the carbohydrate reserves in winter.Further studies can be done on extending the abilities of the model to simulate water and nutrient competition between the grapevine and the understory plants

Figure 9 .
Figure 9. Model verification and validation of the bunch number per vine (A), berry number per bunch (B) and yield per vine (C) for seven Sauvignon blanc sites and five different retained cane number treatments in Marlborough, New Zealand.The solid black line is the one-to-one line between predictions and observations, and the red line is the linear regression between the two.The yield module was calibrated with the four Marlborough long-term sites, namely Upper Brancott, Western Wairau, Seaview Awatere and Central Rapaura, and validated with the remaining sites (Table1).Hawke's Bay Bridge pa and Waipara Valley Glasneven were not included in the yield comparison because retained number of canes per vine changes from year to year and even within the same year owing to management decisions on yield control.

Figure S1 .
Figure S1.Illustration of the set-up of a vertical-shoot-positioned vineyard where our data were collected and definition of different distances.Figure S2.Screenshots of Agricultural Production System sIMulator (APSIM) Next Generation user interface showing the basic simulation setup for a vineyard, list of model simulations and validations and a simulation example of bringing sheep into the vineyard, and sensitivity to water stress.Figure S3.Model verification of the simulated budburst, flowering and véraison dates for Chardonnay, Merlot, Pinot gris and Pinot noir across New Zealand.Figure S4.Model validation of the leaf per shoot, leaf area per vine, leaf area on main shoots and lateral shoots in the MRL_Rapaura_Irrigation trial.Figure S5.Illustration of the dynamics of leaf size in different cohort groups in the Central Rapaura phenology and yield monitoring vineyard in the 2004/05 season.

Figure S2 .
Figure S1.Illustration of the set-up of a vertical-shoot-positioned vineyard where our data were collected and definition of different distances.Figure S2.Screenshots of Agricultural Production System sIMulator (APSIM) Next Generation user interface showing the basic simulation setup for a vineyard, list of model simulations and validations and a simulation example of bringing sheep into the vineyard, and sensitivity to water stress.Figure S3.Model verification of the simulated budburst, flowering and véraison dates for Chardonnay, Merlot, Pinot gris and Pinot noir across New Zealand.Figure S4.Model validation of the leaf per shoot, leaf area per vine, leaf area on main shoots and lateral shoots in the MRL_Rapaura_Irrigation trial.Figure S5.Illustration of the dynamics of leaf size in different cohort groups in the Central Rapaura phenology and yield monitoring vineyard in the 2004/05 season.

Figure S3 .
Figure S1.Illustration of the set-up of a vertical-shoot-positioned vineyard where our data were collected and definition of different distances.Figure S2.Screenshots of Agricultural Production System sIMulator (APSIM) Next Generation user interface showing the basic simulation setup for a vineyard, list of model simulations and validations and a simulation example of bringing sheep into the vineyard, and sensitivity to water stress.Figure S3.Model verification of the simulated budburst, flowering and véraison dates for Chardonnay, Merlot, Pinot gris and Pinot noir across New Zealand.Figure S4.Model validation of the leaf per shoot, leaf area per vine, leaf area on main shoots and lateral shoots in the MRL_Rapaura_Irrigation trial.Figure S5.Illustration of the dynamics of leaf size in different cohort groups in the Central Rapaura phenology and yield monitoring vineyard in the 2004/05 season.

Figure S4 .
Figure S1.Illustration of the set-up of a vertical-shoot-positioned vineyard where our data were collected and definition of different distances.Figure S2.Screenshots of Agricultural Production System sIMulator (APSIM) Next Generation user interface showing the basic simulation setup for a vineyard, list of model simulations and validations and a simulation example of bringing sheep into the vineyard, and sensitivity to water stress.Figure S3.Model verification of the simulated budburst, flowering and véraison dates for Chardonnay, Merlot, Pinot gris and Pinot noir across New Zealand.Figure S4.Model validation of the leaf per shoot, leaf area per vine, leaf area on main shoots and lateral shoots in the MRL_Rapaura_Irrigation trial.Figure S5.Illustration of the dynamics of leaf size in different cohort groups in the Central Rapaura phenology and yield monitoring vineyard in the 2004/05 season.

Figure 10 .
Figure10.Model verification and validation of the single-berry fresh weight (A), total soluble solids (TSS, B) and total titratable acid concentration (TA in g L −1 tartaric acid equivalents, C).Data are for five Sauvignon blanc sites and five different retained cane number treatments in Marlborough, New Zealand.The solid black line is the one-by-one line between predictions and observations, and the red line is the linear regression between the two.The berry properties simulations were calibrated with the four Marlborough long-term sites, namely Upper Brancott, Western Wairau, Seaview Awatere and Central Rapaura, and validated with the remaining sites.

Figure S6 .
Figure S6.Model verification of the dynamics of berry fresh weight for the four long-term monitoring Marlborough (New Zealand) Sauvignon blanc sites in the 2016/17 and 2017/18 seasons.Figure S7.Model verification of the dynamics of berry total soluble solids (°Brix) for the four long-term Marlborough (New Zealand) Sauvignon blanc sites in the 2016/17 and 2017/18 seasons.Figure S8.Model verification of the dynamics of berry titratable acid concentration for the four long-term Marlborough (New Zealand) Sauvignon blanc sites in the 2016/17 and 2017/18 seasons.Figure S9.The effects of sites and year on the simulated dynamics of berry fresh weight and final berry weight.
The Wang-Engel function in APSIM Next Generation included one additional variable T ref (reference temperature) and returned the ratio of dD(T a ) and dD(T ref ) as the final results.When T ref equals to T opt , it has the same result as the original function in

Data sets used for model calibration and validation 2.2.1 Data set 1: regional phenology and yield monitoring.
t is the cumulated development days after véraison and r is the TA degradation rate with t, 0.0923 This data set was used for calibrating parameters related to phenology, bunch number, berry weight, berry TSS and berry TA.Starting in 2004, phenology data (budburst, flowering and fruit development), yield per vine, bunch number, berry weight and berry TSS and TA were collected from Sauvignon blanc vines in four commercial vineyards in Marlborough (Table 1, detailed descriptions of the experiment and yield calibration procedure in Zhu et al. (2020)).The trial was expanded to four other regions and varieties in 2014: Gisborne (two Chardonnay vineyards), Hawke's Bay (two Chardonnay, two Merlot, one Sauvignon blanc, one Pinot gris), Marlborough (two Chardonnay, two Pinot gris, two Pinot noir, eight Sauvignon blanc), Waipara Valley North Canterbury (one Pinot noir, one Sauvignon blanc) and Central Otago (four Pinot noir).In 2014 the regions represented by the data set accounted for 93 % of the New Zealand national vineyard area, and the five varieties in the data set represented also accounted for 93 % of the national varietal mix.The phenology of Sauvignon blanc was calibrated with the four Marlborough long-term sites, namely Upper Brancott, Western Wairau, Seaview Awatere and Central Rapaura, and validated with the remaining sites.For other varieties, all data were used for calibration.

Table 1 .
List of data set used in model calibration and validation.

Table 2 .
A partial list of model parameters and variables.A complete description of model parameters and variables was provided in Supporting Information-Method S1 APSIM Next Generation grapevine model documentation.A reference was provided for some parameter or variable names, while the remainder was calibrated based on the corresponding data set listed in Table1.