Fire-derived phosphorus fertilization of African tropical forests

Study site

The study was carried out in post-agriculture forests at different growth stages near the forest reserve of Yoko (N00°17′; E25°18′; mean elevation 435 m a.s.l.), situated between 29 and 39 km south east of Kisangani, in the Democratic Republic of the Congo. We set up 15 (40 × 40 m) plots, set out in triplicate along five successional stages (15 plots): agriculture and 5, 12, 20, 60 years old secondary forest (respectively, 5 yrs, 12 yrs, 20 yrs, 60 yrs). Additionally, soils were also characterized in three agricultural plots (Ag). We interviewed owners, farmers, and local experts to determine the time-since-disturbance of all plots. Tree height measurements were recorded at the plot level for 20% of individuals of each diameter class. The climax vegetation in the region is classified as semi-deciduous tropical. Climate falls within the Af-type following the Köppen-Geiger classification³³. Annual rainfall ranges from 1418 to 1915 mm with mean monthly temperatures varying from 23.7 to 26.2 °C. Throughout the year, the region is marked by a long and a short rainy season interrupted by two small dry seasons December–January and June–August. Soils in the region are highly weathered Oxisols, being poor in nutrients, with low pH and dominated by sandy texture.

Sampling and sample analysis

Throughfall and bulk precipitation was collected weekly using polyethylene (PE) funnels supported by a wooden pole of 1.5 m height to which a PE tube was attached and draining into 5 L PE container. A nylon mesh was placed in the neck of the funnel to avoid contamination by large particles. The container was buried in the soil and covered by leaves to avoid the growth of algae and to keep the samples cool. We installed eight throughfall collectors in each plot as two rows of four collectors, with approximately 8 m distance between all collectors. On every sampling occasion, the water volume in each collector was measured in the field, and recipients, funnels and mesh were replaced, rinsed with distilled water. A volume-weighted composite sample of the devices per plot was made. All samples were stored in a freezer immediately and sent in batch to Belgium for chemical analysis. The volume-weighted composite samples were first filtered using a nylon membrane filter of 0.45 µm before freezing. Total phosphorus was measured by inductively coupled plasma atomic emission spectroscopy (ICP AES, IRIS interpid II XSP, Thermo Scientific, USA). Although we acknowledge the potential for microbial activity in the collectors during a 1-week, dark, in situ storage of the samples, the use of total phosphorus concentration and lack of algal growth allow for complete phosphorus recovery.

Following analysis, the samples from the replicate field sites per forest stage were pooled into ‘weekly’ forest-type samples, and these were subsequently analyzed for dissolved black carbon (DBC). In short, the pooled water samples were acidified to pH 2 and analyzed for dissolved organic carbon (DOC) concentration via high-temperature catalytic oxidation on Shimadzu TOC-L total organic carbon analyzer following established methodology³⁴. DOC was isolated from the water samples by solid phase extraction (SPE) following Dittmar et al.³⁵. Briefly, SPE cartridges (Varian Bond Elut PPL, 1 g, 6 mL) were conditioned sequentially with methanol, ultrapure water, and ultrapure water acidified to pH 2 using concentrated HCl, then passed through the SPE cartridges by gravity. SPE cartridges were dried under a stream of high-purity N₂ gas. DOC was eluted from the SPE cartridge with methanol (SPE-DOC) and stored at −20 °C until further analysis. DBC was quantified using the benzenepolycarboxylic acid (BPCA) method as detailed in Wagner et al.²⁰. The BPCA approach to quantifying DBC involves chemothermal oxidation of condensed aromatic DOC compounds to benzenehexacarboxylic acid (B6CA) and benzenepentacarboxylic acid (B5CA) products. The B6CA and B5CA oxidation products are robustly measured and derive exclusively from pyrogenic sources³⁶. Condensed aromatic DBC, as measured using the BPCA method, is ubiquitous in aquatic environments globally^21,37,38,39. DBC has also been quantified in throughfall and stemflow in longleaf pine forests that undergo regular prescribed burning⁴⁰. Therefore, we use the BPCA method as a proxy for carbon inputs from biomass burning in the current study. To analyze our samples for BPCAs, aliquots of SPE-DOC (~0.5 mg C equivalents) were combined with concentrated HNO₃ in flame-sealed glass ampoules and heated to 160 °C for 6 h. The resultant BPCA-containing residue was dried and re-dissolved in mobile phase for subsequent analysis. Individual BPCAs were separated and quantified using an HPLC system (UltiMate 3000, Thermo Fisher, Germany) (CV < 5%). Sample DBC concentrations were calculated using the established power relationship between DBC (µM C) and the sum of B6CA and B5CA (nM-BPCA) using the equation from Stubbins et al.⁴¹: [DBC] = 0.0891 × ([B6CA + B5CA])^0.9175; n = 351, R = 0.998, p < 0.0001. We report the concentrations of B5CA and B6CA with DBC to facilitate comparison with other datasets (Table S2).

