No evidence for fractal scaling in canopy surfaces across a diverse range of forest types

Canopy structural complexity is a key emergent property of forest ecosystems that directly relates to their ability to store carbon, cycle water and nutrients, and provide habitat for biodiversity. However, we lack a general framework for quantifying the structural complexity of forest canopies, and it remains to be seen if there are simple rules that explain the huge variation in canopy structure that we observe across the world's forests. A leading candidate for characterizing differences in structural complexity among forest ecosystems are fractal scaling laws, or measures of statistical self‐similarity. If forest canopies were fractal in nature, their structural attributes could be distilled into a single coefficient—their fractal dimension—which could then be used to directly compare structural differences within and across biomes. To test this idea, we used airborne laser scanning (ALS) data acquired across nine landscapes in Australia that span a huge environmental gradient and include everything from dry shrublands, tropical savannas, dense rainforests, to 90 m tall Mountain Ash forests. Using the ALS data, we built high‐resolution 3D canopy height models of each landscape so that we could quantify how closely they followed fractal scaling laws, based on a comprehensive set of fractal dimension estimators. Across all ecosystem types, fractal scaling assumptions were consistently violated, with clear and systematic differences between real‐world forest canopies and simulated fractal surfaces. However, deviations from fractality did vary predictably among sites, with forests in less arid environments, dominated by tall trees with large crowns, exhibiting a higher degree of self‐similarity across scales. Synthesis: Our study conclusively shows that canopy surfaces are not fractal beyond the scale of individual tree crowns. Nevertheless, we still observed ecologically meaningful and generalisable patterns in forest structure across scales, which closely reflect biophysical and physiological constraints on tree size and architecture. These patterns point the way towards a more general framework for characterizing ecological complexity.


