A question of scale: modelling biomass, gain and mortality distributions of a tropical forest

. The quantification of forest carbon budgets is important for understanding the role of forests in the global climate system. Given the variety of different methodologies (inventories, remote sensing, modelling) and spatial resolutions involved, methods for consistent transfer between scales are needed. In this study, the scaling of variables, which drive the carbon budget, 15 was investigated for a tropical forest in Panama. Based on field inventory data from Barro Colorado Island, spanning 50 ha over 30 years, the distributions of aboveground biomass, biomass gain and mortality were derived at different spatial resolutions, ranging from 10 to 100 m. Methods for fitting parametric distribution functions were compared. Further, it was tested under which assumptions about the distributions a simple stochastic simulation forest model could best reproduce observed biomass distributions at all scales. Also, an 20 analytical forest model for calculating biomass distributions at equilibrium and assuming mortality as a white shot noise process is presented. Scaling exponents of about -0.47 were found for the standard deviations of the biomass and gain distributions, while mortality showed a different scaling relationship with an exponent of -0.3. Lognormal and gamma distribution functions fitted with the moments matching estimation method allowed for consistent parameter transfers between scales. Both forest models 25 (stochastic simulation and analytical solution) were able to reproduce observed biomass distributions across scales, when combined with the derived scaling relationships. study provides insights about transferring between scales and its effect on frequency distributions of forest attributes, which is particularly important for the increasing efforts to combine information from sources such as inventories, remote sensing and modelling. Based successive censuses, AGB gains per tree were calculated. The inventory was aggregated for square shaped plots with 10-, 20, 50- and 100-m side Attributes calculated for each of these plots were standing AGB at each census and AGB gain, loss and mortality for each census interval. AGB loss was the sum of AGB of all trees alive in the census at the beginning of an interval and dead in the following census. AGB mortality was AGB loss divided by the standing AGB at the beginning of 100 the interval. Due uncertainties some trees showed negative At individual tree level these negative corresponded positive When single gains, losses and mortalities were 105 downscaled annual through division by 5 years.


