**3.1 The sampling grid design of the French NFI**

To illustrate with concrete examples the use of the interpenetrating grids, the French NFI survey is being presented here. It created its base (1 x 1 km) grid in 2004, as it changed its sampling design from a 10-year periodic to an annual sampling that was first presented in Vidal et al. (2005). It borrows to the design of the FIA forest inventory survey (Roesch and Reams 1999).

As the initial choice for the NFI was to meet the periodicity of 10 years, all the nodes of the grid being sampled within 10 years, 10 interpenetrating grids were devised. There are 8 pairs of solutions to Eq. 1 that meet this objective of fractionation into 10 panels: a= +/- 3 and b = +/- 1, and conversely, a= +/- 1 and b = +/- 3 (Figure 3). The divider that allows to reduce the grid density (i.e., the levelling factor) was set to m = 2: the number of cells is divided by two from one level to another, while the cell area doubles.

Of note, each sequence of 5 successive subgrids (1-5 and 6-10) ensures a regular systematic spatial cover (Figure 3), stressing that this grid also supports the 5-yr periodicity. It thus offers options for forest inference at annual, 5-year or 10-year periods, and brings flexibility in the reporting periodicity. Another interesting property of the square grid emerges when the fractionation is made for 5 panels: as a result of the balanced nesting, any node belonging to a given panel is surrounded by those of the four other panels (Figure 2b). In that situation, the time-periodicity is geographically translated into systematic patterns of five nodes as the central node of that repeated pattern also happens to be the mid-period node.

At the time when the sampling design was conceived, the sampling units obtained from the grid were meant to be measured only once. But in 2010, the choice was made to remeasure all sampling units after 5 years for accurate wood fluxes (e.g., mortality, harvests, recruitment) estimation. The initial grid was consequently split into two translated subgrids having each a periodicity of 5 years. The subsetting levels of the adjacent units from panels t and t+5 were consequently chosen to be identical (Hervé et al. 2016), ensuring that, at any given year and subsetting level, the units of panels t and t+5 would be adjacent, hence reducing time travel for field crews.

**3.2 Using the grid approach for implementing spatial changes in the sampling intensity**

**3.2.1 Optimizing phase sampling**

Double sampling for stratification is a popular and efficient sampling and estimation method, which implies using nested samples of different sizes (Rao 1973, Cochran 1977). Here, interpenetrating grids bring a decisive advantage whereby grid levels offer a direct support to the construction of nested balanced samples.

The French NFI brings an example of such two-phase sampling design based on grid levels, whereby the first phase is a large sample of photo-interpretation plots, while the second is a nested subsample of plots used to produce field measurements (de Vries 1986, Mandallaz 2007). The nested levels of the square grid offer a straightforward means of subsampling the units of the first phase: typically, the first phase is based on the finest grid level (with cells of 10 km2 in the case of the French NFI) while the second phase is a based on higher levels with cells of 20 km2 or 40 km2, depending on the subsampling intensity desired.

The same principles can be used to accommodate further sampling phases (e.g., three phases in the Italian NFI, Tabacchi et al. 2005).

**3.2.2 Stratified sampling with unequal probability**

Grid levelling can also support a stratification-based variation of the sampling intensity within a given geographic domain, for samples drawn with unequal probability. This feature is used in the French NFI where the inclusion probability of field plots (second phase plots) varies according to the strategic significance of the vegetation type observed during the first phase photo-interpretation. Forest plots (forests being a category defined according to UN/FAO’s definitions) are hence subsampled with a probability of ½ (i.e., m=2) while other wooded lands (of lower importance) are sampled with a probability of ¼ (m=4), resulting in unequal probability sampling that preserves the spatial balance and all geometric properties of the samples.

Alongside this stratification, a geographic stratification is performed whereby lower levels are used across two geographic subdomains of the territory (*Landes* and *Mediterranean region*) which are known to have a smaller spatial variability and a lesser importance. The local reduction of the sampling intensity is a common practice in NFI surveys (Tomppo et al. 2010, Vidal et al. 2016, Bouriaud et al. 2020). The grid-based approach however ensures, beyond the spatial balance, the possibility to obtain a spatially systematic coverage of the entire territory at the lowest common density of several zones.

*3.2.3 ***Disturbance-driven variations in the sampling intensity**

