A land surface model combined with a crop growth model for paddy rice (MATCRO-Rice Ver. 1) – Part I: Model description

Crop growth and agricultural management can affect climate at various spatial and temporal scales through the exchange of heat, water, and gases between land and atmosphere. Therefore, accurate simulation of fluxes for heat, water, and gases from agricultural land is important for climate simulations. A land surface model (LSM) combined with a crop growth model (CGM), called LSM-CGM combined model, is a useful tool for simulating these fluxes from agricultural land. Therefore, we developed a new LSM-CGM combined model for paddy rice fields, the MATCRO-Rice model. The main 5 objective of this paper is to present the full description of MATCRO-Rice. The most important feature of MATCRO-Rice is that it can consistently simulate latent and sensible heat fluxes, net carbon flux, and crop yield by exchanging variables between the LSM and CGM. This feature enables us to apply the model to a wide range of integrated issues.


Introduction
In the last 15 years, climate and land surface modelling studies have shown that crop growth and farm management in agri-10 cultural land significantly affect climate via the exchange of heat, water, and gases. For example, applying a regional climate model combined with a crop growth model (CGM) to the United States, Tsvetsinskaya et al. (2001) showed that crop growth can change the surface temperature by 2 to 4 • C. Maruyama and Kuwagata (2010) showed that crop growing season can affect the amount of evapotranspiration by using a land surface model (LSM) combined with a CGM. Levis et al. (2012) incorporated a CGM into an earth system model, and showed that the timing of crop sowing can change the amount of precipitation. Using 15 a dynamic global vegetation model combined with a CGM, Bondeau et al. (2007) showed that the global carbon cycle, which has a significant effect on global warming, is largely modified by crop growth and farm management. Osborne et al. (2009), using a global climate model coupled with a CGM, demonstrated that the crop-climate interaction can affect annual variability in surface temperature. All these studies indicate that crop growth and farm management are key determinants of climate and that climate simulations need to accurately simulate the fluxes of heat, water, and gases in agricultural land. 20 A LSM or dynamic vegetation model (DVM) incorporated with a CGM, called LSM-CGM or DVM-CGM combined models, are a useful tool for simulating the fluxes of heat, water, and gases in agricultural land. Hence, several LSMs and DVMs incorporated with a CGM have been developed (BATS-GF: Tsvetsinskaya et al., 2001;Agro-IBIS: Kucharik, 2003; 3 Land surface model The main outputs of the LSM of MATCRO-Rice are LHF and SHF. The LSM has five modules, which are "energy balance at the canopy and surface water", "within-canopy shortwave radiation", "bulk transfer coefficient for latent and sensible heat", "canopy water balance", and "soil water and heat transfer". Each module is described in detail in the following sections. Before describing each module, we note the following two major modifications from the original LSM, MATSIRO (Takata et al.,5 2003).
1. LAI, crop height, and root depth, which are constant in the original MATSIRO, are dynamically calculated in the CGM and are the inputs to the LSM.
2. Surface water is added above the soil surface to represent a flooded surface in paddy rice fields.
Other minor modifications are described separately in each of the following sections. We note that the photosynthesis model 10 used in MATCRO is described in the CGM section (Section 4).

Energy balance at the canopy and surface water
This module calculates LHF and SHF by solving energy balance at two layers above the soil, canopy and surface water. The module is based on the original MATSIRO (Takata et al., 2003), except for the addition of surface water above the soil and other minor modifications. The energy balance at the canopy and surface water are given as follows: (1) where R nc and R nw are the net radiant flux density at canopy and surface water, H c and H w are the SHF from the canopy and surface water, E c , E t , and E w are the evaporation from wet canopy, transpiration from the canopy, and evaporation from the surface water, respectively, G ws is the heat flux from the surface water to soil, and S tw is the heat flux stored into surface water. 20 It is important to note that the downward flux for R nc , R nw , and G ws indicates a positive flux, whereas downward flux for H c , H w , E c , E t , and E w indicates a negative flux. All variables in the model are listed in Table 3. λ is the physical constant for the latent heat of vaporisation (Table 4). Each radiant, heat, and water flux in Eqs. 1 and 2 are given by the following equations. 10 where R d s (0), R d l (0), and R u s (0) are the downward shortwave radiant flux density, downward longwave radiant flux density, and upward shortwave radiant flux density at the canopy top, respectively, τ cs and τ cl are the canopy transmissivity for shortwave and longwave radiation, respectively, C Hc and C Hw are the bulk transfer coefficients (BTCs) for sensible heat between canopy and atmosphere and between surface water and atmosphere, respectively, C Ec and C Ew are the BTCs for latent heat between canopy and atmosphere and between canopy and atmosphere, respectively, T a , P a , U , and Q are air temperature, air pressure, 15 wind speed, and specific humidity, respectively, f cw is the fraction of wet canopy, T c , T w , and T s (0) are the canopy, surface water, and soil surface temperature, respectively, c pa and c pw are the specific air and water heat, respectively, k w is the water thermal conductivity, ρ w and ρ a are water and air density, respectively, σ is the Boltzmann constant, Q sat is specific humidity at saturation, d w is the depth of surface water, is the longwave emissivity of surface water, and d/dt indicates the time differentiation. The argument of the radiant flux density denotes LAI depth from the canopy top, and the argument of soil 20 temperature denotes soil depth from the soil surface. Therefore, R d s (0), R d l (0), and R u s (0) indicate the radiant flux density at the canopy top, and T s (0) indicates the soil surface temperature.
The variables ρ a and Q sat are physically calculated from the air temperature and air pressure (Appendix A), c pa , c pw , k w , ρ w , 25 and σ are physical constants (Table 4), d w is a simulation setting parameter (Table 5), and is set to 0.96 (Campbell and Norman, 1998). T c and T w are numerically determined to satisfy Eqs. 1 to 11. The numerical method is described in Masutomi et al. (2016).
The original MATSIRO uses C Hc instead of C Ec in Eq. 8 when specific humidity of the air is greater than the saturated specific humidity of the canopy (i.e., Q sat − Q < 0), because dew condensation occurs at canopy of interest. MATCRO does 30 not consider the effect for simplicity. It should be noted that C Hc is used for calculating the evaporation from wet canopy in Eq. 7.

Within-canopy shortwave radiation
The main role of this module is to simulate direct downward photosynthesis active radiation (PAR), scattered downward PAR, and scattered upward PAR at a LAI depth of l from the canopy top by calculating the transmission and reflection of shortwave radiation by leaves within canopies. These PARs are used for calculating carbon assimilation in the CGM (Section 4.1). In addition to the simulation of PARs, transmissivities for shortwave and longwave radiation are simulated in this module. The 5 transmissivities are used for calculating LHF and SHF (Section 3.1).
This module is based on the simple model developed by Watanabe and Ohtani (1995). The model determines radiation within canopies by calculating the transmission and reflection of the radiation within the canopy. In this model, radiation within the canopy is divided into three components (downward direct, downward scattered, and upward scattered) and two wavebands (PAR and near infrared [NIR]). In addition, the following three assumptions are considered in the model for simplicity. 10 1. Leaf orientation is random (i.e., spherical distribution).
2. Leaf reflectivity and transmissivity of the radiation are vertically uniform within a canopy.
3. Scattered radiation income from a zenith angle of 53 • .
It should be noted that the assumption 3 is based on the fact that radiant flux uniformly emitted from a horizontal plane is approximately equal to radiant flux density from a zenith angle of 53 • . From the three assumptions above, we can express 15 analytically the radiant flux density for downward direct (D d i (l)), downward scattered (S d i (l)), and upward scattered (S u i (l)) within canopy for each waveband (i = 1: PAR; i = 2: NIR), as follows: 20 Here, F is a parameter for the distribution of leaf orientation. If we assume spherical distribution for leaf orientation as mentioned above, we have F = 0.5 (Goudriaan and van Laar (1994)). The variable l is a LAI depth from the canopy top. The variable θ is a zenith angle of the sun (Appendix B). The function sec() indicates the secant function. The coefficients, a i , where R d s (0) is the downward shortwave radiant flux density at the canopy top and f df is the fraction of scattered radiation to total radiation. In Eqs. 15 and 16, we assumed that both PAR and NIR are half of R d s (0). According to Goudriaan and van Laar 30 (1994), f df is given as a function of the transmissivity of atmosphere (τ atm ) as follows: where R ex is the extraterrestrial radiation, R sun is the solar constant, and D oy is the number of days from Jan 1. The equations 5 15-19 that calculate D d i (0) are based on formulations by Goudriaan and van Laar (1994), while the original MATSIRO uses different equations.
The transmissivity of canopies for shortwave radiation (τ cs ) is expressed as Here, R u s (0) and R d s (L) are the radiant flux density for upward shortwave at the canopy top and downward shortwave at the 10 bottom of the canopy, respectively. L denotes the LAI, which is calculated in the CGM (Section 4.4). R u s (0) and R d s (L) are represented by where r ij and τ ij are the canopy reflectivity and transmissivity, respectively, i and j represent wavebands (i = 1: PAR; i = 2: 15 NIR) and direct (j = 1) or scattered radiation (j = 2). These are given in Appendix D.
Last, the transmissivity of a canopy for longwave radiation (τ cl ) is expressed as where, d f is the scattered factor. We set d f = sec(2π(53/360)) from the assumption that scattered radiation income is from a zenith angle of 53 • (Watanabe, 1994).

Bulk transfer coefficient for latent and sensible heat
This module calculates BTCs for latent and sensible heat (C Ec , C Ew , C Hc , and C Hw ). The BTCs are used to simulate energy balance (Section 3.1). This module is based on Watanabe (1994), where C Ew , C Ec , C Hw , and C Hc are given by where C E and C H are the BTCs for latent and sensible heat between the entire surface (canopy + surface water) and atmosphere and are given by In Eqs. 24 to 29, κ is the Karman constant, d is the zero-plane displacement height, z a is the reference height at which wind 5 velocity is observed, z M w , z T w , z Qw are the roughness lengths that express the effect of surface water on the profiles of momentum, temperature, and specific humidity, respectively, z M , z T , and z Q are the roughness lengths of an entire surface (canopy + surface water) for the profiles of momentum, temperature, and specific humidity, respectively. z a is a simulation setting parameter (Table 5), and d, z M , z T , z Q , z M w , z T w , and z Qw are the functions of crop height and LAI (Appendix E).
Ψ M , Ψ H , and Ψ E are the diabatic correction factors for momentum, heat, and vapour transport, respectively. The factors are 10 functions of atmospheric stability ζ as follows: The equations above are adopted from Campbell and Norman (1998), whereas the original MATSRIO model employs different 15 equations. The variable ζ is replaced by either the atmospheric stability between the entire surface and atmosphere (ζ) or the atmospheric stability between surface water and atmosphere (ζ w ). These are given by where L M O and L M Ow are the Monin-Obukhov lengths for the exchange between the entire surface and atmosphere and 20 between the surface water and atmosphere, respectively, and are given by where g is the gravitational constant, T w and T c are the temperatures of the surface water and canopy, Θ 0 is the potential temperature, C M and C M w are the BTC for momentum between an entire surface and atmosphere and between water surface MATSIRO uses C M . T w and T c are calculated in Section 3.1. Θ 0 is given by where R dry is the gas constant of dry air. Although the original MATSIRO fixes Θ 0 at 300 K, MATCRO calculates the value according to Campbell and Norman (1998). C M and C M w are given by Now we have six independent equations, Eqs. 24, 25,26,27,37, and 38, for six unknown variables, C Ew , C Ec , C Hw , C Hc , C M , and C M w , respectively. Therefore, we can determine the values of these variables by numerically solving Eqs. 24 to 38.
The numerical method is described in Masutomi et al. (2016).

Canopy water balance 10
The main purpose of this module is to calculate the fraction of wet canopy (f cw ) which is used for simulating energy balance at canopy (Section 3.1). To calculate f cw , this module calculates water balance at canopy. Although the module is based on the original MATSIRO, the amount of water that canopies can hold was replaced by using the method described in Penning de Vries et al. (1989). The variable f cw is given as

15
where w c is the amount of water stored in canopy and w cap is the water capacity of the canopy. The w c is calculated by solving the canopy water balance, which is given by where ρ w is the density of water, I c is the amount of precipitation intercepted by canopy, D g is the amount of water that falls from the canopy onto surface water due to gravity, and E c is the amount of water that evaporates from the canopy (Eq. 7). I c 20 depends on the amount of precipitation (P r ) and LAI (L) and is given by where f int indicates the interception efficiency of precipitation by canopy. According to Rutter et al. (1975) and Penning de Vries et al. (1989), D g and w cap are given as respectively, where D 1 and D 2 are parameters (Rutter et al., 1975), and W sh is the shoot dry weight, which is calculated in the CGM (Eq. 127).

Soil water and heat transfer
This module calculates heat and water transfer in soil. The main role of this module is to determine the temperature at a soil surface (T s (0)), which is used for simulating energy balance of the surface water (Section 3.1). Although this module is based Soil temperature at a soil depth of z from the soil surface (T s (z)) is calculated from the gradient of heat flux in the soil as follows: 10 where c hs is the volumetric heat capacity of the soil and G s (z) is the heat flux at a soil depth of z and is given from the gradient of soil temperature Here, k ts is the soil thermal conductivity. In Eq. 46, we assumed that heat flux at the bottom of the soil layer (z = z max ) is zero.
z max is a simulation setting parameter. When solving Eqs. 45 and 46, the heat flux from surface water to soil (G ws ), calculated 15 in Eq. 10, is used as a boundary condition. The parameter c hs is calculated from the heat capacities of soil components as follows.
where ρ s is the bulk density of soil, c pm is the specific heat of soil minerals, and w s (z) is the volumetric concentration of soil water. ρ s is a soil-type specific parameter determined by soil type at a simulation site, and c pm is given according to Campbell 20 and Norman (1998) . We note that the first term of the right hand side in Eq. 47 indicates the heat capacity of dry soil. Although the original MATSRIO model assigns a default value to the heat capacity of dry soil for all soil types, MATCRO-Rice calculates the value of the heat capacity of dry soil using the bulk density of soil and the heat capacity of soil minerals, as shown in the first term of Eq. 47. It should be noted that the effect of soil organic matter on c hs is not considered in MATCRO. The parameter k ts (z) in Eq. 46 is given by where k ts0 and k tss are the thermal conductivity of dry and saturated soils, respectively, K e is the Kersten number, and w sat is the volumetric soil water concentration at saturation. k ts0 and k tss are parameters. We set k ts0 =0.25 (Campbell and Norman, 1998), and k tss = 1.58 (Best et al., 2011). The parameter w sat is specific to soil type. Equations 48 and 49 for the calculation of k ts (z) are based on the equations developed by Best et al. (2011), while the original MATSIRO employs a different equation.
The variable w s (z) depends on the gradient of water flux and absorption by roots at a soil depth z and is given by where F s (z) and S s (z) are water flux and absorption by roots at a soil depth of z, respectively. For simplicity, the top soil layer is assumed to be saturated, because the surface above soil is flooded. Given the assumption, we do not need to explicitly simulate water flow from a flooded surface into soil. This assumption is not considered in the original MATSIRO. z sat is a simulation setting parameter. F s (z) is calculated from the gradient of water potentials as follows.
where K(z) is the hydraulic conductivity and ψ(z) is the water potential at a soil depth of z. F s (z) in the bottommost layer (z b < z < z max ) represents the base flow, and τ b is the recession constant for base flow. This model uses a simple model for simulating base flow developed by Hanasaki et al. (2008), although the original MATSIRO utilizes a more complicated model (TOPMODEL: Beven and Kirkby (1979)). z b is a simulation setting parameter, and τ b is determined as described in Hanasaki et al. (2008). K(z) and ψ(z) are given by Clapp and Hornberger (1978) as follows. 15 where K s and ψ s are hydraulic conductivity and water potentials at saturation, respectively, and B is a parameter that determines the relationship of hydraulic conductivity or water potentials between saturated and unsaturated soils. K s , ψ s , and B are soil-type specific parameters. S s (z) in Eq. 51 is calculated from the transpiration where E t is the transpiration calculated in Eq. 8 and z rt is a root depth calculated by the CGM (Eq. 131). In Eq. 55, we assumed that S s (z) has no dependency on soil depth.

Crop growth model
The main purpose of the CGM is to simulate rice yield and biomass growth for each organ during a growing period. The CGM 25 has four modules: "net carbon assimilation", "crop development", "crop growth", and "LAI, height, and root depth". Each module is described in detail in the following sections.

Net carbon assimilation
The main role of this module is to calculate net carbon assimilation (A n ) in canopy for simulating crop growth. In addition, the stomatal conductance per unit leaf area for both sides of the leave (g s ) is calculated for simulating roughness length (Appendix E). Although this module is based on the Big-leaf model (Sellers et al., 1992(Sellers et al., , 1996a used in the original MATSIRO, we refined two points in the calculation according to the approach described by de Pury and Farquhar (1997) and Dai et al. (2004). The 5 first refinement is that leaves in a canopy are divided into sunlit and shade leaves. Subsequently, A n per unit leaf area for each the sunlit and shade leaves are calculated. The second refinement is that A n for the entire canopy is calculated considering vertical distribution of nitrogen within the canopy.
A n for the entire canopy is given by 10 where A n,sn and A n,sh are net carbon assimilation per unit leaf area for sunlit and shade leaves, respectively, L sn and L sh are LAI for sunlit and shade leaves, respectively, and overbars represent the amounts per unit leaf area. A n,sn and A n,sh are defined by the difference between gross carbon assimilation and respiration as follows: where A g,x and R d,x are gross carbon assimilation and respiration per unit leaf area, respectively, and the suffix x indicates sn 15 or sh. L sn and L sh are given as follows.
where f sn (l) is the fraction of sunlit leaves at a LAI depth of l and is defined as follows: 20 where F denotes distribution of leaf orientation and θ is a zenith angle of the sun (Appendix B). The effect of photosynthesis down-regulation due to acclimatization to elevated CO 2 is represented as follows: where A g ,x is gross carbon assimilation per unit leaf area for sunlit and shade leaves without photosynthesis down-regulation, 25 f dwn is the factor for photosynthesis down-regulation, γ gd and γ g are parameters that characterize the response to increased CO 2 , and C 0 is the base concentration of CO 2 . The Eqs. 61 and 62 are based on Arora et al. (2009), although the original MAT-SIRO does not consider the effect of photosynthesis down-regulation. We set γ gd = 0.42, γ g = 0.9, and C 0 = 288 according to  Arora et al. (2009). The calculation for A g ,x and R d,x is based on the leaf photosynthesis model developed by Collatz et al. (1991). In their model, A g ,x is determined by three limiting factors: Rubisco, light, and sucrose synthesis, as follows: where ω c,x , ω e,x , and ω s,x are Rubisco-limited, light-limited, and sucrose-limited carbon assimilation per unit leaf area, respectively. To implement smooth transition between each limited state, A g ,x is determined practically by solving the following 5 two equations (Sellers et al., 1996b): where β ce and β pc are the parameters that determine the smoothness of transition between each limited state. β ce is a cropspecific parameter and β pc is a parameter that does not depend on crop type. The variables ω c,x , ω e,x , and ω s,x are given Here, V mc,x and V ms,x are the maximum Rubisco capacity per unit leaf area for ω c,x and ω s,x , respectively, c i,x is the partial 15 pressure of intercellular CO 2 , [O 2 ] is the partial pressure of intercellular O 2 , Q x is the photon flux density for PAR absorbed per unit leaf area by sunlit and shade leaves, e is the quantum efficiency, Γ * is the light compensation point, and K c and K O are the Michaelis constant for CO 2 fixation and oxygen inhibition, respectively. We set [O 2 ] = 20,900 (Collatz et al., 1991). e is a crop specific parameter. V mc,x and V ms,x are given by where V max,x is the reference value for the maximum Rubisco capacity per unit leaf area of sunlit (V max,sn ) and shade (V max,sh ) leaves, s 1 , s 2 , s 3 , and s 4 are parameters that represent temperature dependence of V max,x on V mc,x or V ms,x . The variables s 1 and s 2 are parameterised in Masutomi et al. (2016), whereas s 3 is a parameter that does not depend on crop type and s 4 is a crop-specific parameter. Q t is given by V max,sn and V max,sh are defined by where V max (l) is the reference value for the maximum Rubisco capacity at a LAI depth of l. The vertical distribution of V max (l) depends on that of leaf nitrogen within canopy and is given by where K n is a parameter that represents the vertical distribution of leaf nitrogen, and V max (0) is the reference value for the maximum Rubisco capacity at the canopy top. V max (0) as well as s 1 and s 2 are parameterized in Masutomi et al. (2016), and we set K n = 0.3 (Oleson and Lawrence, 2013). Γ * , K c , and K O are given by where S is the ratio of the partition of RuBP to the caboxylase or oxygenase reactions of Rubisco.
Q x in Eq. (67) is defined by the following equation: 15 Here, Q x is the PAR absorbed by the entire canopy for sunlit (Q sn ) and shade (Q sh ) leaves. Q sn and Q sh consist of direct and scattered components and are given as 20 where Q sn,d , Q sn,s , and Q sh,s are the direct PAR absorbed by sunlit leaves, the scattered PAR absorbed by sunlit leaves, and the scattered PAR absorbed by shade leaves, respectively. These are described by (1 − f sn (l))dl, 25 where D d 1 (l), S d 1 (l), and S u 1 (l) are calculated by the LSM (Eqs. 12 to 14) and k q is a constant that transfers the radiant flux density to photon flux density.
R d,x in Eq. 57 is given by the following equation: where f d is a respiration factor and crop-specific parameter, whereas s 5 and s 6 are parameters that are not crop-dependent. It 5 should be noted that A n,x can be calculated using the equations described in this section (Eqs. 57 to 85) if c i,x is given.
A n,x should be equal to the CO 2 flux between the leaf interior and boundary layer and the CO 2 flux between the leaf boundary layer and the atmosphere. If these requirements are fulfilled the following equation can be derived: where c a is the partial pressure of atmospheric CO 2 , c s,x is the partial pressure of CO 2 at the leaf boundary layer for sunlit 10 and shade leaves, g l is the leaf boundary conductance for vapour per unit leaf area, and g st,x is the stomatal conductance for vapour per unit leaf area for sunlit and shade leaves. From Eq. 86, c i,x and c s,x are defined by The parameters c a and g l are given by 15 c a = (C a * 10 −6 )P a , where w H2O is a constant for the molar weight of vapour, g a is the leaf boundary conductance for heat per unit leaf area (for both sides of the leaf), c h is the leaf transfer coefficient for heat and is a crop specific parameter, U c is the mean wind speed 20 in the canopy (Appendix F). Note that Eqs. 90 and 91 are based on Maruyama and Kuwagata (2008), whereas the original MATSIRO uses C h instead of g a /2 in Eq. 90.
A n,x meets the Ball-Berry relationship (Ball, 1988), which describes the relationship between A n,x , g st,x , and other environmental conditions. The Ball-Berry relationship is given by 25 where m and b are the slope and intercept of the Ball-Berry relationship, and h s,x is the relative humidity at leaf boundary. It is noteworthy that b indicates the stomatal conductance when A n,x is equal to or less than zero (Baldocchi, 1994) and that the effect of water stress on b is not considered in MATCRO-Rice because the surface is flooded. The variables m and b are crop specific parameters, and h s,x is defined by where e s,x is the vapour pressure at leaf boundary and e sat is the saturated vapour pressure. The variable e s,x is expressed as e s,x = (e a g l + e i g st,x )/(g l + g st,x ), where e a and e i are the vapour pressure in the air and leaf, respectively. Eq. 94 is derived from the fact that the water vapour flux from the stomata to leaf surface is equal to the water vapour flux from the leaf surface into the atmosphere, which is shown in the following equation: The parameters e a , e i , and e sat are given by 10 where e i is assumed to be saturated. Last, g s is given by the following equation: where g st is the stomatal conductance for vapour per unit leaf area for both sides of the leaf.

Crop development
The crop development module calculates DV S, which is an index used to quantify developmental stage of crops. DV S is 20 mainly used for determining the timing of transplanting, heading, and harvesting. In addition, DV S is used for partitioning of carbon assimilation into each organ and for estimating LAI and height. This module is based on the formulation by Bouman et al. (2001). DV S is calculated from where GDS is the growing degree seconds at t, mGDS is GDS required until maturation, DV R is the development rate at t, T 0 is the melting temperature of water, and T b , T h , and T o are the minimum temperature, maximum temperature, and optimal temperature for development, respectively. The value of mGDS is parameterized in Masutomi et al. (2016), and T b , T h , and T o are crop-specific parameters. T 0 is a physical constant (Table 4). It should be noted that DV S = 0 represents sowing and DV S = 1 represents maturation. Furthermore, we introduce two parameters that represent the timing of emergence During the transplantation of rice seedling, the seedlings enter transplanting shock, which prevents shoot growth (Bouman et al., 2001). In MATCRO-Rice, the transplanting shock period is defined by DV S, where trDV S is DVS at the time when transplanting shock starts and teDV S is DVS at which transplanting shock ends. Both trDV S and teDV S are parameterized

Crop growth
This module calculates the growth of organs and reserves. The organs considered in MATCRO-Rice include leaf, stem, panicle, and root. In addition, the model considers glucose reserves in leaves and starch reserves in stem. All carbon assimilated in leaves through photosynthesis is first stored in leaf in the form of glucose. Then, the stored glucose is partitioned to each organ and 15 stored in the stem when the amount of the stored glucose exceeds the critical rate to dry weight of leaf. This module is based on MACROS (Penning de Vries et al., 1989).
The dry weights of each organ and reserve are expressed by where W lef , W stm , W pnc , W rot , W stc , W glu are the dry weight of leaves, stems, panicles, roots, starch reserves, and glucose reserves at t, respectively, W lef,0 , W stm,0 , W rot,0 , and W glu,0 represent the initial dry weight at emergence of each organ and reserve, G R,lef , G R,stm , G R,pnc , G R,rot , G R,stc , and G R,glu are the growth rates of the corresponding organ and reserve, L S,lef is the loss rate of leaves due to leaf death, R M,stc is the loss rate of starch reserves in stem due to remobilization, t e is the time at emergence after sowing, and W lef,0 , W stm,0 , W rot,0 , and W glu,0 are simulation setting parameters.

5
The glucose reserve in leaf is supplied through photosynthesis in leaves and remobilization from the stem. Thus, the supply of glucose is given by where, S glu is the supply of glucose to leaf reserve, A n is the net carbon assimilation calculated in Eq. 56, and C CO2,glu and C stc,glu are the conversion factors from CO 2 or starch to glucose, which are chemically determined (Table 4). We assumed 10 that the partition of glucose in leaves to each organ occurs if the following equation is met: where δt is one simulation time step, k glu is the critical ratio at which the partition of glucose happens, and δt is a simulation setting parameter. We set k glu = 0.1 (Penning de Vries et al., 1989). When Eq. 111 is met, the amount of glucose that exceeds the critical ratio is partitioned to each organ and reserve according to the following equation: where G P,glu is the amount of glucose partitioned to each organ and reserve. The growth rate of each organ and reserve is expressed as follows:

20
G R,pnc = G P,glu P R,sh P R,pnc C glu,pnc , where P R,sh is the ratio of glucose partitioned to shoot, P R,lef and P R,pnc are the partition ratios of glucose from shoot to 25 leaf and panicle, f stc is the proportion of glucose allocated to starch reserve in stem, C glu,lef , C glu,stm , C glu,rot , C glu,pnc , and C glu,stc are dry weight of corresponding organs and reserves that are produced from the unit weight of glucose. f stc , C glu,lef , C glu,stm , C glu,rot , and C glu,pnc are crop-specific parameters. f stc is parameterized in Masutomi et al. (2016). We set the values of C glu,lef , C glu,stm , C glu,rot , and C glu,pnc according to Penning de Vries et al. (1989). C glu,stc is a chemical constant. If Eq. 111 is not met, glucose is not partitioned into each organ and reserve, except as the glucose reserve in leaf. Therefore, the growth rate of each organ and reserve are calculated as follows: The partition ratios to each organ are given as where DV S rot1 , DV S rot2 , DV S lef 1 , DV S lef 2 , DV S pnc1 , and DV S pnc2 represent the DV S values at which corresponding 10 partitions change, P rot is the ratio of partitioned glucose to the roots at DV S < DV S rot1 , and P lef is the ratio of glucose partitioned to the leaf and glucose partitioned to shoot at DV S < DV S lef . DV S rot1 , DV S rot2 , DV S lef 1 , DV S lef 2 , DV S pnc1 , DV S pnc2 , P rot , and P lef are crop-specific parameters and are parameterized in Masutomi et al. (2016). In Eq. 121, we assume that no glucose is partitioned to shoot during transplanting shock (teDV S < DV S ≤ teDV S). It is important to note that transplanting shock is considered only when transplanting is conducted.

15
Loss of leaf dry weight due to leaf death (L S,lef ) and remobilization from starch reserve in stem (R M,stm ) occur after heading and they are defined as follows where r dd,lef and r rm,stc represent the ratios of leaf death and remobilization. r dd,lef varies with DV S as follow: where r d1,lef is the ratio of leaf death at harvest (DV S = 1) and it is parameterized in Masutomi et al. (2016). We set r rm,stc = 1.16 * 10 −6 , assuming that all starch stored in stem is remobilized in 10 days after heading (Bouman et al., 2001).
Last, the dry weight of shoot (W sh ), used in Section 3.4, is given by

LAI, crop height, and root depth
Leaf area index (L), crop height (h gt ), and root depth (z rt ) are expressed as where SLW is the specific leaf weight, SLW mx and SLW mn are the maximum and minimum values of specific leaf weight, 15 respectively, k SLW is a parameter that determines the relationship between DV S and specific leaf weight, h gt,aa , h gt,ab , h gt,ba , and h gt,bb are parameters that define the relationship between LAI and crop height, z rt,mx is the maximum root depth, and r rt is the root growth rate. The allometric equations for estimating crop height (Eq. 130) is based on Maruyama and Kuwagata (2010).

Crop yield
Crop yield is calculated from dry weight of the panicle at maturity as follows: where Y ld is the crop yield, W pnc,mt is the dry weight of the panicle at maturity, k yld is the ratio of the crop yield to W pnc,mt .

25
The variable k yld is a crop specific parameter and it is parameterized in Masutomi et al. (2016).
We developed a new LSM-CGM combined model for paddy rice fields called MATCRO-Rice, which is fully described in the present paper. MATCRO-Rice has two features: (i) The model can consistently simulate LHF, SHF, biomass growth for each organ, and crop yield by exchanging variables listed in Table 2; (ii) The model considers water surface in paddy rice fields.
According to our literature survey, MATCRO-Rice is the first LSM-CGM combined model for rice that employs these two 5 features.
The first feature enables us to apply the model to a wide range of integrated issues. For example, by using MATCRO-Rice, we can assess the impacts of paddy rice fields on climate through heat and water fluxes and consistently assess the impacts of climate on rice productivity. Osborne et al. (2009)  The current version (Ver. 1) of MATCRO-Rice has a couple of major limitations. First, nitrogen dynamics is not included in MATCRO-Rice, although it is well known that nitrogen stress significantly affects crop growth, and hence LHF and SHF. This 20 indicates that MATCRO-Rice simulates LHF, SHF, biomass growth, and crop yield with no nitrogen stress. To apply the model to the site with nitrogen stress, it is necessary to include nitrogen dynamics. This feature is an important future challenge.
Second, the impact of water stress on crop growth is not considered in MATCRO. This limitation is not considered a problem in irrigated land but in rain-fed land. If the model is applied in rain-fed lands, the model needs to be improved.
6 Code availability 25 The source code of MATCRO will be distributed at request to the corresponding author (Yuji Masutomi: yuji.masutomi@gmail.com).
The website for MATCRO-Rice will be developed in the near future.
The air density (ρ a ) and the specific humidity at saturation (Q sat ) are calculated physically according to the equation for the state of dry air and the Clausisu-Clapeyron equation, respectively, as follow: where T a is air temperature, P a is air pressure, T x is temperature of the canopy (T c ) or surface water (T w ), T 0 is the melting temperature of the water, R dry and R vap are the gas constants of the dry air and vapour, respectively, e sat (T 0 ) is the vapour pressure at melting temperature of the water, and λ is the latent heat of vaporisation. T a and P a are meteorological inputs (Table 1). T x (T c or T w ) is calculated in Section 3.1. The other parameters are physical constants (Table 4).
Appendix B: Zenith angle θ 10 According to Goudriaan and van Laar (1994), zenith angle of the sun (θ) is calculated as follows.
δ s = − arcsin(sin(23.45(2π/360)) cos(2π(D oy + 10)/365)), where L t is the latitude in radians at the simulation site, δ s is the declination of the sun, h arg is the hour angle from noon 15 (h r = 12), D oy is the number of days from Jan 1 at the simulation site, and h r is the local time at the simulation site.

Appendix C: Coefficients for radiation equations
The coefficients for radiation equations  are calculated as follows: where i indicates the wavebands of radiation (i = 1: PAR; i = 2: NIR), r i and t i are the leaf reflectivity and transmissivity, respectively, F is the distribution of leaf orientation, d f is a scattering factor, A 3,i is a new variable introduced in Eqs. C2 and C3, L is the LAI, r g is the surface albedo for shortwave radiation, D d i (0) and S d i (0) are direct and scattered downward 15 radiant flux density at the canopy top, respectively, and θ is the zenith angle of the sun. r i and t i are crop-specific parameters determined by Sellers et al. (1996b). F is set to 0.5 from the assumption of random leaf orientation (Goudriaan and van Laar, 1994), and d f is sec(2π(53/360)) (Watanabe and Ohtani, 1995). A 3,i is defined in Eq. C8, L is calculated in the CGM (Eq. 128), and r g for surface water is given in Maruyama and Kuwagata (2010). D d i (0) and S d i (0) are given in Eqs. 15 and 16, respectively, and θ is calculated in B1. 20 It should be noted that a i , A 1,i , and A 2,i are not variables determined by constant parameters, while C 1,i , C 2,i , C 3,i , C 4,i , and A 3,i are variables.
Here, z M s , z T s , and z Qs are the roughness lengths of surface water for momentum, temperature, and specific humidity, respectively. In this model, we assume z M s , z T s , and z Qs = 0.001 m (Kimura and Kondo, 1998). c m , c h , and c e are the leaf transfer coefficients for momentum, temperature, and specific humidity, respectively. c m an c h are crop-specific parameters, while c e is calculated in Eq. E15. h gt and L are crop height and LAI, respectively, and are calculated in the CGM (Eqs. 130 and 128).
g s is the stomatal conductance per unit leaf area for both sides of the leaf (Eq. 99). U c is the mean wind speed in the canopy and is calculated in Appendix F. A + , C 0 X , C ∞ X , z + M , z + X , z + * , P 1 * , P 2 * ,P 3X , P 4X , F X are the intermediate variables, and κ is the Karman constant. The symbol " * " indicates "M", "T", or "Q", and the symbol "X" indicates "T" or "Q".

Appendix F: Mean wind speed in the canopy
Mean wind speed in the canopy (U c ) is expressed as where U h is the reference wind speed, and γ m is the coefficient of exponential decrease for wind speed in the canopy.  Ag,x mol(CO2) m −2 (l)s −1 61 gross primary production per unit leaf area of sunlit (Ag,sn) and shade(A g,sh ) leaves A g ,x mol(CO2) m −2 (l)s −1 65 gross primary production without photosynthesis down-regulation per unit leaf area of sunlit (A g ,sn ) and shade(A g ,sh ) leaves