Introduction
Forests are an important pool in the global carbon cycle. It is estimated that they store about 45% of the terrestrial carbon, mainly in trees in the form of living aboveground biomass (AGB) and in the soil (Bonan, 2008). Forests are inherently dynamic systems which are sensitive to climate patterns and a wide range of disturbance regimes (Lewis et al., 2015). Hence, the worldwide quantification of forest biomass stocks and changes are important. Different approaches are being used for this 35 purpose. These include forest inventory, eddy flux measurements, remote sensing and dynamic ecosystem modeling (Mitchard, 2018;Shugart et al., 2015). Forest biomass is variable in space and time, driven by the interplay of productivity and mortality (Rutishauser et al., 2019).
Methods for quantifying forest biomass all operate at different spatial scales, regarding their extents and resolutions.
Inventories measure and map forests at the resolution of individual trees. Direct measurements taken for every tree are usually 40 the stem diameter at breast height (DBH) and the X-and Y-coordinates of the stem foot. Attributes like tree height and tree biomass can be derived from DBH via allometric relationships (Chave et al., 2005). From the individuals, aggregated forest attributes in per area units can be derived. Inventory plots have extents ranging from a few square meters (e.g., national forest inventories) to several hectares (e.g., research megaplots). Eddy flux towers have typical footprint sizes in the range a few hundred meters, for which they acquire net gas exchange between the land surface and atmosphere (Rebmann et al., 2005). 45 Such field measurements are an important information source for deriving tropical forest carbon budgets at continental scales (Brienen et al., 2015;Hubau et al., 2020).
Much larger extents can be covered by remote sensing. It is most common that remote sensing products are generated following area-based approaches. That means they are typically gridded maps consisting of square-shaped pixels with a wide range of side lengths, depending on the sensor platform and type of product. A range of resolutions from centimeters to a few meters 50 are typical for optical imagery from drone, airborne and high resolution spaceborne platforms and gridded products from airborne lidar (e.g., canopy height models). Products on canopy cover and species composition are typically generated at resolutions of 20 to 50 m, owing to the resolutions of frequently used multispectral (Landsat, Sentinel-2) and radar satellites (ALOS-PALSAR, TanDEM-X) (Hansen et al., 2013;Martone et al., 2018). Maps about productivity and disturbances are typically derived from the moderate resolution imaging spectroradiometer (MODIS). They have high temporal resolutions at 55 the trade-off of having spatial resolutions of only 250 to 1000 m (Justice et al., 2002;Running et al., 2004). Remote sensing products about forest structure and biomass are often derived at 100-m resolution (at regional extents) or 500 to 1000-m resolutions (at global extents), because their estimation requires the derivation of aggregated metrics from primary measurements at finer resolutions, such as canopy height or canopy cover (Asner et al., 2013;Baccini et al., 2012;Saatchi et al., 2011). 60 Models which simulate the biomass dynamics and carbon balance of forests can have different complexities. This applies to the included processes and to the resolution. At one end of the spectrum, there are spatially explicit, individual-based models which include physiological processes for modelling the effects of environmental drivers on the carbon balance (Shugart et https://doi.org/10.5194/bg-2022-24 Preprint. Discussion started: 22 February 2022c Author(s) 2022. CC BY 4.0 License. al., 2018. On the other end, there are spatially implicit differential equation models describing the carbon balance in terms of and aggregated biomass stock and fluxes of productivity and mortality (Fisher et al., 2008). While the former type of model is 65 time-demanding to parameterize and computationally expensive, the latter type of model lacks variability and spatial patterns in the state variables. However, spatial variability may be introduced in simple models by simulating them on a grid or by using methods from statistical physics. The spatial resolution of forest simulations can vary and this has effects on the results. Different approaches of how to deal with differences in resolutions of ecological models and data have been used an the past and can broadly be grouped into pre-, in-and post-model scaling (Fritsch et al., 2020). 70 The heterogeneity of a forest landscape can be described by frequency distributions of forest attributes, e.g. biomass. Despite the large differences in measurement areas between the methods, it is common practice to report biomass stocks and carbon fluxes in per hectare units. Having units standardized to one reference scale is important to make values comparable and transferable. Since the actual areas, which these values represent, do often deviate from one hectare, it is crucial to consider the scale-dependence of the shapes of distributions (Chave et al., 2003). Descriptive statistics of distributions and relationships 75 between variables may critically depend on the chosen gridding approach, a phenomenon known as the modifiable areal unit problem (Wong, 2008). In times where field measurements, remote sensing and model simulations are increasingly used in combination, approaches for harmonizing the different datasets with regard to spatial resolution are required.
Hence, with regard to forest biomass distributions and their temporal dynamics the following problems arise: Field measurements provide us with high resolution information, but for limited extends, while remote sensing can cover large 80 extends with limited resolutions. Forest models make use of both kinds of information, either as input or to validate their output against. But it is often unclear how scale affects observed and simulated distributions.
In this study, it was investigated how the distributions of several relevant forest variables change with scale. Different ways of describing them with distribution functions across scales from 10 to 100 m were compared. It was tested whether a simple carbon balance model can reproduce the biomass distribution from productivity and mortality distributions at all scales. 85 Additionally, a method for direct calculation of the biomass distribution is presented. The three main questions of the study were: 1) How do forest ecosystem attributes such as standing biomass and biomass gains, losses and mortality vary with spatial scale? 2) Can a stochastic simulation forest model reproduce realistic biomass distributions from gain and mortality information at different scales? 3) Can a stochastic analytical forest model reproduce these biomass distributions as well?

