Predicting mangrove forest dynamics across a soil salinity gradient using an individual-based vegetation model linked with plant hydraulics

In mangrove forests, soil salinity is one of the most significant environmental factors determining forest distribution and productivity as it limits plant water uptake and carbon gain. However, salinity control on mangrove productivity through plant hydraulics has not been investigated by existing mangrove models. Here we present a new individual-based model linked with plant hydraulics to incorporate physiological characterization of mangrove growth under salt stress. Plant hydraulics was associated with mangroves' nutrient uptake and biomass allocation apart from water flux and carbon gain. The developed model was performed for two coexisting species – Rhizophora stylosa and Bruguiera gymnorrhiza – in a subtropical mangrove forest in Japan. The model predicted that the productivity of both species was affected by soil salinity through downregulation of stomatal conductance. Under low-soil-salinity conditions (< 28 ‰), B. gymnorrhiza trees grew faster and suppressed the growth of R. stylosa trees by shading that resulted in a B. gymnorrhiza-dominated forest. As soil salinity increased, the productivity of B. gymnorrhiza was significantly reduced compared to R. stylosa, which led to an increase in biomass of R. stylosa despite the enhanced salt stress (> 30 ‰). These predicted patterns in forest structures across the soil salinity gradient remarkably agreed with field data, highlighting the control of salinity on productivity and tree competition as factors that shape the mangrove forest structures. The model reproducibility of forest structures was also supported by the predicted self-thinning processes, which likewise agreed with field data. Aside from soil salinity, seasonal dynamics in atmospheric variables (solar radiation and temperature) were highlighted as factors that influence mangrove productivity in a subtropical region. This physiological principle-based improved model has the potential to be extended to other mangrove forests in various environmental settings, thus contributing to a better understanding of mangrove dynamics under future global climate change.

1 Introduction

Mangrove forests grow in intertidal zones in tropical and subtropical regions (Giri et al., 2011) and store a large amount of carbon (C) especially in their soil, commonly referred to as “blue carbon”. It has roughly 4 times higher ecosystem-scale carbon stock than other forest ecosystems (Donato et al., 2011), characterizing them as globally important C sinks (Mcleod et al., 2011; Alongi, 2014; Taillardat et al., 2018; Sharma et al. 2020), therefore playing an important role in climate change mitigation. However, mangrove forests have declined worldwide; at least 35 % of the mangrove forests had disappeared in the 1980s and 1990s predominantly because of deforestation due to conversion to aquaculture ponds, rice fields, urban development, and palm oil plantations (Friess et al., 2019). Deforestation has been continuing until now particularly in Southeast Asia, with a recent estimate of mangrove loss rates between 0.11 %–0.70 % (Friess et al., 2019, 2020). The loss of mangrove soil C through mineralization following deforestation has been of concern as a source of carbon emission to the atmosphere in addition to the loss of C sequestration capacity (Atwood et al., 2017; Sharma et al. 2020; Adame et al., 2021). To facilitate effective mangrove conservation, management, and restoration, a better understanding of C sequestration rates and the soil C dynamics, hence mangrove blue C dynamics, under different environmental conditions and climate change is urgently needed.

While the mangrove soil C dynamics are complex and involve physical, biogeochemical, and ecological processes (Kristensen et al., 2008; Alongi, 2014; Bukoski et al. 2020) that still remain poorly understood, one of the most important variables determining soil C dynamics may be related to mangrove productivity. Mangroves supply their products, such as leaf litter and dead roots, to the soil C pool (Kristensen et al., 2008; Alongi, 2014; Ouyang et al., 2017), which are closely related to forest structural variables such as canopy height and above-ground biomass (AGB) (Saenger and Snedaker, 1993; Komiyama et al., 2008). Such autochthonous C accounts for a significant amount of total soil C in mangrove forests (Xiong et al., 2018; Sasmito et al., 2020). Therefore, the aim of this study is to successfully quantify and predict the biomass dynamics and growth processes of mangroves in different environmental conditions. These results would take a step forward in our understanding of mangrove C sequestration rate and soil C dynamics.

Although data and insights on mangrove AGB distributions in relation to environmental variables have recently increased (Simard et al., 2019; Roval et al., 2016, 2021), there is still no established way to predict the dynamics of mangrove AGB in the changing environmental conditions. Generally, the ecosystem's response to environmental variables is nonlinear, and biomass dynamics is cumulatively affected by a nonlinear response. Therefore, predicting the effect of one environmental variable on mangrove biomass dynamics is difficult if based only from the monitoring data on mangroves' biomass, which are exposed to the effects of multiple environmental variables. This makes the assessment of environmental impacts on mangrove biomass dynamics challenging if datasets from only the field-based monitoring approach are used.

The dynamic vegetation model (DVM) simulates vegetation or forest growth based on physiological principles that include processes such as tree competition, establishment, and mortality (Fisher et al., 2018). This model could be a way to overcome the limitation of field-based approach and predict mangrove biomass dynamics under multiple environmental variables. Various DVMs (e.g., big-leaf, cohort-based, individual-based) have been developed mainly for terrestrial ecosystems and have successfully reproduced the dynamics of various forests in the temperate, tropical, and boreal regions (Fisher et al., 2018). Recently, DVMs have advanced in physiological expression of stomatal conductance under water stress by incorporating a plant hydraulic model that explicitly solves plant water flux (Bonan et al., 2014; Xu et al., 2016; Li et al., 2021). Recent studies also identified plant hydraulics as a critical factor that determines the plants' biomass allocation pattern to leaves, stem, and roots (Magnai et al., 2000; Trugman et al., 2019b; Portkay et al., 2021), the variations of which could drive a significant variation in plant productivity (Trugman et al., 2019a).

In mangrove forests, the salt in soil porewater (soil salinity) is one of the significant environmental factors that determine the mangroves' distribution, productivity, structure, and zonation pattern (Ball and Farquhar, 1984; Clough and Sim, 1989; Sobrado, 2000; Ball, 2002; Suarez and Medina, 2005; Suwa et al., 2009; Barr et al., 2013; Nguyen et al., 2015). Therefore, it is essential to properly represent the effects of soil salinity on mangrove growth considering species differences in the tolerance of salinity in order to accurately predict the mangrove biomass dynamics. Soil salinity imposes highly negative water potential in the substrate, making the water acquisition energetically challenging for plants, which acts in a similar way to water stress (Reef and Lovelock, 2015). With this perspective, the theoretical works of Perri et al. (2017, 2019) demonstrated the importance of considering the plant hydraulics for predicting the photosynthetic and transpiration rates under salt stress. However, although there are several individual-based DVMs for mangroves (e.g., FORMAN by Chen and Twilley, 1998; Kiwi by Berger and Hildenbrandt, 2000; mesoFON by Grueters et al., 2014; and BETTINA by Peters et al., 2014), no model has yet considered salinity control role in photosynthesis and transpiration through plant hydraulics, suggesting room for improvement in the physiological representation of the mangrove biomass dynamics under the impacts of soil salinity. It is expected that the nutrient uptake rate is also affected by soil salinity through the regulated transpiration rate (Simunek and Hopmans, 2009), making nutrient availability one of the key factors controlling mangrove growth especially under high-soil-salinity conditions (Lovelock et al., 2004, 2006a, b; Feller et al., 2007; Reef et al., 2010). Nonetheless, the modeling studies have not explicitly considered the role of nutrient uptake in mangrove growth.