Canopy 3D model

UAV surveys were carried out in Yoko on February 9, 2020. The Yoko survey consisted of five flights due to the dispersed locations of the 15 inventory plots and in total covered an area of 302 ha. Two UAV platforms were used in the surveys: (i) a consumer-grade DJI Mavic 2 Pro. This UAV was equipped with a Hasselblad L1D-20c camera (20 megapixels, 5184 × 3456 pixels, ca. 77° FOV). The onboard GNSS supports GPS and GLONASS. (ii) a customized DJI Phantom 3 Advance. We removed the DJI camera-gimbal system and mounted a GoPro Hero 3 camera (12 megapixels, 4000 × 3000 pixels, with 2.92 mm F/2.8 123° HFOV lens) and connected the camera to a RTK/PPK (real-time kinematic and post-processing kinematic) enabled GNSS receiver to determine camera exposure position at centimeter level. In brief, the Mavic camera provides higher image quality but lower GNSS accuracy, while the GoPro camera provides accurate positioning with PPK solution but had lower resolution of images. The combination of both platforms provides high positional accuracy as well as high-quality images⁴². A 90% forward image overlap and 80% side overlap were used during the flights for both UAV/camera systems. The flight height was set at 180 m from the ground level, providing an average ground sampling distance (GSD) of ca. 0.04 m px-1 for Mavic and 0.10 m px-1 for GoPro. During the flights, a Reach RS (Emlid Ltd.) base station was mounted on a tripod placed in an open area within the surveyed region to provide correction data for PPK georeferencing. The coordinate of the base station was determined using the ca. 8 h average value of the single solution throughout the survey. The images of the GoPro were georeferenced using the RTKLib^© software in PPK mode. This setup provides meter-level absolute accuracy but centimeter-level precision (relative accuracy). The images were then processed using the Pix4D Mapper software (https://www.pix4d.com/) where both image types were combined in a single SfM-workflow. The outputs (i.e., 3D point cloud, DSM and RGB mosaics had a centimetric precision (RMSE estimated at ca. 0.04 m)⁴².

Canopy roughness and canopy rugosity calculation using canopy 3D model

To show a potential canopy trapping effect, we related the canopy 3D properties to the dry deposition loads in the different plots, using the UAV structure-from-motion 3D product. It has been shown that dry deposition velocity (v_ds) of particles is linearly related to the surface roughness length (z₀)^27,28, and a simple parameterization for various surfaces—across various biomes—has rendered the following relation²⁷:

Hence, comparing across the successional stages of our forest succession, the canopy trapping hypothesis should manifest through a consistently different v_ds along the chronosequence, if the canopy complexity varies in z₀. As such, we aimed to determine z₀ along the chronosequence, in order to determine the v_ds, to see if this related to the differences in dry deposition loads. Hence, z₀ was determined following Lettau⁴³ and De Vries et al.⁴⁴:

where C is a constant—assumed 0.5⁴⁴, H is the average obstacle height and λ is the obstacle density of the roughness elements. H was derived based on field inventories and field-based measurements of tree height. In turn, λ was determined following de Vries et al. and Li et al.^44,45 as:

In other words, λ integrates the positive Δy divided by the Δx over a certain transect of the canopy height model of a certain experimental plot. This was done on 32 transects in the north–south direction of the plots, as well as 32 transects in the east–west direction for every plot. However, the calculation of λ over height transects is prone to instrument noise in the canopy height models. To filter out this noise, the transect data can be smoothed by applying a rolling means of a certain size, i.e., including a certain number of neighbor points along the transect. To determine the best size of the rolling mean, the calculated λ value can be plotted against the number of measurements that are included in the moving average. The resulting slope of λ versus the size of the moving average has two components; a first steep descent, followed by a less steep stable slope. It is accepted that the first steep descent is caused by filtering out unwanted noise in the data, while the second part is the actual λ (Fig. S3). Instead of detecting a breakpoint in this curve per plot, we qualitatively inspected the different curves, and selected a general number of measurements to be included in the rolling average (19 points in total, i.e., 9 points on both sides of the point itself) to determine λ, to avoid that an inconsistent rolling mean would drive differences in the λ estimation across plots. The drone flights were performed 6 months after the end of the monitoring period. We do not expect this to affect the forest canopy structure, but the agricultural field were heavily overgrown already by then, compared to the bare fields they were at the beginning of the monitoring period. For this reason we divided the standard deviation of the agricultural plots by three in Fig. 2, for visualization purposes, assuming a linear increase of variability over time, and targeting the average (after 6 months of monitoring) situation in the monitoring period.

Additionally, we calculated canopy rugosity⁴⁶. In essence, this is a much simpler proxy that estimates the standard deviation of the vertical variation of the upper canopy layer. For this, we used transects that were generated based on the CHM of the plots, as described above, to calculate the standard deviation of canopy height. The final ‘rugosity’ parameter is an average of the standard deviation over all the different transects per plot.

Data analysis

We calculated weekly throughfall fluxes of P by multiplying the concentration of P with the weekly recorded rainfall volume per plot. We aggregated this data to annual fluxes by calculating a weekly average and multiplying by 52 weeks per year, to account for some missing values—in case a throughfall collector was stolen, or sampling was logistically impossible that week (on one occasion). Up-to-date, it has proven difficult to quantify net nutrient deposition in forest ecosystems, since most studies are based on throughfall collection. This method uses rainfall collectors set up in open-field, paired with collectors under the forest canopy (i.e., throughfall collection). The measured deposition load in the open field is typically considered ‘wet’ deposition, while the throughfall collectors are considered to be the sum of wet deposition, dry deposition and canopy leaching, minus canopy uptake of nutrients. Contrary to other direct quantification methods of dry deposition^8,47, this method is attractive for its simplicity and also the potential to fully include canopy complexity effects on the dry deposition budgets. However, the downside of the method is that it is complex to disentangle dry deposition, from a canopy process such as leaching and uptake. There have been a variety of theoretical models that attempt disentangling both exogenous and endogenous effects in throughfall deposition loads, but these show varying success and results (for an overview see ref. ⁴⁸). In this study, we approached this differently by combining the deposition loads of a proxy for biomass burning (DBC) with the P deposition loads. Theoretically, subtracting the weekly wet deposition, i.e., the deposition loads as recorded in the open field, from the throughfall deposition load in the different forest types, results in the weekly net canopy effect, i.e., either canopy uptake/leaching or the added dry deposition effect of the canopy. We did this for both DBC deposition loads as well as P deposition loads, and subsequently assessed if the P deposition could be consistently predicted using only DBC deposition loads across the plots. We performed standardized major axis (SMA) regression between DBC deposition load and P deposition load, and between rainfall volume and P deposition. Where assumptions of the model fit were violated, we logtransformed both variables and performed a log–log SMA regression, but these were backtransformed for visualization purposes. SMA regressions were fitted using the ‘smatr’ package⁴⁹, using the robust estimation method. All analyses were carried out with the R software⁵⁰.

Wind trajectories and satellite data

To show how biomass burning aerosols might be sourced throughout the year across the African continent, we crossed a burnt area dataset with backward wind trajectories ending at the site, similar to previous work at the same site¹⁹. The fire pixel dataset covering the entire African continent was obtained from NASA’s MODIS Collection 6 NRT⁵¹. The backward wind trajectories were generated using the Hybrid Single Particle Lagrangian Integrated Trajectories (HYSPLIT) model provided by the National Oceanic and Atmospheric Administration Air Resources Laboratory (NOAA ARL)^52,53,54,55 with the reanalysis 2.5 degree archive meteorological dataset. For the entire study period, we generated one trajectory at noon for every day, ending above the study site pixel at 500 m above ground level, going back 1 week. Instead of displaying all the trajectories separately, we constructed a convex hull around the trajectories to effectively show the potential source areas to be crossed by winds (Fig. S3). Additionally, we extracted the hourly simulated BCSMC from above our study site from the Modern-era Retrospective Analysis for Research and Applications version 2 (MERRA-2) from the Goddard Earth Observing System Model, version 5⁵⁶.

Uncertainties

This study reports some of the first in situ measured P deposition loads for central Africa. However, due to logistical constraints and local conditions, we had to design the monitoring setup in a robust and low-tech manner. These adaptations for the long-term monitoring might bring about some additional uncertainties we want to report here. First, the rainfall and throughfall collectors we used in the field did not have cooling (no grid power available, and safety issues in maintaining technological setups in the field), nor did we add preservation chemicals to the collection bottles (due to the environmental risks involved with the disposal of mercury. Nevertheless, we judge that the consequences of this are minor for total phosphorus concentration in the samples. The DBC quantification, by definition, was performed on filtered samples and thus did not include particulate BC (PBC). However, we use DBC merely as a tracer for (dissolved) total P sourced by biomass burning. Given the pre-dominance of biomass burning as the source of combustion-derived aerosols in the Congo (consistent source), the long-range transport of aerosols to our study site, and the sub 0.2 μm size fraction of these aerosols, we are confident that DBC is a reliable tracer for combustion-derived inputs. Finally, we use ‘wet’ and ‘dry’ deposition throughout this manuscript, while in fact our wet deposition was not collected with ‘wet-only’ collectors. In practice, the open-field depositional loads will also have dry deposition collecting on the collector surface. Hence, our ‘wet’ deposition estimates are slightly overestimating wet deposition, while the ‘dry’ deposition estimates are slightly underestimating actual dry deposition.

Source: Ecology - nature.com