| INTRODUC TI ON
Forests vary widely in their structural make-up and complexity.
Sometimes, only a few species coexist and form stands of great regularity, such as the cathedral-like beech forests of the temperate zone.Underneath their canopies time seems to stand still, with juvenile trees waiting up to a century for release into the upper layers (Petrovska et al., 2022).In other places, such as the dynamic rainforests of the tropics, hundreds of species grow next to and into each other.Lianas crawl up the trunks and across crowns, saplings compete with each other for sunspots and frequent treefall gaps let new generations burst through into the sunlight (Bazzaz & Pickett, 1980).
Such structural diversity is a key emergent property of forest canopies that directly relates to how they store carbon, cycle water and nutrients, and provide habitat for biodiversity (Fotis et al., 2018;Gough et al., 2019;Malhi et al., 2022;Shugart et al., 2010).And yet, despite its importance, we continue to lack a simple, quantitative framework that identifies generalities in 3D forest structure across data sources, scales and forest types (Atkins et al., 2023;Ehbrecht et al., 2021;Storch et al., 2018).
One of the most promising candidates for general patterns in 3D forest structure is fractal scaling laws, or descriptors of selfsimilarity, which have a long history in the study of tree architecture (cf.da Vinci's rule, Minamino & Tateno, 2014).As one quickly realizes when drawing a tree, the large branches of a tree tend to look like the tree itself, smaller branches tend to look like the large branches, and so do even smaller branches or twigs.Each part of the tree is, again, a little tree-sometimes, as in resprouting or clonal trees, literally so (cf. the Japanese daisugi, Toda, 1974).In practice, such self-similarity across scales allows us to reduce the complexity of tree crowns and express it through much simpler coefficients that summarize how structural properties change with each other (West et al., 1999;Zeide, 1991).Within limits (Coomes et al., 2011;Jucker et al., 2022;Lines et al., 2012;Muller-Landau et al., 2006;Owen et al., 2021), this allows us to predict tree properties from their parts, and has a wide range of applications in ecology, from understanding how trees optimize their crowns in response to environmental constraints or anthropogenic impacts (Arseniou & MacFarlane, 2021;Niklas, 1994), to predicting vegetation growth (Niklas, 2006), to greatly improving estimates of global forest biomass (Chave et al., 2014;Jucker et al., 2017;Sullivan et al., 2018).
Given the utility of fractal geometry for understanding tree architecture and allometric growth, it is natural to look for self-similarity at scales beyond individual trees (Zeide, 1991).Could, for example, space-filling create similar constraints at stand level as within trees and thus predict tree numbers and sizes (West et al., 2009)?And could the canopy surface be regarded as a fractal-like 'green carpet', which, during succession, is repeatedly stretched out and folded to maximize light interception (Birnbaum, 2001)?Expressing differences in forest structure via simple coefficients has obvious practical value, but the implications of scaling go beyond that: Scaling rules tell us something about the processes that generated them, typically multiplicative in nature (Mitzenmacher, 2003), provide a direct link to spatial patterns via autocorrelation (Mark & Aronson, 1984), and, in being dimensionless, allow for comparisons across data sources and ecosystem types.Ultimately, they could allow us to derive a general framework for ecological complexity that applies not just to forests but also to other ecosystems (Torres-Pulliza et al., 2020).
However, summarizing forest structure through simple metrics such as fractal dimension makes sense only if and when the underlying objects are indeed fractal by nature and exhibit some form of self-similarity over several orders of magnitude (Halley et al., 2004).
Since most properties of ecosystems are scale-dependent, and not in a consistent way (Atkins et al., 2023), this is controversial (Loke & Chisholm, 2022;Madin et al., 2023).If we follow, for example, metabolic scaling theory, forests can be understood as bundles of water-conducting tubes, with the number of twigs and branches in trees following similar scaling laws as the number of trees observed in a stand (West et al., 2009).But how far does this self-similarity extend?Does what is true for an individual tree-branches growing upwards and sideways in ways reminiscent of the trunk-also hold for an entire forest or forest landscape (Brown et al., 2002)?What would such a 'bundle of bundles' look like: small trees growing on top of medium-sized trees, medium-sized trees growing on large trees, and large trees held together by a mythological world tree (Figure 1)?
In reality, forests are not truly self-similar.When seeing a picture of a forest from above, we rarely struggle to tell at what scale it was taken, or confuse a grove of trees with an individual tree.
Instead, forests assemble horizontally, and how they develop depends a lot on external, landscape-scale factors, such as soil substrates and topography (Clark & Clark, 2000;Cushman et al., 2022;Jucker et al., 2018;Kane et al., 2015), as well as disturbance regimes, climate and biogeographic differences between sites (Banin et al., 2012;Muller-Landau et al., 2021).Whether forests show fractal behaviour thus depends a lot on whether the underlying drivers are fractal-like.This may be the case for certain disturbance types, such as wildfires (Malamud et al., 1998).But in many systems, local heterogeneity in soil composition and topography will introduce irregular patterns in species composition that disrupt scaling rules (e.g.Jucker et al., 2018).Climate, and in particular aridity, should be a key determinant of the scale at which deviations from fractality occur.
Both dry and humid environments can exhibit scaling behaviour (Jucker, 2022;Kéfi et al., 2007;Staver et al., 2019), but even in the absence of strong external drivers of vegetation structure, forest canopies have been observed to display a distinct break in self-similarity at scales that approximate the size of large canopy trees (10-100 m; airborne laser scanning, aridity, canopy height, forest structure, fractal scaling, self-similarity, structural complexity Shugart et al., 2010).Since hydraulic risks impose strong limitations both on the maximum size of trees (Lines et al., 2012;Niklas, 2007;Olson et al., 2018;Owen et al., 2021;Shenkin et al., 2020) and overall canopy height and density (Ehbrecht et al., 2021;Tao et al., 2016), aridity should be a good predictor of the transition from small-scale dynamics to landscape-level processes and the deviations from fractality it entails.
Fractal scaling can, of course, still be meaningful when forests are statistically self-similar only in some properties or under rescaling of their vertical dimension (the technical term for this is self-affinity, but for simplicity will be subsumed under self-similarity throughout the article).The coastline of Britain, for example-which played a key role in the development of fractal geometry (Mandelbrot, 1967)does not consist of many small coastlines of Britain.And yet it can still be described by a fractal dimension, as it is similar across scales in the way it is folded in 2D space.The empirical evidence for statistical self-similarity in forest landscapes is, however, ambiguous.
Fractal-based indices have been found to be useful for characterizing variation in forest structure at small scales (Ehbrecht et al., 2021) and hold potential for biomass inference (Guzmán Q. et al., 2020).
At the same time, as for many measures of structural complexity, it is possible that indices contain the same information as simpler forest properties, such as canopy height (Atkins et al., 2022;Gough et al., 2019), or simply summarize the properties of a few representative trees (Liu et al., 2022;Seidel et al., 2019).At landscape level, studies have found some evidence for fractal scaling in patchy forest types, fragmented by flooding, fires or human disturbance (Hastings et al., 1982;Krummel et al., 1987;Malamud et al., 1998), or in gap patterns in tropical forests (Solé & Manrubia, 1995).However, fractal scaling exponents are often instable across scales (Leduc et al., 1994;Palmer, 1988), and, given that comparisons across environmental gradients or with theoretical fractals are rare, it is unclear whether this is merely an artefact of methodology (Huang et al., 1994) or whether it reflects a genuine property of forest landscapes.

F I G U R E 1
A hypothetical self-similar forest stand, where the way trees form a canopy resembles the way individual trees grow their crowns.The image was generated using OpenAI's DALL-E 2 (https:// openai.com/ produ ct/ dall-e-2).The prompt to generate this picture was: 'Lots of trees growing on top of another tree, which itself grows on a tree, in the style of Albrecht Dürer, white background'.A visual representation of simulated fractal surfaces, and how they compare with empirically derived canopy surfaces, can be found in Figure S8 of the Supporting Information.
It thus remains to be tested (i) if there are similarities in scaling relationships across diverse forest landscapes, (ii) if these patterns can be described as fractal or fractal-like and (iii) if violations of scaling assumptions can be linked to aridity and the size of large canopy trees.To tackle these three major knowledge gaps, we compiled airborne laser scanning (ALS) data from nine treescapes in Australia (each 5 × 5 km in size) that span a wide range of canopy structuresfrom short-statured dry shrublands to tropical savannas and rainforest, all the way to 90 m tall Mountain Ash forests-and used them to generate detailed 3D canopy height models (CHMs) of each site.
Using the CHMs, we assessed the scale dependence of canopy structure with a variety of methods and tested whether any of the forest ecosystems exhibited fractal-like properties in their canopy surfaces.We examined how forest types differed in their structural patterns and whether differences in deviations from fractality could be related to aridity.Finally, for all sites, we manually delineated the tallest trees in a 1-km 2 area (225 trees in total) and related their vertical and horizontal dimensions to the scaling properties observed at the canopy scale.This allowed us to examine whether links between aridity and fractal scaling deviations could be mediated by hydraulic limitations on tree size.

| Airborne laser scanning data
We compiled ALS data acquired at nine of Australia's Terrestrial Ecosystem Research Network (TERN) SuperSites (Librarian, 2021;Karan et al., 2016).Each site comprises a representative 5 × 5 km 2 patch of continuous, relatively homogeneous forest cover and was scanned with the same RIEGL Q560 airborne laser scanner between 2012 and 2015 at a nominal resolution of 0.3 m or average pulse densities of ≥20 m −2 and a maximum scan angle of 20°.Together, the TERN SuperSites span a wide range of ecosystems, including arid and semi-arid shrublands and woodlands (Alice Mulga, Chowilla), frequently burned tropical savannas (Litchfield) and tropical rainforests (Robson Creek), temperate woodlands (Credo) and forests (Zigzag Creek), dry sclerophyll forests (Rushworth Forest), as well as tall wet sclerophyll forests (Warra, Watts Creek).

| Terrain and canopy surface modelling
For each of the nine sites, we downloaded the raw point clouds from the TERN website (cf.Data Availability section), processed them with a standardized pipeline (Fischer et al., 2023) and constructed high-resolution (0.5 m) gridded digital terrain models (DTMs) and digital surface models (DSMs) of each landscape.We then obtained canopy height models (CHMs) by subtracting the DTMs from the DSMs (Figure 2).All ALS processing was done in R 4.2, using the lidR (Roussel et al., 2020) and terra (Hijmans, 2023) packages, and with the LAStools software suite (Isenburg, 2023).Results have been stored on Zenodo (Fischer & Jucker, 2023).Additionally, we relied on the data.table(Dowle & Srinivasan, 2023), ggplot2 (Wickham, 2016), viridis (Garnier et al., 2022) and gridExtra (Auguie, 2017) packages for data manipulation and visualization.
Various algorithms can be used to construct CHMs, and the way they fill or do not fill small canopy gaps will impact the assessments of forest structure, particularly when variation across scales is concerned.Based on fractal theory, for example, an even and smooth canopy surface will look like a 2D plane, but if the same surface is speckled with tiny gaps that extend down to the ground, it will look like a 3D volume.To account for algorithm-induced biases and uncertainty in our results, we created four different CHMs, including a simple rasterization of the highest returns per 0.5 m grid cell (hereafter CHM highest ), a Delaunay triangulation of the ALS first returns (CHM tin ), a surface interpolated via the so-called 'spikefree' algorithm (Khosravipour et al., 2016) that adjusts canopy reconstruction to the scan's pulse density (CHM aspikefree ) and a surface with a customized 'local spikefree' algorithm (CHM lspikefree ) that uses the same adaptive smoothing algorithm, but adjusts its parameters to within-scan pulse density variation (Fischer et al., 2023).To remove extreme outliers and potential noise from the CHMs, for each site, we calculated the 99.99th percentile of canopy height and set all pixels that were 10% higher than that threshold to NA.Note that all four CHM types were used throughout, with metrics averaged across them, but, for simplicity, we used only CHM lspikefree for visualizations.

| Assessing general scaling properties of forest canopies
To characterize how forest structure changes across spatial scales, we computed variation in canopy height at four different scales (2, 20, 200 and 2000 m), as this is one of the most common methods to assess how surfaces change across scales (Torres-Pulliza et al., 2020).
To quantify height variation at each of these scales we used the standard deviation of height (sd height ), as this was more robust to outliers, discretization effects and missing data than the full height range (cf. Figure S2 for discretization effects).To put scale changes into a wider ecological context, we computed each site's maximum canopy height (99th percentile) as well as a height-normalized rumple index rumple norm10 .The latter was defined as canopy surface area divided by planar area, where the site's maximum canopy height was rescaled to a standard 10 m.Rescaling removes the dependence on canopy height, and we chose 10 m (instead of, e.g. 1 m) to obtain an easily interpretable ratio.Finally, we linked sites to estimates of annual precipitation and site water balance (SWB, a measure of aridity) obtained from the 1-km gridded CHELSA climatology (1981-2010;Brun et al., 2022;Karger et al., 2017).Throughout our analysis, sites will be shown in ascending order of SWB (i.e. from driest to wettest, Table 1).As a convenient shorthand, we will refer to the four wettest and closed-canopy sites as 'wet' forest types (Robson Creek, Zigzag Creek, Watts Creek and Warra), and to the others as 'dry' forest types (Alice Mulga, Chowilla, Credo, Litchfield, Rushworth Forest).

| Assessing self-similarity of forest canopies across scales
Whether objects are self-similar or self-affine is rarely tested, and no formal criteria seem to exist.In practice, methods to calculate fractal dimension involve the fitting of linear regression lines to log scaled quantities (Zhou & Lam, 2005), and the existence of fractal-like behaviour should be testable by how close the data are to a linear relationship on log-log scales.A few studies, for example, report the explained variance (Guzmán Q. et al., 2020), or whether fitted regression slopes match the slopes between the endpoints of the data (Torres-Pulliza et al., 2020).However, log scale transformations mask substantial variation in the original, typically biologically more relevant quantities and neither high R 2 values (≥0.9) nor comparable slopes between endpoints and regression estimates rule out deviations from linearity (Loehle & Li, 1996).Crucially, R 2 will depend on  Dubuc et al., 1989).Therefore, we here employed two custom heuristics that did not suffer from this problem.
First, we compared observed changes in height variation at 20 and 200 m (sd 20m and sd 200m ) to expectations under fractal scaling (assuming ideal log-linear relationships between sd 2m and sd 2000m ).
Second, we obtained fractal dimension (D) estimates from a number of methods and, for each of them, calculated an index of deviation from fractal scaling (ΔD).ΔD is simply the difference between D estimated from the lower half of the range (D lower ) and D estimated from the full range (D), with the midpoint of the range determined on log scales (Figure 3).Note that, by construction, this is equivalent to the difference between D and D upper .As fractal dimension D has lower and upper bounds of 2 and 3, respectively, ΔD would be correlated with D, with the largest deviations possible at D = 2.5, and the smallest at the upper and lower bounds.To remove this effect, we normalized ΔD to a value between −1 and 1 by dividing by the maximum possible deviation.The maximum TA B L E 1 Structural, floristic and climatic differences between the nine TERN SuperSites, ordered from driest to wettest.Maximum canopy height (99% percentile), canopy cover (estimated at 2 m above-ground) and rugosity (normalized rumple index rumple norm10 where maximum canopy height is rescaled to 10 m) were calculated from canopy height models derived from ALS data.
Mean elevation above sea level was calculated from the digital terrain models.Mean annual precipitation (MAP) and site water balance (SWB) were obtained from the 1-km gridded CHELSA climatology 1981-2010 (Brun et al., 2022;Karger et al., 2017).Note that canopy height and cover are strongly related to SWB, with moist forests generally much taller and denser than dry forests and shrublands.F I G U R E 3 Illustration of how the deviation from fractal dimension D is calculated, with an example based on a variant of the height variation method (D sd ).For a given fractal or non-fractal surface (continuous red line), we construct the ideal fractal scaling line from the endpoints (continuous orange line), estimate D and then repeat the process for the lower half of the range (from the lowest point to the midpoint on log scales, dashed red line).ΔD is calculated as the difference between D lower and D, divided by the maximum possible difference ΔD max .Upper and lower bounds both for the overall fractal scaling exponent and the scaling exponents of the lower and upper half (2 and 3 in both cases) mean that ΔD max can be simply calculated as D - 2 (D < 2.5) or 3 - D (D ≥ 2.5).The fractal dimension estimators we compared are visualized in Figure S3 and included:

Site
a. D hvar : The original variation method (e.g.Dubuc et al., 1989;Torres-Pulliza et al., 2020).This method calculates the average range (max-min) in canopy height for grid cells of increasing cell lengths, and then obtains D hvar = 3 - S hvar , where S hvar is the (positive) slope of height variation plotted against grid cell length on log scales.
b. D sd : An analogue of the variation method, but using the standard deviation in height per grid cell instead of the full height range, with D sd = 3 - S sd , where S sd is the (positive) slope of height variation plotted against grid cell length (both log-transformed).
c. D cube : The cube counting method (e.g.Zawada et al., 2019) is a 3D extension of the box counting method, with D = -S cube , where S cube is the (negative) slope of the number of 3D cubes that are needed to cover the 3D surface at a specific cube length, plotted against cube length on log scales.To detect self-affinity (and not only self-similarity), the vertical extent was rescaled to match the horizontal extent.
d. D surfarea : The surface area method, a custom version of the prism method (Clarke, 1986).This is a method that calculates the canopy surface area at various scales, as if the variation between the measurement points did not exist, and estimates D surfarea = 2 - S surfarea where S surfarea is the (negative) log-log slope of surface area against measurement scale.The original method uses four corner points to calculate the surface area, but presents several difficulties in accurately characterizing fractal dimension (De Santis et al., 1997).Here, we employed a method that uses a centre grid cell plus eight adjacent grid cells, as implemented in the sp R package (Bivand et al., 2013;Jenness, 2004).As in D cube , we rescaled the surface to match the horizontal extent.
e or falls below 50%.We also note that this method, while expected to work well for fractal surfaces (Loke & Chisholm, 2022), discards a lot of 3D information present in surface models and thus may not be able to fully characterize deviations from fractality.
To compute fractal scaling estimates, it is easiest to proceed in dyadic instead of decadic steps (usually requiring approximate integer steps in between), so we computed all statistics at the following step sizes: 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 m, and calculated D both across the full range and from 2 to 64 m (D lower ).
Although the exact scale range was 2-2048 m, for simplicity, we will refer to this as a scale range of 2-2000 m throughout the article.
To make methods comparable, we had to resolve ambiguities in how scale was defined.For instance, the step size of the variation method (D hvar or D sd ) is typically defined via ε or half the tile size at which it is measured, whereas the step size for the cube counting method (D cube ) is typically defined via the side length of the cube (i.e. the whole tile size).Further ambiguity comes from the fact that scale on a raster grid can be taken as the whole extent in grid cells, or as the distance between grid cell centre points (whole extent - 1).For true fractals, such differences matter little, but making estimators consistent is crucial for the estimation of non-linearities and breakpoints.In particular, for D surfarea , we defined the baseline scale such that surface area at the smallest scale (2 × 2 m 2 ) was computed using all 16 grid cells (0.5 × 0.5 m 2 each).We deemed this the closest analogue to how height variation is computed across all 16 grid cells for D hvar or how a cube is fitted around all 16 grid cells at 2 m resolution.
Correspondingly, at each larger scale, surface area was then also To test the validity of this approach, we also fit OLS and SMA regression models with the lmodel2 package (Legendre, 2018), and assessed consistency between regression slopes and endpointbased slope estimates via correlation.