Figure 1 (a) Location of Ishigaki Island; (b) location and (c) aerial photo of the study site – Fukido mangrove forest. The white line in panel (c) indicates the boundary of mangroves and other land covers where mangroves are assumed to inhabit the areas of elevation < 1.0 m + mean sea level, which was delineated based on a lidar-derived digital elevation model (DEM). The blue lines indicate small creeks. The circular makers indicate survey plots' locations along with four transects (T–A to T–D), while the pie charts indicate species composition in each plot. The red arrows indicate outlets of rivers flowing into the mangrove forest (R1 to R4). The aerial photo and DEM products were obtained from Asia Air Survey Co. Ltd., Japan. Shorelines are from Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG).

Here we hypothesized that the individual-based DVM incorporating plant hydraulic traits can reasonably predict mangrove biomass, structure, and species zonation pattern across a soil salinity gradient without empirical expression of the soil salinity influence on mangrove productivity. Such model would advance the understanding of mangrove biomass dynamics under multiple environmental stresses, which ultimately influence the mangrove soil carbon dynamics. To test the hypothesis and contribute to the improvement of the physiological representation of mangrove growth specifically under soil salinity impacts, we developed a new individual-based DVM for the mangrove forest. The developed model is based on a terrestrial individual-based DVM – the SEIB-DGVM (Spatially-Explicit Individual-Based Dynamic Global Vegetation Model, Sato et al., 2007). We added a plant hydraulic model to SEIB-DGVM and coupled it with the photosynthetic model to consider the impacts of soil salinity on the mangrove water uptake and carbon gain. We also explicitly considered the role of nutrient uptake on biomass dynamics. Furthermore, a novel biomass allocation scheme linked with plant hydraulics and resource uptake rate was introduced as the mangroves' strategy to cope with salt stress and enhance the rate of production. We tested the developed model and determined the reproducibility of forest structures (e.g., species composition, biomass) in a subtropical mangrove forest in Japan with two coexisting species (Rhizophora stylosa and Bruguiera gymnorrhiza).

2 Materials and methods

2.1 Study sites

Our study site for the model application is an estuarine mangrove of the Fukido River (Fukido mangrove forest) in Ishigaki Island, Japan (Fig. 1, 24^∘20^′ S, 124^∘15^′ E). The site is characterized as a subtropical region. According to the climatological normal data obtained by the Japan Meteorological Agency, the annual-mean air temperature is 24.5 ^∘C, with a monthly average of 29.6 ^∘C in July and 18.9 ^∘C in January (see also Fig. 4). The mean monthly precipitation is 142 mm in July and 135 mm in January. Four small rivers (R1–R4) flow into the Fukido mangrove forest, while the river R2 has two outlets (Fig. 1c). The mean discharge rates of the rivers in October 2012 were less than 0.03 m³ s⁻¹ for R1, R3, and R4 and around 0.05 m³ s⁻¹ for R2 (Mori et al., unpublished data). The tide is semi-diurnal with the highest and lowest amplitude of 1.8 and 0.8 m, respectively (Egawa et al., 2021).

The site is vegetated by two species, R. stylosa and B. gymnorrhiza. The trees of R. stylosa dominated the seaward zone, especially areas close to the river mouth (Fig. 1c), while B. gymnorrhiza dominated the landward zone. The species R. stylosa is classified as a relatively salt-tolerant species, while B. gymnorrhiza is classified as a less salt-tolerant but shade-tolerant species (Putz and Chan, 1986; Sharma et al. 2012; Reef et al., 2015). According to Ohtsuka et al. (2019), the Fukido mangrove forest is a mature and intact mangrove forest designated as natural protection area by Ishigaki city, where distinct disturbances to the mangroves have not occurred since at least 1977 based on aerial photograph analysis.

2.2 Field data collection

We used the tree census data of the Fukido mangrove forest shown in Suwa et al. (2021) to assess model performance. The tree census data were collected from the survey plots established along four transects (T–A, T–B, T–C, and T–D), shown in Fig. 1c. The details of the survey protocol are described in Suwa et al. (2021). The stem biomass of individual trees (M_S, g) was estimated from a common mangrove allometric equation proposed by Komiyama et al. (2005), which was validated with various mangrove species:

$$\begin{array}{}\text{(1)}& M_{\mathrm{S}}=\mathrm{70}\mathit{\rho }{\left[{\left(\text{DBH}/\mathrm{100}\right)}^{\mathrm{2}}H\right]}^{\mathrm{0.931}},\end{array}$$

where ρ is the wood density (g cm⁻³), DBH is the stem diameter at breast height (m) divided by 100 for the unit conversion from meter to centimeter, and H is the tree height (m). However, tree height data were occasionally absent at some plots, especially along T–C and T–D, and in such cases, the tree height was estimated using a DBH–H allometric relationship (Fig. S1a and b in the Supplement). The AGB at each plot (Mg ha⁻¹) was then calculated from the estimated stem biomass.

The crown diameter was also measured for some selected trees, besides the data shown in Suwa et al. (2021). The trees for crown measurement were randomly selected at each transect. The diameters parallel and perpendicular to the transect line were measured for each tree, and the crown diameter (D_crown, m) was represented by the average of the values from the two directions. In total, crowns of 81 trees of R. stylosa and 103 trees of B. gymnorrhiza were measured (Supplement Fig. S1c and d).

Soil salinity and porewater dissolved inorganic nitrogen concentration (DIN) were also measured at each plot as environmental drivers of mangrove production. Soil samples were collected by inserting a PVC pipe into the soil at each plot, and soil porewater was extracted from the surface 10 cm soil sample. The porewater samples were kept frozen and brought to the laboratory for analysis. Salinity of the porewater (soil salinity) was measured using a salinity meter (PAL-SALT, ATAGO Co., Ltd., Japan), while DIN concentrations were measured using a QuAAtro 2-HR (SEAL Analytical Ltd., Germany, and BLTEC K.K., Japan). These measurements were conducted from August to September 2013. The summary of the environmental and vegetation variables at each plot is provided in Table S1 in the Supplement.

2.3 Model description

The mangrove growth model was formulated based on an individual-based model, SEIB-DGVM (Sato et al., 2007). The forest dynamics was represented by a 30 m × 30 m computational domain. In this domain, the irradiance distribution, tree establishment, death, and changes in plant morphology subsequent to growth were simulated (Sato et al., 2007). A feature of SEIB-DGVM is that it explicitly solves the effects of shading by neighboring trees on the light acquisition. The SEIB-DGVM thus provides the advantage in describing tree competition for light more than the other types of DVMs such as big-leaf or cohort-based models (Fisher et al., 2018). In SEIB-DGVM, the crown of each tree is represented by a cylindrical-shaped object divided by 0.1 m thick crown layers to account for the within-crown vertical variability in irradiance distribution. It is assumed that leaf biomass is evenly distributed in the crown layers.

Originally, the SEIB-DGVM defines four biomass pools – leaf, trunk, fine root, and stock (non-structural storage pool); the trunk includes both the above-ground stem and the below-ground coarse root (Sato et al., 2007). In this study, we considered the stem and coarse root separately to explicitly consider the role of coarse root turnover in the biomass dynamics (Castaneda-Moya et al., 2011; Adame et al., 2014). Additionally, we added a new biomass pool – the above-ground root, especially for Rhizophora species whose above-ground root, or “prop root”, could account for nearly 60 % of their AGB (Nishino et al., 2015; Vinh et al., 2019).

Figure 2 The model framework newly added to SEIB-DGVM for describing mangrove growth. The red box and arrows indicate the substrate conditions given in the model. The black boxes and arrows indicate processes computed in the hourly time steps, while the blue ones are for the daily time step.