Unanticipated deviations from a desired sampling may be dictated by the need to address a local survey of a forest windthrow or sanitary event (IFN 2009). Locally increasing the spatial sampling intensity while preserving the link with the initial sampling design, and some spatial balance, is possible by using the base grid at its original density. In this situation, levels can be used to construct the nested samples at varying sampling intensity. Of note, this forms an alternative to explicit disturbance-based designs (van Deusen 2000) which still has not received the attention it deserves.

In all surveys, the maximum possible sampling effort is however limited in practice by resources rather than statistical-related constraints. Increasing the intensity of the sampling in a given zone therefore needs to be compensated for by a decrease in other parts of the territory. The use of levels brings practical solutions that ensure the spatial and temporal coherence of the samples. The known location and number of potential units per grid panel and level enables an efficient planning, and allocation of resources.

**3.3 Temporal variations in the sampling intensity: an original development for trading systematicity with adaptability**

In annual surveys, unavoidable events (e.g., COVID-19 crisis of 2020 causing interruptions of the work) or unexpected budget-driven variations of the sampling effort across time may request changes in the sampling effort, challenging the initial sampling design. In the French NFI survey for instance, incremental reductions of the sampling effort have been implemented since 2012, with the second phase sampling effort being gradually decreased by 20% as compared to the initial one. These incremental reductions cannot be obtained from grid-levelling since switching from one level to the next reduces by 2 the number of nodes.

Nevertheless, and this is an original use of the present design, grid nesting also enables the design of a spatially balanced sample with a controlled and more continuous reduction of the sampling effort, by using high-level grids (low spatial density) as a means to exclude sampling units from lower levels (higher spatial density). The principle is to exclude units belonging to both a low level (ex. 40 km2) and a high level of the grid (ex. 320 km2). For instance, in 2020 the French NFI reduced by 14% its sampling effort for forests, which was obtained by subtracting the nodes also belonging to grid levels 5 and 8 from the level 2 grid, yielding a reduction of the sampling effort by: 1/2^(5-2) + 1/2^(8-2) = 9/64 ~ 0.14. This is illustrated in Figure 4.

While reducing the sampling intensity, the intersection of different levels of a grid leads to a spatially balanced pattern of units’ location since the location of the units not sampled is systematic. But, when multiple panels with different reduction rates are merged, the systematicity is however altered, and the regularity is more or less reduced (Figure 4).

To test the impact of sampling effort reduction effect on the spatial balance, a reduction of the sampling intensity was simulated over a virtual portion of grid of 150 x150 nodes (see details in Supplementary Information 2). To this end, the average distance to the four nearest neighboring points was computed for the inner 100x100 cells only, to avoid border effects.

Using higher levels to suppress nodes (Figure 5a) generated only a limited number of possible distances to the nearest neighbor remaining nodes (as also suggested in Fig. 4 low right, and Supp. Fig 1), with a distribution of distances across the simulation area that varies slowly according to the level of reduction. While reducing the number of nodes, some neighboring nodes remain at the original distance of the panel grid, but increasingly less as nodes are being removed, down to the moment when no more node is at that original distance (Supp Fig. 1). Once reached this limit, which corresponds to removing one every two nodes, all nodes are at distance in order of √2 larger than on the original panel grid, hence the upper limit at √2 visible in Figure 5a.

The standard deviation of the mean distance between neighboring points revealed a maximum spatial heterogeneity at sampling reduction intensity of 31%, reaching around +10% (Figure 5b). This level of reduction is critical, as one cannot expect to produce it with alternative options easily (sub-griding in 4 panels would lead to reduction by 20%, sub-sampling on a next homothetic level would yield a 50% reduction). Also, this standard deviation was found to locally saturate along this gradient (e.g., 7-12%, or 20-25% reduction intensity) suggesting that small reductions (in the range of what they are currently implemented to the yearly adjustments of the sampling effort) would not cause substantial distortions to the spatial balance of the samples. For greater reductions rates of the sampling intensity (25-45%), greater values of the standard deviation of the distance among points suggested that the variability in distances among points was critical and does not support the hypothesis of a spatial balance. Between 40% and 50%, the evolution is less regular (smooth). This is due to the limited number of sub-grids with high levels that can be removed. A simulation based on a bigger initial grid would lead to a more detailed and precise evolution without modifying the illustrated trends. In such configurations, resizing the initial base grid may form a logical option.