Field data analysis
The 50-ha forest inventory plot (ForestGEO) on Barro Colorado Island (BCI), Panama, served as a study site (Condit, 1998;Condit et al., 1995Condit et al., , 2019Hubbell et al., 1999). The inventory comprised single tree measurements of all trees with a DBH ≥ 1 cm (ca. 250,000 trees), including the spatial position of each individual (Condit et al., 2012). Aboveground biomass values of each tree were calculated based on a common allometric equation for tropical trees (Chave et al., 2005). Censuses used in this study comprised data from 1985 to 2015 collected in 5-year intervals. Based on unique tree identifiers and AGB values in successive censuses, AGB gains per tree were calculated. The inventory was aggregated for square shaped plots with 10-, 20-, 50-and 100-m side lengths. Attributes calculated for each of these plots were standing AGB at each census and AGB gain, loss and mortality for each census interval. AGB loss was the sum of AGB of all trees alive in the census at the beginning of an interval and dead in the following census. AGB mortality was AGB loss divided by the standing AGB at the beginning of 100 the interval. Due to measurement uncertainties some trees showed negative AGB gains. At the individual tree level these negative AGB gains corresponded to 29.6% of the positive AGB gains. When aggregating AGB gain as the sum of single tree AGB gains the negative proportion would decrease to 12.3% (10-m), 7.2% (20-m), 1.5% (50-m) and 0% (100-m). To obtain unbiased, positive AGB gains a correction was done at the single tree level: All negative AGB gains were set to zero and all positive AGB gains were reduced by 29.6% to maintain the overall AGB gain. AGB gains, losses and mortalities were 105 downscaled to annual values through division by 5 years.
For investigating the scaling behavior of AGB beyond the 50-ha plot, an airborne lidar dataset from 2009 covering the whole 15-km² island (Lobo and Dalling, 2014) was used. AGB was predicted at 100-m resolution from mean top-of-canopy height of a 1-m resolution canopy height model using a power law regression model (Knapp et al., 2020) and aggregated at the 200-, 500-and 1000-m scale. Pixels intersecting the coast line were excluded to have pure forest pixels only. 110

Frequency distributions
The frequency distributions of all four variables were fitted using the R package 'fitdistrplus' (Delignette-Muller and Dutang, 2015; R Development Core Team, 2014) and using two alternative probability distribution functions (PDF): the 1) gamma and 2) lognormal distribution function. Both distribution types can be described with two parameters, respectively. If the parameters are known, the descriptive statistics of the PDFs (arithmetic mean and standard deviation SD) can be calculated 115 from established equations and vice versa (Table 1). Two different methods were used to fit the PDFs to the empirical data using: 1) maximum likelihood estimation (MLE) and 2) moments matching estimation (MME). The optimization criterion for MLE is to find distribution parameters which maximize the likelihood for all observations. MME aims at finding the parameters for which the calculated moments match the empirical moments.