The original SEIB-DGVM does not have a plant hydraulic module, and the effects of soil water on stomatal conductance were empirically parameterized. It also does not account for plant nutrient uptake; thus, the plant growth depends solely on photosynthesis. The biomass allocation is modeled based on scaling law (Trugman et al., 2019a). In this study, these processes that control plant growth were almost entirely modified to describe mangrove growth under salt stress (Fig. 2). The following sections explain the modification of the SEIB-DGVM for this study related to plant hydraulics. Other modifications to the SEIB-DGVM are summarized in Notes S3 and S4.

2.3.1 Inclusion of plant hydraulic module

The plant hydraulic module implemented in this study is primarily based on the model developed by Xu et al. (2016) in a soil–plant–atmosphere continuum scheme. Here we describe essential processes in the plant hydraulic module which will be related to the new biomass allocation model in the next section.

Table 1 Parameters constraining plant morphology, biomass proportion, and stoichiometry.

R. s: R. stylosa; B. g: B. gymnorrhiza.
^a Derived from DBH–H_max relationship. See Note S1 and Fig. S1 for details. ^b Computed in the model. See Fig. 3c and d. ^c Derived from DBH–$D_{\text{crown}}^{\ast }$ relationship. See Note S1 and Fig. S1 for details. ^d Average of values reported in Tamooh et al. (2008), Castañeda-Moya et al. (2011), Adame et al. (2014), Robertson and Alongi et al. (2003), and Muhammad-Nor et al. (2019). ^e Estimated from prop root allometry in the Fukido mangrove forest. See Fig. S3.

The plant water uptake rate (≈ sap flow rate; J_sap, kg H₂O per tree per second) is calculated as

$$\begin{array}{}\text{(2)}& J_{\text{sap}}={\displaystyle \frac{{\mathrm{\Psi }}_{\mathrm{s}}-{\mathrm{\Psi }}_{\mathrm{l}}-{\mathrm{\Psi }}_{\mathrm{h}}}{R_{\text{whole}}}},\end{array}$$

where R_whole is the whole-plant hydraulic resistance (MPa s tree kg⁻¹ H₂O), and Ψ_s and Ψ_l are the soil and leaf water potential (MPa), respectively; the Ψ_h = ρ_w g H 10⁻⁶, which is the gravitational water potential drop from the ground to the crown (MPa), where ρ_w is the water density (kg m⁻³) and g is the gravitational acceleration (m s⁻²). The parameter Ψ_s can be expressed as a sum of the matric potential and osmotic potential (Ψ_π, MPa). The parameter Ψ_π can be expressed as the difference in the osmotic potential between the soil and plant, which is linearly related to soil salinity and the partial uptake of the salt by mangroves represented by the salt filtration efficiency, ε (fraction) (Perri et al., 2017). Here, a constant water temperature value (25 ^∘C) was used to compute Ψ_π; however, note that sensitivity of Ψ_π to change in temperature is significantly small compared to salinity. Alternatively, the matric potential is negligibly small compared to Ψ_π in mangrove forests where the soil is usually water-saturated due to frequent tidal flooding (Perri et al., 2017). The parameter R_whole can be expressed as the sum of the root to stem hydraulic resistance (R_root) and the stem to leaf hydraulic resistance (R_stem), both expressed in (MPa s tree kg⁻¹ H₂O). The parameter R_root is given by

$$\begin{array}{}\text{(3)}& R_{\text{root}}={\displaystyle \frac{R_{\mathrm{r}}}{M_{\text{FR}}}},\end{array}$$

where R_r is the fine root hydraulic resistance per unit biomass (MPa s g kg⁻¹ H₂O) and M_FR is the fine root biomass (g per tree). The parameter R_stem is given by

$$\begin{array}{}\text{(4)}& R_{\text{stem}}={\displaystyle \frac{a_{\mathrm{1}}H}{K_{\text{sap}}{A}_{\text{sap}}}},\end{array}$$

where a₁ is the correction factor for tree height (H) to water path length, K_sap is the stem hydraulic conductivity per unit sapwood area ($\mathrm{kg}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{MPa}}^{-\mathrm{1}}$), and A_sap is the sapwood area of a tree (m² sapwood area per tree), which is calculated from the DBH and diameter ratio of the heartwood relative to the entire stem (β_heart, Table 1; Trugman et al., 2019b). The parameter K_sap can be expressed as a product of saturated xylem conductivity (K_sap,sat) and a factor representing the effect of xylem cavitation (Xu et al., 2016):

$$\begin{array}{}\text{(5)}& K_{\text{sap}}=K_{\text{sap,sat}}{\left(\mathrm{1}+{\left({\displaystyle \frac{{\mathrm{\Psi }}_{\mathrm{l}}}{P_{\mathrm{50}}}}\right)}^{a_{\mathrm{2}}}\right)}^{-\mathrm{1}},\end{array}$$

where P₅₀ (MPa) is the water potential at which 50 % of the xylem conductivity is lost and a₂ is an empirical parameter (dimensionless). The change in leaf water potential is governed by the equation

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_{\mathrm{l}}}{\mathrm{d}t}}={\displaystyle \frac{J-T_{\text{whole}}}{C_{\mathrm{p}}\text{LA}}},\end{array}$$

where T_whole is the whole-plant transpiration rate (kg H₂O per tree per second), LA is the whole-plant leaf area (m² leaf area per tree), and C_p is the plant capacitance ($\mathrm{kg}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{MPa}}^{-\mathrm{1}}$). The parameter T_whole is calculated by vertically integrating the product of the leaf-level transpiration rate and the leaf area in each crown layer. The leaf-level transpiration and photosynthetic rates and stomatal conductance are calculated using a leaf flux model of Bonan et al. (2014), where the stomatal conductance is estimated from an optimization approach of Cowan and Farquhar (1977) using the marginal water use efficiency ($\mathit{\lambda }=\mathrm{\Delta }{A}_{\mathrm{n}}/\mathrm{\Delta }E$, where λ is the optimal water use efficiency (WUE), and A_n and E are the leaf net photosynthetic rate and the transpiration rate, respectively) and regulated by Ψ_l. See Note S4 in the Supplement for the detailed calculations of A_n and E, as well as the linkage to Ψ_l.

The processes for transpiration, photosynthesis, plant water uptake, and change in leaf water potential were computed in hourly time step (Fig. 2). Overall, high salinity increases sensitivity of the leaf water potential to plant transpiration (Eqs. 2 and 6), which in turn may cause stomatal closure even with a low transpiration rate. It also increases the optimal WUE value leading to lower stomatal conductance (Ball and Farquhar, 1984; Clough and Sim, 1989; Barr et al., 2014; Perri et al., 2019), thereby lowering the photosynthetic and transpiration rates.

Figure 3 Schematics of (a, b) allometric and (c, d) physical constraints on tree height (H_max, H_con) and crown diameter ($D_{\text{crown}}^{\ast }$, D_crown,con), where the H_con and D_crown,con in panels (c) and (d) are for the tree with crown filled by yellow color and (e) the newly added biomass allocation scheme to SEIB-DGVM. See Note S1 for the derivation of allometric constraints from field data.

2.3.2 Inclusion of hydraulics and growth optimality-based biomass allocation

The biomass allocation occurs at the daily time step in the new biomass allocation scheme introduced in this study (Fig. 2). At each time step, four variables were considered for biomass allocation of individual trees – the daily C (C_grow, g C per tree per day) and N (N_grow, g N per tree per day) resources that can be used for tree growth, the daily minimum leaf water potential (Ψ_l,daymin, MPa), and the midday photosynthetically active radiation (PAR) at the crown top (PAR_top, $\mathrm{\mu }\mathrm{mol}\mathrm{photon}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). The C_grow and N_grow were computed from the daily C and N uptake rates, where the N uptake rate was calculated by multiplying the porewater DIN concentration and plant water uptake rate (see Note S5 in the Supplement for the detail). Biomass was allocated according to these variables to optimize the plant hydraulics and enhance the uptake rate of growth-limiting resource (C or N) under the constraints summarized in Table 1. Allometric and physical constraints were considered for H and D_crown (Fig. 3a–d; see Note S1 in the Supplement for the derivation of the allometric constraints).

