Electromagnetic Induction Is a Fast and Non-Destructive Approach to Estimate the Inﬂuence of Subsurface Heterogeneity on Forest Canopy Structure

: The spatial forest structure that drives the functioning of these ecosystems and their response to global change is closely linked to edaphic conditions. However, the latter properties are particularly difﬁcult to characterize in forest areas developed on karst, where soil is highly rocky and heterogeneous. In this work, we investigated whether geophysics, and more speciﬁcally electromagnetic induction (EMI), can provide a better understanding of forest structure. We use EMI (EM31, Geonics Limited, Ontario, Canada) to study the spatial variability of ground properties in two different Mediterranean forests. A naturally post-ﬁre regenerated forest composed of Aleppo pines and Holm oaks and a monospeciﬁc plantation of Altlas cedar. To better interpret EMI results, we used electrical resistivity tomography (ERT), soil depth surveys, and ﬁeld observations. Vegetation was also characterized using hemispherical photographs that allowed to calculate plant area index (PAI). Our results show that the variability of ground properties contribute to explaining the variability in the vegetation cover development (plant area index). Vegetation density is higher in areas where the soil is deeper. We showed a signiﬁcant correlation between edaphic conditions and tree development in the naturally regenerated forest, but this relationship is clearly weaker in the cedar plantation. We hypothesized that regular planting after subsoiling, as well as sylvicultural practices (thinning and pruning) inﬂuenced the expected relationship between vegetation structure and soil conditions measured by EMI. This work opens up new research avenues to better understand the interplay between soil and subsoil variability and forest response to climate change.