| Comparing canopy height models and ideal fractal surfaces
In theory, any systematic deviation from ΔD = 0 would indicate that the underlying surface is non-fractal.However, data are noisy, and, for fractals, there may be additional biases and errors due to estimation methods as well as discretization and sampling extent (Loke & Chisholm, 2022).To separate meaningful deviations in ΔD from expected deviations, we simulated fractal Brownian motion and generated corresponding 2D landscapes on an 8192 × 8192 grid cell extent.We used both spectral synthesis and random midpoint displacement to simulate fractal surfaces (Fournier et al., 1982;Saupe, 1988), the former implemented in GRASS GIS (GRASS Development Team, 2017) and called via rgrass (Bivand, 2023), the latter in a custom implementation based on Saupe (1988).We chose the implementation where random variation is added not only to the midpoints but also to all grid cells at each step, and extended this procedure by adding random variation for three more iterations without creating new midpoints.This has little impact on the overall surface, but reduces the drop-off in variation across scales as the resolution of the grid is approached.To assess bias and uncertainty in estimators, we simulated fractal surfaces from D = 2.1 to D = 2.9 in steps of 0.1 for both surface generation methods.To compare estimators with CHMs, we calculated additional surfaces with D between 2.725 and 2.975 in steps of 0.25, as preliminary results suggested this as the most likely range of canopy surface D.