The parameters C_grow and N_grow are allocated to the respective biomass pools in a scheme shown in Fig. 3. We applied the concept that plants keep their favorable hydraulic conditions throughout the growth periods by adjusting the morphological structures (Magnai et al., 2000). In this regard, we introduced a parameter Ψ_lk – the critical leaf water potential (MPa) – at which plants aim to maintain their leaf water potential (note that Ψ_lk is different from Ψ_l,min at which plants close the stomata). It was then considered that when Ψ_l,daymin falls below Ψ_lk, the plant tries to reduce R_whole by allocating biomass to either the fine root or stem, which reduces R_whole more effectively (Case 1 and 2 in Fig. 3; note that decreases in R_stem and R_root were expressed by a negative value):

$$\begin{array}{}\text{(7a)}& {\displaystyle }{M}_{\mathrm{S},t}=M_{\mathrm{S},t-\mathrm{1}}+{\text{dM}}_{\mathrm{S}}\mathrm{\text{ if }}{\displaystyle \frac{{\text{dR}}_{\text{root}}}{{\text{dM}}_{\text{FR}}}}>{\displaystyle \frac{{\text{dR}}_{\text{stem}}}{{\text{dM}}_{\mathrm{S}}}},\text{(7b)}& {\displaystyle }{M}_{\text{FR},t}=M_{\text{FR},t-\mathrm{1}}+{\text{dM}}_{\text{FR}}\mathrm{\text{ if }}{\displaystyle \frac{{\text{dR}}_{\text{root}}}{{\text{dM}}_{\text{FR}}}}<{\displaystyle \frac{{\text{dR}}_{\text{stem}}}{{\text{dM}}_{\mathrm{S}}}},\end{array}$$

where M_FR,t and M_S,t are the fine root and stem biomass (g per tree) at time step t (day), and dM_FR and dM_S are the daily biomass increment potential of fine root and stem (g per tree per day), respectively, which are limited by either of C_grow and N_grow and represented as

$$\begin{array}{}\text{(8a)}& {\displaystyle }\begin{array}{rl}{\text{dM}}_{\text{FR}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{C}}_{\mathrm{M}}}}\times min[& {\mathrm{C}}_{\mathrm{grow}}(\mathrm{1}-F_{\text{gr}})(\mathrm{1}-F_{\text{CR,C}}),\\ & {\mathrm{N}}_{\mathrm{grow}}(\mathrm{1}-F_{\text{CR,N}}){\mathrm{CN}}_{\mathrm{r}}],\end{array}\text{(8b)}& {\displaystyle }\begin{array}{rl}{\text{dM}}_{\mathrm{S}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{C}}_{\mathrm{M}}}}\times min[& {\mathrm{C}}_{\mathrm{grow}}(\mathrm{1}-F_{\text{gr}})(\mathrm{1}-F_{\text{AR}}),\\ & {\mathrm{N}}_{\mathrm{grow}}(\mathrm{1}-F_{\text{AR}}){\mathrm{CN}}_{\mathrm{w}}],\end{array}\end{array}$$

where C_M is the carbon mass per unit dry weight in plant tissue (g C g⁻¹ DW), F_gr is the growth respiration fraction, F_CR,C and F_CR,N are the fractions of C_grow and N_grow, respectively, to be allocated to the coarse root to realize β_FR (target fine root biomass relative to coarse root; Table 1), F_AR is the fraction of the resources to be allocated to the above-ground root to realize β_AR (target prop root biomass relative to stem; Table 1, also see Fig. S3 in the Supplement) which was determined from an allometric model using DBH obtained in our study site by Yoshikai et al. (2021), and CN_r and CN_w are the CN ratios in fine root and woody tissue (g C g⁻¹ N), respectively, that convert the unit of N_grow to C_grow. In Eq. (7), the ${\text{dR}}_{\text{root}}/{\text{dM}}_{\text{FR}}$ is calculated from Eq. (3), while the dR_stem/dM_S is calculated from

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{{\text{dR}}_{\text{stem}}}{{\text{dM}}_{\mathrm{S}}}}={\displaystyle \frac{{\text{dA}}_{\text{sap}}}{{\text{dM}}_{\mathrm{S}}}}\times {\displaystyle \frac{{\text{dR}}_{\text{stem}}}{{\text{dA}}_{\text{sap}}}},\end{array}$$

where ${\text{dA}}_{\text{sap}}/{\text{dM}}_{\mathrm{S}}$ is obtained from Eq. (1) by calculating the increase in DBH with stem biomass increment dM_S without height growth, and ${\text{dR}}_{\text{stem}}/{\text{dA}}_{\text{sap}}$ is given from Eqs. (4) and (5), where Ψ_l,daymin is used in Eq. (5). It should be noted that the variables Ψ_l,daymin, C_grow, and N_grow change with various factors including atmospheric and substrate variables and tree competition, and no absolute optimal biomass proportion achieves the condition ${\text{dR}}_{\text{root}}/{\text{dM}}_{\text{FR}}$ = ${\text{dR}}_{\text{stem}}/{\text{dM}}_{\mathrm{S}}$ throughout the computational period. Also, due to the different CN ratios in fine root and woody tissues, the increment in stem biomass (dM_S) with a unit N resource is greater than that of the fine root biomass (dM_FR) under N-limited conditions (Eq. 8, Table 1).

Alternatively, if plants are not stressed by the lowered leaf water potential (Ψ_l,daymin > Ψ_lk), the resources are allocated to a plant organ that effectively increases the uptake rate of either C or N, limiting the growth rate. Under N-limited conditions, plants allocate biomass to the leaves to increase whole-plant transpiration capacity, which increases N uptake rate nearly proportionally (as suggested by Eq. S22 in the Supplement) (Case 3 in Fig. 3); this is considering that the limited uptake of N is due to the small transpiration rate rather than water uptake regulation by hydraulic resistance. The increase in leaf biomass increases D_crown and dLAI (leaf area index per 1 m vertical height) depending on the crown diameter constraints ($D_{\text{crown}}^{\ast }$ and D_crown,con; see Note S6 in the Supplement for the details). However, if the increase in leaf biomass is inhibited by dLAI_max (maximum dLAI; Table 1) and crown diameter constraints, the resources are allocated to the stem for height growth, which in turn will make a new crown layer and eventually allow further leaf accommodation (Case 4 in Fig. 3). Under a C-limited condition, the limited C uptake rate may be attributed to low light availability or small whole-plant leaf area. In this regard, we introduced a criterion PAR_k, where the photosynthetic rate is reduced by half of the light-saturated photosynthetic rate, allowing the assumption that the limited C uptake rate is due to low light availability if PAR_top is lower than PAR_k. In this case, the resources are allocated to the stem for height growth to acquire better light conditions under tree competition (Case 5 or 6 in Fig. 3); otherwise, the resources are allocated for an increase in leaf area (Case 3 or 4 in Fig. 3). Lastly, the residual C_grow or N_grow after the biomass allocation is allocated to the stock pool.

2.4 Simulation configuration