Scaling equations
In the following, a general equation for scaling of probability distributions is formulated (based on considerations by Dubayah 125 et al. (1997)). Let θ be a parameter of a probability distribution (e.g., SD, shape, rate, etc.). Let λ and Λ be the side lengths of plots of two different sizes and κ be a scaling exponent characterizing how the parameter θ changes between scales. The values of θ (characterizing the frequency distributions; see Table 1 for examples) at scales λ and Λ shall be called θλ and θΛ, respectively. Then θΛ can be calculated from θλ in the following way: If values for the parameter θ are known at different scales, κ can be derived via log-log-linearization.
(2) log 10 ( Λ ) = • log 10 ( log 10 ( Λ ) = 2 • log 10 ( Λ ) Hence, in log-log-space, the ratio of θ at two different scales is a linear function of the ratio of areas with κ as slope. If the ratio 135 of side length instead of areas is used, the slope of the relationship becomes 2κ. In the further analyses, the area ratio was used (Eq. (3)). The focus of this study was to investigate how the scaling equation can be used to calculate the SD of forest attribute https://doi.org/10.5194/bg-2022-24 Preprint. Discussion started: 22 February 2022 c Author(s) 2022. CC BY 4.0 License.
distributions at a certain target scale from a known SD at a reference scale, which will be called rescaling further on. A graphical example is given in Fig. 1 for the SD of the biomass distribution using 10 m as the reference scale.

Figure 1. Scale dependence of the standard deviation (SD) of the distribution of aboveground biomass (AGB). The left graphic shows the SD values at each scale (Λ) plotted over the respective plot side lengths Λ. The right graphic shows the scaling relationship using plot areas (Λ² instead of Λ) and after standardizing with 10 m as the reference scale (λ), and plotting the values in a log-log-linear
way. The slope of the grey regression line is the scaling exponent κ.

Comparing frequency distributions 145
Two metrics were used to quantify the agreement between different PDFs. The first metric was the relative overlap between PDFs (OVL; Inman and Bradley, 1989), i.e., the intersection of the areas under the curves (AUC) divided by the union of the areas under the curves of two PDFs. The R package 'overlapping' was used for computation of OVL (Pastore and Calcagnì, 2019).
The second metric was the relative error of standard deviations (RESD), which we defined as the absolute difference between the SD of an estimated (e.g., rescaled or simulated) PDF1 and an observed (empirical) PDF2, normalized by the SD of the latter.
Hence, two PDFs are in close agreement if OVL is close to one and RESD is close to zero. OVL is a good quantitative measure 155 for visual agreement of the main bodies of the distributions. RESD, in contrast, is sensitive to the influence of the extreme values in the tails of the distributions, which is often not apparent from visual inspection of the PDFs.

Forest simulation model
A simple grid-based forest growth model was used to simulate forest dynamics over time and obtain the resulting AGB frequency distribution at mature stage. The model is based on the model suggested by Fisher et al. (2008), in which the change 160 of AGB (B) is described as a differential equation involving an AGB growth parameter (G) and an AGB mortality parameter (M).
The analytical solution of this equation provides the average trajectory of forest biomass (over a large area) for given G and M and the average mature AGB is G/M. 165 For simulating small scale heterogeneity including patch dynamics, Fisher et al. (2008) used a grid-based approach. Each grid cell represented 10 m × 10 m of forest area. At each simulation time step, they applied Eq. (7) to calculate the change of AGB in each cell, based on a constant AGB gain (G) and a stochastic mortality (M). A mortality event in their case meant that B was set back to zero in a patch and growth restarted from bare ground. 170 In this study, the parameters G and M were assumed to be distributed like the AGB gain and AGB mortality values derived from the field inventory. Variability in gain G represents differences in site conditions or growth rates of trees and was assumed to vary in space but stay constant in time. Hence, in the simulation each patch was assigned a random G drawn from the AGB gain distribution and this G was kept constant for the patch throughout the simulation run. Variability in mortality rate M reflects the stochastic nature of mortality. In a common year, only small amounts of AGB are lost per patch, due to small trees 175 dying from competition. Large AGB losses due to dying large trees are rare events. Hence, mortality rate M was drawn randomly from the AGB mortality distribution for each patch at each time step. This simulation approach was applied at different grid cell sizes (10-, 20-, 50-, 100-m). The total simulation area was set to 25 km² (5 km × 5 km) to obtain smooth output distributions. The simulation time step was one year and 300 years were simulated to reach a mature forest. This stochastic simulation forest model will be called SFM further on. 180 The AGB distributions in equilibrium state after 300 years of simulation time were analyzed by comparing them to the AGB distributions from the field inventory. The simulations were conducted with G and M coming from different PDFs. The PDFs could either be gamma or lognormal, fitted with either MLE or MME and be derived from the field at either 10-, 20-, 50-or 100-m scale (reference scale). The scaling equation was used to rescale SD and PDF parameters from one reference scale to the other three scales, to conduct simulations at all four scales for each of the possible input constellations (pre-model scaling). 185 The resulting biomass distributions at the end of the simulations were aggregated from the finer to the coarser scales (postmodel scaling). For evaluating the goodness-of-fit for the AGB distributions, OVL and RESD were calculated pairwise https://doi.org/10.5194/bg-2022-24 Preprint. Discussion started: 22 February 2022 c Author(s) 2022. CC BY 4.0 License. between each simulated (and aggregated) distribution and the empirical distribution at the respective scale. For evaluating the simulation run as a whole, the mean OVL and mean RESD over all single distributions was calculated.

Derivation of biomass distributions from theorywhite shot noise 190
As a second approach, the AGB distribution can also be derived analytically as a function of G and M, assuming M to be white shot noise (WSN). This type of noise is used in environmental modeling to describe pulsed processes that occur at irregular time intervals (e.g., rainfall at daily resolution). WSN assumes pulse intensities and interval lengths both to be exponentially distributed and is, therefore, characterized by two parameters: the mean intensity and the mean occurrence probability.
It was assumed that the AGB mortality at small spatial resolution (10 m × 10 m) can be described by WSN. Since observed 195 AGB mortality is never exactly zero, due to small trees dying every year, the mean occurrence probability can be assumed as one. The mean intensity then directly corresponds to the mean AGB mortality M ̅ . It can be shown that under this WSN assumption the AGB for a given G and M ̅ value follows a gamma distribution of the following form (Ridolfi et al., 2011): For a range of k different G values, which themselves follow a gamma or lognormal distribution, AGB of a forest landscape 200 can be described by a mixture of gamma distributions.
The exponential distribution of mortality intensity, as assumed by the theory, has a monotonically decreasing PDF. Such PDFs are usually only observed for mortality at small scales, while mortality follows unimodal distributions at larger scales (see e.g. Fig. 4). Hence, for deriving AGB distributions from WSN theory at larger scales, there are different options, of which two 205 were investigated here. 1) The simplest one is to apply the WSN theory directly at all scales, ignoring the actual shape of the mortality distribution. 2) Alternatively, WSN theory was applied to derive the AGB distribution only at 10-m scale, and the SD of this distribution was rescaled in a post-model scaling using Eq. (1) and an empirically derived scaling exponent κ. The rescaled SD was used to approximate the AGB distribution with a lognormal distribution. In the following, this stochastic analytical forest model, involving WSN theory, will be called AFM. The main steps in the analysis, from inventory data to 210 frequency distributions to scaling relationships and forest model applications at different spatial scales are summarized in