| Comparison of scale changes with limits on tree size
To assess how limitations on tree growth and allometry affect structural changes across scales, we divided the CHMs at all nine TERN sites into 25 patches of 1 × 1 km 2 and, in each of them, manually delineated the tallest tree per km 2 (25 per site, 225 in total).In each patch, we found the tallest trees by first finding the 2.5 × 2.5 m 2 window with the highest mean canopy elevation and then used the draw() function in the R terra package (Hijmans, 2023) to carry out delineations and saved them to polygon shapefiles.Since the tallest trees were often emergent trees, they could be delineated with high confidence.However, we note that at some sites such as Alice Mulga, Watts Creek or Zigzag Creek, canopies were sometimes highly connected, so over-or undersegmentation cannot be ruled out.To derive robust estimates of tree dimensions, we first excluded all pixels below 50% of the 25th height percentile (e.g.within-crown gaps, or mismatches between delineation and crown perimeter) and then calculated tree height (95th percentile), crown area (ca) and its derivative crown radius (cr) = ( √ ca ∕ ), and an index of crown smoothness (ca divided by the interquartile range of height, dimensionless).This is motivated by the idea that smoother crowns would be expected to have either less variation in height or spread their variation over a larger crown area.
We compared the correlation of the three summary tree properties (height, crown radius, crown smoothness) of the tallest trees per km 2 with statistics obtained at the landscape level (height variation across scales and deviation from fractal scaling) and expected clear correlations for tree size due to the limits it imposes on horizontal canopy variation.We also expected a positive link between crown smoothness and losses in height variation, as increased crown smoothness means reduced height variation at the smallest scale.To put results into context, we computed all potential permutations of how canopy metrics could be linked to tree metrics by shifting trees between sites (9! = 362,880 permutations in total) and compared average absolute correlations across the permutations to the actually observed ones.

| How does forest structure change across scales?
Canopy height variation displayed highly non-linear changes across scales, both untransformed and under logarithmic transformation (Figure 4a,c, and Table 2).From 2000 down to 2 m and across all CHM algorithms, forest canopies lost 51%-79% of their height variation.However, only 6% of maximum height variation was lost between 2000 and 200 m (95% range 0%-16% across sites), and 12% between 200 and 20 m (6%-20%).Instead, the majority of height variation-47% (40%-55%)-was lost over just 18 m (from 20 to 2 m).Height variation loss across the entire range of scales (from 2000 to 2 m) showed a strong negative correlation with rugosity across sites and CHM algorithms (rumple norm10 , r = −0.92, Figure S4) and occurred earlier and more gradually at taller, wet sites than at shorter and more open dry sites (mean loss of 11% in height variation from 2000 to 200 m, compared to 2% at dry sites, Table 2).

| Do forest canopies follow fractal scaling rules?
We found no evidence for fractal self-similarity across any of the nine forest sites.First, ideal fractal scaling assumptions (i.e.linear changes in variation on log scales) were clearly and systematically violated across all sites and CHMs (Figure 4b,d, Table 2).Under fractal scaling, the majority of the height variation loss would have been  S1, Figure S5).Based on these best estimators, fractal surfaces had an average ΔD = 0.04 (±0.19), clearly distinct from values obtained for the CHMs (ΔD = 0.83 ± 0.12; Figure 5).

F I G U R E 4
Observed and theoretical changes (i.e.expected under ideal fractal scaling) in canopy height variation across scales for the nine TERN SuperSites.The left row shows how height variation (standard deviation of canopy height in m) changes across three orders of magnitude from 2 to 2000 m, on original scales (a) and on log-transformed scales (c).The right column (panels b and d) shows the corresponding plots for how height variation would change across scales for an idealized fractal surface, that is, when a straight line is fitted to the endpoints on log scales.Sites are coloured by aridity from driest to wettest, and the underlying canopy height model algorithm is CHM lspikefree .Note that the observed curves on log scales (c) are also non-linear across the lower half of the range (2-64 m) and an intermediate scale range (e.g.10-100 m).This gives already a visual indication that overall deviations cannot be fully explained by artefacts due to resolution or sampling extent (Loke & Chisholm, 2022).

| How do changes across scale relate to biological limitations of tree growth?
The size of the tallest trees varied widely across sites, from heights of 8.1 m (Alice Mulga) to 90.2 m (Warra) and crown radii from 5.3 (Alice Mulga) to 31.8 m (Robson Creek).As expected, maximum height was tightly coupled to cross-site shifts in water balance (r = 0.91), with a corresponding, but weaker increase in crown size (r = 0.53).Tree crown properties were also clearly correlated with scaling properties (Figure 6), showing moderate to good correlations with height variation loss across spatial scales (r = 0.49 and 0.78) and stronger correlations with deviations from fractal scaling (r = − 0.89 and − 0.7).Specifically, the strongest losses in height variation and smallest deviations from fractal scaling occurred in wet forests dominated by tall trees with large crowns.In contrast, crown smoothness was only linked to height variation loss (r = 0.73) and not ΔD (r = 0.12).Nevertheless, overall correlations were clear (mean absolute r = 0.62) and unlikely to be due to chance: When comparing all potential permutations of the data, more than 99.95% had lower average correlation values (Figure S7).

| One tree does not a forest make…
Quantifying structural complexity in forests is a major research objective in ecology (Atkins et al., 2018;Ehbrecht et al., 2021;Seidel et al., 2019) and we finally have enough high-resolution data covering broad spatial scales to make this possible (Lines et al., 2022).Fractal theory provides an attractive approach to summarize structural complexity, as it promises to reveal hidden generalities in seemingly heterogeneous systems and to reduce complexity to a single number, a fractal dimension (Brown et al., 2002;Torres-Pulliza et al., 2020;Zeide, 1991).Our results show, however, that assumptions about fractal self-similarity of forest canopies are not warranted.Across all nine TERN SuperSites-spanning dry shrublands, tropical rain forests and 90 m tall Mountain Ash stands-canopy surfaces consistently and systematically deviated from the assumptions of fractal scaling.If forest canopies had fractal properties, height variation would gradually decrease across scales, but this is not what we observed.Height variation remained largely stable from 2000 down to 20 m across forest types and then decreased abruptly, irrespective of estimation technique (Figure 4a).This means that scaling assumptions were violated both at the landscape (2-2000 m) and stand scale (2-200 m) and self-similarity of forest canopies cannot be assumed (cf.also Figure S8 for a visualization of self-similarity in fractal surfaces and its lack in forest canopies).This is consistent with breakpoints for many forest structure metrics around 10-100 m (Atkins et al., 2023;Kane et al., 2011).It should also generalize to other sites, since the TERN SuperSites were explicitly chosen to be representative of their surroundings, with a wide coverage of ecosystem types and representative within-site variation (Karan et al., 2016).
As deviations from fractal scaling in canopies exceeded those from simulated fractal surfaces, we also expect our results to be robust to fractal estimation techniques as well as potential biases due to image discretization and sampling extent (Loke & Chisholm, 2022).There were exceptions, most notably for the two estimators D hvar and D box , which suggested some overlap between canopies and simulated fractal surfaces, but these estimators still showed the same qualitative patterns and were among the worst-performing in identifying fractality (Table S1).We can therefore confidently conclude that the tree is not the forest, and the stand is not the landscape.