The model was applied to the Fukido mangrove forest to test its performance in reproducing the forest structural variables (species composition, mean DBH, and AGB). The model was forced with atmospheric variables (air temperature, relative humidity, atmospheric pressure, wind speed, and cloud fraction) and substrate conditions (soil salinity and porewater DIN). Direct and diffused solar radiation and longwave radiation were calculated in SEIB-DGVM from the given variables such as cloud fraction, air temperature, and latitude (Sato et al., 2007). The atmospheric variables for the Fukido mangrove forest given to the model were derived from a global reanalysis product JRA-55 (Kobayashi et al., 2015). For long-term simulation (i.e., more than 100 years), the yearly atmospheric variation in 2013, a year when the field data collection was conducted, was repeatedly given in the simulation.

Simulations with different soil salinity, or the “salinity gradient simulation”, which varied from 18 ‰ to 36 ‰ with 2 ‰ intervals, were conducted to reproduce the forest structural variables across a soil salinity gradient. For the porewater DIN, a spatially averaged DIN (average of DIN measured at the survey plots: 200 µmol L⁻¹) was given to the model as the representative value of the porewater DIN in this forest. In each simulation, soil salinity and the porewater DIN were set as constant due to lack of data and model on the temporal variations in substrate conditions. We also conducted “plot-wise simulation”, or the simulation for each survey plot, by giving the measured soil salinity and porewater DIN at each plot. Note that the results shown in this paper are from the salinity gradient simulation; the results of the plot-wise simulation are provided in Fig. S5 in the Supplement and discussed later.

The initial condition was set as bare land (no vegetation) for all simulations. Tree establishment occurs at 1 m × 1 m grid cells at a yearly time step according to light condition at the forest floor and a parameter of establishment probability (P_establish, ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) prescribed for each species (Sato et al., 2007). The species that will establish at a grid cell is determined according to a fraction of total biomass of each species in the computational domain such that a species occupying a larger fraction has a higher probability of establishment. On the other hand, it is sometimes randomly determined by a probability Est_random, where the value of Est_random was set to 0.05 in this study. This corresponds to Scenario 4 in the tree establishment scheme in SEIB-DGVM (see Sato, 2015, for the details). We followed Sato et al. (2007) for the initial conditions (tree morphology and biomass proportion) of the established trees.

The SEIB-DGVM uses stochastic models for the processes of tree establishment and mortality, and for this reason the result of a simulation varies every time. In this regard, we conducted ensemble simulations (20 runs) for each soil salinity in the salinity gradient simulation and extracted the general trends.

Table 2 Model parameters related to plant hydraulics and productivity.

R. s: R. stylosa; B. g: B. gymnorrhiza.
^a Value for Rhizophora mangle. ^b The value used for terrestrial forest ecosystem was applied due to lack of information. ^c The minimum of the reported values was adopted. ^d Value for Rhizophora apiculata. ^e The average value of data in Castañeda-Moya et al. (2011) was adopted. ^f Value for Rhizophora stylosa.

The model parameter settings related to plant hydraulics and productivity are summarized in Table 2. Other minor model parameters are summarized in Table S2 in the Supplement. The parameter values for the two species in the Fukido mangrove forest, R. stylosa and B. gymnorrhiza, were determined based on literature. If the data for a focal species were unavailable from the literature, the data from the genus or family were applied. Some parameter values were adapted from other mangrove genus or terrestrial ecosystems, and in this case, the same value was given to the two species (Table 2). The values of two plant hydraulic trait parameters – Ψ_lk (critical leaf water potential) and β₀ (sensitivity of marginal WUE to leaf water potential in Eq. S21 in the Supplement; see Note S3) – were calibrated to reproduce the AGB and mean DBH of each species across the soil salinity gradient.

The Fukido mangrove forest's age is unknown, which makes the comparison between the model and field data difficult. However, considering that it is an old and mature forest intact at least since 1977 (Ohtsuka et al., 2019), we assumed that the forest structural variables of the Fukido mangrove forest are in steady states. We conducted long-term simulations for 450 years with this assumption, and we extracted the modeled DBH and AGB in steady states (> 300 years) and compared them with the field data.

Lastly, we performed sensitivity analysis of the plant hydraulic trait parameters (ε, P₅₀, Ψ_lk, and β₀) to see the relative importance of each parameter in reproducing the observed pattern of the forest structure, specifically AGB, across the soil salinity gradient. We changed the value of a target parameter of one species (either R. stylosa or B. gymnorrhiza) to the one determined for the other species, which is shown in Table 2, and ran the salinity gradient simulation. Note that to examine the sensitivity to Ψ_lk, we changed the values of Ψ_lk and Ψ_l,min to keep the buffer between the two parameter values. Also, to save on computational cost, we run only one simulation for each sensitivity test instead of the ensemble approach described above. Model sensitivities are shown in Fig. S6 in the Supplement.

Figure 4 Seasonal variations in atmospheric forcing variables: (a) solar radiation, (b) air temperature, and (c) vapor pressure deficit (VPD). Modeled seasonal dynamics: (d) monthly mean and standard deviation of species-specific gross photosynthetic rate (P_g, g C per square meter ground per day), (e) transpiration (T, mm d⁻¹) normalized by leaf layer index (LAI, m² m⁻²) of the respective species, and (f) midday (Ψ_l,midday) and predawn (Ψ_l,predawn) leaf water potential of each species. R. s: R. stylosa; B. g: B. gymnorrhiza. Solar radiation is expressed as daily sum, while air temperature and VPD are expressed as daily mean. The leaf water potential shown is the median value of individuals. Here, the modeled dynamics are from a simulation of soil salinity set as 30 ‰, and the results of a year when LAI reached 1.55 are shown. At this time, the LAIs of R. stylosa and B. gymnorrhiza are 0.87 and 0.68, respectively. In panel (d), seasonal variations in $P_{\mathrm{g}}/\mathrm{LAI}$ measured by Okiomoto et al. (2007) are also shown as reference, the data of which are from an area with LAI = 1.55 in the Fukido mangrove forest in 2000–2001.

3 Results

3.1 Modeled seasonal and diurnal dynamics

Seasonal variations in atmospheric forcing variables and modeled species-specific gross photosynthetic rate (P_g) and transpiration (T) normalized by the leaf area index (LAI) and midday (Ψ_l,midday) and predawn (Ψ_l,predawn) leaf water potential are shown in Fig. 4. The modeled variables are from one of the ensemble simulations with soil salinity set as 30 ‰. The model demonstrated strong seasonality in photosynthesis and transpiration primarily due to seasonality in solar radiation and air temperature. The model predicted the peak of $P_{\mathrm{g}}/\mathrm{LAI}$ in June, with values ∼ 5.1 and 4.9 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for R. stylosa and B. gymnorrhiza, respectively, and the peak of $T/\mathrm{LAI}$ in July–September, with values ∼ 1.07 and 0.85 mm d⁻¹ for each species, respectively. The $P_{\mathrm{g}}/\mathrm{LAI}$ and $T/\mathrm{LAI}$ were predicted to be depressed during winter (December–February), with values ∼ 3.0 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for both species and ∼ 0.43 and 0.36 mm d⁻¹ for each species, respectively. We compared the modeled leaf-level P_g with the field-estimated values in the Fukido mangrove forest by Okimoto et al. (2007). Their measurements were conducted in an area where the LAI is 1.55, the same LAI as the one shown in Fig. 4d; thus, the effects of LAI on leaf-level P_g could be eliminated for comparison. Although the modeled $P_{\mathrm{g}}/\text{LAI}$ of both species is slightly lower than the one obtained by Okimoto et al. (2007) (∼ 1.0 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), especially from June to August, overall, the model agreed well with their results.