Scaling relationships of standard deviations of forest attribute distributions
The relationships between spatial scale (plot area) and standard deviations of the frequency distributions were derived for 220 AGB, gain, loss and mortality within the 50-ha forest plot for the scale range from 10 to 100 m. The scaling exponents for AGB (κ = -0.461), AGB gain (κ = -0.468) and AGB loss (κ = -0.467) were similar, whereas it was different for AGB mortality (κ = -0.301). The log-log-linear relationships are shown in Fig. 3.  . (1-3)). Plotted are the log10transformed ratios of the SDs over the log10-transformed ratios of plot areas, with 10 m being the reference scale (λ). The top axis shows the plot side lengths (Λ) and the right axis shows the values of the SD for each scale, respectively.
Beyond the 50-ha plot, AGB was analyzed via a lidar-derived AGB map (Fig. 4). The lidar analysis demonstrates that at 100m scale the inventory-based SD of the AGB distribution within the 50-ha plot (59 t ha -1 ) matches the lidar-based SD of the 230 whole island (59 t ha -1 ), while the lidar-based SD within the 50-ha plot was somewhat smaller (49 t ha -1 ). At the 200-m scale, https://doi.org/10.5194/bg-2022-24 Preprint. Discussion started: 22 February 2022 c Author(s) 2022. CC BY 4.0 License. the lidar-based SD within the 50-ha plot (29 t ha -1 ) conforms the scaling relationship. At the 500-m scale, it is somewhat higher than expected (20 t ha -1 ), but at this large pixel size a strict spatial overlap with the 50-ha plot is not possible and the area contributing to this value is 150 ha (6 pixels). The island-wide SDs at 200-m (47 t ha -1 ) and 500-m (37 t ha -1 ) resolution are higher than predicted by the scaling relationship. This additional variability can be explained by the larger extent and possible 235 autocorrelation in the spatial distribution of biomass across the island. At 1000-m scale, the island-wide SD (10 t ha -1 ) is close to the scaling relationship, which can be explained by the fact that only five square kilometers in the center of the island, around the 50-ha plot, were analyzed while all square kilometers intersecting the shore line were excluded.

Fitted distribution functions 245
For each of the four variables of interest (AGB, gain, loss and mortality) the best describing distribution function was identified, by fitting either lognormal or gamma distribution functions with the MLE or MME method. Selected fits for all four spatial  From each fit at one scale (reference scale) the expected SDs at all other scales (target scales) were calculated using the scaling equation (Eq. (1)) with the empirically derived scaling exponents κ, respectively. From these rescaled SDs the function 255 parameters and probability density curves were derived using the equations in Table 1. An example for downscaling AGB from 100-m to the three smaller scales is shown in Fig. 6.  To identify the best fitting distribution and method for each variable, i.e., the most consistent across scales, two criteria were used: 1) the mean overlaps (OVL) between empirical and fitted probability densities and 2) the mean relative error of the SD (RESD). Mean in both cases refers to mean over all four spatial scales. Table 2 lists the values for these criteria for all distributions and methods and Fig. S1 shows the best matching distribution function for each variable at all spatial scales. 265 Table 2. Goodness-of-fit criteria values for the different distribution functions and fitting methods for all four variables. Best matching distributions, i.e., those with the largest mean overlap (OVL) and smallest mean relative error of standard deviation (RESD), are highlighted in bold font. Mean hereby refers to the mean over all four spatial scales (10, 20, 50, 100 m)