| … but we can see the forest from its many trees
Even though there was little evidence for fractal scaling, we did find several generalities in how forest structure changed with scale across  the nine TERN SuperSites.Qualitatively, all forest canopies exhibited the same saturation in height variation at large scales, irrespective of forest type, with great consistency in deviation from fractal scaling rules.Furthermore, differences in variation across sites were strongly linked to water availability, with taller and wetter forests changing more across scales-and more gradually-than shorter and drier forests.To give an example, at the driest site (Acacia woodland at Alice Mulga), height variation decreased by only 6% between measurement scales of 2000 and 20 m, and then dropped by 45% from 20 to 2 m.In contrast, at the wettest site (Mountain Ash forest at Warra, Tasmania), variation decreased by 30% from 2000 to 20 m, and then by an additional 40% when going from 20 to 2 m (Table 2).
It is not surprising that there should be an aridity-linked pattern both in how height variation changes across scales and in deviations from fractal scaling given the strong control that water availability exerts on canopy structure across the globe (Ehbrecht et al., 2021;Tao et al., 2016).What was surprising was that deviations from fractal scaling increased with aridity and did not reverse even in the driest ecosystems or under known fractal disturbance agents.Arid shrublands and woodlands are strongly shaped by patchy local nutrient pools, large-scale disturbances such as fire, and facilitative process that are usually multiplicative in nature, all of which should be present at our sites and make fractal-like clumping more likely (Green, 1989;Hastings et al., 1982;Kéfi et al., 2007).In part, this deviation from expectations in dry ecosystems could be due to lower species diversity and more conservative trait strategies (Anderegg et al., 2021;Peters et al., 2021) and the more homogeneous canopy units this entails.But even the presence of frequent fires, which are usually associated with fractal behaviour (Malamud et al., 1998), was not linked to fractal-like canopy surfaces in a seasonally wet savanna site (Litchfield, Figure 5).
Conversely, it seems that fragmentation due to human disturbance at the wettest site (Warra, Tasmania)-as seen in large, cleared areas in the canopy height model (Figure 2)-mirrored the fractal-like patterns of individual tree crowns at smaller scales and thus rendered the canopy more similar to simulated fractals (Figure S8).
However, as strong correlations with maximum canopy height and crown radius show, the overriding effects are more likely due to tree-level factors, particularly the physiological limits on how tall and wide trees can grow (Jucker et al., 2022;Niklas, 2007).For one, limits on tree height already cause deviations from expected scaling rules at the tree level (Lines et al., 2012) and mean that the canopy surface quickly saturates in height and cannot fold freely across scales like a true fractal would (Birnbaum, 2001).Second, crown radius determines the extent of within-crown variation and thus implicitly sets the scale at which disturbance effects become visible.Above the scale of individual trees (>0.1-1 ha), we start to measure the equilibrium forest structure (i.e.variation between different successional phases), while below this scale, the canopy surface mosaic disintegrates into large canopy trees, the gaps they leave behind when they die and the clusters of small, regenerating trees that follow (Shugart et al., 2010).All of these are comparatively homogeneous entities whose structure mostly depends on tree size.
Knowing the structure of individual trees can thus provide a lot of information on the overall surface structure of forest canopies, which is in agreement with studies that have found little information in fractal dimension beyond large trees or groups of a few trees (Liu et al., 2022;Seidel et al., 2019) and underlines the general importance of large individual trees for the assessment of forest structure (Lutz et al., 2018;Meyer et al., 2018).D cube (ds) , also showed good predictivity of fractal dimension (Figure S6, Table S1), but was systematically shifted towards negative deviations from fractality and thus more likely to misclassify fractals (Figure S5, Table S1), so we did not include it among the best estimators.Palmer, 1988;Solé & Manrubia, 1995), and digital surface models are routinely used to infer scaling relationships both in canopies (Jucker, 2022;Staver et al., 2019) and ecosystems more generally (Loke & Chisholm, 2022;Torres-Pulliza et al., 2020).Compared to laser-derived point clouds, canopy height models also have the advantage of avoiding occlusion effects, which should make them more robust to acquisition details such as flight altitude, scan angles or pulse densities (Asner et al., 2013).However, by not representing understorey vegetation, they do not fully account for the complex 3D space-filling of trees and differ from how most forest-dwelling organisms will experience their environment.This could affect our results, particularly at the smallest scales.Replicating our approach using high-resolution 3D data from terrestrial laser scanning would thus strongly complement our analysis.Crucially, this would require extents of at least 1 ha, and great care should be taken to account for occlusion and pulse density variation or other instrument-induced properties (Ehbrecht et al., 2021).

| A multidimensional problem
It is possible that such a study would reveal more fractal-like patterns at subcanopy scales, but we would expect many of our findings at large scales to also hold under a full 3D assessment.First, the basic arguments hold irrespective of how fractality is assessed: Most forest landscapes do not look like the trees that grow in them (Figure 1), and their structure is shaped by a wide range of external factors many of which are not fractal or do not combine to produce fractality.Second, as evidenced in this study, fundamental limitations on tree dimensions will also impose limitations on 3D space filling that are not comparable to how theoretical fractals fold in space.
Third, fractal self-similarity is a property that can often be found even when reprojecting to a lower dimension (e.g. by flattening a 3D tree into a 2D representation or by cutting a slice out of a 3D surface; Loke & Chisholm, 2023).The consistent lack of self-similarity in canopy surfaces across landscapes should thus be a good indicator of potential violations in full 3D canopy structure at the same scales.