The midday leaf water potential showed seasonal variations as with photosynthesis and transpiration (Fig. 4f). Due to the partial salt uptake of R. stylosa (as indicated by the lower ε value of this species, Table 2) that alleviates the osmotic potential difference between the soil and plant, the predawn leaf water potential of R. stylosa was constantly higher than that of B. gymnorrhiza (Fig. 4f). Rhizophora stylosa also showed larger magnitude of leaf water potential reduction at midday during summer compared to B. gymnorrhiza and higher leaf-level transpiration rate (Fig. 4e and f). During winter, due to the lowered transpiration rate, the leaf water potential reduction at midday was resultantly alleviated compared to summer.

Diurnal variations of the simulated photosynthesis, transpiration, and leaf water potential of the two species during summer and winter under two different salinity conditions (30 ‰ and 24 ‰) are shown in Fig. 5. Compared to 24 ‰ salinity, both species showed significantly lowered leaf-level transpiration rates under 30 ‰ salinity especially during summer (Fig. 5b and e), suggesting downregulation of stomatal conductance under high-soil-salinity conditions. On the other hand, the decrease in leaf-level photosynthetic rates was not significant (Fig. 5a and d). The leaf water potential during nighttime was lower when soil salinity was 30 ‰ compared to conditions when salinity was 24 ‰, due to the different osmotic potential in soil porewater. The leaf water potential, however, showed almost the same levels at midday during summer, which were close to the values of Ψ_lk determined for each species (Fig. 5c, Table 2). The reduction in leaf water potential to the level of Ψ_lk suggests the role of dynamic biomass allocation, which adjusts the whole-tree transpiration demands and hydraulic conductivity, in constraining the leaf water potential dynamics (Fig. 3). In contrast, the diurnal dynamics in leaf water potential during winter showed similar magnitude of reduction of the water potential at midday between the two soil salinity conditions (Fig. 5f).

Figure 5 Simulated averaged diurnal variations in (a, d) photosynthesis and (b, e) transpiration of R. stylosa (R. s) and B. gymnorrhiza (B. g) normalized with LAI of the respective species, and (c, f) leaf water potential of the two species for summer (June–August) and winter (December–February) under two soil salinity conditions (30 ‰ and 24 ‰). The variations under 30 ‰ soil salinity correspond to the results shown in Fig. 4. The variations under 24 ‰ soil salinity are from the results of a year that showed the same LAI (1.55). The diurnal variations in leaf water potential were derived based on the median value of individuals.

Figure 6 Visualization of forest structures over 200 years under different soil salinity (sal), 20 ‰, 24 ‰, 30 ‰, and 34 ‰, taken from one of the ensemble simulations. The brown-colored objects represent the stem while the yellow- and green-colored objects represent the crowns of R. stylosa and B. gymnorrhiza, respectively. The forest floor shown is the 30 m × 30 m wide computational domain.

Figure 7 Temporal dynamics in above-ground biomass (AGB), leaf area index (LAI), and mean diameter at breast height (DBH) of R. stylosa (R. s) and B. gymnorrhiza (B. g) in four soil salinity conditions: (a) 20 ‰, (b) 24 ‰, (c) 30 ‰, and (d) 34 ‰. Note that trees with DBH < 0.05 m were not included in the calculation of mean DBH. Solid lines show median and shading the 90th percentile from ensemble simulations.

3.2 Modeled biomass dynamics under different soil salinity

Figure 6 shows the changes in the forest structures for over 200 years under different soil salinity conditions, 20 ‰, 24 ‰, 30 ‰, and 34 ‰, from one of the ensemble simulations (the present-day average soil salinity of the survey plots is 28 ‰). The time-series results of AGB, LAI, and mean DBH of the two species are shown in Fig. 7. Trees with DBH < 0.05 m were not accounted for in the calculation of the mean DBH because it is sensitive to the presence of small trees. Overall, the model demonstrated the significant influence of soil salinity on species composition and forest structure.

The model predicted that B. gymnorrhiza dominates over R. stylosa when soil salinity is 20 ‰ or 24 ‰ (Fig. 7a and b). Under soil salinity of 20 ‰, the AGB of B. gymnorrhiza exponentially increased up to 200 Mg ha⁻¹ after 60 years since the initial condition. It slightly decreased after that and was kept almost constant at 175 Mg ha⁻¹ after 150 years. The LAI of this species showed almost the same trend with AGB, while the mean DBH showed fluctuation especially in the first 200 years (Fig. 7a). The sudden decrease in the mean DBH is attributed to the onset of formation of forest gaps resulting from deaths of large B. gymnorrhiza trees that promoted the establishment of small trees (Fig. 6). After the decrease in the mean DBH, it gradually increased again and saturated at 0.17 m (Fig. 6a). Alternatively, the AGB and LAI of R. stylosa were significantly lower than B. gymnorrhiza, with its peak at only 25 Mg ha⁻¹ and 1 m² m⁻², respectively. This can also be seen in the decreasing number of R. stylosa trees subsequent to forest growth (Fig. 6). In contrast to AGB and LAI, the mean DBH of R. stylosa reached around 0.2 m after 75 years, as large as that of B. gymnorrhiza in steady state (Fig. 7a). This suggest that some R. stylosa trees can grow until mature conditions (see also Fig. 6), while trees of this species with DBH > 0.05 m disappeared in all ensemble simulations after 300 years (Fig. 7a). The trees of R. stylosa sometimes emerge due to the random factor in the establishment process, but most of the trees did not grow more than DBH of 0.05 m in the canopy of B. gymnorrhiza.

The trend in forest growth under 24 ‰ salinity was similar to that of 20 ‰ (Figs. 6 and 7b) but showed a slightly lower and higher peak for B. gymnorrhiza and R. stylosa, respectively, of the AGB, LAI, and mean DBH. This suggests decreased productivity of B. gymnorrhiza compared to 20 ‰ soil salinity and increased productivity of R. stylosa albeit with the increase in soil salinity. The survival rate of R. stylosa was higher than the results for 20 ‰, resulting in the high mean DBH of this species throughout the simulation period (Fig. 7b).

When the soil salinity was 30 ‰, the AGB of B. gymnorrhiza significantly decreased compared to the results for 20 ‰ and 24 ‰ salinities, becoming equivalent to those of R. stylosa (Fig. 7c). The LAI and mean DBH also showed a significant decrease, suggesting significantly lowered productivity of B. gymnorrhiza. The AGB and LAI of R. stylosa significantly increased compared to the results for 20 ‰ and 24 ‰, but the mean DBH significantly decreased.

The model predicted that B. gymnorrhiza cannot grow well at 34 ‰ soil salinity and that R. stylosa dominates under this salinity condition (Figs. 6 and 7d). Despite the further decrease in AGB, LAI, and mean DBH of B. gymnorrhiza, those of R. stylosa showed almost the same level for these parameters at 30 ‰ soil salinity.

3.3 Comparison between modeled and field-measured forest structural variables

Figure 8 shows the field-measured and modeled mean DBH and AGB of R. stylosa and B. gymnorrhiza across the soil salinity gradient. The field data clearly showed the effects of soil salinity on forest structural variables – decrease in mean DBH for both species and decrease in AGB of B. gymnorrhiza but increase in AGB of R. stylosa with increasing soil salinity. The model reproduced well the said patterns across the soil salinity gradient, and the values are within or close to the field data variations (Fig. 8). The change in species composition is also well reproduced, suggesting that the model can reproduce the forest structural variables across the soil salinity gradient.

Figure 8 Comparison of field-measured and modeled (a) mean DBH and (b) AGB of R. stylosa and B. gymnorrhiza along with soil salinity gradient. From each ensemble simulation, modeled mean DBH and AGB in steady states (> 300 years) were extracted and pooled for all ensembles, and the median (solid line) and the 90th percentile (shading) of the pooled samples are shown. Note that trees with DBH < 0.05 m were not included in the calculation of the mean DBH.

