The mathematical growth model elaborated by Salmo and Juanico (2015) considered the following species-specific factors: α for the biological growth parameter, β for the biological leaf area index, the Malthusian growth rate Ω , and maximum diameter at breast height (DBH) D _max . This model is implemented numerically with Equation (1) by setting Δt equal to one day with the factors: D _t ∈ℝ for the plant’s DBH at time t with D _t ∈[0.5,15) cm and Ω,α,β,D _max ∈ℝ with α≠−2 . The functions σ(x,y) , η(x,y) , and K _t (x,y) represent the stressor responses to salinity, inundation, and time-dependent competition, respectively. The input (x,y) is the location of the plant within an area of size L×L . The types of landscapes and stressors associated with this area are described in the next section.

The equation is based on the Malthusian (exponential) growth model and modulated by the stressor responses. The landscape for the σ , η , and the K competition field was changed.

Figure 1A illustrates the landscape being modeled with slope parameter sl=0.52 , while Figure 1B displays the corresponding inundation stress field of this landscape patch. The landscape was a L×L patch with the diagonal representing the coastline. By default, this diagonal is the mean sea level (MSL). Both the functions σ(x,y) and η(x,y) have the input parameter sl for the slope of the terrain. The slope describes the reach of the mean high tide (MHT). For this setup, the MHT can reach approximately halfway inland as shown by the orange dashes in Figure 1A . Since none of the trees present are within reach of the MHT, the corresponding inundation stress at their location is effectively zero, as indicated in Figure 1B . Consider a landscape with sl=0 as illustrated in Figure 1C where the MHT reaches the left corner of the patch at (0,0) , while Figure 1D shows the inundation stress field as an increasing linear function starting with a value of 0 from the origin (0,0) to 0.5 at the MSL line (diagonal). For this case, the whole patch is submerged during high tide.

Figures 2A, B elaborate further what sl=0 means in terms of inclination. Let L=40.0 meters, then side b=28.3 meters. If the difference between the MSL and MHT is 1 meter, then the inclination of the landscape, which is the angle between c and b in Figure 2B is about 2.05 degrees. Hence, the landscape modeled by sl=0 has the terrain underwater by a depth of 1 meter at the coordinate (40,40) . Consider the case represented by Figures 2C, D . If d=14.1 such that the MHT reaches halfway inland, then the inclination is 4.05 degrees. A case where the MHT fails to penetrate inland is instantiated by Figure 3 , resembling a coastal cliff with slope value sl≈1 .

The models of the salinity field S(x,y) and inundation field I(x,y) are embedded on the landscape. Stressor fields ( Supplementary Material ) are treated as piecewise functions of the slope parameter sl . The salinity field consists of values between 0 and 72. This field then influences the salinity response according to the following:

For a given slope, a plant located at (x,y) will encounter inundation stress according to another piecewise model. The inundation response is determined by the following:

For the individual plant competition ( Supplementary Material ), the model used was the FON approach, which sums up the single field intensities of neighboring trees into an aggregate field strength F(x,y) (Berger and Hildenbrandt, 2000). In FON, it is assumed that individual growth is impossible if the quantity referred as “strength of neighborhood” F _A exceeds 0.5. Therefore, the growth of the plant is stopped when F _A >0.5 , i.e., when the function C(F _A )=1−2F _A is less than zero.

The change in the number of individual plants at different life stages is described by a stochastic compartmental model ( Supplementary Material ). The demographic events consist of recruitment, mortality (Schaal and Leverich, 1982) and growth (Fulton, 1993), which may be expressed as state-transition equations.

• Recruitment: Tree →           r             Tree + Seedling

• Seedling death: Seedling →         m seed             Dead

• Sapling death: Sapling   →         m sap             Dead

• Tree death: Tree →         m tree             Dead

• Seedling to sapling growth: Seedling →           g seed       Sapling

• Sapling to tree growth: Sapling →       g sap     Tree

The growth rates g _seed and g _sap of seedling and sapling, respectively, can be estimated by the DBH growth rate across the defined size at the transition: 2.5 cm and 5.0 cm. In this manner, the practical measurement of the growth rates can be made from field surveys.

The state-transition equations are implemented with a stochastic simulation algorithm formulated by Gillespie (1976). Recruitment was established as seeds take root outside the crown of any existing tree. The dispersal rate multiplied by the time for the next event gives the maximum distance a seed can travel from the parent tree (given that seeds float in seawater).

With linear stability analysis ( Supplementary Material ), the average, long-term dynamical behavior of the stochastic model can be described with the expansion method by Van Kampen (1992). From this average (sometimes referred as “mean-field”) dynamics, a system of ordinary differential equations can be examined further for its bifurcation properties. A dimensionless parameter ξ , reminiscent of the basic reproduction number in epidemic models, can be derived from the transition rates of the stochastic model. The result is the following index:

By determining the value of the transition rates over one or two years from a trial plantation in a particular site, then ξ could be estimated. The restoration index (Equation 4) has a critical value equal to one, which is the value that separates the average, long-term behavior of the model into two: (1) endemic equilibrium if ξ>1; and (2) extinction equilibrium if ξ<1. The endemic equilibrium is the state described by the survival of the plant population long after it was established. The extinction equilibrium is the opposite state in which the population died out eventually. The endemic equilibrium is the preferred outcome of restoration.

The sensitivity analysis aims to get a good measure of how a time series diverges given an elementary effect. Particularly, the sensitivity of the following outcomes were analyzed:

The behavior of the stochastic model can be affected by the numerical value of the factors given in Equation 1 (namely, α, β, Ω ) and the slope sl , which are implicit in Equation 2 and 3. Although considered as constants, U _i and ds in Equation 3 are likewise included in the sensitivity analysis. The ranges of the parameter values are:

The above sample ranges were restricted to values within reasonable exploration. For example, Ω>0.25 would model mangroves capable of maturing to tree status in less than three months. Using the framework of elementary effects analysis ( Supplementary Material ), the partial difference with respect to a single factor of the model must be averaged (Morris, 1991). This averaging yields an effect with mean μ (taken as an absolute value) and standard deviation σ . A high μ suggests that the concerned factor generally shifts the outcome of the model by a large degree. A low σ implies that the concerned factor indeed affects the outcome of the model.

The effect of coastal mangroves to the inland propagation of storm surge was simulated in quasi-3D, which builds on earlier quasi-2D simulation results for wave-vegetation dynamics (Zhang et al., 2012). The simulations involved a combination of several models of the relevant physics and biology of the coastal vegetation system ( Figure 4 ). The physical part accounted for the water flow and the spatial profile of the sea floor (bathymetry) and land (topography). The biological aspect accounted for the fully developed mangrove along the coastline, representing its interaction with incoming sea waves by the tree density, average height, and average DBH (Suzuki et al., 2012). This parametrization effectively interfaced the result of the mangrove restoration simulation with the coastal wave model using XBeach version 1.22.4714:4905M (Roelvink et al., 2015). The model considered a simplified three-section vertical structure for the individual Rhizophora tree with a crown, trunk, and prop roots ( Supplementary Figure 1 ). The tree density (per ha) from the restoration simulations were taken as vegetation input for the coastal wave model. For Tacloban, Leyte, due to the thin existing mangroves, the assumed tree density was 125/ha. For another site in Iloilo, the assumed tree density was 1,250/ha based on field observations. Both assumed tree densities were lower than the values projected from the restoration simulations. This underestimate was intended to buffer for any overestimation arising from the restoration model.

To describe the storm surge development, propagation and impact, three factors of the wave dynamics of Haiyan ( Supplementary Material ) were considered: wave conditions, including wind direction; bathymetry and topography, including the coastline orientation; and demographic landscape. The conditions of Haiyan were simulated from the available wind data. The bathymetry and topography were reconstructed from available NASA SRTM dataset (jpl.nasa.gov). The model bathymetry was generated using the following steps with Delft3D. First, a grid was defined using spherical coordinates spanning the area of interest. Second, with RGFGRID the grid was orthogonalized and refined to obtain a median pixel size of about 40m×40m . Then, with QUICKIN, a triangular interpolation was used to generate the model bathymetry from SRTM15+ data with an intrinsic resolution of 500m×500m at the equator (Tozer et al., 2019). The SRTM15+ dataset has also been utilized for inundation models with reasonably satisfactory results in tsunami simulations (Serra et al., 2021; Qiu et al., 2022) and sea-level rise (Wang and Marsooli, 2021). Lastly, the population census statistics during the period of the Haiyan disaster in Tacloban, Leyte were used to estimate the spatial density of human settlements in the urban area.

Due to the dearth of mangroves in Tacloban, Leyte (Carlos et al., 2015), a different location with a dense coastal mangrove and along Haiyan’s path was sought. A site in the northeastern Iloilo, Pan de Azucar Island (also known locally as Tambaliza) that belonged to a town called Concepcion, fits the description ( Figure 5 ). The mangrove of the island was mapped and encoded in the bio-shield simulations. Particularly, the estimate of the tree density, estimated at around 1,250/ha, was a crucial input to the simulations. The simulation results for Pan de Azucar was used as a baseline to verify if endemic mangroves exerted any appreciable coastal protection effect against the category-5 disturbance brought by Haiyan.