Introduction
Forests play major roles in biogeochemical cycles (water, carbon, nitrogen) with strong repercussions for our societies including water availability, carbon sequestration, energy, and biodiversity. However, rising temperatures and change in rainfall distribution due to climate change lead to increasing drought intensity and frequency, higher potential evapotranspiration, and therefore impose greater and more frequent water stress for trees. This can affect forest structure and functioning [1] and in turn major forests' services [2]. A better understanding of the impacts of the edaphic environment on forest structure and functioning is a major scientific challenge to better anticipate the effects of climate change and impacts on our societies. However, how climatic conditions determine water stress depends strongly on edaphic conditions [3].
Vegetation transpiration (T) is a major component in the water cycle accounting for 40% of continental rainfall at the global scale [4][5][6][7]. According to Li et al. [8], T accounts for nearly 50% of precipitation in forest areas. Anticipating the evolution of T is complex because T is linked to several parameters including climatic demand, vegetation development, and water availability and accessibility from the ground. The latter issue of water resources availability is still particularly poorly understood in forest soils despite its critical importance for forest development and survival [3,9]. Both soil and vegetation are spatially structured, but the difficulties to acquire spatialized data of the soil and subsoil properties have limited the development of knowledge regarding the interrelation between forest structure and functioning with local variability of soil and subsoil properties. In this article, we define soil in the pedological sense (soil horizons above the weathering rock) and the subsoil represents what lies beneath the soil, i.e., the intact and weathered rock (fractures, faults, conduits). We use in this article the notion of soil/subsoil because geophysical measurements, as implemented in this work, do not allow to clearly distinguish the two entities.
On the one hand, seedlings located on shallower soils are less likely to grow, survive, and develop, and tree density should thus be locally lower as the forest stand matures. On the other hand, on such poor soils, a lower tree density implicate a greater light availability for regeneration, that it is less subjected to adult competition. In other words, seedlings are more likely to survive because there is less competition for light, but also more likely to die due to greater water stress.
For mature trees, the effects of locally shallower soil are also contradictory. Increasing water stress negatively impacts growth, but decreases also the competition. In addition, it has been shown that trees adapt their structures (i.e., leaf area and rooting depth) and functions to the constraints of different local environments [10,11]. Indeed, the eco-hydrological equilibrium theory posits that vegetation develops according to constraints imposed by available water reserves [12,13]. However, when stress regimes change, these equilibria turn into imbalances, and trees that have benefited from significant water resources may become the most sensitive to new and more intense water stresses [14].
In order to better understand the relationships between the spatial structure of the soil and the spatial structure of vegetation, we need to improve our means of investigating spatial variations in the soil and subsoil. Several approaches are used to characterize the water resource availability to the vegetation, including pedological pits [15] and punctual soil water estimates using capacitive or neutron probes [16]. However, the "deep" water resource (beyond 1 or 2 m) is still questioned within the forest ecologists community because it remains difficult to characterize deeper horizons and the solution of "deeper" drilling (until 4 or 6 m) is costly, destructive, and has a limited spatial representativeness [17,18].
Indeed, most approaches classically used in forest ecology (e.g., soil pits, humidity sensors) have a limited investigation depth and a limited spatial representativeness. Surface based geophysics is an alternative that offers a wide range of tools based on different physical principles (electrical, electromagnetic, seismic, gravimetric, and magnetic), [19]. In agronomy and soil sciences, geophysics has entered the list of commonly used tools [20][21][22]. Garré et al. [23] define the term agrogeophysics. Such non-destructive approaches make it possible to obtain spatialized information on the physical properties of the underground and to carry out temporal monitoring. Geophysical approaches offer larger spatial representativeness than classical approaches (e.g., soil pits, time domain reflectometry (TDR) probes). Despite few forest ecology studies using geophysics [10,[24][25][26][27], this approach is still under-exploited in relation to its potential. In addition, all these studies use the electrical resistivity tomography (ERT) technique, which allows to produce geoelectrical cross-sections. However, other approaches such as electromagnetic induction (EMI), which allows to produce geoelectrical maps, may have relevant and complementary strengths with ERT in describing the spatial variability of forest soils/subsoils properties. Indeed, the EMI allows to quickly cover large areas to produce geoelectric maps over several hectares in one day. While at the same time, ERT hardly exceeds four cross-sections, but provides a more detailed view of the geoelectric properties' distribution in depth. The combination of ERT and EMI is widely used in other fields [28][29][30][31], because EMI allows to investigate large areas and to spatially expand some of the interpretations from ERT.
In this paper, we investigate whether an EMI approach can identify the spatial variability of soil/subsoil properties that could provide insights into the spatial variability of forest stands' development (plant area index, PAI). We investigated two complementary forest sites, one is a naturally regenerated forest composed of Aleppo pines and Holm oaks and the other is a monospecific plantation of Altlas cedar. We also performed several ERT cross-sections on each site to better characterize soil thickness and helped interpreting the EMI signal.

Experimental Sites
We carried out geophysical experiments at two nearby sites in south-eastern France ( Figure 1). Both sites have a Mediterranean climate with annual rainfall of about 700 mm in Font-Blanche and about 770 mm in Valliguières and mean temperatures of 14.5 and 14 • C, respectively. The two sites are located on similar karstified cretaceous limestone with urgonian facies including rudist. We assume that the water resource available to plants is contained in the soil and the weathered rock (fracture, fault, conduits), also called epikarst. They differ by the level of management: one is a naturally regenerated forest (Font-Blanche) and the other one is a plantation (Valliguières). The general objective of this experimental design was to evaluate the effects of silvicultural treatments (thinning and pruning) on tree growth [32,33]. Due to the stony and karstic nature of the ground, a subsoiling was carried out at 50 cm depth before planting. Seven soil pits distributed over the experimental site to highlight the spatial heterogeneity of the site revealed a stony rendzina (calcisol) whose thickness ranged between 15 and 70 cm according to pedologic pit exploration. The rock fraction is rather homogenous and high: about 75-80%. A low density understory has naturally grown after the planting on some plots and includes various species: Quercus Ilex L, Phillyrea L, and Juniperus oxycedrus.
The experimental treatments were set up in 16 plots. The treatments were defined as follows: Four thinning treatments (density after thinning): A = 1200 trees/ha (i.e., no thinning), B = 800 trees/ha, C = 600 trees/ha, D = 400 trees/ha. Thinning was applied in 1992 (one replication). Between 400 and 450 trees per hectares were then selected in each plot and are termed crop trees. The soil is a stony rendzina, whose thickness is highly variable over the study site. Four soil pits distributed over the experimental site to highlight the spatial heterogeneity of the site revealed soil depths ranging from 0 to 60 cm with a rock fraction ranging between 50% in the near surface and 90% below 50 cm.
This experimental site is dedicated to the long-term monitoring of forest biogeochemical cycle as part of the ICOS (Integrated Carbon Observation System) network and to the study of the effects of rainfall pattern changes on forest structure and functions. It includes a main control plot (80 m × 80 m) and different sub-plots (25 m × 25 m) with more detailed measurements. The M30 sub-plot (25 m × 25 m) was equipped in 2008 with a rainfall exclusion system using PVC gutters installed 2 m above the ground and that excludes 30% of rainfall. The P30 plot received approximately 30% additional rainfall during the years 2010 and 2011. The PM30 plot is a control plot equipped with gutters positioned upside-down to account for the potential effects of gutter shade. The TD plot is a control plot without gutters. The site is equipped with an eddy-covariance flux tower, rising 17 m above the ground, where meteorological variables, carbon, water, and energy fluxes are continuously monitored. More details are available in Moreno et al. [34].

Electromagnetic Induction
Electromagnetic Induction (EMI) and in particular Frequency Domain EM (FDEM) is a technique based on the propagation of electromagnetic waves that allows to quickly map the apparent electrical conductivity (ECa) of the underground (up to several hectares per day). We used the EM31-MK2 apparatus developed by Geonics Limited. The transmitter and receiver coils are spaced 3.6 m apart and the operating frequency of the instrument is 9.8 kHz. The field conditions (low electrical conductivity) make the device used in low induction number (LIN) conditions [35,36]. This configuration allows to reach an investigation depth around 3 m in horizontal dipole (HD) and 4 m in vertical dipole (VD). These investigation depths are consistent with the expected development of root zones for the three species. The investigation depth varies slightly according to the underground Eca [35,37]. The surveys were done in spring (March 2015 in Font-Blanche and May 2017 in Valliguières) when the soil and weathered rock is filled with water in order to maximize the geophysical contrast with massive limestone. The measurements positioning was quite challenging in the two forest contexts where the trees are evergreen. On the site under artificial regeneration (Valliguières), we georeferenced the extremities of each inter-row and we acquired continuously in HD mode only. The points were realigned along a straight line connecting these two georeferenced extremities. On the site under natural regeneration (Font-Blanche), we made measurements every 5 m above georeferenced stakes, which are used for other monitoring. We thus carried out the measurements in HD mode only. We acquired 22,213 points in Valliguières. At Font-Blanche, we carried out the measurements in HD and VD mode. EMI measurements did not include M30, PM30, P30 plots because of iron structures (support for gutters) that can interfere with EMI. Similarly, the area close to the flux tower has been excluded as well. We removed 57 points out of 256 due to these metal objects, which disturb EMI signal. The ECa maps shown in Figures 2 and 3 were calculated with Surfer software (Golden Software ® , Golden, CO, USA) using the kriging method. We did not apply any temperature correction on ECa because the temperature was quite stable during the acquisition. Moreover, we checked that there was no device drift by returning to the starting point at the end of each survey.

Electrical Resistivity Tomography
Electrical Resistivity Tomography (ERT) is a commonly used geophysical method in environmental studies consisting of estimating the spatial variability (lateral and vertical) of electrical resistivity along a transect. In Font-Blanche, ERT measurements were performed along two 126 m long transects (Figure 2). In Valliguières, four ERT transects were installed all along plots 11 and 12 (600 stems/ha), 21 and 22 (1200 stems/ha), 41 and 42 (800 stems/ha), and 91 and 92 (400 stems/ha), for a length of 110 m. The choice of transects was positioned in accordance with both the EMI map and the forest structure. We used an ABEM Terrameter SAS 4000 [38] with 64 regularly spaced electrodes every 2 m. The surveys were done in spring (March 2017 in Font-Blanche and April 2018 in Valliguières) when the soil and weathered rock is filled with water in order to be consistent with EMI surveys. The ERT and EMI acquisitions were not performed in the same year for logistical reasons, but we assume that this time lag is not limiting for a relative comparison of geophysical properties between the two techniques. Indeed, in these forest contexts where the topography is very flat, most of the variations in geophysical properties are related to seasonal variations (i.e., temperature and water content). The positioning of the end of the profiles was done with a differential GPS (Trimble Geo7X, Sunnyvale, CA, USA). Resistivity data were acquired using three different complementary arrays including Wenner-Schlumberger, Dipole-Dipole, and Gradient [39]. The three arrays were merged before inversion. We only kept data with a repetition error of less than 1% for processing. Up to ten additional points could be removed per section using Res2Dinv's "exterminate bad data points" function which identifies incoherent resistivity values [40]. Apparent resistivities were processed using the Res2Dinv 3.57.36. The models presented below were all obtained after 3 iterations. This limited number of iterations is a good compromise that allows to reach acceptable RMS without generating extreme resistivity values to reach a mathematical optimum that would risk to move away from the field reality. The mesh was refined and the resulting ERT models have a 1 m lateral resolution and a vertical resolution that ranges gradually from 0.4 m near the surface to 2.5 m at a depth of 20 m.

Ecophysiological Measurements
The commonly used Leaf Area Index (LAI) is defined as the projected area of green foliage per unit horizontal ground surface area. It corresponds to the main surface exchange of carbon and water fluxes between the forest canopy and the atmosphere. LAI is most commonly estimated with indirect optical methods based on the light interception properties of the foliage. With these methods it is often not possible to distinguish the foliage from other plant parts (i.e., stem, branches, cones, trunk) and the term Plant Area Index (PAI) as the sum of all plant surfaces is often more appropriate.
We estimated PAI using hemispherical photographs. Photos were taken when the sun was low (mostly at sunset) or, very rarely, when the sky was uniformly overcast. Photos were made with a Nikon D3200 (Tokyo, Japan) digital camera with a Sigma 4.5 mm EX DC HSM fisheye lens. The camera was positioned skywards and the top of the camera was directed to the north. Camera wettings were ISO 400, automatic white balance, aperture F5.6 to F8, and shutter speed adjusted manually so as to expose for the sky. The RAW files were developed with the DCRAW software, with no gamma correction. The resulting Tiff files were then processed and analyzed using a set of macros written for the ImageJ software [41]: extraction of the blue channel, automatic contrast adjustment, manual double thresholding, and calculation of the gap fraction. These macros replicate the thresholding method proposed by Leblanc et al. [42]. Gap fraction data were then converted to PAI using Miller's formula, as approximated by Welles and Norman [43].
At Font-Blanche, photos positions followed a regular 5 m × 5 m (TD plot) and 10 m × 10 m grid (main plot) at Font-Blanche ( Figure 2C), and were taken at 40 cm above the ground on a self-levelling platform (manufactured at INRAE-URFM, Avignon, France). PAI is routinely monitored at Font-Blanche every two months. For this study, we used data from the 8 April 2015, the closest date to the EMI survey.
At Valliguières, photos were taken at 1.3 m above ground using a tripod in plots 11,12,21,22,41,42,91,92, and 102. Positions followed two transects per plot, with photos taken every 5 m for the transect closest to the ERT transect, and every 10 m otherwise ( Figure 3A). Photos were taken on the 11 October 2017.

Font-Blanche: A Naturally Regenerated Forest
The EMI mapping ( Figure 2B) shows a strong spatial variability of soil/subsoil conditions at the scale of the experimental plot. It can be observed that, in general, areas further south are less conductive (red) than the north. Unfortunately, there are gaps in data acquisition in several areas that are related to the presence of metals (i.e., structures carrying gutters (M30 and PM30), flux tower, and other equipment). These gaps limit the spatial extension of our interpretation.
However, we note a satisfactory consistency between the EMI and ERT results (Figure 2A,D). The coherency between both approaches is well-known and documented by numerous authors in various geological contexts [28][29][30]. The ERT cross-section ( Figure 2A) shows considerable variations in resistivity that are consistent with field observations and EMI in Font-Blanche. Indeed, at the distance 0 m (Figure 2A) a very resistant zone was observed on the surface agreed with a limestone outcrops. Also, such findings are consistent with recent papers in karst environments with similar geological and pedological contexts [10,26] that have shown that the most resistive areas of ERT profiles concur with most stony areas made of massive limestone, whereas the most conductive areas in the near surface concur with soil and weathered limestone. We will extend this interpretation initially developed for ERT to the EMI results in the same idea as the work of Grellier et al. [44].
In addition, we confirmed the empirical link between EMI signal and soil thickness estimates with a miner's bar across the experimental plot. We found a significant correlation (p-value = 3.10 −4 ) between ECa (EMI in horizontal dipole) and estimated soil depth ( Figure S1). Altogether, this information allow us to interpret the EMI signal, assuming that the most resistant areas (red) of the EMI map ( Figure 2B) are those where the soil is less developed.
PAI estimates ( Figure 2C) revealed an heterogeneous distribution of the vegetation cover. The southern part of the plot presents lower PAI values than the northern part. Agreement between the ECa map ( Figure 2B) and the PAI map ( Figure 2C) appears visually, with a tendency of greater PAI associated to higher conductivity zones.

The Valliguières: Atlas Cedar Experimental Plantation
We investigated using EMI in all 16 plots, totaling about 8 hectares spread over four separated zones in the Valliguières forest. We note that the ECa is quite contrasted between the plots (Figure 3). Some plots are fairly homogeneous ( Figure 3A), such as plots 51 and 52 (northernmost plot, see position on Figure 1B), while others are more heterogeneous, such as plots 21 and 22 (southernmost plot, see position on Figure 1B). EMI allows to identify the general geological structure directed east/west, which agreed with the syncline general orientation [45]. The four ERT cross-sections were carried out on four contrasting areas of the EMI map. The resistivities of the first three meters of the ERT cross-sections are consistent with the EMI map. Indeed, plots 21 and 22 are located in the most conductive zones and the ERT cross-sections show significantly more conductive materials in the near surface. Contrary to Font-Blanche, the deep parts (more than 5 m) of the ERT cross-sections show quite contrasted resistivities in the cross-sections (particularly between plots 22 and 92). We suppose that this variability is related to the variability of limestone layers inherent to the Urgonian facies [46]. Carrière et al. [47] showed on a nearby site that this variability in petrophysical properties between limestone layers induces a variability in geophysical signature detectable by ERT.

Combined Interpretation
Geophysical data helped us to better understand the variability of soil/subsoil conditions on the two sites. Firstly, we found that ERT and EMI were consistent together. For example, the ERT-11/12 cross-section ( Figure 3B), performed in the most conductive area of the EMI map ( Figure 3A), reveals a significantly higher conductivity than the ERT-41/42 cross-section, which is performed in a more resistant area. At the Font-Blanche site, the EMI signal was closely related to thickness estimates ( Figure S1). This concurs with the results of Nourtier et al. [26] and Carrière et al. [10], who found a close relationship between ERT signal and soil water reserves estimated through pedological pits.
Secondly, we seek to relate the distribution of soil/subsoil properties to the vegetation cover structure. The results shown in Figures 2 and 4 highlight that, for the Font-Blanche site, the least conductive zones that we interpret as the rockiest match the zones where vegetation is the least developed (low PAI). On the contrary, areas with high ECa exhibit higher PAI. Indeed, when the data are binned per PAI classes, a significant relationship appears (p-value = 5.10 −6 ) in the forest of Font-Blanche ( Figure 4A), which is a natural regeneration, not disturbed by human activities (no thinning or pruning). By contrast, this relationship was not significant (p-value 0.34) for the Valliguières forest, when all the data were considered together ( Figure 4B). We explain this weak correlation by anthropogenic interventions (subsoilling, regular plantation, pruning, and thinning) and geological variability that makes the relationship more complex. In Valliguières, anthropization and geological effects appear as confounding factors to our interpretation at the whole site scale. This is why we carried out the same analysis on subgroups of plots that underwent similar managements ( Figure 5), to assess whether PAI variations can still be related to soil/subsoil heterogeneity when the anthropic component is the same. Moreover, at the plot scale we assume less variability in bedrock properties because changes in limestone facies [46] are progressive along NNW/SSE axis (i.e., perpendicular to the pluri-kilometric scale anticline [45]). In these reduced areas, we observe trends that are consistent and complementary with our observations made at Font-Blanche ( Figure 5), with higher PAI on the most conductive areas. We assume that vegetation is less developed on the less conductive zones because they are rockier and therefore less favorable to vegetation development. However, the limited number of observation points (n = 20 and n = 9) and the smaller range of ECa ( Figure 5) lead us to interpret plot scale result with caution.

Discussion
The spatial variability of the EMI signal is influenced by two main factors on the two sites: (i) the thickness of soil and weathered rock layer; (ii) the bedrock petrophysical properties. While the former is expected to relate to the tree functioning and thus the structure of the forest, the latter is supposed to have less influence on ecological properties. Indeed, the experimental plot of Font-Blanche is quite small (<1 ha) and the bedrock properties are quite homogeneous at the plot scale, despite the intrinsic karst heterogeneity, allowing to interpret EMI spatial variations mainly as variations in soil/weathered rock depth variations. Here, a relationship between variations in PAI and EMI conductivity emerges ( Figure 4A). However, the experimental site of Valliguières was investigated over a larger area (~20 ha). Here, the variability of ERT cross-sections at depth reveals a greater variability in the bedrock petrophysical properties ( Figure 3). Therefore, the variability of the EMI signal cannot be compared directly between plots where the ERT shows contrasting electrical resistivities at depth (i.e., plots 21 and 91, Figure 1B). Indeed, including plots with contrasting bedrock conditions, will noise the relation we seek to highlight (soil and weathered rock thickness), which might no longer be the main variable that influences the EMI signal.
In addition, forest management activities affect soil/subsoil properties (subsoiling) and pruning and thinning affect tree growth and PAI. Such management can have significant effects on tree development [32] by exerting a stronger control than natural soil/subsoil properties. In the managed forest of Valliguières, even if conductivities exhibit a large spatial variation, management and treatments, combined with bedrock variability, hide the possible link between PAI and spatial variation of conductivity ( Figure 4B). However, considering only more homogeneous plots with the same treatment, this link seems to hold again ( Figure 5). Clearly, the bedrock variabilities and forest management activities appear as confounding factors and should be taken into account in interpretation.

Conclusions and Perspectives
Throughout this paper, we have shown that geophysical approaches such as ERT and EMI can provide relevant information to better understand forest structure at the experimental site scale. We interpreted the variability of EMI signal as a variability of the soil/weathered rock thickness. This approach makes it possible to detect areas more or less favorable for vegetation development. Indeed, in areas of high conductivity, we observed higher PAI indicating a greater biomass. Our results are consistent with the eco-hydrological equilibrium hypotheses [12,13] and with different empirical observations [11,48]. This theory describes how water-limited natural vegetation systems are in balance with their climatic and soil environments. The theory assumes that biomass productivity is proportional to available water resources. More recently, Nourtier et al. [26] and Carrière et al. [10] have taken these assumptions further and shown that vegetation in areas with higher biomass productivity (a priori favorable soil condition) is more vulnerable to severe water stress. We verify this assumption with another geophysical technique (EMI) in this paper. Moreover, we also show that the relationship between geophysical properties and biomass productivity is also sensitive to the petrophysical bedrock properties and the sylvicultural practices. Consequently, linking geophysical and ecophysiological data must be done carefully. Indeed, the link between these different types of data cannot be done directly. A detailed interpretation phase must be respected for the geophysical results in order to identify the various factors that may influence the signal and potentially disturb the analysis. Moreover, it would be relevant to characterize a forest site at various depths by combining several EMI devices. The exploration of the seasonal variation of EMI signal could also be relevant to characterize the spatial variability of forest soils and subsoils. Considering these factors, using geophysics in forest ecology opens new research avenues to better understand the structure and functioning of forest stands and their evolution in the global change context.