Estimating Biomass and Carbon Sequestration Capacity of Phragmites australis Using Remote Sensing and Growth Dynamics Modeling: A Case Study in Beijing Hanshiqiao Wetland Nature Reserve, China

Estimating the biomass of Phragmites australis (Cav.) Trin. ex Steud., i.e., a common wetland macrophyte, and the associated carbon sequestration capacity has attracted increasing attention. Hanshiqiao Wetland Nature Reserve (HWNR) is a large P. australis wetland in Beijing, China, and provides an ideal case study site for such purpose in an urban setting. In this study, an existing P. australis growth dynamics model was adapted to estimate the plant biomass, which was in turn converted to the associated carbon sequestration capacity in the HWNR throughout a typical year. To account for local differences, the modeling parameters were calibrated against the above-ground biomass (AGB) of P. australis retrieved from hyperspectral images of the study site. We also analyzed the sensitivity of the modeling parameters and the influence of environmental factors, particularly the nutrient availability, on the growth dynamics and carbon sequestration capacity of P. australis. Our results show that the maximum AGB and below-ground biomass (BGB) of P. australis in the HWNR are 2.93 × 103 and 2.49 × 103 g m−2, respectively, which are higher than the reported level from nearby sites with similar latitudes, presumably due to the relatively high nutrient availability and more suitable inundation conditions in the HWNR. The annual carbon sequestration capacity of P. australis in the HWNR was estimated to be 2040.73 gC m−2 yr−1, which was also found to be highly dependent on nutrient availability, with a 50% increase (decrease) in the constant of the nutrient availability KNP, resulting in a 12% increase (23% decrease) in the annual carbon sequestration capacity. This implies that a comprehensive management of urban wetlands that often encounter eutrophication problems to synergize the effects of nutrient control and carbon sequestration is worth considering in future practices.


Introduction
Phragmites australis is a typical emergent wetland plant and a globally distributed species [1]. P. australis plays a crucial role in maintaining the ecological function of wetlands and providing a variety of ecosystem services, such as water purification, habitat provision, biodiversity conservation and carbon sequestration [2][3][4][5][6]. Among these critical functions, its potential use as a major carbon sink for the city can effectively support the carbon management of the city. To this end, we adapted the growth dynamics model of P. australis originally developed by Asaeda and Karunaratne [14] to estimate the biomass (including the AGB and below-ground biomass (BGB)) and carbon sequestration capacity of P. australis in the core area of the HWNR in this study. Hyperspectral remote sensing images, together with an empirical retrieval equation proposed by Li et al. [33] for the HWNR, were used to obtain the plant biomass data for model calibration. The temporal variation in biomass and carbon sequestration capacity of P. australis throughout its growth cycle in the HWNR were simulated and analyzed, and the effects of the key environmental factors on the phenological points and maximum AGB were further discussed. The specific objectives of this study were as follows: (1) to establish a growth dynamics model that is calibrated against the hyperspectral image data to estimate the AGB and BGB of P. australis in the HWNR; (2) to illustrate the variation trend throughout the growth cycle of P. australis in the HWNR; and (3) to analyze the influence of environmental factors on the biomass and carbon sequestration capacity of P. australis in the HWNR. To accomplish these objectives, data and methods used in the present study, including the introduction of the study area are documented in Section 2. The detailed modelling results are presented in Section 3, and then Section 4 discusses the effects of the key environmental factors on plant growth and carbon sequestration capacity as well as the implications for wetland management. Finally, the main conclusions are drawn in Section 5.

Study Site
The Hanshiqiao Wetland Nature Reserve (HWNR) (40 • 06 N, 116 • 48 E) is located in Shunyi District, Beijing, China. The Reserve, which was established by the Beijing Municipal Government in April 2005, is a typical urban wetland in the North China Plain. The HWNR has a total area of 1900 ha and comprises a core area of 163.5 ha, which was chosen as our study area ( Figure 1). growth trend and carbon sequestration capacity of P. australis in the HWNR to help assess its potential use as a major carbon sink for the city can effectively support the carbon management of the city. To this end, we adapted the growth dynamics model of P. australis originally developed by Asaeda and Karunaratne [14] to estimate the biomass (including the AGB and below-ground biomass (BGB)) and carbon sequestration capacity of P. australis in the core area of the HWNR in this study. Hyperspectral remote sensing images, together with an empirical retrieval equation proposed by Li et al. [33] for the HWNR, were used to obtain the plant biomass data for model calibration. The temporal variation in biomass and carbon sequestration capacity of P. australis throughout its growth cycle in the HWNR were simulated and analyzed, and the effects of the key environmental factors on the phenological points and maximum AGB were further discussed. The specific objectives of this study were as follows: (1) to establish a growth dynamics model that is calibrated against the hyperspectral image data to estimate the AGB and BGB of P. australis in the HWNR; (2) to illustrate the variation trend throughout the growth cycle of P. australis in the HWNR; and (3) to analyze the influence of environmental factors on the biomass and carbon sequestration capacity of P. australis in the HWNR. To accomplish these objectives, data and methods used in the present study, including the introduction of the study area are documented in Section 2. The detailed modelling results are presented in Section 3, and then Section 4 discusses the effects of the key environmental factors on plant growth and carbon sequestration capacity as well as the implications for wetland management. Finally, the main conclusions are drawn in Section 5.

