Liu Huan (2021). Modeling of annual net primary production of a forest in the Taramakau Valley, Westland, New Zealand. Journal of Environment and Health Science (ISSN 23141628). 2021 (02). https://doi.org/10.58473/JEHS0007
2023. Copyrights Certificate Registered by Brock Chain Technology: Brock Chain ID. (ae8ca7ba8666537a8c2d23fae459836289fff107fa292e1f62b1923f46c45c7c) ;
2016. Copyrights Register Information: The majority of these materials are registered as book '著作权人：刘焕；作品：《研究生文凭进展（第三版）》' 2016, which can be cataloged in National Copyright Database: http://qgzpdj.ccopyright.com.cn/
2016. 版权注册信息：本文大多数内容已经以图书形式登记注册在全国版权数据库，登记入库信息：著作权人：刘焕；作品：《研究生文凭进展（第三版）》 2016;可在全国版权登记数据库检索 http://qgzpdj.ccopyright.com.cn/
The formally published serials is the printing <Journal of Environment and Health Science (ISSN 23141628)>, and the serials NO. is the month/year when the materials is accessible on this website, authorized by publisher；
正式发表的期刊是印刷版《环境与卫生科学杂志（ISSN 23141628）》，期刊期号为文章内容在本网站上网年/月，出版人许可自行正式发表。
Originality Certificate: The originality of English content is 87% tested by Turnitin (International).
All Copyrights Reserved.
Latest revised on 31/05/2023.
Article Download (Google Drive) Open Access RepositoryZenodo
论文下载（提取码: 4qns，百度网盘链接与提取码仅仅允许本香港主机网站用户使用，禁止转载）
Scilit Database Retrieval URL : Link Researchgate Retrieval URL: Link
Web of Science Retrieval URL: https://publons.com/wosop/publon/66151053/
Bielefeld Academic Search Engine: Link
Article 6. Modeling of annual net primary production of a forest in the Taramakau Valley, Westland, New Zealand
Advisor: George Perry, School of Environment, The University of Auckland
Abstract Model provides a useful tool for the estimation of annual net primary production (NPP) in the terrestrial ecosystem. This project used Vensim software to design a numerical model of annual net primary production in a forest in the Taramakau Valley, Westland, New Zealand, based on the empirically derived correlations between forest growth and available resources. Atmospheric transmissivity, light use efficiency (LUE), light saturation point and water limitation factor were parameterized in this model. Light use efficiency (LUE) is considered as a function of atmospheric transmissivity in this model, which well indicates the daily variation in LUE. Sensitivity analysis was conducted to validate the parameters in this model. The results of this model have been validated by comparing with field studies or other models. The estimated NPP is 938.335g C/m2 in 2007, which is a reasonable value. Finally, the limitations and recommendations for future model development have been discussed, with application on the estimation of production in carbon sink forest.
Key words: Numerical Model, Net Primary Production, Light Use Efficiency, Atmospheric Transmissivity, Light Saturation Point, Water Limitation Factor, Carbon Sink Forest.
论文题目：新西兰西部森林年净初级生产力（NPP）的数值模拟（英文）
摘要：模型预测法为估算陆地生态系统年净初级生产力（NPP）提供了一个有效的工具。此课题采用 Vensim 软件对新西兰西部地区taramakau山谷森林的年净初级生产力进行了数值模拟。模型输入与输出关系是森林生长量和可利用资源之间的经验数值关系。该模型的四个主要参数包括大气透过率、光能利用效率、光饱和点和可利用水资源。其中光能利用效率（LUE）在数值模型中被定义为大气透过率的函数，大气透过率可以有效反映光能利用效率的日变化。此课题应用了灵敏度分析方法对 9 个模型参数的有效性进行了验证。通过与现场测定值或与其他模型预测值的比较，此课题模型的预测结果已得到充分验证。模型预测的陆地生态系统年净初级生产力（NPP）在 2007 年为 938.335g C /平方米，在合理的估值区间内。模型发展的局限性和趋势也在最后进行了探讨，并且应用于碳汇林的生产力估测。
关键词：数值模拟、年净初级生产力、光能利用效率、大气透过率、光饱和点、干旱因子、灵敏度分析、碳汇林
1. Introduction Over the past 30 years, a considerable number of global vegetation models have been developed to investigate the land carbon sink, which is one of the major concerns on the implementation of Kyoto Protocol. Usually, an essential part of these models is the submodel of simulating net primary production (NPP). Net primary production is the mainly (even only) source of carbon sink in the terrestrial ecosystem, which explains most annual carbon fluxes between the atmosphere and biosphere. Most NPP simulation models are based on the interactions between vegetation growth and its available resources, such as radiation, available soil water, soil nutrient, and ambient temperature.
Among these resources, radiation is the only energy source for vegetation growth, and hence understanding the conversion efficiency of irradiance into biomass (light use efficiency) becomes essential for the simulation of NPP. Light use efficiency (LUE) is often assumed to be constant in the previous models (Whitehead et al., 2001). However, the daily variations in LUE did occur in New Zealand forest due to environmental limitations (Dugan & Whitehead, 2006). This becomes an indication for the improvement of NPP simulation model in this project.
This project uses Vensim software to design a simulation model of annual net primary production in a Wineberry and Fuschia dominated forest in the Taramakau Valley, Westland, New Zealand, which employs a bottomup approach on the basis of the interactions between the forest growth and its available resources (including radiation, precipitation, and temperature).
2. Review of terrestrial net primary production modeling Generally, net primary production models can be divided into two groups: empirical modeling and theoretically modeling. The first group models are based on the empiricallyderived correlations of net primary productivity with available resources, and the other group of models are based on biochemicalmechanistic processes (Mulligan & Wainwright, 2004). In this section, both types of these models have been critically reviewed, with emphasis on their representative characteristics and limitations. Finally the choice of modeling types is discussed critically.
The Miami model, which empirically derives the associations between NPP and the mean annual temperature and precipitation, distinguishes plant types in the derivation of the correlation functions (Leith, 1975). This model can make reasonable estimations only for current NPP rates and their distribution, but not for the future NPP. This is mainly due to two reasons. Firstly, the Miami model does not consider the change of vegetation density. Secondly, the model assumption that vegetation distribution responses immediately to the climate change is sometimes not reasonable (Mulligan & Wainwright, 2004).
King et al., (1997) and Post et al., (1997) subsequently modified the Miami model. They quantified the changes in ecosystem carbon density by using a model of combining the empirical with the βfactor which added a carbon dioxide fertilization term to the Miami model. This βfactor was derived from the Farquhar photosynthesis formulation. A parameter (k1) was introduced as a scaling factor, which accounted for the efficiency losses in the actual increase of vegetation carbon density as a result of carbon dioxide fertilization. Particularly, this model distinguished relatively big amount of vegetation types in the global scope. However, the submodel of βfactor did not take into account of the C4 photosynthesis pathway (King et al., 1997; Post et al., 1997).
Polglase and Wang (1992) also developed a biochemical model based on the βfactor. They distinguished ten types of biomes in this model. Baseline productivity was generated for each biome, which was based on the averages of the characteristic climate and productivity. Then βfactor modified the baseline productivity values. However, the parameters used in this βfactor were derived by using the average temperature. As the authors themselves suggested, significant inaccuracy would occur when there was significant variation in seasonal temperature (Polglase & Wang, 1992).
Lenton (2000) designed a quasibiochemical model, in which photosynthesis was a function of temperature and atmospheric carbon dioxide, and respiration was depended on temperature and the vegetation carbon density. However, a particular problem can be caused by its assumption that the photosynthetic area is constant even though vegetation density may change (Mulligan & Wainwright, 2004).
Svirezhev and von Bloh (1997) developed another type of quasibiogeochemical model. In this model, the maximum net primary productivity was generated as a parameter, which was modified by factors indicating the temperature response, water stress, carbon dioxide fertilization, and intraspecific competition. However, compared with other biochemical models, the functions interpreting the carbon dioxide fertilization and competition by Svirezhev and von Bloh (1997) were simply integrated into this model, which made it difficult to be understood (Mulligan & Wainwright, 2004).
Leuning et al., (1995) presented a processbased model. The model is based on the leaflevel measurements of photosynthesis and respiration, which is then scaled up to the canopy level. Submodels of radiative transfer, energy balance, evaporation and photosynthesis are integrated. Stomatal conductance and photosynthesis is modified by a site waterbalance submodel. This model successfully predicts the forest canopy phenology (Leuning, et al., 1995). However, this model requires speciesspecific data for parameterization, which is not feasible for the application at a regional or global scale.
TRIFFID model calculates NPP by simulating photosynthesis and respiration. In this model, photosynthesis is calculated using the full Farquhar formulation, and respiration is determined by the rate of Rubisco carboxylation, both of which were modified by a water stress factor and the canopy level factor using Beer’s law. TRIFFID model is considered to be a relatively comprehensive vegetation model, incorporating intensively detailed information of biochemistry, physiognomy and forest dynamics. However, its complication leads it to be somewhat opaque and mathematically intractable (Cox, 2001; Mulligan & Wainwright, 2004).
Huntingford et al., (2000) subsequently presented a simplified version of the TRIFFID model, which was designed to investigate the qualitative response of the global vegetation to climate change. However, in this model, the description of a basic respiration required very detailed biochemical expression, and a single generic vegetation type relied on a broad parameterization, which still led this model to be inconvenient (Huntingford et al., 2000).
BIOME3 is one of the most comprehensive and complicated models for the estimation of NPP. It is also impenetrable and computationally intensive. In this model, photosynthesis is calculated by two simultaneous equations. The first relies on the Farquhar model modified by light and Rubisco limitations. An optimized daytime productivity and Rubisco carboxylation rates are given. The second is based on a stomatal conductance model that incorporates the effects of water stress. Then the NPP and leaf area index are adjusted by calculating an ‘extensive range’ of values from experiments. This model distinguishes different vegetation types by their maximum transpiration rates. Nevertheless, this complicated model requires a relatively high cost to apply (Haxeltine and Prentice, 1996; Mulligan & Wainwright, 2004).
The majority of these above models aim to simulate NPP in the global scale. However, because the environmental processes are usually scaledependant, it is still not clear how these models interpret the data changes with the change of scales (Zhang & Wainwright, 2004).
Discussion of choice between empirical modeling and theoretical modeling: For the empirical models, the repeatable accuracy of modeling results, the availability of model input variables (such as convenience and economic cost accessing the modeling input data), the sensitivity of model parameters that influences the inaccuracy of modeling results (the more sensitive for the model parameters, the higher inaccuracy of modeling results, because all the parameters are estimated according to the empirical data, if the parameter is very sensitive, a tiny inaccuracy in parameter estimation will lead to a big error in final results) and the economic tradeoff between the fieldspecific parameter determined by field sampling and estimation only of parameters are the main consideration factors to build a model; For the theoretical modeling, not only the empirical modeling results are assessed according to the above factors, but also the underlying theories are made more understandable. However, for the modeling of NPP at macroscale, this article does not choose the theoretical biochemistry modeling as the first choice because biochemical modeling can only give better understanding at molecular scale rather than macroscale, but these biochemical parameters increase the sensitivity of model parameters which may increase the inaccuracy of final modeling results.
3. Model description and implementation This section will describe the equations used and their corresponding implementations in the Vensim system. Generally, the selection of these subequations facilitates the availability of input data. The final output of ‘Net primary production’ is calculated by an integrated equation. Among the five input variables, I employ the ‘Lookup’ function for the daily information input, while the ‘IF THEN ELSE’ function is used to facilitate the input process of quarterly data.
This model is based on a daily time scale with a spatial scale of m2. In the initial settings of this model, the initial time, final time, time step, and time unit is ‘1’, ‘365’, ‘1’, and ‘day’, respectively.
3.1. Absorbed Photosynthetically Active Radiation (PAR)
Fig 1. The interactions in submodel of ‘Absorbed PAR’. (See PDF Article).
The daily absorbed PAR is considered as a function of ‘Daily radiation’ above canopy and ‘Daily fraction of absorbed PAR’. In the equation editing box of both ‘Daily radiation’ and ‘Daily fraction of absorbed PAR’, a ‘Lookup’ function is used for the daily data input. The unit of ‘Daily radiation’ is MJ/m2, while the ‘Daily fraction of absorbed PAR’ is without unit.
Fig 2. The editing box of ‘Lookup’ function for Daily radiation.(See PDF Article).
In the [Graph Lookup] box of ‘Daily radiation’, the days (x) axis is scale from 1 to 365, and the (y) axis is scaled by entering a value of 35 for the [Ymax]. Then the daily radiation data can be input in the left box (Fig 2).
Fig 3. The editing box of ‘Lookup’ function for ‘Daily fraction of absorbed PAR’.(See PDF Article).
Similarly, in the [Graph Lookup] box of ‘Daily fraction of absorbed PAR’, I scale the days (x) axis from 1 to 365, and (y) axis from 0 to 1. Then the data of daily fraction of absorbed PAR can be input in this box (Fig 3).
Fig 4. Equation editing for ‘Day’ variable.(See PDF Article).
Then a shadow variable <Time> is added, as a determinant to ‘Day’ (Fig 1). In the equation editing box of ‘Day’ variable, an equation of ‘Day = MODULO(Time, 365)+1’ is typed (Fig 4).
The variable of ‘Radiation’ is determined by both ‘Daily radiation’ and ‘Day’. Then, the variable of ‘Radiation’ can be interpreted as ‘Radiation = Daily radiation(Day)’ in the equation editing box (Fig 5). Similarly, for the variable ‘Fraction of PAR’, an equation of ‘Fraction of PAR = Daily fraction of absorbed PAR(Day)’ is added.
Fig 5. Equation editing for ‘Radiation’ variable.(See PDF Article).
However, not all the radiation can be utilized by forest canopy. To reflect this fact, I use ‘Light utilization’ linked with ‘Light saturation point’. The variable of ‘light saturation point’ is set to be a constant value (Q) of 22.04 MJ/m2. Then in the variable of ‘Light utilization’, an equation of‘IF THEN ELSE(Radiation<Light saturation point, Radiation,Light saturation point)’ is added.
Finally, the variable of ‘Absorbed PAR’ is determined by both ‘Fraction of PAR’ and ‘Light utilization’. The PAR is estimated to be 45% of total solar radiation. In the equation editing box, the equation is ‘Absorbed PAR = Fraction of PAR*Light utilization*0.45’.
3.2. Water limitation factor Water deficit is one of the major factors limiting forest productivity. Hence, there will be a limitation factor of available water h(X) (Nicholas et al., (2007)) in this model, where X is the proportion of cumulative precipitation (∑P) to cumulative evapotranspiration (∑E) over a certain period. Parameter s is the slope of the curve relating saturation water vapor pressure to temperature, and γ is the psychrometric constant. Parameter a is empirical coefficient. In this model, the sitespecific values of parameter a, s, and γ is estimated to be 1.8, 1.12, and 0.066, respectively, which are different from the value published in Nicholas et al., (2007) article.
‘Water limitation factor’ is determined by the input variables of ‘Radiation (3 months)’ and ‘precipitation (3 months)’ with a unit of ‘MJ/m2’ and ‘mm’, respectively. A shadow variable <Time> is added to determine these two input variables (Fig 6), which require quarterly climate data. For example, in 2007, the quarterly radiation values are 1500.24MJ, 511.28MJ, 730.27MJ, and 1598.99MJ from January to December at this site, and the 3month precipitations are 282mm, 503.6mm, 369mm, and 532.4mm from January to December. Data are derived from NIWA.
Fig 6. The interactions in submodel of water limitation factor.(See PDF Article).
I employ the IF THEN ELSE function to input these quarterly data. In the variable of ‘Radiation (3 months)’, an equation of‘IF THEN ELSE(Time<91, 1500.24, IF THEN ELSE(Time<182, 511.28, IF THEN ELSE(Time<274, 730.27, 1598.99)))’ is added (Fig 7). Similarly, the ‘Precipitation (3 months)’ is interpreted as ‘IF THEN ELSE(Time<91, 282, IF THEN ELSE(Time<182, 503.6, IF THEN ELSE(Time<274, 369, 532.4)))’. These two input variables determine the ‘Water balance’
Fig 7. The equation editing box for ‘Radiation (3 months)’.(See PDF Article).
Nicholas et al., (2007) equations contain two variables, the ‘Water balance’ and ‘Water limitation factor’ (Fig 6). In the equation editing box of ‘Water balance’, the equation of ‘1.8*"Precipitation (3 months)"/(1.12*0.45*"Radiation (3months)"/(1.12+0.066))’ is typed. I also use IF THEN ELSE function to edit ‘Water limitation factor’: ‘IF THEN ELSE(Water balance<1, Water balance, 1)’.
3.3. Light use efficiency Dungan and Whitehead (2006) reported that atmospheric transmissivity explained most of daily variability in light use efficiency (ranging from 75% to 90%) for Fuchsia and Wineberry species in 1999 in this proposed forest. Therefore, in this project, light use efficiency (LUE) will be mainly considered as a function of atmospheric transimissivity (ε):
LUE = y0 + a × exp(0.5×(ln(ε/x0)/b)2) equation 3
Based on their experiment results, the estimated values of parameter y0, a, x0, and b, are 0.28, 0.795, 0.18, and 0.78, respectively.
The atmospheric transmissivity is expressed as the fraction of irradiance reaching the top of the canopy after absorption and scattering by clouds and atmospheric turbidity (H/H0) (Bristow & Campbell, 1984). Bristow and Campbell (1984) proposed a model of estimating the atmospheric transmissivity (ε) which considered ε as a function of the difference between daily maximum and minimum temperature (ΔT). In this model, the parameters of aB, bB, cB are empirical coefficients with an estimated sitespecific value of 0.66, 0.23, and 0.8, respectively, which is different from Bristow and Campbell (1984) article.
Fig 8. The interactions in submodel of ‘Light use efficiency’.(See PDF Article).
A method which is similar to the ‘Daily radiation’ is employed for the input of ‘Daily temperature variation’. A ‘Lookup’ function is used for the ‘Daily T variation’ input. The shadow variable <Time> and a MODULO function are employed for the interpretation of variable ‘Days’. Then the ‘Temperature variation’ is determined by the ‘Daily T variation’ and ‘Days’ (Fig 8).
For the variable of ‘Atmospheric transmissivity’, the BristowCampbell equation ‘0.66*(1EXP(0.23*(Temperature variation^0.8)))’ is input in the equation editing box.
Then the ‘Light use efficiency’ is determined by ‘Atmospheric transmissivity’. The Dugan – Whitehead (2006) equation is typed in the equation editing box, which is ‘0.28+0.795*EXP(0.5*((LN(Atmospheric transmissivity/0.18)/0.78)^2))’.
3.4. Net primary production Finally, the ‘Net primary production’ will be calculated by an integrated equation:
NPP = PAR * FPAR * LUE * h(X) equation 5
Fig 9. The interaction between ‘Net primary production’ and its determinants.(See PDF Article).
In the equation box of ‘Net primary production’, a level variable is selected with an initial value of zero (Fig 10). Then the ‘Net primary production’ = INTEG [‘Absorbed PAR*Light use efficiency*Water limitation factor’] is typed in the below box.
Fig 10. The equation editing box of NPP.(See PDF Article).
4. Model test and validation The model equations can be tested using the ‘check model’ function in the Vensim system. In this section, the output of each submodel will be validated independently by comparing with either field studies or other model results.
4.1. Absorbed PAR In the submodel of absorbed PAR, the total estimated PAR is 1953.351MJ/m2 in 2007 (Fig 11). The mean fraction of absorbed PAR is 86.12%. Dugan and Whitehead (2006) did a field study in 1999 at the modeling site, and reported that the total PAR and the fraction of absorbed PAR were 3463MJ/m2 and 93%, respectively.
Fig 11. The absorbed PAR in 2007 in the modeling forest.(See PDF Article).
In this model, the PAR is estimated base on the assumption that the PAR is a constant proportion (45%) of total solar radiation. The total solar radiation data are measured by the meteorological station of Reefton Ews, which is approximately 6km away from the modeling forest. Surprisingly, the total estimated PAR in 2007 (1953.351MJ/m2) is significantly lower than the field measurements at the modeling site in 1999 (3463MJ/m2). In addition to the temporal and spatial variability in PAR, the different measurement methods used by Dugan and Whitehead (2006) may be also an important reason to explain this difference. In their field study, the PAR is directly measured by a quantum sensor (Li190SA, LiCor, Lincoln, Nebraska), which is different from the measurement equipment for total solar radiation used in the meteorological station.
The mean fraction of absorbed PAR in 2007 (86.12%) measured by Earth Resources Observation System Data Center (EROS) is also lower than the field measurement made by Dugan and Whitehead (2006) in 1999 (93%). Apart from the temporal variation in this fraction value, another reason was because the methods used by EROS calculated the absorbed PAR only and remove the PAR reflected by forest canopy, while Dugan and Whitehead (2006) calculated the intercepted PAR, which included both of them.
4.2. Light use efficiency
Fig 12. The estimated light use efficiency in 2007 in the modeling site.(See PDF Article).
The mean estimated light use efficiency by this model is approximately 0.65g C/MJ in 2007 (Fig 12). This result was consistent with Dugan and Whitehead (2006), who estimated that light use efficiency was approximately 0.60g C/MJ in 1999 in the modeling forest.
4.3. Water limitation factor
Fig 13. The estimated water limitation factor in 2007 in modeling site.(See PDF Article).
As can be seen in Fig 13, the productivity in 2007 was reduced by the water limitation factor only in the first 3 months. According to NIWA Climate database, the quarterly precipitation in 2007 was 282mm, 503.6mm, 369mm, and 532.4mm from January to December, which means that the first 3 months was the driest season in 2007. Therefore, the estimated value of water limitation factor matches to the precipitation records.
4.4. Net primary production
Fig 14. The estimated net primary production in 2007 in the modeling site.(See PDF Article).
The final output of net primary production in 2007 is 938.335g C/m2 (Fig 14). Dugan et al., (2004) conducted a simulation study using a processbased, biochemicalmechanistic model in the same forest in 2003, which showed that the estimated net primary production in 2003 was closely consistent with the estimation made by this project.
In addition, the net primary production of the modeling forest can be also obtained from the Earth Resources Observation System Data Center (EROS). According to their estimation, the average NPP was 806.7g C/m2 in 2007, which was slightly lower but still close to the result of this project. Therefore, the final output of 938.335g C/m2 in this project should be a reasonable value.
Compared with other global models, the estimated value of this project is lower than the simple TRIFFID, DEMETER, and TRIFFID models, but higher than BIOME3 model. Especially, the cumulative absorbed PAR in this project is 58 W/m2. This value reaches the light saturation point in TRIFFID model, above which little NPP increases with the increase of PAR in TRIFFID model.
5. Parameterization There are nine parameters in total in this model. For some submodels which have their own patent without publication of their parameterization methods, I used the reference values which are closest to the scenario of the modeling forest.
5.1. Atmospheric transmissivity There are three empirical coefficients (aB, bB, cB) used in this submodel. For the parameter aB, I use a reference value (0.66) calculated at a highland site which has similar topography to the proposed forest. For the bB and cB, Guillerdo et al., (2004) gave the equations for parameterization, which were calculated on the basis of the ΔT and φ, the daily difference between maximum and minimum temperature, and latitude of the proposed forest (42° 46' S), respectively. The ΔT data can be obtained from the NIWA Climate Database. I also built a simple model to calculate the bB, and cB values using Vensim software. The average calculated value of parameter bB and cB is 0.227 and 0.790, respectively. However, in this parameterization process, I used the daily ΔT data in 2007 only, which might not be sufficient to precisely determine these empirical coefficients (ideally, the data should cover the past ten years). Therefore, these are only approximate values.
5.2. Light use efficiency Dugan and Whitehead (2006) proposed an equation (equation 3) to quantify the atmospheric transmissivity effect on the light use efficiency (LUE), without publication of parameter values. Nevertheless, they reported their relationship for each species (Fuchsia and Wineberry) (Fig 16).
Fig 16. The relationship between daily values of LUE and atmospheric transmissivity for fuchsiaand wineberry.Sources from Dugan and Whitehead (2006).(See PDF Article).
In the parameterization process, I firstly converted these speciesspecific results into mixedspecies value, according to each species’ proportion of basal area in the modeling forest. This conversion equation is:
LUE (mixed species) = LUE (Fuchsia)*0.063 + LUE (Wineberry)*0.937
Then four paired results, which are (0.1, 0.43), (0.2, 0.62), (0.4, 0.49), and (0.6, 0.38), are obtained for the calculation of these parameter values in equation 3. The estimated parameter values of y0, a, x0, and b are 0.28, 0.795, 0.18, and 0.78, respectively.
5.3. Light saturation point The light saturation point (Q) used in this model was estimated from the global model TRIFFID (Fig 13). This model result showed a limitation value yearly for light utilization.However, I converted this yearly value into a daily value of 20.02MJ/m2.
However, this light saturation point which is based on a daily temporal scale sometimes may not properly explain the photosynthesis process (e.g. light saturation may occur when the daily radiation is below this limitation value if the hourly radiation is high during that day). Most physiological studies employ a time scale of seconds to determine the light saturation value. In practice, the NPP shows an asymptotic relationship with the absorbed PAR. However, this parameterization of Q assumes that the relation between NPP and PAR is linear.
5.4. Water limitation factor In equation 1, parameter ‘s’ is the slope of the curve relating saturation water vapor pressure to temperature, and ‘γ’ is the psychrometric constant. These two parameters are constant, which can be found in the references (“HyperPhysics: Thermodynamics, and The Asce Standardized Reference Evapotranspiration Equation”). However, the empirical coefficient ‘a’ is sitespecific. I used a reference value of 1.8 which was derived from a temperate Douglasfir forest in Australia.
6. Sensitivity analysis Sensitivity analysis not only contributes to the model validation, but also becomes essential for the future model development. There are a number of reasons to do sensitivity analysis for the modeling. Firstly, sensitivity analysis helps to better understand the behavior and interactions between the output and parameters; in the multicomponent model, the submodels can verified by sensitivity analysis; it helps to improve the model parsimony by rejecting the insensitive parameters; it contributes to identify which parameters require additional field research to strengthen the knowledge base; sensitivity analysis also examines which inputs are highly correlated to the final outputs (Mulligan & Wainwright, 2004).
Hamby, (1994) summarized a large number of sensitivity analysis methods for model parameters. In this project, the sensitivity analysis was implemented using Pearson’r correlation methods. Ten tests in total were employed for this correlation test. In each test, all the parameters were randomly selected within a reasonable range which was determined by historical records. Then the correlation (R) between each parameter and final output was calculated using IFA Online Statistic Service (Table 1). The higher the correlation (R), the more sensitive is the parameter.
Table 1. Sensitivity analysis for nine parameters in this model.(See PDF Article).
However, there are some drawbacks in this correlation test used for sensitivity analysis: the method is based on the assumption that the input/output relationship is linear (β is an exponential parameter in this model); and some input parameters strongly correlated to one another (such as aB and β) may result in significant input/output correlations (Hamby, 1994).
In each test, parameter aB, bB, cB, y0, β, x0, b, a, and Q are random selected within a range of [0.5, 0.9], [0.15, 0.30], [0.70, 1.00], [0.15, 0.40], [0.40, 0.70], [0.10, 0.30], [0.70, 1.00], [1.0, 4], and [18.00, 25.00], respectively. As can be seen from the results (Table 1), parameter aB, β, x0, and a are highly correlated with the final output, which show a high sensitivity in this model.
Particularly, the parameter Q is insensitive in this model. As discussed above, this light saturation point (Q) does have limitation in this model. Therefore, this parameter Q may be rejected to ensure the parsimony of the model, or requires further studies in the future development.
7. Conclusion 7.1. Application and Advantages This simulation model well estimated the net primary production in 2007, which has been validated in section 4.4. Among these five input variables, daily precipitation, daily radiation, 3month radiation, 3month precipitation, and daily variation of temperature can be downloaded from NIWA Climate Database. The fraction of absorbed PAR can be derived from Earth Resources Observation System Data Center (EROS), with a spatial scale of 1 km2. These fraction data are processed per 8 days, which can be downloaded in HDF format, and then can be transferred into excel for the calculation of mean values. Therefore, all of the data used in this model is readily and freely available for forest managers. Another advantage of this model is that it well indicates the daily variations in LUE, which improves the accuracy of simulating forest growth.
Another advantage of this numerical modeling is that the parameters can be estimated on the basis of empirical values to reduce the modeling cost. For example, the parameters of water limitation factor or light saturation point can be estimated according to the past study on the relevant species, which have been effectively done in my study. There are four steps in parameterization: estimation of parameters according to the empirical values; conducting sensitivity analysis for evaluating the accuracy and credibility of estimated parameters; validation of the estimated parameters by comparing the final modeling results with the field observations at experimental scale; final promotion of this modeling in larger scale. Generally, the more sitespecific or speciesspecific parameters required by field measurement, the higher cost for building a model. Then there is a tradeoff between the field measurement and the estimation based on the empirical value, which have been deduced in my modeling study.
7.2. Limitations and future development 7.2.1. Lack of considerations on the available soil nutrient and pest impact In practice, soil available nutrients, such as nitrogen, potassium, and phosphate, usually become a limitation for the forest growth. However, this model is based on the assumption that there are always sufficient soil nutrients in this proposed forest. In addition, this model also assumes that the proposed forest is pestfree, which contradicts the fact that New Zealand forest usually suffers from a number of disease, insects, and invasive mammals.
7.2.2. Temporal and spatial scale This model employed a daily time scale, which matches to the meteorological data from NIWA Climate Database. It may not be feasible to predict the future annual NPP using this model, because it seems to be impossible to predict the future longterm climate data on the basis of a daily scale.
Lots of environmental process and patterns are scaledependent. In this model, the m2 of spatial scale is used, which is also a common spatial scale in other NPP estimation models. The empirical data used for parameterization in this model is based on the local scale (the whole forest range), and then converted into 1 m2. Nevertheless, the relatively small scale makes it impossible to indicate the species composition and distribution. Further more, this model cannot be applied to other forests in NZ with different species composition, which requires additional field studies for the empirical parameterization.
7.2.3. Future model development This model is nonspatial, which does not consider the ‘species flow’. However, the ‘Light use efficiency’ used in this model can be significantly influenced by the change of species composition, because variations in LUE do occur among different species. Further more, the application of 3S technology on the biodiversity monitoring of vegetation ecosystem, novelly presented by the article 3 of this journal, facilitates the integration of species biodiversity on this carbon sink model in global scale.
In the parameterization process, I integrate the speciesspecific LUE of two species into a ‘mixedspecies’ value according to their proportion of basal areas in the proposed forest. However, this proportion may change with the change of species composition over a longer period. The future model development should integrate the species composition change into the determination of ‘mixedspecies’ LUE.
This is the revised materials in book “Proceedings for Degree of Postgraduate Diploma in Environmental Science (3rd Edition).” Published in 2016. The ‘chapter’ content mentioned in this article is in previous book. Revised on 25/01/2021. Secondly Revised on 01/09/2021. Thirdly Revised on 09/01/2022. Fourthly Revised on 29/11/2022. Fifthly revised on 24/04/2023. Latest revised on 31/05/2023.
References   Akihiko, I. A., & Takehisa, O. (2002). A simulation model of the carbon cycle in land ecosystems (SimCYCLE): a description based on drymatter production theory and plotscale validation. Ecological Modelling [Ecol. Model.]. 151, 143176.
https://doi.org/10.1016/S03043800(01)004732    Cleugh, H. A., Leuning, R., Mu, Q., & Running, S. W. (2007). Regional evaporation estimates from flux tower and MODIS satellite data. Remote Sensing of Environment [Remote Sens. Environ.]. 106(3), 285304.
https://doi.org/10.1016/j.rse.2006.07.007    Conover, W. J. (1980). Practical Nonparametric Statistics. New York: John Wiley & Sons.    Dungan, R. J., & Whitehead, D. (2006). Modelling environmental limits to light use efficiency for a canopy of two broadleaved tree species with contrasting leaf habit. New Zealand Journal of Ecology [N. Z. J. Ecol.]. 30(2), 251259.    Dungan, R. J., Whitehead, D., Mcglong, R., Allen, B., & Ducan, R. (2004). Simulated carbon uptake for a canopy of two broadleaved tree species with contrasting leaf habit. Functional Ecology [Funct. Ecol.]. 18, 3442.
https://doi.org/10.1111/j.13652435.2004.00809.x    Garcia, J. V. (1994). Principios F'ısicos de la Climatolog'ıa.: Ediciones UNALM.    Guillermo, A., Baiggrio, A., Esequiel, B., Villegas, C., Irene, T., Jose, F., et al. (2004). Atmospheric transmissivity: distribution and empirical estimation around the central andes. International Journal of Climatology, 24, 11241136.
https://doi.org/10.1002/joc.1060    Hall, F. G., Huemmerich, K. F., & Goward, S. N. (1990). Use of narrowband spectra to estimate the fraction of absorbed photosynthetically radiation. Remote Sensing and Environment, 15, 271282.
https://doi.org/10.1016/00344257(90)900976    Hamby, D. M. (1994). A review of methodology for parameter sensitivity analysis in environmental modelling. Environmental monitoring and assessment., 32, 135154.
https://doi.org/10.1007/BF00547132       Helton, J. C., Garner, J. W., Marietta, M. G., Rechard, R. E., Rudeen, D. K., & Swift, E. (1993). Uncertainty and Sensitivity Analysis Results Obtained in a Preliminary Performance assessmentfor the Waste Isolation Pilot Plant. Nuc. Sci. and Eng, 114, 286331.
https://doi.org/10.13182/NSE93A24041    Hoffman, E., & Gardner, R. H. (1983). 'Evaluation of Uncertainties in Environmental RadiologicalAssessment Models'. In H. R. Meyer (Ed.), Radiological Assessments: a Textbook on Environmental Dose Assessment. Washington, DC: Nuclear Regulatory Commission Report No. NUREG/CR3332.       Leuning, R., Helen, A., Cleugh, S. J., & Zegelin, D. H. (2005). Carbon and water fluxes over a temperate Eucalyptus forest and a tropical wet/dry savanna in Australia: measurements and comparison with MODIS remote sensing estimates. Agricultural and Forest Meteorology [Agric. For. Meteorol.]. 129, 151173.
https://doi.org/10.1016/j.agrformet.2004.12.004    Manzul, K. H., & Ysuoka, Y. (2005). Estimation of net primary productivity by integrating remote sensing data with an ecosystem model. Remote Sensing of Environment [Remote Sens. Environ.]. 94, 298310.
https://doi.org/10.1016/j.rse.2004.10.004    Mu, Q., Heinsch, F. A., Running, S. W., Cleugh, H. A., & Leuning, R. (2006). Evaluation of Satellitebased Evapotranspiration Using Fluxnet Data. 2000 Florida Ave., N.W. Washington DC: American Geophysical Union.    Mulligan, M., & Wainwright, J. (2004). Modelling and Modelling Building. In J. Wainwright & M. Mulligan (Eds.), Environmental Modelling: Finding Simplicity in Complexity: Wiley.    Myneni, P., Hoffman, B., Knyazikhin, A., Privette J L, J, G., Tian Y, et al. (2002). Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment [Remote Sens. Environ.]. 83, 214231.
https://doi.org/10.1016/S00344257(02)000743    Nicholas, C., Coops, A., Rachhpal, S., Jassal, B., Leuning R C, Andy, T., et al. (2007). Incorporation of a soil water modifier into MODIS predictions of temperate Douglasfir gross primary productivity: Initial model development. Agricultural and Forest Meteorology [Agric. For. Meteorol.]. 147, 99109.
https://doi.org/10.1016/j.agrformet.2007.07.001    Polglase, P. J., & Wang, Y. P. (1992). Potential CO2enhanced Carbon Storage by the Terrestrial Biosphere. Aust. J. Bot., 40, 641656.
https://doi.org/10.1071/BT9920641    Post, W. M., King, A. W., & Wullshleger, S. D. (1997). Historical variations in terrestrial biospheric carbon storage. Global Biogeochem. Cycles, 11(1), 99109.
https://doi.org/10.1029/96GB03942    Rosati, A., Metcalf, S. D., & Lampine, B. D. (2004). A Simple Method to Estimate Photosynthetic Radiation Use Efficiency of Canopies. Annals of Botany, 93, 567574.
https://doi.org/10.1093/aob/mch081    Smith, B., Knorr, W., Widlowski, J. L., Pinty, B., & Gobron, N. (2008). Combining remote sensing data with process modelling to monitor boreal conifer forest carbon balances. Forest Ecology and Management [For. Ecol. Manage.]. 255(12), 39853994.
https://doi.org/10.1016/j.foreco.2008.03.056    Turner, D. P., Urbanski, S., Bremer, D., Wofsy, S. C., Meyers, T., Gower, S. T., et al. (2003). A crossbiome comparison of daily light use efficiency for gross primary production. Global Change Biology, 9, 383395.
https://doi.org/10.1046/j.13652486.2003.00573.x     