An effective bio-shield must offer protection to a populated town from storm surge damage, while minimizing the cost of establishing this bio-shield. Tacloban was the town in the Philippines that suffered the most damage by Haiyan in 2013. The economic viability of establishing the bio-shield needed an estimate of the value at risk for Tacloban. A household survey was deployed across the city between December 2016 and February 2017, which gathered data from a convenience sample of 5,000 households scattered across Tacloban’s metropolitan area of which 55.7% responded ( Supplementary Table ). The survey sought to estimate how much economic value was at risk of inundation damage by the Haiyan storm surge with the following index:

Equation 5 tags a value of every member as equivalent to the GDP per capita in Tacloban. This value represented the average income or potential income that the household member would contribute to Tacloban’s economy. The tangible assets included the detachable possessions such as household appliances and other valuable items. HV , thus, estimated the casualty value that Tacloban would lose due to a sudden catastrophic disturbance. The same estimate was applied to Pan de Azucar.

The HV index in Equation 5 did not attempt to place a value on the real property for data privacy reasons. The economic valuation of the bio-shield also did not consider the carbon storage benefits of mangroves. The spatial resolution of the SRTM15+ bathymetry is not the highest available. However, due to computational limitations, this dataset was chosen. For further simplicity, the bio-shield simulations did not include sea level rise due to the melting of the polar ice caps. Studies have shown that a reasonable rate of sea level rise in Philippine coasts is about 15 mm/y, which is about nine times the global average (Holden and Marshall, 2018). At this rate, MSL would recede by about 38 cm in 25 years, which the present study considers small enough to assume a quasi-steady coastline over the course of the simulations. Lastly, the mangrove restoration and bio-shield simulations only assumed Rhizophora mangroves to simplify the estimates of growth rates and drag coefficients. The simulated patch is large enough (L=40 m) to generate an extended mangrove, but sufficiently small to manage the computational memory requirements.

Based on multiple runs of the stochastic simulations, the Malthusian parameter Ω , and salinity response parameters dsalt and U _i exerted the most influence on the long-term tree population sizes ( Figures 6A, B ). Figure 6C also shows that Ω exerted substantial influence on the plant’s time to maturity, which is essential for establishment and survivability. A separate analysis for the time derivative of Equation 1 was made for the factors Ω , α , β , and D∈[0.5,15] cm. The increments for Ω , α , and β were maintained including their respective ranges with the exception of D , where it was scaled by a factor of 1.45. The domain of the random samples from the functions σ , η , and K was [0,1] . With the factors dsalt and U _i omitted, the curvature of σ did not play a role in analyzing the elementary effects of the factors in consideration. Figure 6D shows the result of a million samples taken from the factor space, with Ω also exerting the strongest influence among the factors even though the scale was 10⁻³ . The slope relatively exerted a moderate effect on the long-term outcomes, while the model was least sensitive to α and β suggesting its robustness to species-specific variations (e.g., among different Rhizophora species) and foliage cover.

Due to stochastic dynamics, the result of each simulation is a time series of tree, sapling, and seedling population taken as the average of ten runs. A typical outcome of a restored mangrove after a run time of 25 years consists of a distributed population across a flat coastal landscape ( Figure 7 ). The initial seedling plantation ( Figure 7A ) generates the distribution of mangroves 25 years post-planting ( Figure 7B ). The landscape-rendered distribution is shown in Figure 7C . The youngest plants are mostly situated at the frontline along the MSL, while a few can be found in spaces between mature trees.

The time evolution of the population size based on tree, sapling, and seedling compartments are shown in Figure 8 for the two extreme slope cases, sl=0 (tidal flat) and sl=1 (coastal cliff). Stressor gradient effects are apparent on for sl=0 , in which the maturity is delayed on some batches because of the seaward increase of stress due to inundation and salinity. For the case sl=1 , the majority of initial seedlings reach their sapling and tree stages at the same time in the absence of salinity and inundation stressor gradients. The seedling population dropped quickly but not to zero, as seedlings planted in sites where stressors are high go through stunted growth. This gap is consistent with the gap between the number of seedlings before the drop and sapling population immediately after. The spike in the tree population near year 5 corresponds to a batch situated near the edge of the cliff. The peak tree density (considering the area approximately 0.08 ha) are 9,500/ha and 11,900/ha on the tidal flat and coastal cliff, respectively. At 25 years post-planting, the tree densities are about the same, at around 3,500/ha. The consistent decline in tree density after an initial peak indicates population thinning due to competition. It involves the DBH growth facilitated by the opening of space due to the death of some trees.

An extended simulation of a particular restoration effort ( Figure 9A ) reveals that the average tree density ( Figure 9B ) tends to increase in the long term after it peaks then declines in the medium term ( Figure 9C ). This average result appears to be guaranteed if ξ>1 in the first few years post-planting ( Figure 9D ). Furthermore, moving the restoration zone inland can improve long-term outcomes as can be illustrated with a time series ( Supplementary Figure 2 ) or a phase portrait ( Supplementary Figure 3 ). The restoration index, like the basic reproduction number in epidemic models, provides a short-term gauge of the effort’s long-term success in obtaining dense mangroves. Considering coastal protection, it is necessary that the coastal mangroves have a high tree density over the long term to serve as bio-shield against storm surges with high, but less common, intensity.