270
For three of the four variables, there was a clear best fit according to the criteria (Table 2): AGB (lognormal MME), AGB loss (lognormal, MME) and AGB mortality (gamma, MLE or MME equally good). For AGB gain, however, the best fitting lognormal distribution functions derived from the MLE and MME methods were considerably different from each other (Fig.  275   5). While the MLE fits had the higher overlap with the field data, they produced large errors when used for calculating the SD at the other scales (large mean RESD). Consistent rescaling of the SD of AGB gain was only achieved when using the MME fit, albeit the weaker overlap of these distributions with the field data. Comparing Fig. S1b with Fig. S2a illustrates the differences of the methods for AGB gain in detail and Fig S2b+c illustrates the RESD for both.

Simulation results 280
Simulations with the stochastic forest model (SFM) were conducted at the scales of 10, 20, 50 and 100 m. The resulting biomass distributions at the end of the simulations were aggregated to the coarser scales, respectively. By testing a range of 64 different combinations of input distribution functions (G and M) and reference scales (at which these distributions were derived), well performing combinations were identified based on the goodness-of-fit criteria OVL and RESD (Table S1, Fig.   S3 to S6). It was found that simulated AGB distributions matched the field distributions particularly well when G and M 285 distributions were derived at 50-m scale (reference scale). The distribution functions identified as best for gain and mortality, based on the field data analyses (Table 2), were also among the best for G and M, based on the simulation analyses (Table S1).
If G was modelled as a lognormal (MLE fit at 50-m scale) and M as gamma (MME fit at 50-m scale) distribution, rescaling and subsequent stochastic simulations resulted in consistent AGB distributions at all scales (Fig. 7).

295
It was found that several versions of the SFM simulation could produce realistic AGB distributions across scales. The best matching AGB distributions were observed for simulations using 50-m as reference scale (mean OVL = 77%, mean RESD = 9%, Table S2). At the 20-m reference scale, the largest observed mean OVL was 75% and the smallest observed mean RESD https://doi.org/10.5194/bg-2022-24 Preprint. Discussion started: 22 February 2022 c Author(s) 2022. CC BY 4.0 License. was 12% (albeit, the two were from different simulations, Table S1). At the 10-m reference scale, the best mean OVL was the smallest among all reference scales (65%), while the best mean RESD was good (13% , Table S1). At the 100-m reference 300 scale, the best mean OVL was better (70%), but the best mean RESD was far higher than for any other reference scale (34% ,   Table S2).
One consistent finding across scales was that simulations using lognormal distributions and the MLE fitting method were more abundant among the best simulations than the ones using gamma distributions and the MME method. In fact, at any reference scale, the approach which used only lognormal and MLE for modeling G and M ranked among the best simulations. 305 Approaches using gamma or MME for more than one distribution (G and M) never ranked among the best (Table S1 and S2).
Examples for simulated AGB distributions at all scales (with OVL and RESD values), when starting from different reference scales, are shown in Fig. S3 to S6.

White shot noise calculations
As an alternative to the simulations with the stochastic forest model (SFM), the calculation of AGB distributions from an 310 analytical forest model (AFM) was tested. The mortality was assumed as white shot noise, i.e., only its total mean (M ̅ = 0.0285) was required. Since the assumption of WSN mortality was considered most appropriate at the finest scale, the PDF (gamma mixture distribution) was derived at 10-m scale and upscaled with the scaling relationship, using a lognormal distribution as approximation for the mixture distribution. The distributions resulting from the AFM (blue lines in Fig. 7) are close to the field data and the SFM outputs. A quantitative analysis based on OVL and RESD is presented in Fig. S7, which also shows how a 315 direct application of WSN at larger scales leads to deviations in the AGB distributions.