| Implications for the analysis of ecological complexity
Our results have a number of implications for studies of complexity, both in forest canopies and beyond.For one, they add to the F I G U R E 6 Correlations between changes in forest canopy structure across scales and maximum tree size, based on the crown characteristics of the tallest trees per km 2 across the nine TERN SuperSites.Tree crown outlines were derived via manual delineations and then used to extract canopy height values, ignoring small canopy gaps within crowns and low pixels along crown edges.Shown are height variation loss between 2 and 2000 m (%) and deviations from fractal scaling, each plotted against maximum tree height, measured as the 95th percentile of height pixels within the tallest tree crowns per km 2 (panels a and d), maximum crown radius of the same trees (b and e) and crown smoothness (crown area divided by interquartile range of height variation, c and f).Each point represents one site and the average across 25 trees.ΔD was computed across the four best fractal dimension estimators, as in Figure 5. Sites are coloured in order of aridity, from driest to wettest.Dashed lines show linear regression lines, grey areas the corresponding 95% confidence intervals.Crown smoothness evidence that fractal dimension, taken as a scale-invariant indicator of complexity, is a problematic measure, both in terms of how it is computed and how it is interpreted (De Santis et al., 1997;Halley et al., 2004;Huang et al., 1994;Loke & Chisholm, 2022;Madin et al., 2023).First, our results show that the theoretical foundation for using fractality in forest canopy assessments is weak, since its core assumption (scale-invariance) is violated across a wide range of forest types.Second, D also seems problematic as a heuristic.Across surfaces, biases in estimates due to scale, CHM derivation and estimation method were quantitively similar or exceeded between-site variation in D, and even occurred for supposedly ideal fractal surfaces.This echoes Huang et al. ''s (1994) statement that 'information about the "correct" D is somewhere, if only we knew where to look'.Finally, all deviations from fractality were surprisingly hard to diagnose, even when visualizing the underlying data, as the data appeared near-linear for most estimators, despite highly variable D (this is a particular problem with the often-used box or cube counting approach, as seen in Figure S3).The reason is likely that many common estimators, such as cube and box counting, are based on mathematical arguments and log scale transformations that are correct only for true fractals.For non-fractal objects, such transformations implicitly also transform the accompanying error structures, an issue that does not seem to have been adequately addressed so far.
All of this does not mean that fractal dimension or general scaling assessments should be abandoned.D may very well be a useful measure at smaller scales or for forest properties other than the canopy surface (Hastings et al., 1982;Malamud et al., 1998).
However, our results do mean that fractality or scaling relationships should not be pre-supposed and great care should be taken to verify assumptions when linearity on log scales is expected.
This harks back to reviews on fractality (Halley et al., 2004;Loke & Chisholm, 2022) and similar debates about log-inferred scaling relationships and how to correctly measure them (White et al., 2008), be it in size distributions, fractal scaling or allometries.We here trialled two heuristics that could be extended to other use cases, such as the fractal deviation index ΔD, which can be applied across estimators, and a simple predictive test, which predicts data from ideal scaling lines fitted on log scales and then assesses the overlap between predictions and observed data on the original scales.In both cases, we calculated the heuristics from lines fitted between endpoints, which was justified by the lack of noise in our scaling curves, but the concept should not be difficult to extend to regression-based estimates.
As we expect assumptions of fractal scaling to be violated more often than not, it may also be indicated to complement or replace fractal dimension by a conceptually simpler and more easily understandable quantity.To assess canopy roughness, this could, for example, take the form of height-normalized rugosity estimates, such as the rumple norm10 index proposed in this study, which correlated highly with D, but is simply the canopy surface area divided by planar area for canopies rescaled to a standard 10 m maximum height.
Alternatively, if we use the height variation method as inspiration, D quantifies the (log-transformed) ratio between height variation at the largest scale and height variation at the smallest scale.As we showed in this study, this can be easily reformulated on original scales as percentage change in average height variation between 2 and 2000 m (and, equivalently, for other scale ranges).Such a metric offers the same information as fractal dimension and deviations from it, but without requiring interpretation of log scale relationships and without the many implications regarding cross-site comparability or scale-free properties.

| Towards a comparison across systems and scales
The value of fractal scaling lies in the mechanisms that can be inferred from it (Brown et al., 2002) and in the simplicity of providing a single number that could be compared across systems (Torres-Pulliza et al., 2020).But even in its absence, we do not have to give up on a more general framework of ecosystem complexity, and our results point a way towards it.Most biological and biologically relevant surfaces should experience saturation effects similar to canopy height variation, because, ultimately, physical and physiological factors will impose size limits.This is true for surfaces made up of organisms such as coral reefs or forest canopies due to limitations on organism size and light availability, but likely also holds true for more fractal-like surfaces such as mountainous terrain, which from a far enough distance starts to look like a flat plane.As a result, if we measure systems at sufficiently large scales, we can define the scale at which saturation occurs-this is known as the 'sill' of a variogram in geostatistics-and then use this as a system-intrinsic reference threshold from which changes across scales are measured.This can take the form of either measuring scale shifts across two to three orders of magnitude down from the saturation threshold or by calculating the scale at which a given proportion (e.g.50% and 75%) of the variation measured at the saturation threshold is lost.Both measures can be computed for any surface-provided that saturation occurs-and should provide relevant biological information without requiring complex assumptions about fractal scaling.They also provide a guideline on the scale range that should be sampled.In this regard, our study shows that while forest canopies deviate from the assumptions of fractal scaling, they do so in a way that is linked to environmental determinants and tree architecture and thus provide a promising way towards a more general framework for measuring and comparing the ecological complexity of ecosystems.

AUTH O R CO NTR I B UTI O N S
Fabian Jörg Fischer and Tommaso Jucker conceived of the idea.
Fabian Jörg Fischer carried out the analysis and led the writing of the paper.Tommaso Jucker contributed substantially to revisions.

