Modeling mangrove responses to multi-decadal climate change and anthropogenic impacts using a long-term time series of satellite imagery

Determining the effects of simultaneously operating climate change and anthropogenic impacts on mangroves remains a challenge for the Persian Gulf and the Gulf of Oman. The objective of this study was to spatially quantify the spatial extents and biomass of dwarf and tall mangroves in response to drought intensities (Standardized Precipitation Index, SPI), sea-level rise (EC), land use/land cover change (LULC), and surface runoff and freshwater inflow to upstream mangrove catchments along the coasts of southern Iran over a 35-year period (1986–2020). The study drew on a long-term time series of 105 Landsat satellite images, maritime data, rainfall data, field surveys, and models to quantify independent variables (i.e., the WetSpass-M model to quantify surface runoff) and mixed models analysis to determine the important drivers of change. Although mixed model analysis indicated that sea-level rise was the main climate change driver of the development of spatial extents and biomass of both tall and dwarf mangroves, its influence was embedded in the larger temporal climate context of the region. Our findings show that spatial extents and biomass development are neither tightly coupled nor developed temporally or directionally synchronously in dwarf and tall mangroves. This study also shows a precipitous decline in rainfall amounts, SPI, and freshwater runoff volumes starting in 1998. The ensuing longterm drought that is still ongoing but decreased slightly in intensity in the last few years shaped the overall spatiotemporal development pattern of the mangrove structures over the entire study period. The correlation between biomass and the independent variables was similar for dwarf and tall mangroves and was positive with SPI and runoff amounts and negative with EC in all three study sites.