Discussion
It was shown how the distributions of the forest attributes, which are critical for the carbon balance, vary with spatial scale.
Methods for fitting the distributions best across all scales and transferring their standard deviations between scales were identified. The obtained distribution functions for biomass gain and mortality were further tested with a stochastic simulation 320 forest model and with a stochastic analytical forest model at multiple scales for their ability to reproduce observed biomass distributions.

Scaling of forest attributes
For transferring between the distributions at the different spatial scales, empirical scaling coefficients were derived. The exponents for scaling the SDs of the distributions were very similar for AGB stocks, gains and losses (all between -0.46 and -325 0.47) while the one for mortality was considerably different (-0.3). Theory states that the SD of an independent identically distributed random variable changes with the square root of the area ratio between scales, i.e., a scaling exponent of -0.5. The slight deviation from -0.5 is an indication of spatial structure in the data, i.e., that grid cells are not fully independent. Further analyses should look into the spatial patterns and mechanisms at the individual level, i.e., tree positions and their size and biomass distribution, by using methods such as point process models (Lister and Leites, 2018) or network analysis (Schmid et 330 al., 2020), to understand how they shape the scaling at aggregated levels.
The very different scaling exponent of mortality is likely caused by the fact that mortality is a ratio of AGB loss and stock.
This characteristic of being a ratio makes mortality non-additive. When aggregating mortality values from smaller to larger scales, the weighted mean has to be used, with AGB as a weighting variable. It further has the consequence that the mean mortality over the whole 50-ha plot is not the same when calculated at different resolutions. The mean mortalities of all 10-335 and 20-m patches were 0.035 and 0.031, respectively and only stabilized above the 50-m scale at 0.029 (Fig. S1d). All the other attributes in this analysis were additive and had a stable mean across all scales.
The lidar-based scaling of AGB beyond the 50-ha plot and across the whole island showed that the SD of the distribution decreases further, but does no longer strictly follow the scaling relationship from within the plot. With increasing scale, the spatial pattern of AGB distribution is possibly increasingly dominated by environmental covariates, e.g., topography, which 340 may alter the slopes of scaling relationships. Additionally, the spatial coverages of the AGB maps were different, due to exclusion of pixels intersecting the coast line, which led a very small sample size of only five pixels at the largest scale (1000m). Further analyses could investigate the scaling in heterogenous landscapes, e.g., based on synthetically generated landscapes with spatial autocorrelation or based on regional or global remote sensing products. With such approaches, an invariabilityarea-relationship has been suggested, which can describe the temporal dynamics of vegetation productivity and connects local 345 to global scales (Wang et al., 2017).

The choice of the proper distribution fit
It was found that for the majority of forest attributes (AGB, loss, mortality) a best fitting distribution function, which was also consistent across scales, could be identified. The example of AGB gain, however, demonstrated that different methods for distribution fitting can have different outcomes. The target of the MLE method is to find the distribution function which 350 maximizes to likelihood of all observed values. Therefore, it usually results in high overlaps of empirical and fitted distributions. However, it does not necessarily lead to a good match of the moments of empirical and fitted distribution.
Because of this, SDs of the fitted PDFs deviated from the empirical SDs, and they deviated differently at each scale. This led to inconsistent PDFs when applying the scaling relationship for rescaling of the SDs. By definition moments matching is the target of the MME method. With MME, SDs were correctly estimated at all scales and could consistently be rescaled. The 355 price for the matching moments were the lower overlaps between empirical and fitted distributions obtained from MME.