ACK N O WLE D G E M ENTS
We acknowledge TERN and Airborne Research Australia (ARA), who collected the airborne laser scanning data and made them openly available under a CC-BY 4 licence, as well as the

F
Canopy height models of the nine TERN SuperSites (each 5 × 5 km 2 ), mapped at 0.5 m resolution and ordered from driest to wettest.Colours are rescaled to show the full extent in height at each site.Tick marks are not rescaled and correspond to 8 m increases in canopy height.Thus, trees from the first two sites (Alice Mulga and Chowilla, below 16 m) would fit into the lower two rungs at the last two sites (Warra and Watts Creek, both reaching beyond 80 m).See Figure S1 in the Supporting Information for a zoomed-in image of the largest trees at each of the sites.fractal dimension estimators, with some estimators creating the impression of strong linearity by default and only revealing non-linearity after transformation (e.g., the height variation estimator in Torres-Pulliza et al., 2020 makes deviations more obvious than the original estimator in possible deviation is 3 -D for surfaces with D greater or equal 2.5, and D -2 for surfaces with D below 2.5.A perfectly fractal surface would thus have ΔD = 0, while a maximally deviating surfaces would have either ΔD = 1 (highest D at small scales) or ΔD = −1 (highest D at large scales).
computed from 16 equally spaced grid cells (spaced at 1 m resolution for a scale of 4 m, at 2 m resolution for a scale of 8 m etc.).This definition intuitively deals with the issue of finite resolution in raster imagery(Huang et al., 1994;Loke & Chisholm, 2022): Even if the scale gets coarser, the resolution remains the same.It also suggested a general improvement for fractal estimators: In analogy to D sd , D hvar and D cube , we calculated D sd (ds) , D hvar (ds) and D cube (ds) (downsampled estimators), where we used the same methods as before, but applied them only to 16 equally spaced points at each scale, which should intuitively account for discretization effects.We used a fixed-grid approach for D sd , D hvar , D cube and D box , and a moving grid for the downsampled estimators and D surfarea , since we expected the latter ones to be more sensitive to grid position.The number of grid positions was all possible ones up to 1024 (32 * 32) positions at the 64 m scale.Beyond this scale, to reduce computational costs, we kept the number of grid positions at a constant 1024, with the grid position shifting in equal steps across tiles.For each landscape, we focused on the 4048 × 4048 m 2 central square of the CHM (corresponding to 8192 × 8192 grid cells of 0.5 m resolution).Any cells outside this central square were discarded to avoid edge effects, as were any patches with more than 5% NA cells per tile.Methods that rely on exhaustive sampling of the grid (e.g.D cube or D cube (ds) ) were rescaled accordingly.For all methods, we report slope estimates and D based on the endpoints, mirroring recent work on surface complexity of coral reefs (Torres-Pulliza et al., 2020).
expected at large scales, with 30% between 2000 and 200 m (range of 22%-41% across sites) and only 14% at 20-2 m (13%-15%)-the exact opposite of what we found.Fractal scaling assumptions were also violated across more restricted ranges: If we assumed linearity only across 2-200 m, then larger scale losses (200-20 m) were again much smaller (13%) than expected under fractal scaling (47%).Violations of scaling assumptions were observed regardless of the estimator used to quantify fractal dimension (Figure5, FigureS5).Across the eight estimators we tested, the mean D calculated from the CHMs of the nine study sites was 2.81, varying from 2.73 (Credo) to 2.87 (Alice Mulga).However, without exception, D was smaller across the lower half of the range (measured on log scales from 2 to 64 m), with D lower = 2.69, varying from 2.53 at Credo to 2.78 at Alice Mulga.Normalizing by the maximum possible deviation, this yielded a mean ΔD of 0.70 (±0.24SD), with shorter forests at dry sites showing larger deviations (ΔD = 0.75) than taller forests at wet sites (ΔD = 0.64).By contrast, when the same fractal dimension estimators were applied to simulated fractal surfaces, we found that ΔD was much closer to 0 on average (ΔD = 0.05 ± 0.36 for surfaces generated from midpoint displacement, ΔD = 0.28 ± 0.24 for surfaces generated from spectral synthesis), indicating that the deviations from fractality in the CHMs were not methodological artefacts.This was even more evident when we only retained the best-performing estimators.D sd , D sd (ds) , D hvar (ds) and D surfarea reliably classified fractal surfaces as fractal, with the IQR of ΔD covering 0 (FigureS5) and noticeably lower RMSE than D hvar , D cube and D box (FigureS5) and lower rates of misidentification of fractal surfaces than D cube (sd) (Table

TA B L E 2
Canopy structural indices across the nine TERN SuperSites.Shown are: average standard deviation of canopy height across tiles of 2-2000 m (in m), predictions for 20 and 200 m based on log-linear scaling between 2 and 2000 m (grey, in brackets), fractal dimension estimates D sd (based on standard deviation of height) across the full, lower and upper half of the range, and the corresponding deviation ΔD sd , which is the absolute difference between D sd and D sd (lower) , normalized by the maximum difference 3.0 - D sd (all dimensionless).All values are averages across four different types of canopy height models.
Loke and Chisholm (2022) a specific box length, plotted against box length on log scales.Here, we employed the exact algorithm fromLoke and Chisholm (2022), but replaced the median by the mean in slicing the canopy.The median is problematic for vegetation characterization in open forests where canopy cover approaches

Site sd 2m sd 20m observed (fractal) sd 200m observed (fractal) sd 2000m D sd full D sd lower D sd upper ΔD sd
Expected deviations ΔD from fractal scaling assumptions compared to observed ΔD for the nine TERN SuperSites.ΔD was calculated as difference between slope estimates over the entire scale range (2000 m down to 2 m) and slope estimates over the lower half of the scaling range (on log scales, 64 to 2 m), based on the best-performing D estimators (D sd , D sd (ds) , D hvar (ds) and D surfarea ) and normalized by the maximum possible deviation for a given D estimate.