Introduction
Climate change and land use/land cover changes (LULC) can disturb hydrological and sedimentary regimes that largely control mangrove structure and biomass (Jia et al., 2018). Climate change, reflected in increased air temperatures, and intensities of storm events as well as altered patterns in rainfall, drought occurrences, and ocean circulation, has led to a decline in the area, structure, production potential, health, survival, and functioning of mangrove ecosystems around the world (Gilman et al., 2008;Mafi-Gholami et al., 2020a, b). Rising sea levels associated with climate change pose a particularly grave hazard to the long-term survival of mangroves if they cannot expand landward (Abdollahi et al., 2017;Mafi-Gholami et al., 2019, 2020a. LULC changes can affect the timing and balance of monthly and seasonal surface and subsurface flow regimes and volumes of runoff entering coastal areas, amplifying climate change effects on mangroves (Eslami-Andargoli et al., 2010). Synergies among simultaneously occurring multiple stressors and disturbances (Gilman et al., 2008;Lewis et al., 2011) can severely decrease the viability of mangrove seedlings and the structure, growth, and net production of mature mangroves, resulting in diminished spatial extents and ecological functions (Eslami-Andargoli et al., 2010;Zhang et al., 2012;Hutchison et al., 2014;Asbridge et al., 2016;Godoy et al., 2018;Mafi-Gholami et al., 2020b).
Quantifying the effects of reduced runoff volumes on the structure and production of mangroves requires detailed long-term regional records of river discharge flows into coastal environments, which often do not exist in much of the developing and underdeveloped world (Eslami-Andargoli et al., 2010;Asbridge et al., 2016;Hu et al., 2018;Gao et al., 2019). To derive a time series of surface runoff from long-term rainfall data and LULC changes in upstream mangrove catchments, remote sensing may be a promising approach, offering comprehensive regional coverage and high temporal resolution for monitoring both mangroves and climate change over a long time (Hu et al., 2018;Pirasteh et al., 2020;Finger Dennis et al., 2021;Wang et al., 2021;Thomas et al., 2021). Satellite and aerial imagery are more accurate, cost-effective, and efficient than in-situ measurements for detecting changes in the extent and biomass of mangroves and LULC changes that occur at large scales and over long periods of time (Zhang et al., 2017;Bihamta et al., 2020;Lucas et al., 2020;Lymburner et al., 2020;De Jong et al., 2021;Finger Dennis et al., 2021).
In this study, we used long-term rainfall and sea-level data and longterm satellite imagery combined with extensive field surveys to prepare a 35-year time series (1986-2020) of drought intensity, sea level rise, surface runoff values, and LULC of upstream catchments to model the development of large-scale spatial extents and biomass of some of the remaining dwarf and tall mangroves in the Middle East. Previous studies documented a significant decline in rainfall amounts in 1998, switching from wet to dry conditions in the region, and average annual rates of sealevel rise of 3.4-4.8 mm yr − 1 on the Persian Gulf and 10.8 mm yr − 1 on the Gulf of Oman with further increases in sea-levels and temperatures and decreases in rainfall forecast over the coming decades (Mafi-Gholami et al., 2017;Mafi-Gholami et al., 2020a, b). It is unclear, however, to what extent recent LULC changes due upstream development of infrastructures and agricultural lands have further reduced freshwater volume inputs to the coasts, and thus exacerbate adverse climate change effects on the remaining mangroves (Etemadi et al., 2016).
In support to address adaptation planning for climate change impacts on mangroves worldwide (Ellison, 2015), we present a novel, largescale, spatially explicit approach for simultaneously quantifying climate change and anthropogenic impacts on a regional scale that can be applied in developed and developing countries. The specific objective of this study was to quantify the relationship between the spatial extents and biomass of dwarf and tall Avicenna marina mangroves to rising sea levels, drought occurrences and surface runoff volumes in the Persian Gulf and the Gulf of Oman region over the past 35 years. We also quantified to what extent rainfall amounts and human-induced LULC changes in upstream catchments have impacted surface runoff volumes. Because pure mangroves form dwarf and tall structures that occupy different locations along the land-sea continuum, we hypothesized that the combined climatic and anthropogenic effects differentially affected the two mangrove structures.

Study area
The study area includes three pure mangroves sites on the northern coasts of the Persian Gulf from 25 • 34 ′ 13 ′′ N to 27 • 10 ′ 54 ′′ N and the Gulf of Oman from 58 • 34 ′ 07 ′′ E to 55 • 22 ′ 06 ′′ E, encompassing a total area of 10025.55 ha (Fig. 1). Long-term climatic data of the study area show an average annual rainfall of 90 mm, an average annual temperature of 27.2 • C with minimum and maximum temperatures of 8.9 • C to 22.4 • C (January) and 30.1 • C to 48 • C (July), and an average annual relative humidity of 35%.

Methodology framework and data sources
This study combined different data sources to relate mangrove responses to climate change and anthropogenic impacts over the most recent 35-year (1986-2020) period (Supplemental Fig. 1). Field data sampled in 2017 provided the basis for estimating mangrove biomass that was regressed against NDVI values derived from Landsat images. Estimating the biomass change point that separates sample plots of dwarf from tall mangroves and applying the corresponding NDVI threshold to the time series of Landsat images produced maps depicting spatial extents and biomass of both mangrove structures. The accuracies of the constructed maps were evaluated and checked against aerial photographs.

Elevation capital (EC), drought intensity and surface runoff in each catchment area
Sea-level rise is related to the EC and can be computed from (Eq. (1)) where R sl is the absolute sea-level rise rate (3.3 mm yr − 1 ), R sub is the subsidence/uplift rate (-4 Annual drought intensity over time was quantified using September rainfall amounts recorded in 31 synoptic rainfall stations distributed throughout the study sites (Iranian Meteorological Organization, IRMO). Drought intensity is typically expressed by the Standardized Precipitation Index (SPI) that reflects long-term rainfall amounts, drought occurrences, and their impacts on natural ecosystems (Wang et al., 2019). An interpolated spatially explicit SPI raster map was produced for each upstream catchment area applying the inverse distance weighting (IDW) model in the ArcGIS software (version 10.7.1) to annual SPI values of all stations.
Due to the limited availability (7 years) of river discharge data into upstream catchments, no direct long-term estimate of freshwater volumes entering the coastal areas could be obtained. Instead, the amount of annual total surface runoff (m 3 / year) in each upstream catchment was computed as Eq. (2) (Abdollahi et al., 2017) using the WetSpass-M model (Hydrological and Hydraulic Engineering Organization, University of Brussels): where ∑ is the sum over 12 months, S rm is monthly surface runoff (mm), C sr is an actual runoff coefficient, P m is monthly rainfall (mm), I m is monthly interception (mm), and C h is a dimensionless coefficient that represents soil moisture conditions scaled between 0 and 1.
Computing the components of Eq. (2) relied on a time series of monthly rainfall data and values of various coefficients that were adopted from previous studies (De Groen and Savenije, 2006;Pistocchi et al., 2008;FRWMO, 2017;IWRMC, 2019;IRIMO, 2019;LP DAAC, 2020) and then optimized for the study area during the model calibration stage (see supplement for more details). LULC changes were quantified as changes in the runoff coefficient of each upstream catchment were based on satellite images.
The Z-statistic of the Mann-Kendall (MK) test that is embedded in the MAKESENSE 1.0 software (Zhang et al., 2012) was used to determine the change point year in the in the 35-year time series of SPI and runoff values in each study site. The MK test statistic is calculated from Eq. (3): and n is the sample size.
The statistic S is approximately normally distributed when n ≥ 8, with mean (E) and variance (V) (Eq. (4) and Eq. (5)): where t i is the number of ties of the event i. The standardized statistic (Z) for a one-tailed test is (Eq. (6)): At the 5% significance level, the null hypothesis of no trend is rejected if |Z| > 1.96. Positive values of Z indicate increasing and negative values denote a decreasing trend.

Field surveys and estimation of amounts of mangrove biomass
Field surveys were done in 2017 in 385 sample plots (190 in Khamir,110 in Tiab, 85 in Jask) established using the Global Positioning System (GPS) at a distance of 150 m from one another along 10 (Khamir) and 6 linear transects (Tiab and Jask) that extended from the landward margins (location of dwarf mangroves) to the seaward edges (location of tall mangroves) (Supplemental Fig. 2). Plots sizes were 30 × 30 m size and corresponded to the 30 m spatial resolution of Landsat images. In each plot, tree diameters were recorded by mangrove structure.
Without known allometric relationships of above-(AGB) and belowground biomass (BGB) of mangroves to DBH for this study area, we used the allometric equations of Comley and McGuinness (2005) (Eqs. (7) and (8)): where AGB is the above-ground biomass, BGB is the below-ground biomass, and DBH is the diameter of each mangrove stem. AGB and BGB values of all stems in each plot were summed to arrive at plot-level biomass estimates.

Separating dwarf and tall mangrove structures
Examining dynamics of dwarf and tall mangroves required spatial separation and the determination of exact boundaries between structures, which is imprecise when using simple classification of vegetation index maps or through on-screen separation of Landsat images with medium resolution. We thus applied a two-step process that first led to the determination of the change point, i.e., the sample plot with the greatest probability of all points in the biomass frequency distribution (Lubes-Niel et al., 1998) to separate tall from dwarf mangroves. Change points were identified using the non-parametric Pettitt-Mann-Whitney test (α = 0.05) and the Cumulative sum (CUSUM) method in the Change Point Analyzer (CPA) software and confirmed by a t-test of the difference between the mean biomass values of tall and dwarf mangroves at the identified change point.
The second step involved finding the threshold NDVI value derived from Landsat images of the sample plots that corresponded to the biomass change point. Three cloud-free Landsat images (paths/row # 158/042, 159/041, 160/041) that coincided with the field survey dates in August and September 2017 and corresponded spatially to each sample plot were selected to extract NDVI values that were then related to the estimated biomass values through linear regression. The regression model was based on 70% of the sample plots of each study area; the remining 30% were used to validate the model. Following the preparation of NDVI maps and computation of the NDVI threshold value, study sites were classified into dwarf and tall mangrove areas and their spatial extents and average biomass values were computed for the year 2017.

Constructing the time series of spatial extent and biomass of tall and dwarf mangroves
The 35-year time series was constructed using 105 Landsat images (paths/row # 158/042, 159/041 and 160/041) that were acquired from the USGS Earth Explorer (Earth Explorer, 2019) portal and taken in summer months between 1986 and 2020 to avoid seasonal phenological differences in mangroves. Because dwarf mangroves (<1 m height) are often inundated during high tide, Landsat images taken at low tide were selected with the aid of daily recorded wave height data from 22 tide gauges located along the coasts to determine the spatial extent of dwarf mangroves and the location of the seaward edges of the tall mangroves (Mafi-Gholami et al., 2019;Gao et al., 2019). Following the required geometric and radiometric image corrections and using the NDVI threshold value, annual NDVI maps with boundaries around dwarf and tall mangrove areas as well as maps of LULC of upstream catchments were generated in the ENVI software and the biomass and spatial extents of both mangrove structures were computed for the entire time series for each study site.

Validation of the maps
The accuracy of the constructed maps for years 2016 and 2017 was assessed by field sampling of 190 plots during the summer of 2017. In addition, 200 plots (i.e., 90 in Khamir, 65 in Tiab, and 45 in Jask) were established in 2020 to assess the accuracy of the 2018-2020 maps. The accuracy of the 1986-2016 maps was validated by checking whether control points indicated the presence of dwarf or tall structures using texture, pattern, and lightness/darkness information derived from aerial photographs and Quickbird images of the years 1985, 1992, 1995, 2001, 2004, 2009, 2012 and 2014 (National Geographical Organization of

Fig. 2.
One-year SPI values (calculated using monthly rainfall data) over the 35-year period (1986-2020) in the upstream catchment area. Iran). Stratified random sampling was used to assess user accuracy, producer accuracy and overall accuracy of the produced maps (Eslami-Andargoli et al., 2010), which determined a consistent accuracy of > 95% of the classified images.

Relationship of mangroves to the independent variables
Pearson correlation coefficients of EC, SPI, surface runoff, and area and biomass of dwarf and tall mangroves were computed for each study site separately for the wet (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997) and dry (1998-2020) periods. Repeated measurement mixed modelling with a first-order autoregressive temporal correlation structure was applied to examine the statistical significance of the fixed effects of EC, SPI and runoff values and their two-and three-way interactions on the responses by study site for each period. Statistical analyses were performed in SAS 9.4 (SAS Inc., Carey NC) (https://www.sas.com/en_us/company-information/profile. html).

Drought occurrences in the upstream catchments between 1986 and 2020
Average annual SPI values showed that the region experienced a wet period (positive SPI values) between 1986 and 1997 followed by a dry period (negative SPI values) that started in 1998 and was most severe in 2006. Droughts lessened slightly thereafter in the region, particularly on the Gulf of Oman, but continued into 2020 (Fig. 2).
Computing runoff volumes over time with observed annual rainfall values and LULC values held constant at 1986 baseline values indicated only small contributions of LULC changes to runoff volume decreases (e. g., 1.1% in Kamir, 1.3% in Tiab, and 1.6% in Jask in 2020).
The NDVI threshold value of 0.57 had the highest probability (>0.98 in all transects) separating dwarf and tall mangroves (Table 1, Supplemental Fig. 4).

Development of spatial extents and biomass of dwarf and tall mangroves between 1986 and 2020
Spatial extents of mangroves developed uniquely over time in each study site and differed by mangrove structure (Fig. 4). The greatest area changes occurred during the wet period for both structures in Khamir and Tiab and for tall mangroves in Jask. Areas of dwarf mangroves increased in Tiab (268 to 830 ha, +46.82 ha yr − 1 ) and decreased in Table 1 Changes of mean total biomass in NDVI values higher and lower than 0.57 in the different study habitats along the northern coasts of the Persian Gulf and the Gulf of Oman.  Khamir (5757 to 4287 ha, -122.55 ha yr − 1 ) and Jask (346 to 312 ha, -2.82 ha yr − 1 ) while those of tall mangroves increased strongly in all three study areas (Khamir: 1485 to 3423 ha, +161.6 ha yr − 1 ; Tiab: 88 to 338 ha, +20.85 ha yr − 1 ; Jask: 109 to 257 ha, +12.34 ha yr − 1 ).
The development of biomass of tall (but not dwarf) mangroves followed different patterns in the three study regions during the wet period but followed a more synchronized decline (both structures) during the dry period (Supplemental Fig. 8).

Correlation structure among EC, SPI, and runoff
EC values increased by 4.8 mm yr − 1 (Khamir), 3.4 mm yr − 1 (Tiab), and 10.8 mm yr − 1 (Jask), with cumulative increases between 1986 and 2020 of 16.4 cm in Khamir (from 135.1 cm to 151.5 cm), 11.2 cm in Tiab (from 104.5 cm to 116 cm), and 36.7 cm (from 26.5 cm to 63.2 cm). The correlation structure was highly variable before and after the change year, with differences among sites and between the two periods (Supplemental Table 1). Of the few statistically significant correlations, most changed the sign/direction before and after the change year. The correlation pattern in both periods was the same in all three sites (with the exception of Khamir in 1986-1997) only for EC and surface runoff, which also had the greatest correlation coefficient by period in the three sites.

Correlation structure between spatial extent and biomass of dwarf and tall mangroves
The correlation structure differed between the pre and post-changeyear periods and was unique to each study site (Supplemental Table 2). During the wet period, both areas and biomass of dwarf and tall mangroves were significantly negatively correlated in Khamir and positively correlated in Tiab. While dwarf and tall mangroves were significantly negatively correlated, biomass was positively correlated in Jask. During the dry period, area and biomass of dwarf mangroves were significantly negatively correlated in Jask, whereas the areas of dwarf and tall mangroves were significantly negatively correlated in Tiab. The biomass of dwarf and tall mangroves was significantly positively correlated in Khamir and Tiab. The area and biomass of tall mangroves was significantly positively correlated in Khamir.

Relationship of EC, SPI, and runoff with spatial extents and biomass of dwarf and tall mangroves
The correlation of the independent and response variables differed among the three study sites and between the pre-and post-change-year periods in each site (Supplemental Table 3). Consistently significant correlations occurred between EC and all four response variables (negative) in Khamir during the dry period and for all three independent variables with each dependent variable in Tiab during the wet period. Runoff was consistent (but not always significantly) and positively correlated with the area of tall mangroves during both periods in all three study sites. Overall, greater correlation coefficients (positive or negative) with all response variables except biomass of tall mangroves were observed for EC and runoff than for SPI.
The best-fitting repeated measures mixed models for areas and biomass of dwarf and tall mangroves differed between the pre-and postchange-year time periods and among the three study sites (Supplemental Table 4). All models indicated a strong linear or curvilinear influence of sea-level rise (EC) on the responses of dwarf and tall mangroves that was only slightly moderated or accentuated by SPI and runoff during both periods. During the wet period, EC was positively related to the area of dwarf mangroves in Tiab, the biomass of dwarf and the areas of tall mangroves in all three sites, and the biomass of tall mangroves in Tiab and Jask; increasing rainfall (SPI) and runoff typically, but not always, led to an increase in area and biomass of both mangrove structures. During the dry period, increases in EC were negatively related to the areas and biomass of both mangrove structures except for the area of dwarf mangroves in Tiab and the biomass of tall mangroves in Jask that both showed a U-shaped parabolic shape. Similar to the wet period, increasing rainfall (SPI) and runoff typically, but not always, led to an increase in area and biomass of both mangrove structures.

Table 2
Least squares regression (LSR) results using 70% of the data. Modeling the ABG and BGB of mangroves in tall and dwarf mangroves in different habitats using NDVI as a predictor.

Discussion
The use of satellites images was instrumental for elucidating the regional pattern of spatial extents and biomass responses of dwarf and tall mangroves to climate change that developed dynamically and asynchronously in the three study sites over the past 35 years. Consistent with global results (Ellison, 2015), sea level rise was the main climate change driver responsible for the declining areas and biomass of both mangrove structures in this study region. In this study, the effect of sea level rise was embedded in the local and regional temporal climate context of a precipitous decline in rainfall amounts, SPI, and freshwater runoff volumes in 1998 that ushered in a still-ongoing long-term drought. Exposure to oil pollution following major oil spills into coastal waters (ICZM, 2017) can further stress mangroves (Lewis et al., 2016) as may local infrastructure projects (Godoy et al., 2018) such as dam construction or the development of a water supply network to farmlands in upstream catchments (e.g., 2002 in Jask) that can drastically reduce river discharges and sediment volumes (Mafi-Gholami et al., 2017).
Although sea level rises over of the past 35 years may not have been linear as assumed in this study, increasing areas of tall mangroves during the wet period despite rising sea levels (Fig. 5) that caused some modest retreat at seaward margins (Mafi-Gholami et al. 2020a) highlight the importance of rainfall for mangrove dynamics. Abundant rainfall during the wet period appears to have lessened local sea level rises and the magnitude of inundation at seaward margins through increased runoff and sediment loads transported from upstream catchments to coastal areas that can raise mangrove beds (Lovelock et al., 2017;Woodroffe et al., 2016). In addition, the conversion of dwarf into tall mangroves at their intersection due to improved conditions for height growth and biomass production and the migration, establishment, and expansion of dwarf mangroves into upland saltmarshes, appear to have contributed to a net expansion of mangrove area.
Immediate and clearly detectable biomass reductions of both structures with the onset of dry conditions (Fig. 6), indicate that freshwater deficiency rapidly enhances stress, causing immediate growth reductions and mortality of sensitive seedlings and saplings, including a reduction of the competitiveness of dwarf mangroves versus saltmarshes that can halt or even reverse further upland incursions (Lovelock et al., 2017;Gilman et al., 2008;Eslami-Andargoli et al., 2010;Mafi-Gholami et al., 2017;Osland et al., 2017). During droughts, local sea level rise is enhanced due to lower runoff volumes, sediment loads, and sedimentation rates that can increase the potential for erosion (Ellison, 2008), influence temporal and spatial changes of mangrove boundaries/ shorelines and reduce mangrove areas and production potentials (Lovelock et al., 2017;Mafi-Gholami et al., 2020a). Although tall mangroves typically have lower sensitivity and vulnerability to environmental stresses and disturbances than dwarf mangroves (Ellison, 2015;Gao et al., 2019), intensified erosion and regression of seaward tall mangrove margins may also be influenced by the large increases in human activities and numbers of fishing ports and vessels in recent years that can lower sedimentation rates (Ellison, 2008;Mafi-Gholami et al., 2020b).
After accounting for the effects of rising sea levels, the effects of both SPI and surface runoff on spatial extents and biomass of both mangrove structures were mostly positive within each period. However, the lack of statistical significance of SPI and runoff during the wet period indicates that area and biomass dynamics were not tightly coupled to rainfall when SPI values were positive. In contrast, the statistical significance of either SPI or runoff during the dry period indicates slightly greater biomass under less droughty conditions. Compared to the effects of rising sea levels (EC), however, the effect sizes of SPI and runoff were small and those of anthropogenic LULC changes on runoff or mangroves were almost negligible.

Conclusions
This study extends previous results (e.g., Rovai et al., 2016;Osland et al. al., 2017;Gabler et al., 2017;Navarro et al., 2020;Lucas et al., 2020;Lymburner et al., 2020) and documents the benefits of a more detailed investigation of the temporal development of each mangrove structure. Asynchronous spatial dynamics of dwarf and tall mangroves within the region underscore the need for spatially explicit, local climate change adaptation planning that is embedded into the larger regional scale. The loss of biomass of tall mangroves and the loss of area of dwarf mangroves should be seen as an early warning system of a threat from drought, coastal erosion, inundation or competition from saltmarshes. The decline of tall mangrove biomass in Khamir and Jask that was out of synch with rainfall and available freshwater runoff during the wet period points to the other environmental factors (e.g., geomorphological characteristics, local tidal ranges, the degree of water salinity) or anthropogenic influences (e.g., oil pollution, mangrove exploitation, LULC changes prior to 1986) that may have heightened their vulnerability in this region. Sufficient rainfall that results in positive SPI values and greater surface runoff appears to increase growth and biomass and encourage migration of mangroves into adjacent upland saltmarshes that compensates or even reverse area losses due to inundation by sea water. Because accelerating sea-level rise and predicted increases in temperatures and drought intensities in the region in the coming decades will likely subject mangroves to further stresses, future studies should investigate the differences in sensitivity and vulnerability of dwarf and tall mangroves to the various aspects of climate change and how this may inform climate adaptation strategies. High-resolution hyperspectral optical and radar data could enhance the precision of separating tall and dwarf mangroves for more reliable assessments of the impact of climate change in future studies.

Authors statement
We declare that we have no known competing financial interests or personal relationships that could have influenced the work reported in this paper. This paper has not been submitted or published previously.

Declaration of Competing Interest
The authors appreciate the GeoAI Smarter Map and LiDAR Lab of the Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University (SWJTU), for the technical research support. This work is the outcome of one of the joint research studies with international University collaborators and has not been published or submitted previously.
The authors declare no conflict of interest. This research received funds to publish this research outcome from the start-up funds with Project Number: A1920502051907-6, from the Faculty of Geosciences and Environmental Engineering (FGEE), Southwest Jiaotong University (SWJTU), China. Data and code are available upon request.