Scaling of the forest models
It was found that both, the stochastic simulation forest model and the analytical stochastic forest model could reproduce the biomass distributions observed in the field at different spatial scales. For the SFM, it proved to be best, to chose an intermediate reference scale (50-m) for deriving the information about gain and mortality and using a pre-model scaling to drive simulations at the different other scales. In this way, the SFM could reproduce the biomass distributions at the four simulation resolutions and also their aggregations obtained in post-model scaling. Choosing the lower (10-m) or upper (100-m) end of the scale range as reference scale led to less agreement of the simulated and the field distributions. Using the small scales as reference scale is probably problematic, because mean mortality is positively biased at these scales, due to the non-additivity of mortality discussed earlier. The rather weak performance of simulations using 100-m reference scale needs further investigation. 365 The SFM produced the most consistent results when the input distributions were fitted using MLE. This was unexpected, considering the earlier finding that MME allows for more consistent scaling of SDs based on the field data analysis. An explanation could be that, while MME fits are more accurate with regard to the SD, MLE fits are overall "closer" to the real distributions (higher OVL). The deviation of MME fits from the real distributions (smaller OVL) may pose a problem in simulations where these distributions are used for drawing random numbers many times. MLE fits might be more appropriate 370 in this case, especially if they were obtained at intermediate scales, from where distribution parameters, like SD, can be rescaled in both directions without causing much bias.
A fast solution for approximating the AGB distribution of the forest from the AGB gain and mortality distributions was given with the AFM using white shot noise. This approach does not require iterative simulation over time. However, it requires the assumption that mortality follows an exponential distribution, which is not the case at any scale. At the 10-m scale the mortality 375 distribution is at least more similar to an exponential distribution than at the larger scales. Indeed, assuming white shot noise directly at large scales did not approximate the AGB distributions well. But, assuming white shot noise only at the 10-m scale and applying a post-model upscaling, enabled a good approximation of the AGB distribution across all scales without running a simulation.

Implications for applications 380
Quantifying the carbon budget of forests is one of the most important applications of forest models and remote sensing.
Increasingly the two techniques are being combined for this purpose, either using remote sensing information for model parameterization, initialization, calibration or validation or using model information as synthetic ground truth for remote sensing interpretation (Plummer, 2000). In this context it is very important that remote sensing products and models describe variables at the same spatial resolution. In case they are operating at different resolutions appropriate scaling of the variables 385 should be conducted.
A detailed model-based understanding of how AGB distributions emerge from growth and mortality could, e.g., be used to estimate disturbance intensities from remote sensing (Williams et al., 2013). A detailed understanding how canopy height changes at a certain scale relate to biomass changes, allows for a direct quantification of the net forest carbon changes fro m remote sensing (Knapp et al., 2018). Knowledge about the scaling of forest attribute distributions is also required for 390 downscaling of gridded maps (i.e., super-resolution), for purposes such as data assimilation (Hill et al., 2011;Rödig et al., 2017) or pixel-to-point comparisons (Rammig et al., 2018). Hence, with the advancement of forest monitoring and modeling the demand for consistent scaling methods is increasing. https://doi.org/10.5194/bg-2022-24 Preprint. Discussion started: 22 February 2022 c Author(s) 2022. CC BY 4.0 License.

Conclusions
The study has shown how the distributions of variables, which are important for the carbon budgets of forests, vary with spatial 395 scale. Biomass and its gain, loss and mortality could all be described with parametric distribution functions (gamma or lognormal) of varying spread at different spatial resolutions. The spread of these distributions, in the form of standard deviations, was described as a function of scale using power law relationships. These scaling relationships allowed for successful up-and downscaling of the respective distribution functions. The application of up-and downscaling for forest models was demonstrated. It was shown that the scaling relationships can be used to reproduce biomass frequency distribution 400 with a stochastic simulation forest model across a range of scales using measured parameters about gain and mortality from a single reference scale as input. It was further shown, that under the assumption of mortality being a white shot noise at the small 10-m scale, the equilibrium biomass distributions of the forest could be calculated using a stochastic analytical forest model. Such scaling approaches will hopefully facilitate the trans-scale integration of data about forest dynamics from various sources of information, such as forest inventory, remote sensing and modeling. This is important for enabling the processing 405 of information at the scale appropriate for the prevailing questions and feasible with given computational resources.

Code and data availability
A publicly available forest inventory dataset was analysed in this study. It can be found on Dryad: https://datadryad.org/stash/dataset/doi:10.15146/5xcp-0d46. The BCI lidar dataset is publicly available from the Office of Bioinformatics, Smithsonian Tropical Research Institute. The R code for the analyses and simulations is available upon request 410 from the corresponding author.