Figure 9 shows the field-measured and modeled relationship of tree density and mean individual stem biomass. Although there are some discrepancies between the model and field data especially for soil salinity conditions > 30 ‰, the model reproduced the overall pattern of the field data.

Figure 9 The relationship of tree density and mean individual stem biomass (M_S). Triangles show field data while circles show modeled values from one of the ensemble simulations with different soil salinity settings (from 18 ‰ to 34 ‰ with 2 ‰ increments) plotted from 300–450 years (with interval of 50 years), which are in steady states in terms of forest structural variables (see Fig. 7). Note that trees with DBH < 0.05 m were not counted in calculating tree density and mean M_S. The line represents the full density curve proposed by Tabuchi et al. (2013): y = 20 389x^−1.567.

4 Discussion

4.1 Model performance

Forest growth is influenced by leaf-level and whole-plant CO₂, water and nutrient fluxes, and forest-scale tree competition, which are all interconnected. The leaf-level fluxes were simulated using a well-established stomatal optimization scheme with the marginal WUE linked with leaf water potential (Bonan et al., 2014; Xu et al., 2016). The model predicted the distinct seasonal dynamics in photosynthesis and transpiration as well as leaf water potential in the Fukido mangrove forest (Figs. 4 and 5). The modeled seasonal variations in leaf-level photosynthesis ($P_{\mathrm{g}}/\mathrm{LAI}$) agreed well with the one measured by Okimoto et al. (2007) in this forest (Fig. 4d). Although there are no data on the seasonal variations in transpiration in this forest, studies on other subtropical mangrove forests, such as the Everglades National Park, Florida (Barr et al., 2014), and China (Liang et al., 2019), that incorporated the eddy-covariance approach also showed strong seasonality in transpiration, similar to the one predicted for the Fukido mangrove forest in this study (Fig. 4e). The evapotranspiration rate normalized by LAI in the Everglades measured by Barr et al. (2014) was 0.4–1.2 mm d⁻¹, which is close to the variation of the modeled $T/\mathrm{LAI}$ in the Fukido mangrove forest (Fig. 4e). These results suggest that the model produced realistic seasonal dynamics for transpiration in the Fukido mangrove forest.

Tree growth was driven by C and N uptake rates in the developed model resulting from the leaf-level and the whole-plant CO₂ and water fluxes. The modeled growth rates at a soil salinity condition where B. gymnorrhiza is the dominant species (sal < 28 ‰) showed close values to the ones measured by Ohtsuka et al. (2019) at a B. gymnorrhiza-dominated site in the Fukido mangrove forest (Fig. S4 in the Supplement). This suggests that the model also reasonably predicted the growth rate of each species in addition to the leaf-level CO₂ and water fluxes.

This model also showed reasonable reproducibility of the self-thinning process arising from tree competition. This was inferred from the decrease in tree density with the increase in individual tree biomass patterns based on the agreement of the measured and modeled tree density-mean M_S relationship, except for those with soil salinity > 30 ‰ (Fig. 9). An exponent value close to $-\mathrm{3}/\mathrm{2}$ was obtained, similar to what is observed in the Fukido mangrove forest (Suwa et al., 2021; Fig. 9) and in many mangrove forests as well (e.g., Analuddin et al., 2009; Deshar et al., 2012; Khan et al., 2013; Tabuchi et al., 2013; Azman et al. 2021). This was achieved by implementing the species-specific morphological traits especially the DBH–$D_{\text{crown}}^{\ast }$ relationship (Fig. 3b; see also Note S1 and Fig. S1). The underestimation trend of modeled tree density at high soil salinity (> 30 ‰) where R. stylosa starts to dominate (Fig. 9) may be attributed to the inaccurate representation of the crown morphological trait of this species, which generally gives larger $D_{\text{crown}}^{\ast }$ compared to observed values (see Note S1 and Fig. S1c). Basically, the crown diameters of individuals determine the tree accommodation spaces, and therefore the overestimated crown diameter may have resulted in the underestimation of the tree density. Crown size representation could be a factor that drives a large part of the uncertainty in DVMs (Meunier et al., 2021). Nevertheless, the data are remarkably scarce in the case of mangroves. The morphological traits of crown size should be investigated in future studies for a more realistic representation of mangroves' tree competition and forest dynamics in the model. It is also important to note that some variables of model prediction such as LAI, shoot / root biomass ratio, morphological plasticity in accordance to changes in environmental variables (as shown in Fig. S7 in the Supplement), and leaf water potential dynamics have not been validated due to lack of data, and future research is needed to address these aspects.

Overall, this is the first modeling study to introduce detailed physiological and mechanistic representations of the mangrove forest growth controlled by photosynthesis, water and nutrient (represented by DIN) uptake, and tree competition, and it is a realistic or accurate reproduction of mangrove growth processes. The remarkable agreement of modeled forest structures with field data across a soil salinity gradient validated our hypothesis – individual-based DVM incorporating plant hydraulic traits can reasonably predict mangrove growth processes under salt stress without empirical expression of the soil salinity influence on mangrove productivity. However, the model still does not account for the plant-to-soil feedback through water uptake, which has been identified by a mangrove-growth–groundwater-flow coupled model (Bathmann et al., 2021) as an important factor affecting both mangrove and substrate conditions (soil salinity). Alternatively, the said model also demonstrated that the forest structural variable and soil salinity dynamics could reach steady states after some time from the initial condition, a setting that is considered to describe the Fukido mangrove forest (Ohtsuka et al., 2019). Our modeling results, which did not include the plant-to-soil feedback, therefore may be valid only for the steady states and still hold uncertainty in the developmental stage. This further implies that model application may be limited only to mature mangrove forests, and further model improvement is needed for its application to forests during the developmental stage (after plantation) or during the recovery stage (after disturbances such as typhoons and deforestation).

4.2 Soil salinity and interspecific competition shaping the forest structural variables

Overall, the model explained that the changes in mean DBH and AGB of the two coexisting species with change in soil salinity are due to the difference in their salt tolerance and interspecific competition (Figs. 7 and 8). While the model predicted the contrasting changes of AGB of the two species, both species showed a decrease in productivity with an increase in soil salinity as seen in the monotonic decrease in DBH (Figs. 7 and 8). This decrease in productivity can be partly explained by the downregulation of stomatal conductance under high-soil-salinity conditions (Fig. 5). In addition, the changes in the biomass allocation pattern that increased the allocation to the stem and roots relative to leaves with an increase in soil salinity have influenced productivity (Fig. S7) – a pattern that reduced the whole-plant photosynthesis and transpiration and increased the carbon (through the stem and root respiration and root turnover) and nitrogen (through the root turnover) cost relative to the unit leaf area. It should be noted that such morphological plasticity with changes in soil salinity predicted by the model qualitatively agrees with the implications by other studies (e.g., Suwa et al., 2008, 2009; Vovides et al., 2014; Nguyen et al., 2015; Chatting et al., 2020), but future studies are needed for quantitative and systematic validation.