Study Site
The Hanshiqiao Wetland Nature Reserve (HWNR) (40°06′ N, 116°48′ E) is located in Shunyi District, Beijing, China. The Reserve, which was established by the Beijing Municipal Government in April 2005, is a typical urban wetland in the North China Plain. The HWNR has a total area of 1900 ha and comprises a core area of 163.5 ha, which was chosen as our study area ( Figure 1).  The Reserve is located in the Chaobai River catchment, and the Caijia River flows through it and becomes its main natural water source. This region has a semi-humid continental monsoon climate, with an annual average temperature of 11.8 • C and annual precipitation of 644 mm, mainly in July and August [28]. There are 182 species of birds and 222 species of plants belonging to 144 genera and 55 families recorded in the Reserve at present. Macrophytes such as P. australis and Typha orientalis are widely distributed, and P. australis is the dominant species [28].

Imagery Dataset
Hyperspectral images were acquired for P. australis interpretation in the study area, retrieving its AGB. P. australis usually starts growing in March and reaches its maximum biomass in July [34] and undergoes a controlled harvest by the administration agency of the HWNR in late November. Thus, images of 1 November 2020, 17 June 2021, and 31 July 2021 were obtained from the Orbita hyperspectral satellite constellation (OHS) (https://www.obtdata.com/, accessed on 12 October 2021) to reflect the conditions of P. australis at the peak and the end of the growth cycle.
Developed by Zhuhai Orbita Aerospace Science & Technology Co. Ltd. (Zhuhai, China), the OHS is the first hyperspectral satellite constellation in the world with surface coating technology for sensors [35]. The constellation consists of eight hyperspectral satellites, and the first batch of four satellites was successfully launched on 26 April 2018, followed by the second batch of four satellites launched on 19 September 2019. Satellites of OHS have an on-orbit life span of more than 5 years, and the revisit cycle is 5 days. A dualfrequency global navigation satellite system (GNSS) receiver (GPS and Beidou) is used to measure the position of OHS in its orbit. The sensor altitude of OHS is around 500 km above the ground, resulting in images of a spatial resolution of 10 m. Each satellite of the OHS has 256 bands with a wavelength ranging from 400 nm to 1000 nm at approximately 2.5 nm intervals. However, only 32 bands can be selected due to limitations in data transmission and storage. The main parameters of OHS are listed in Table 1 [36]. High-resolution images (pixel size < 10 m) such as SPOT images with a pixel size of 2.5 m and QuickBird images with a pixel size of 2.4 m have been widely used to form a dataset of training samples for P. australis interpretation, as they contain abundant spatial information [37,38]. In this study, a high-resolution BJ-2 image of June 2021 was used for selecting training samples for interpretation. The BJ-2 images are provided by the BJ-2 satellite constellation, which was developed by Surrey Satellite Technology Ltd. (https://www.sstl.co.uk/, accessed on 2 August 2021) and launched in July 2015. The constellation has a revisit cycle of one day and provides sub-meter resolution images of red, green, and blue bands at a pixel size of 0.8 m and infrared images at a pixel size of 3.2 m. With a relatively high spatial resolution, BJ-2 images have been utilized for various applications such as monitoring landslides, extracting hydraulic engineering constructions, classifying land-cover types, etc. [39][40][41]. In this study, the 0.8 m resolution BJ-2 image that contained red, green, and blue bands was acquired, geo-rectified to the WGS 84 datum (50N), and further used as a reference image for the registration of the OHS images.

Model Formulation
The present model was adapted from the P. australis growth model originally developed by Asaeda and Karunaratne [14]. Similarly, the biomass of P. australis (g DM m −2 ) was divided into above-ground and below-ground organs. The above-ground organs were further subdivided into flowering shoots including fertile leaves and inflorescences, as well as non-flowering shoots including (sterile) leaves and non-flowering secondary shoots. The below-ground organs were subdivided into roots and new and old rhizomes. The changes in biomass of each organ throughout the growth cycle were estimated by incorporating the net growth in each organ as a function of photosynthesis, respiration, mortality, and inter-organ translocation. Briefly, the biomass of each organ was calculated as follows: where B, P, R, M, and T r denote biomass, gross photosynthesis production, respiration loss, mortality biomass, and inter-organ translocation, respectively. The specific framework of the model is shown in Figure 2. The model stratified the above-ground part into 1 cm thick horizontal layers, where the net growth budget and shoot elongation were calculated separately. Based on the above organ division, we selected five state variables-i.e., biomass of shoots, inflorescences, roots, old rhizomes and new rhizomes-to characterize plant growth. The details of the growth equations for each state variable are documented in the Supplementary Materials. The growth equations for each layer and organ were solved simultaneously using the fourth-order Runge-Kutta method, and the time step of the iteration was set to one day.
classifying land-cover types, etc. [39][40][41]. In this study, the 0.8 m resolution BJ-2 image that contained red, green, and blue bands was acquired, geo-rectified to the WGS 84 datum (50N), and further used as a reference image for the registration of the OHS images.

Model Formulation
The present model was adapted from the P. australis growth model originally developed by Asaeda and Karunaratne [14]. Similarly, the biomass of P. australis (g DM m −2 ) was divided into above-ground and below-ground organs. The above-ground organs were further subdivided into flowering shoots including fertile leaves and inflorescences, as well as non-flowering shoots including (sterile) leaves and non-flowering secondary shoots. The below-ground organs were subdivided into roots and new and old rhizomes. The changes in biomass of each organ throughout the growth cycle were estimated by incorporating the net growth in each organ as a function of photosynthesis, respiration, mortality, and inter-organ translocation. Briefly, the biomass of each organ was calculated as follows: where B, P, R, M, and Tr denote biomass, gross photosynthesis production, respiration loss, mortality biomass, and inter-organ translocation, respectively. The specific framework of the model is shown in Figure 2. The model stratified the above-ground part into 1 cm thick horizontal layers, where the net growth budget and shoot elongation were calculated separately. Based on the above organ division, we selected five state variables-i.e., biomass of shoots, inflorescences, roots, old rhizomes and new rhizomes-to characterize plant growth. The details of the growth equations for each state variable are documented in the Supplementary Materials. The growth equations for each layer and organ were solved simultaneously using the fourth-order Runge-Kutta method, and the time step of the iteration was set to one day. The growth dynamics model reflects the biomass change during a typical growth cycle of P. australis and neglects interannual variation. The major phenological points and The growth dynamics model reflects the biomass change during a typical growth cycle of P. australis and neglects interannual variation. The major phenological points and the corresponding phenological events in the present model are summarized in Table 2. As the phenological points are largely correlated with latitude and further modulated by local conditions such as meteorology and terrain, the empirical regression equations that describe the relationship among the phenological points of P. australis derived by Asaeda and Karunaratne [14] through reported case studies from multiple countries (see Supplementary Materials) were used as the basis for further calibration, as shown in the following. All dates used to describe phenological points were Julian days. The formation of new roots and fertile leaves, as well as the beginning of translocation of dry matter from old rhizomes to fertile leaves and roots The formation of sterile leaves and non-flowering secondary shoots, as well as the beginning of translocation of dry matter from old rhizomes to sterile leaves and non-flowering secondary shoots The formation of new rhizomes and peduncles, together with the beginning of translocation of dry matter from old rhizomes to peduncles, as well as the beginning of translocation of photosynthesized material to below-ground plant organs The appearance of panicles as well as the beginning of translocation of dry matter and photosynthesized material from peduncles to panicles t e The ending of mobilization of dry matter from rhizomes to shoots and roots t s The commencement of shoot senescence, together with the beginning of translocation of accumulated shoot dry matter to below-ground organs

Model Calibration
Following Asaeda and Karunaratne [14], the key parameters (see Supplementary Materials) in the P. australis growth dynamics model include parameters that are related to the generic ecophysiology of P. australis (e.g., half saturation constant of age for shoot photosynthesis) and parameters that are more dependent on local conditions and are thus highly site-specific (e.g., the phenological points). Generic parameters adopted values from the original model and needed no calibration. Site-specific parameters needed to be further calibrated using field data.
In general, observation of the occurrence of most phenological events is difficult or sometimes impossible. Among those phenological events, the start of plant growth is relatively easier to detect [42]. Once the date when growth began as denoted by t b is specified, the date when panicles start forming (t f ), the date when the upward translocation from old rhizomes terminates (t e ) and the date that marks the onset of senescence (t s ) can be preliminarily derived using the respective empirical relationships of the various phenological points as functions of t b . Further calibration and the resulting phenological points of P. australis in the HWNR are shown in Section 3.2.
To determine the order of calibration of the remaining site-specific parameters, their influence on the model predictions of the AGB and BGB were assessed by carrying out a sensitivity analysis, in which a shifting range of ±50% for each parameter was implemented one at a time. The most sensitive modeling parameters thus identified were then calibrated one at a time through comparison against remote sensing retrieved biomass data by trial and error. The values of the parameters associated with the least sums of squared residuals between model predictions and retrieved data were adopted as the calibrated values, as presented in Section 3.

Model Execution
To carry out the simulation, input data such as the latitude were specified. Regarding the initial values, we assumed the initial leaf area index (LAI) to be zero as a normal practice in previous studies [17,33,43,44]. Since P. australis in the study area are harvested at the end of each year, the initial AGB was assumed to be zero as well. A series of simulations assuming different initial BGB values were conducted, and the value associated with the least mean square error between the predicted AGB in a year and the retrieved AGB was adopted. The model was implemented in MATLAB R2020b (MathWorks, Natick, MA, USA).

Preprocessing of Hyperspectral Images
Radiometric calibration and atmospheric correction were carried out on the OHS images to obtain surface reflectance required by subsequent analyses. Radiometric calibration was performed to convert the original digital number (DN) values into apparent radiance using the following formula [45]: where Rad is the apparent radiance, gain refers to the gain coefficient of the absolute radiometric calibration, TDI Stage is a parameter that describes how active pixels are placed in a time-delay-integration (TDI) image that is provided by the instrument supplier, and offset is the absolute radiometric calibration offset coefficient. Atmospheric correction was performed using the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) module [46] to obtain surface reflectance with Rad. Images with surface reflectance were registered to the corresponding BJ-2 image with the registration root mean square error (RMSE) less than 0.5 pixel. The registered OHS images were resampled to their original pixel size (10 m) by the nearest-neighbor resampling technique [47]. Radiometric calibration, atmospheric correction, and registration were performed using ENVI 5.3 (L3Harris Technologies, Broomfield, CO, USA), and resampling of the registered OHS images was carried out using ArcGIS 10.3 (ESRI, Redlands, CA, USA).

Interpretation of P. australis
In this study, as the plant growth model is a generic zero-dimensional plot-scale model, the potential inter-annual variation in the spatial distribution of P. australis in the study area was neglected and P. australis was interpretated using the OHS hyperspectral image of June 2021, as the plant was matured at that time. A total of 344 polygon samples of water body and vegetation were randomly selected from the high-resolution BJ-2 image. As Nymphaea tetragona is the dominant floating plant and T. orientalis is the most common coexistent species of P. australis in the HWNR [48,49], vegetation samples were further categorized into P. australis, N. tetragona, T. orientalis, and other vegetation ( Figure S1), resulting in five classes in total. The polygon samples were divided into two groups. One group consisting of 237 training samples was used for interpretation, and the other group consisting of 107 validation samples was used for evaluating the accuracy of interpretation (Table 3). Separability metrics including Jeffries-Matusita distance (JM) and the transformed divergence (TD) were adopted to evaluate the efficiency of the selected training samples [50]. The values of these two metrics generally range from 0 to 2, and the corresponding separability is considered to be good if the metric value is greater than 1.9. Support vector machine (SVM) was used to extract P. australis in the study area. SVM is an advanced non-parametric machine learning classifier that has been proved to be robust and reliable for P. australis interpretation using hyperspectral images [51]. This method follows the idea of inputting samples into a high-dimensional space and constructing a hyperplane to separate samples into different classes, and a separation margin defined by the hyperplane indicates the difference between two classes of samples, with the optimal hyperplane that can maximize the margin selected to generate the ultimate separation result [52]. With the advantage that fewer training samples are generally required, SVM also demonstrates better adaptability and faster training speed than other machine learning methods [53].
In the present study, the radial basis function (RBF) kernel was selected for SVM, as it has been proven to be more accurate than other kernels [54]. Following relevant previous studies [55], the two key parameters of SVM-i.e., the value of gamma in RBF kernel and penalty parameter-were set as 0.01 and 100, respectively. Overall accuracy (OA) and Kappa coefficients were calculated to assess the accuracy of interpretation results. Producer's accuracy (PA) and user's accuracy (UA) were also calculated to estimate the degree of misclassification for each class, as PA and UA represent the percentage of correct interpretation results in the validation samples and the training samples pixel by pixel [56]. Based on the feature information revealed by the high-resolution BJ-2 image, artificial visual interpretation was also performed on the results of SVM to correct the misclassified pixels and obtain the final interpretation result. All the steps described in this section were carried out in ENVI 5.3 (L3Harris Technologies, Broomfield, CO, USA).

Retrieval of P. australis AGB
The AGB of P. australis was retrieved using an empirical model constructed by Li et al. [33]. The model was developed based on the observed AGB data and a canopy spectrum of P. australis in the HWNR. The coefficient of determination (R 2 ) of the model was 0.72, which indicates that the retrieval model was reasonably reliable. The specific equation is as follows: where x represents the maximum surface reflectance in the wavelength range of 510-560 nm. Selection of the maximum surface reflectance and the calculation were performed in ArcGIS 10.3 (ESRI, Redlands, CA, USA).

Evaluation of Carbon Sequestration Capacity
The carbon sequestration capacity of P. australis in the study area was calculated based on the photosynthesis equation [57]: The ratio between the molecular weight of the organic material (starch) and CO 2 sequestered by photosynthesis was 1:1.62, which means per 1 g of organic material fixated approximately 1.62 g CO 2 (0.44 g C). The net primary production (NPP) was calculated to represent the net accumulation of organic material generated via photosynthesis. NPP is defined as the amount of total organic material minus the loss of material caused by plant respiration [58]. As shown in Section 2.3.1, the net growth (i.e., living biomass) in the model was defined as the production of photosynthesis minus the loss due to respiration and mortality. Hence, in the present study, daily NPP (g m −2 ) was estimated as follows: where ∆B t represents the variation in simulated living biomass between two adjacent days (i.e., B t and B t−1 ), and M t refers to the simulated mortality biomass of day t. Thereafter, daily carbon sequestration capacity (gC m −2 ) was estimated using the calculated NPP and the above ratio. Monthly and annual carbon sequestration capacities were further obtained by adding up daily carbon sequestration capacities throughout each month and the entire year, respectively. The total annual carbon sequestration capacity of P. australis in the HWNR was calculated by multiplying the annual carbon sequestration capacity per square meter with the distribution area of P. australis in the HWNR.

Distribution and Retrieved AGB of P. australis
As shown in Table 4, the training samples used in the present study were reliable, as the JM and TD values between P. australis and other classes were consistently higher than 1.9. Table 5 shows the accuracy assessment results of the interpretation. With an overall accuracy of 92.21% and a Kappa coefficient of 0.88, the accuracy of the interpretation result was satisfactory. According to the producer's accuracy and user's accuracy, the error of interpretation mainly arose from the difficulty in identifying N. tetragona and T. orientalis. Producer's accuracy showed that pixels of N. tetragona and T. orientalis were mainly misclassified as water body (165/895) and P. australis (194/694), respectively. Fairly good results had been achieved for P. australis interpretation, with the producer's accuracy and user's accuracy attaining 94.71% and 96.95%, respectively. The misclassified pixels were manually corrected, resulting in the distribution of P. australis shown in Figure 3, with a total distribution area of 86.58 ha.  Figure 3 shows the retrieved results and the distribution pattern of the AGB of P. australis in the study area. Among the three representative months of June, July, and November, P. australis in July had the highest averaged AGB (2863.59 ± 131.86 g m −2 ), followed by June (1156.72 ± 247.09 g m −2 ) and November (498.38 ± 145.04 g m −2 ). A statistical analysis on the pixel numbers was further performed, and we defined the pixels ranged between lower and upper quartiles as the concentrated range of the AGB. Results show that AGB for June, July, and November mainly concentrated between 991.33 and 1319.43 g m −2 , 2774.30 and 2951.85 g m −2 , and 405.57 and 597.21 g m −2 , respectively. The above results show a consistent trend that P. australis kept growing until it experienced mortality in November. Areas with high AGB values were mainly concentrated in the southern part of the study area throughout the three months. On the contrary, low AGB values were prevalent in the northern part.
interpretation mainly arose from the difficulty in identifying N. tetragona and T. orientalis. Producer's accuracy showed that pixels of N. tetragona and T. orientalis were mainly misclassified as water body (165/895) and P. australis (194/694), respectively. Fairly good results had been achieved for P. australis interpretation, with the producer's accuracy and user's accuracy attaining 94.71% and 96.95%, respectively. The misclassified pixels were manually corrected, resulting in the distribution of P. australis shown in Figure 3, with a total distribution area of 86.58 ha.

Calibration and Sensitivity Analysis of Growth Dynamics Model
According to the retrieved AGB, upward translocation of biomass to the above-ground organs was observed at the end of April, then the AGB reached the peak value around August and finally declined to the minimum around November. To reflect the typical growth pattern and make the model more suitable for simulating the growth cycle of P. australis in the HWNR, the following calibration of the key parameters was performed: (1) The dates of the phenological points were calibrated based on the preliminary value obtained by the empirical relationships in Section 2.3.2. Most phenological points were moved backward to match the actual growth cycle of P. australis in the HWNR. The specific time of some representative phenological events during the growth cycle of P. australis is shown in Figure 4.
australis in the study area. Among the three representative months of June, July, and November, P. australis in July had the highest averaged AGB (2863.59 ± 131.86 g m −2 ), followed by June (1156.72 ± 247.09 g m −2 ) and November (498.38 ± 145.04 g m −2 ). A statistical analysis on the pixel numbers was further performed, and we defined the pixels ranged between lower and upper quartiles as the concentrated range of the AGB. Results show that AGB for June, July, and November mainly concentrated between 991.33 and 1319.43 g m −2 , 2774.30 and 2951.85 g m −2 , and 405.57 and 597.21 g m −2 , respectively. The above results show a consistent trend that P. australis kept growing until it experienced mortality in November. Areas with high AGB values were mainly concentrated in the southern part of the study area throughout the three months. On the contrary, low AGB values were prevalent in the northern part.

Calibration and Sensitivity Analysis of Growth Dynamics Model
According to the retrieved AGB, upward translocation of biomass to the aboveground organs was observed at the end of April, then the AGB reached the peak value around August and finally declined to the minimum around November. To reflect the typical growth pattern and make the model more suitable for simulating the growth cycle of P. australis in the HWNR, the following calibration of the key parameters was performed: (1) The dates of the phenological points were calibrated based on the preliminary value obtained by the empirical relationships in Section 2.3.2. Most phenological points were moved backward to match the actual growth cycle of P. australis in the HWNR. The specific time of some representative phenological events during the growth cycle of P. australis is shown in Figure 4.  (2) The constant of the nutrient availability and the biomass transfer rate from fertile leaf, sterile leaf, and peduncle to new rhizomes were slightly increased, and the respiration rate of above-ground organs was reduced in our model, to achieve higher peak biomass, as shown by the retrieved results. (3) The mortality rates of the above-ground organs at the onset of senescence, especially those of shoots and flowers, were increased after t s to match the occurrence of abrupt biomass loss at that stage.
The seasonal variation in the predicted AGB of P. australis was compared with the average value of the retrieved data across the study area (Figure 5a). For all retrieved data points, the model predictions agreed satisfactorily (R 2 = 0.9678) and captured the general growth pattern. peak biomass, etc., were all successfully reproduced by our model prediction, which was consistent with the growth pattern of P. australis reported in the literature [14]. Figure 5b shows the seasonal variation in the BGB of P. australis, which was in line with the typical growth pattern of the below-ground organs-e.g., the reduction in the BGB during the early growing season due to the supply of stored material to shoot growth-and the increase in the BGB after tp (new rhizomes start forming) because of the downward translocation of material from photosynthesis and shoot dry matter, etc. [14]. Additionally, the BGB at the end of the growth season was slightly higher than that at the beginning, which can be attributed to the need of nutrient storage and preparation for growth in the following year [17]. Through sensitivity analysis of the key modeling parameters, the growth dynamics model was found to be most sensitive to maximum specific net daily photosynthesis rate, specific mortality rate of shoots, fraction of shoot transfer to rhizome, constant of the nutrient availability, and fraction of shoot biomass for elongation. Table 6 shows the sensitivity of the model output on key modelling parameters for the study site. For instance, a 50% increase in the maximum specific net daily photosynthesis rate resulted in a 56.35% increase in the AGB at ts (the date when shoots start senescence), and a 50% decrease resulted in a 58.01% decrease in the AGB at ts. By contrast, a 50% increase or decrease in the other four parameters led to much lesser effects on the AGB at ts. In addition, the outcome of the model was also affected marginally by some other parameters such as respiration rates, mortality rates, and fractions of transfer, etc. For instance, a mildly increasing biomass growth rate in the early-growing season, followed by a drastically increasing biomass growth rate during the mid-growing season, the decline in biomass due to senescence at the end of the growing season, the time of peak biomass, etc., were all successfully reproduced by our model prediction, which was consistent with the growth pattern of P. australis reported in the literature [14]. Figure 5b shows the seasonal variation in the BGB of P. australis, which was in line with the typical growth pattern of the below-ground organs-e.g., the reduction in the BGB during the early growing season due to the supply of stored material to shoot growth-and the increase in the BGB after t p (new rhizomes start forming) because of the downward translocation of material from photosynthesis and shoot dry matter, etc. [14]. Additionally, the BGB at the end of the growth season was slightly higher than that at the beginning, which can be attributed to the need of nutrient storage and preparation for growth in the following year [17].
Through sensitivity analysis of the key modeling parameters, the growth dynamics model was found to be most sensitive to maximum specific net daily photosynthesis rate, specific mortality rate of shoots, fraction of shoot transfer to rhizome, constant of the nutrient availability, and fraction of shoot biomass for elongation. Table 6 shows the sensitivity of the model output on key modelling parameters for the study site. For instance, a 50% increase in the maximum specific net daily photosynthesis rate resulted in a 56.35% increase in the AGB at t s (the date when shoots start senescence), and a 50% decrease resulted in a 58.01% decrease in the AGB at t s . By contrast, a 50% increase or decrease in the other four parameters led to much lesser effects on the AGB at t s . In addition, the outcome of the model was also affected marginally by some other parameters such as respiration rates, mortality rates, and fractions of transfer, etc.

Estimation of Carbon Sequestration Capacity
As shown in Figure 5, living AGB and BGB of P. australis ranged from 1.87 to 2.93 × 10 3 g m −2 and 1.75 × 10 3 to 2.49 × 10 3 g m −2 , respectively. In addition, mortality biomass ranged from 0 to 203.95 g m −2 and 0.07 to 0.62 g m −2 for the above-and below-ground part, respectively. Therefore, daily carbon sequestration capacity for the above-and below-ground part varied in the range of −0.05 to 19.71 gC m −2 and −2.58 to 6.67 gC m −2 , respectively. Maximum daily carbon sequestration capacity was attained on the 207th day. Annual carbon sequestration capacity was 2040.73 gC m −2 yr −1 , which translated to a total amount of 1766.86 tC being sequestrated annually by P. australis in the HWNR.
As shown in Figure 6, monthly whole-plant carbon sequestration capacity was found to slowly increase from January to April, which was consistent with the phenological pattern that the above-ground part of P. australis started growing at the end of April (120th day). Growth of sterile shoot in May led to a rapid increase in the AGB and monthly whole-plant carbon sequestration capacity, i.e., 509% higher than that in April. The monthly whole-plant carbon sequestration capacity reached its peak in July with the value of 528 gC m −2 , which was mainly contributed by the largest monthly carbon sequestration capacity of the above-ground part in July of 577.19 gC m −2 .

Estimation of Carbon Sequestration Capacity
As shown in Figure 5, living AGB and BGB of P. australis ranged from 1.87 to 2.93 × 10 3 g m −2 and 1.75 × 10 3 to 2.49 × 10 3 g m −2 , respectively. In addition, mortality biomass ranged from 0 to 203.95 g m −2 and 0.07 to 0.62 g m −2 for the above-and below-ground part, respectively. Therefore, daily carbon sequestration capacity for the above-and belowground part varied in the range of −0.05 to 19.71 gC m −2 and −2.58 to 6.67 gC m −2 , respectively. Maximum daily carbon sequestration capacity was attained on the 207th day. Annual carbon sequestration capacity was 2040.73 gC m −2 yr −1 , which translated to a total amount of 1766.86 tC being sequestrated annually by P. australis in the HWNR.
As shown in Figure 6, monthly whole-plant carbon sequestration capacity was found to slowly increase from January to April, which was consistent with the phenological pattern that the above-ground part of P. australis started growing at the end of April (120th day). Growth of sterile shoot in May led to a rapid increase in the AGB and monthly whole-plant carbon sequestration capacity, i.e., 509% higher than that in April. The monthly whole-plant carbon sequestration capacity reached its peak in July with the value of 528 gC m −2 , which was mainly contributed by the largest monthly carbon sequestration capacity of the above-ground part in July of 577.19 gC m −2 . As new peduncles and roots started to form from the end of July (205th day), belowground carbon sequestration capacity increased afterwards. Daily below-ground carbon As new peduncles and roots started to form from the end of July (205th day), belowground carbon sequestration capacity increased afterwards. Daily below-ground carbon sequestration capacity became positive from the 207th day, and monthly carbon seques-tration capacity of the below-ground part became higher than that of the above-ground part in September. As both the above-and below-ground part continued senesce from the middle of September (260th day), monthly whole-plant carbon sequestration capacity kept decreasing and became negative in November and December.

Comparison in Phenological Points and Biomass of P. australis from Other Sites
The biomass and phenological points of P. australis from four additional sites were collected from the literature to compare with the results in the HWNR (Table 7). Compared with P. australis from higher-latitude area (e.g., Nesyt fish pond, Czech Republic, and Loch of Balgavies, Scotland), the maximum AGB of P. australis in the HWNR was larger, and its start dates of shoot growth and panicle formation were relatively later, which was consistent with the general pattern of latitudinal variation reported by Karunaratne et al. [59] and Clevering et al. [60] (Figure 7). However, significant variability can also be found between biomass of P. australis in the HWNR and those in similar latitudes, such as Baiyangdian Lake (Hebei, China) [61] and Yeya Lake wetland (Beijing, China) [62], which indicated that biomass is also affected by other confounding local factors, such as nutrient level [63,64], inundation depth, and hydroperiod (i.e., inundation frequency and duration) [65,66]. This variability can be found in other latitudes as well, e.g., P. australis in Loch of Balgavies, Scotland, and Vejlerne Nature Reserve, Denmark-two sites with similar latitudes that also have significant differences in biomass [67,68]. in the latter (1.59 mg/L and 0.067 mg/L) [73], which can also lead to the higher maximum ABG of P. australis in the HWNR.

Influence of Nutrient Availability on Growth and Carbon Sequestration of P. australis
Our results show that the constant of the nutrient availability (KNP) was one of the key environmental factors of the growth dynamics model (Table 6), which can further influence the carbon sequestration capacity. As eutrophication poses potential threat in the HWNR [74], we further explored how it influences the growth and carbon sequestra-  Figure 7 shows that the maximum AGB of P. australis in the HWNR was considerably higher than that in Baiyangdian Lake (Hebei, China) [61,65] and Yeya Lake wetland (Beijing, China) [62], which can be attributed to the influence of nutrients and inundation depth. Wang et al. [70] studied the optimum inundation depth for P. australis in different growing stages and found that a relatively low inundation depth (≤0.10 m) can facilitate the germination and establishment of P. australis at the early stage (from the end of April to May) and that higher inundation depth is more suitable at the rapid growth stage (from the end of June to the end of August). Due to the annual precipitation and evaporation pattern [28], a prominent seasonal variation trend with lower inundation depth in spring and higher inundation depth in summer was exhibited in the HWNR, which was consistent with the general requirement of suitable inundation depth for P. australis [70] and contributed to the relatively high biomass of P. australis in the HWNR.
In addition, previous studies have shown that higher nutrient concentration facilitates the growth of P. australis [71,72]. Comparing the nutrient concentration between the HWNR and Yeya Lake wetland, the annual average values of total nitrogen and total phosphorus in the former (1.93 mg/L and 0.33 mg/L) were significantly higher than those in the latter (1.59 mg/L and 0.067 mg/L) [73], which can also lead to the higher maximum ABG of P. australis in the HWNR.

Influence of Nutrient Availability on Growth and Carbon Sequestration of P. australis
Our results show that the constant of the nutrient availability (K NP ) was one of the key environmental factors of the growth dynamics model (Table 6), which can further influence the carbon sequestration capacity. As eutrophication poses potential threat in the HWNR [74], we further explored how it influences the growth and carbon sequestration capacity of P. australis by setting five scenarios with different constant of the nutrient availability. Assuming the calibrated K NP in Section 3.2 as the base value, Scenario 1 was set to the base scenario, while Scenarios 2-5 were set by varying K NP at intervals (−50%, −25%, +25%, +50%). Figure 8a,b shows simulated results of the AGB and BGB in all five scenarios. in the latter (1.59 mg/L and 0.067 mg/L) [73], which can also lead to the higher maximum ABG of P. australis in the HWNR.

Influence of Nutrient Availability on Growth and Carbon Sequestration of P. australis
Our results show that the constant of the nutrient availability (KNP) was one of the key environmental factors of the growth dynamics model (Table 6), which can further influence the carbon sequestration capacity. As eutrophication poses potential threat in the HWNR [74], we further explored how it influences the growth and carbon sequestration capacity of P. australis by setting five scenarios with different constant of the nutrient availability. Assuming the calibrated KNP in Section 3.2 as the base value, Scenario 1 was set to the base scenario, while Scenarios 2-5 were set by varying KNP at intervals (−50%, −25%, +25%, +50%). Figure 8a,b shows simulated results of the AGB and BGB in all five scenarios.
(a) (b) Figure 8. Biomass variation in P. australis subject to scenarios with varying nutrient availability: (a) in the above-ground organs; (b) in the below-ground organs. ts: the date when shoots start senescence; tp: the date when peduncles and new rhizomes start forming. With higher values of K NP , both the maximum AGB and BGB in Scenarios 4-5 were significantly higher than in Scenario 1. Conversely, they were lower in Scenarios 2-3 than in Scenario 1. This pattern indicated the facilitating effect of K NP on the growth of P. australis, which was consistent with the findings in several other cases (e.g., Burullus Lake, Egypt; Swazaro, Poland) [71,72]. Subject to varying K NP , new rhizomes began to grow at different rates after t p (the date when peduncles and new rhizomes start forming) because of different amounts of current photosynthesized material transferred to the below-ground organs, leading to contrasting biomass in the below-ground part. After t s (the date when shoots start senescence), the senescence of the above-ground organs led to abrupt decline in the AGB and also diminished the contrast among the various scenarios. At the same time, the transport of photosynthetic materials from the above-ground part to the below-ground part widened the gap in the BGB. This is also consistent with the finding reported by Karunaratne et al. [59] that increased nutrient availability leads to higher living biomass as well as more photosynthetic materials transported to the below-ground part.
As a result of the biomass change, our results show that the annual carbon sequestration capacity increases by 12% when K NP increases by 50%, and decreases by 23% when K NP decreases by 50% ( Figure 9); i.e., annual carbon sequestration capacity of P. australis is enhanced by a higher K NP [71,72,75,76]. Hence, comprehensive analysis of the ecological consequences of nutrient removal and carbon sequestration capacity should be performed in future wetland management to optimize the trade-off between eutrophication control and carbon sequestration. Additionally, as shown in Section 2.5, carbon sequestration was estimated based on the contribution of both living parts (i.e., above-and below-ground part) and the mortality part in the present study. Our results show that the magnitude of carbon sequestration capacity of each part had a positive correlation with increasing K NP (Figure 9), which was consistent with reported findings from other studies [77]. However, the effect of the decomposition of mortality part was not considered in our carbon sequestration analysis. In fact, this process can trigger the carbon release and cause loss of carbon sequestration capacity of P. australis [77,78]. As the K NP increases, enhanced carbon sequestration capacity from the mortality part can cause larger loss in carbon sequestration capacity if not treated properly. As such, we recommend that more frequent harvest and regular collecting of mortality organs should be carried out [77,79] after the senescence begins so that carbon sequestration loss due to decomposition can be reduced.
With higher values of KNP, both the maximum AGB and BGB in Scenarios 4-5 were significantly higher than in Scenario 1. Conversely, they were lower in Scenarios 2-3 than in Scenario 1. This pattern indicated the facilitating effect of KNP on the growth of P. australis, which was consistent with the findings in several other cases (e.g., Burullus Lake, Egypt; Swazaro, Poland) [71,72]. Subject to varying KNP, new rhizomes began to grow at different rates after tp (the date when peduncles and new rhizomes start forming) because of different amounts of current photosynthesized material transferred to the belowground organs, leading to contrasting biomass in the below-ground part. After ts (the date when shoots start senescence), the senescence of the above-ground organs led to abrupt decline in the AGB and also diminished the contrast among the various scenarios. At the same time, the transport of photosynthetic materials from the above-ground part to the below-ground part widened the gap in the BGB. This is also consistent with the finding reported by Karunaratne et al. [59] that increased nutrient availability leads to higher living biomass as well as more photosynthetic materials transported to the below-ground part.
As a result of the biomass change, our results show that the annual carbon sequestration capacity increases by 12% when KNP increases by 50%, and decreases by 23% when KNP decreases by 50% ( Figure 9); i.e., annual carbon sequestration capacity of P. australis is enhanced by a higher KNP [71,72,75,76]. Hence, comprehensive analysis of the ecological consequences of nutrient removal and carbon sequestration capacity should be performed in future wetland management to optimize the trade-off between eutrophication control and carbon sequestration. Additionally, as shown in Section 2.5, carbon sequestration was estimated based on the contribution of both living parts (i.e., above-and below-ground part) and the mortality part in the present study. Our results show that the magnitude of carbon sequestration capacity of each part had a positive correlation with increasing KNP (Figure 9), which was consistent with reported findings from other studies [77]. However, the effect of the decomposition of mortality part was not considered in our carbon sequestration analysis. In fact, this process can trigger the carbon release and cause loss of carbon sequestration capacity of P. australis [77,78]. As the KNP increases, enhanced carbon sequestration capacity from the mortality part can cause larger loss in carbon sequestration capacity if not treated properly. As such, we recommend that more frequent harvest and regular collecting of mortality organs should be carried out [77,79] after the senescence begins so that carbon sequestration loss due to decomposition can be reduced. Figure 9. Annual carbon sequestration capacity in above-and below-ground and mortality parts of P. australis growing with different nutrient availability. Figure 9. Annual carbon sequestration capacity in above-and below-ground and mortality parts of P. australis growing with different nutrient availability.

Conclusions
In this study, an existing growth dynamics model was adapted to estimate biomass and carbon sequestration capacity of P. australis during a typical growth cycle in the Hanshiqiao Wetland Nature Reserve (HWNR) in Beijing, China. The AGB retrieved from hyperspectral images of the study area was used for model calibration. The results show that the AGB reached its peak of 2930 g/m −2 in late August and that the maximum BGB of 2486 g m −2 is observed in early November. The maximum monthly carbon sequestration capacity is attained in July, and the annual carbon sequestration capacity is 2040 gC m −2 yr −1 , which translates to 1766 tC to be sequestrated by P. australis growing in the HWNR. In addition to latitude, inundation depth and nutrient level are found to have significant effects on the growth of P. australis in the HWNR and account for its relatively high biomass compared with nearby sites with similar latitudes. Scenarios with varying nutrient availability were also simulated, and the results show that nutrient availability has a facilitating effect on the growth and carbon sequestration capacity of P. australis in the HWNR. Based on the modelling results, in future practice, we recommend the comprehensive management of urban wetlands that often encounter eutrophication problems to synergize the effects of nutrient control and carbon sequestration. In addition, it is notable that relative to living parts (above-and below-ground part), the contribution of mortality part to carbon sequestration was disproportionately enhanced by increasing nutrient availability, leading to enlarged potential loss of carbon sequestration capacity due to decomposition of mortality organs. Therefore, more frequent harvest and regular collecting of the mortality part of P. australis should be carried out after the senescence begins to reduce the carbon sequestration loss.
The zero-dimensional plot-scale plant growth dynamics model is useful for a preliminary estimation of the plant biomass and carbon sequestration capacity throughout the vegetated area in a wetland. A spatially explicit model to account for spatial heterogeneity in plant growth can be developed and applied in the future. The empirical retrieval model established by Li et al. [33] was based on field biomass data gathered in 2017, and interannual variation may lead to some uncertainty in our application. We only consider the contribution of plant photosynthesis in our carbon sequestration estimation, and other relevant processes in the wetland carbon cycle such as microbial decomposition and respiration in soil [80] can be incorporated in future ecosystem carbon sequestration estimation.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/s22093141/s1, Figure S1: Ground-truth samples of ground objects; Table S1: Parameters used in the growth dynamics model; Table S2: Total annual biomass of mortality of P. australis growing with different nutrient availability; Equations (S1)-(S17): Equations used to simulate the dynamic growth of P. australis in study site.