The coastal inundation results for Tacloban corroborated with field measurements conducted about two weeks after the storm (Lagmay et al., 2015). The water level results from the simulations of the actual scenario for Tacloban in the San Jose peninsula and downtown area ( Figure 10 ) agreed with field data. The water level due to the surge was between 4-5 m in San Jose, whereas it ranged from 5 to 6 m in the downtown area (Lagmay et al., 2015). The simulation results at the peak of the surge ( Figure 11A ) indicated water levels within similar ranges in San Jose and downtown. This corroboration was sufficient to assume that the model represented coastal inundation to an appreciable degree.

The primary purpose of the study was to determine if mangroves could have protected Tacloban, Leyte when Haiyan generated intense storm surges. For this purpose, it is important to consider that the dissipative effect of vegetation assumed in XBeach has been validated from extensive laboratory and field tests (Roelvink et al., 2015). Tacloban, by default, did not have mangroves that Pan de Azucar had during Haiyan. The test scenario for Tacloban, thus, considered hypothetical mangroves. The second question was in which parts of the coast would restoration efforts be situated. Given the complexity of interaction between the wind profile and Tacloban’s bathymetry and topography, this question was not trivial to address. The solution was to provide an optimal and economically viable reason for coastal mangrove restoration. The restoration strategy was to minimize the total area occupied by the mangroves, which should minimize the relocation of existing households. The intended outcome was to maximize the area that the storm surge could not inundate, which translated to maximizing the coastal protection capacity.

The hypothetical mangroves, which directly face Haiyan’s wind direction, are located around the San Jose Peninsula toward the southeastern coast ( Figure 10 ). A thin mangrove exists on the northern tip of the San Jose Peninsula. Simulating the Haiyan storm surges over Tacloban in November 2013 ( Figure 11A ), the hypothetical mangroves do not completely shield Tacloban totally from the inundation ( Figure 11B ). The incoming storm surge could still wash over San Jose peninsula, although the mangroves slowed down the transmitted water flow moving toward the southwestern portion of Tacloban ( Supplementary Movie 1 ). Consequently, areas in the southeastern sector of Tacloban are notably drier. This southeastern sector coincides with the location of the San Jose, Caibaan, and Marasbaras villages, which is home to about 5,500 households.

The presence of mangroves in the Pan de Azucar Island of Concepcion, Iloilo enabled testing the protective bio-shield effect ( Figure 12A ) in comparison to a hypothetical absence of mangroves ( Figure 12B ). This comparison is the exact opposite to the one applied to Tacloban—a converse hypothesis test. The existing mangroves were facing an oblique direction to Haiyan’s wind velocity. Simulating Haiyan conditions revealed that the inland inundation would have been wider and deeper in the absence of mangroves, possibly penetrating into Sitio Proper (See and Wilmsen, 2022). The 1,250/ha mangroves in the surveyed area lessened the momentum of the reflected waves that penetrated inland, acting like shock absorbers that limited the energy of water flow ( Supplementary Movie 2 ). This physical effect may also be a consequence of the oblique “angle of attack” of the winds relative to the mangrove orientation, although this aspect may require further confirmation. However, this result would suggest that the bio-shielding impact of mangroves may be site specific. Thus, the protection capacity of mangroves, which ground-truth reports in the area corroborated, were a factor in the relatively lower degree of damage and per-capita casualty rates in Pan de Azucar relative to Tacloban during the Haiyan storm surges in November 2013.

The comparative simulations indicate that mangroves along selected coasts of Tacloban ( Figure 10 ) could lead to patches of dry sites when Haiyan occurred ( Figure 11B ). The dry patches correspond to some of the most populous villages of the city. Based on the 2010 population census, the dry patches would have been the sight of 5,500 households. With the data from the household survey, the estimated value attributable to those dry patches is HV= Php 2.9 billion (5,500 households x 5 members/household x $2,600 GDP per capita multiplied by Php 40.00/US$1.00 (based on the foreign exchange rate for 2013). The total area of the new mangroves ( Figure 10 ) is about 158 ha, determined through the following considerations. Each cell on the simulated patch is approximately a square with about 36-m length on all sides. The total number of planted cells is 1,216, whereas the total area is 158 ha.

Considering Php100,000 as the minimal spend of fully vegetating a cell, which includes the possible average cost of displacing an existing household, then the total cost of planting is only around Php120 million. The tradeoff between the potential benefit of saved value (Php2.90B) and cost of bio-shielding efforts is staggering—a factor of 24x favoring the benefit. Thus, the rational economic decision would be to situate the restoration efforts along selected areas of the Tacloban coastline.