The sensitivity analysis of the plant hydraulics trait parameters provided some insights into the different salt tolerance of the two species that shaped the forest structures along the soil salinity gradient (Fig. S6). For example, it showed the substantial contribution of the partial salt uptake of R. stylosa, represented by the lower ε, to the salt tolerance of this species (Fig. S6a) at the possible expense of higher P₅₀ value (Table 2, Fig. S6c and d), which is considered as the coordinated functional traits (Jiang et al., 2017). The model showed the highest sensitivity to the parameters Ψ_lk and Ψ_l,min (Fig. S6e and f), suggesting that the mangroves' capacity to reduce the leaf water potential is one of the most important functional traits characterizing their salt tolerance, as suggested by Reef and Lovelock (2015). The response of AGB to changes in Ψ_lk, a parameter controlling biomass allocation pattern, also indicates the substantial impact of biomass allocation dynamics influenced by salinity on plant productivity. On the other hand, the model showed minimal sensitivity to β₀ (Fig. S6g and h). While the higher stomatal conductance of R. stylosa than B. gymnorrhiza (as shown in Figs. 4 and 5) qualitatively agreed with the implications by Clough and Sim (1989) and Reef and Lovelock (2015), the model results suggested that the choice of −0.6 for β₀ already leads to efficient stomatal openings for photosynthesis compared to the case of −0.4 for β₀ (Table 2). This may explain the small variations in the simulated leaf-level photosynthetic rates between the two species and among the different soil salinity levels (Figs. 4 and 5). Understanding the mangroves' stomatal behavior relative to soil salinity and covariation of leaf water potential and photosynthesis have not been well established from field data (Perri et al., 2019). Further field-based research and data implementation to the model are needed for better and more reliable representation of mangroves' stomatal conductance and associated regulation of photosynthesis under salt stress.

The model specifically predicted that B. gymnorrhiza competes over R. stylosa when soil salinity is favorably low for the growth of B. gymnorrhiza (sal < 28 ‰), an observation that is consistent with our field data and the data from other mangrove forests (Putz and Chan, 1986; Enoki et al., 2014). This result may be attributed to the following model parameter settings based on literature – higher wood density (ρ), smaller specific leaf area (SLA), and higher leaf turnover rate (TO_l) of R. stylosa than B. gymnorrhiza (Table 2). Higher ρ indicates the requirement of higher biomass increase for the height or radial growth of the stem. Smaller SLA and higher TO_l indicate the higher requirement of C and N to produce new leaf tissues or to keep the same amount of leaves, i.e., the need of R. stylosa for more C and N resources for growth compared to B. gymnorrhiza. The biomass requirement of prop roots, which lowers the biomass allocation to the stem (Fig. S3), and the smaller $D_{\text{crown}}^{\ast }$ of R. stylosa compared to B. gymnorrhiza (Fig. S1c and d) may also have contributed to the former's lower growth rate. Consequently, B. gymnorrhiza grew faster and suppressed the growth of R. stylosa by severe shading (Figs. 6 and 7). The higher growth rate of B. gymnorrhiza compared to R. stylosa at relatively low-salinity conditions agrees with the study by Jiang et al. (2019).

Interestingly, our model was able to simulate unique conditions not previously reported by other modeling works. For instance, the model predicted that R. stylosa trees could grow until the mature conditions under the canopy of B. gymnorrhiza-dominated forest provided the chance of favorable light conditions, resulting in the high mean DBH but low AGB of this species at relatively low soil salinity (∼ 24 ‰) (Figs. 6 and 7). Simulating this kind of process may only be possible through the individual-based approach with calculations of detailed irradiance distribution as done by the SEIB-DGVM in this study. Alternatively, the model predicted the significantly lowered growth rate of B. gymnorrhiza at high-soil-salinity conditions (sal > 30 ‰) where B. gymnorrhiza cannot grow until mature conditions, which resulted in the low AGB and small mean DBH of this species. This reduced the suppression of B. gymnorrhiza on R. stylosa and generated the Rhizophora stylosa-dominated forest (Figs. 6 and 7). Despite the abundant population of R. stylosa, the sizes of individuals were relatively small due to high salt stress and resulted in the high AGB but small mean DBH of this species.

4.3 Effects of other factors and further model improvement

Besides soil salinity, this study highlighted the importance of atmospheric variables as important drivers controlling mangrove production. This is seen in the photosynthesis–transpiration seasonal dynamics with peak during summer (June–September) and depression during winter (November–March) (Fig. 4) that none of the previous mangrove modeling studies has examined yet. The model predicted winter depression primarily due to low solar radiation and air temperature. Specifically, low air temperature (< 20 ^∘C) significantly reduced photosynthetic capacity – the maximum carboxylation rate and the maximum electron transport rate (Aspinwall et al., 2021); this, in turn, decreased the marginal WUE ($\mathrm{\Delta }{A}_{\mathrm{n}}/\mathrm{\Delta }E$), leading to the downregulation of stomatal conductance, a behavior of mangroves' stomata observed under low-temperature conditions (Akaji et al., 2019; Aspinwall et al., 2021). This resulted in the depression of photosynthesis and transpiration during this season. The significance of atmospheric control on stomatal conductance and associated dynamics in winter compared to salinity control is also highlighted in the similar magnitude of reduction of the leaf water potential at midday between the different soil salinity conditions in this season compared to summer (Fig. 5c and f). Such winter depression lowers the production of mangroves in subtropical regions and may be differentiated from tropical mangroves in terms of productivity. This could be a key factor in explaining and predicting the latitudinal gradients in mangroves' structural variables such as canopy height and AGB with the highest values at the equatorial region (Saenger and Snedaker, 1993; Simard et al., 2019; Rovai et al., 2021).

The model gave significantly better prediction of the AGB spatial distribution when the spatially averaged DIN concentrations were applied to the substrate condition compared to plot-wise DIN values (Fig. S5, plot-wise simulation). This suggests that N availability was better represented by the spatially averaged value in this study. Porewater DIN in mangrove forests is highly heterogeneous horizontally (Inoue et al., 2011) and vertically (Kristensen et al., 1998; Lee et al., 2008) even in very small scales such as 10 cm. The DIN measured from one soil core sample might not have captured the representative value at each plot due to such heterogeneity. Differences in the predicted AGB between the two cases highlight nutrient availability in affecting mangrove production and biomass dynamics in this forest. Therefore, an appropriate representation of nutrient availability is critical for accurate prediction of mangrove production. More detailed measurement of porewater nutrient concentrations in space and time is needed for a more reliable model prediction, and future works will account for this aspect. Similarly, future works should consider biogeochemical processes which control nutrient dynamics in the substrate. For example, the porewater of the Fukido mangrove forest is rich in ammonia compared to nitrate (Table S1), contrary to the groundwater flowing into this forest, which is rich in nitrate (Mori et al., unpublished data). This suggests that biogeochemical processes, such as mineralization of organic matter, N fixation, and denitrification (Reef et al., 2010), are important drivers controlling nutrient dynamics in the forest, which ultimately affects soil organic matter dynamics. These factors should therefore be taken into consideration in future works as one of the plant-to-soil feedbacks in addition to water uptake processes.

5 Concluding remarks

This paper presents a new individual-based model modified from SEIB-DGVM for a better physiological representation of mangrove growth under the impact of soil salinity. The plant hydraulics was incorporated and linked with the plant production process (C and N uptake) and biomass allocation. The developed model showed high reproducibility of the complex nonlinear patterns in species composition and forest structural variables in a subtropical mangrove forest shaped across a soil salinity gradient without empirical parameterizations of soil salinity influence on mangrove productivity. While there are still some important processes to be accounted for to further improve the model (e.g., plant-to-soil feedback and soil biogeochemical processes), the physiologically improved model predicted the various key ecological processes such as seasonal dynamics in photosynthesis and transpiration, interspecific competition, and self-thinning process, together with forest structure. Thus, including plant hydraulic traits that incorporate species differences in the ability to deal with salinity is critical and adequate for predicting the dominant dynamics in mangrove forests. Although the model has been tested using only two species in one site, owing to its physiological principles that do not hold empirical expressions of influences of environmental variables on mangrove productivity, it can be potentially extended to other mangrove species in various environmental settings. Therefore, it may contribute to predicting how the mangrove biomass dynamics will respond to future changes in global climate.