Vegetation shadow casts impact remotely sensed reflectance from permafrost thaw ponds in the subarctic forest-tundra zone

Thermokarst lakes and ponds are a common landscape feature resulting from permafrost thaw, but their intense greenhouse gas emissions are still poorly constrained as a feedback mechanism for global warming because of their diversity, abundance, and remoteness. Thermokarst waterbodies may be small and optically diverse, posing specific challenges for optical remote sensing regarding detection, classification, and monitoring. This is especially relevant when accounting for external factors that affect water reflectance, such as scattering and vegetation shadow casts. In this study, we evaluated the effects of shadowing across optically diverse waterbodies located in the forest–tundra zone of northern Canada. We used ultra-high spatial resolution multispectral data and digital surface models obtained from unmanned aerial systems for modeling and analyzing shadow effects on water reflectance at Earth Observation satellite overpass time. Our results show that shadowing causes variations in reflectance, reducing the usable area of remotely sensed pixels for waterbody analysis in small lakes and ponds. The effects were greater on brighter and turbid inorganic thermokarst lakes embedded in post-glacial silt–clay marine deposits and littoral sands, where the mean reflectance decrease was from -51 to -70%, depending on the wavelength. These effects were also dependent on lake shape and vegetation height and were amplified in the cold season due to low solar elevations. Remote sensing will increasingly play a key role in assessing thermokarst lake responses and feedbacks to global change, and this study shows the magnitude and sources of optical variations caused by shading that need to be considered in future analyses.


Introduction
Abrupt permafrost thaw is known to accelerate greenhouse gas emissions relative to active layer deepening. However, most of its mechanisms have not yet been incorporated in Earth system and global climate models (Kuhn et al. 2018). This is due to the small relative area of the disturbances, despite their biogeochemical significance (Walter Anthony et al. 2018;Heslop et al. 2020;Turetsky et al. 2020). The most frequent and widespread mechanism of abrupt permafrost thaw is thermokarst lake and pond formation (Bouchard et al. 2017;Walter Anthony et al. 2018;Wauthy et al. 2018). These aquatic environments are less than 10,000 m 2 and are generally less than 5 m deep (Bégin & Vincent 2017). Several studies have shown that the majority of the thermokarst ponds and lakes are biogeochemical hotspots, releasing carbon dioxide (CO 2 ) and especially methane (CH 4 ) (Vonk et al. 2015;Matveev et al. 2016;Kuhn et al. 2018;Zandt et al. 2020).
Optical satellite, airborne and unmanned aircraft systems (UAS) remote sensing are essential tools to map and understand the dynamics of important constituents of waterbodies; e.g., cyanobacterial blooms (CYAN), chlorophyll a (Chl), fluorescent dissolved organic matter (fDOM), colored dissolved organic matter (CDOM), dissolved oxygen (DO), total suspended solids (TSS) (Toming et al. 2016;Peterson et al. 2020;Sagan et al. 2020). However, there are many challenges for the remote sensing mapping and optical monitoring of small waterbodies, from spatial resolution to coregistration errors (Pekel et al. 2016;Muster et al. 2019;Olefeldt et al. 2021). In addition, vegetation surrounding waterbodies may cast shadows, causing difficulties for the analysis of the water spectral characteristics.
Ecologically, shadowing limits the amount of incoming visible and ultraviolet radiation, making it an important factor for aquatic ecosystems. The reduction of incident solar radiation may also affect the water thermal regime, stratification, presence, concentration and behavior of certain chemical species, photosynthetic, photochemical and photobiological transformations, aquatic ecosystem structure and productivity, primary production, microbial communities, and greenhouse gas fluxes (Magnuson et al. 1997;Vincent 2009;Przytulska et al. 2016;Williamson et al. 2020).
In optical remote sensing applications, shadows changing at-sensor solar radiance, are mostly considered as noise that should be removed. Shadows cause biases in spectral indexes (e.g., NDVI, Bowen ratio), indicators (e.g., crop productivity), and spectral properties (Stagakis et al. 2012;Aboutalebi et al. 2019). They can also be problematic for image classification, causing loss of information (Mora et al. 2015;Movia et al. 2016;Milas et al. 2017;Al-Najjar et al. 2019). In optical remote sensing of water, shadows can cause falsepositive detections (Feyisa et al. 2014;Fisher et al. 2016;Pekel et al. 2016;Xie et al. 2016). Distinguishing the small differences between shadows and water, and resolving the effects on spectral reflectance that vary through space and time, remain important challenges (Tian et al. 2017;Guo et al. 2020;Pickens et al. 2020;Yan et al. 2020).
Shadowing problems in remote sensing have been addressed by techniques such as thresholding (Parmes et al. 2017), 3D reconstruction based on 2D objects morphometrics (Zhu et al. 2015), reduction and information recovery (e.g., gamma correction, multisource data fusion) (Movia et al. 2016), and predictive models (Hung et al. 2012). Although the constraints imposed by shadows and potential solutions have received increasing attention in the remote sensing literature, research is still lacking about the effects of shadows on optically contrasting waterbodies (Cordeiro et al. 2021).
Our goal was to develop an experimental study to evaluate the spectral, spatial, and temporal effects of shadowing on optically diverse permafrost thaw lakes surrounded by variable height and dense vegetation. We studied thaw lakes of the boreal forest-tundra transition zone in Nunavik, subarctic Canada. This region provides an ideal model system for assessing and improving remote sensing products for permafrost thaw waters. It is undergoing fast ecosystem changes (e.g., shrubification) and its thermokarst waterbodies are optically diverse, spanning a wide range of turbidities, dissolved organic carbon concentrations and colors, even across small distances (Watanabe et al. 2011). Spectrally calibrated data obtained with a UAS were used to assess the reflectance bias promoted by shadows, aiming at bridging the gap to coarser resolution multispectral satellite data (Lu et al. 2020). Cast shadowing was analyzed by assessing the area of shady and sunlit pixels throughout the year, and the reflectance bias across the range of watercolors.

Study area
The KWAK valley study area is located in Subarctic Canada, Northern Quebec, in the transition zone between sporadic and discontinuous permafrost, showing boreal forest-tundra (Bhiry et al. 2011) (Fig. 1). Since the end of the Little Ice Age, warming and increasing precipitation have led to rapid changes in the landscape. The thawing of ice-rich permafrost mounds has created numerous thermokarst lakes and caused a successional vegetation cover shift characterized by shrubification and terrestrialization processes, as well as the expansion of tree vegetation (Payette et al. 2004;Vallée and Payette 2007;Bouchard et al. 2014).
The UAS survey area of 0.3 km 2 contains 92 thermokarst lakes and ponds varying from 2 to 2000 m 2 . The waterbodies show diverse optical properties, from black to light brown and with contrasting color differences even over small distances, with blue, green, brown, black, and white waters (Laurion et al. 2010;Bouchard et al. 2011;Watanabe et al. 2011;Freitas et al. 2019). Since limnology and water geochemistry are mostly controlled by autogenic processes, varying concentrations of organic (natural dissolved organic matter) and inorganic (post-glacial silt-clay marine sediments) contents are observed, resulting in different colors and greenhouse gas flux characteristics (Laurion et al. 2010;Bégin and Vincent 2017;Folhas et al. 2020).
The Eastern Hudson Bay region has been showing an increasing greening trend over recent decades (McManus et al. 2012;Beck et al. 2015;Ju and Masek 2016). In KWAK, the dense and high vegetation surrounding the shallow (1-3 m) thermokarst lakes, is characterized by trees (Picea mariana, Picea glauca and Larix laricina) and erect shrubs (Betula glandulosa, Salix sp., Alnus crispa and Myrica gale) (Bégin and Vincent 2017). Some lakes show aquatic grass meadows and are surrounded by mosses (Sphagnum spp.) and herbaceous plants such as Carex (Bouchard et al. 2014) (Fig. 2).

General methodological framework
KWAK is a remote area and it is only accessible by helicopter in the summer. It shows large sectors of continuous surficial water logging, numerous lakes and high and dense shrubs. This imposes serious logistical constraints, both for conducting visits to the site, as well as for direct measurements in the multiple ponds and lakes, since walking and equipment transportation from lake to lake is extremely difficult. Hence, the acquisition of ultra-high resolution aerial imagery with UAS is the best way to obtain high multispectral data from a large number of lakes, representative of the spatial diversity of the area.
The data were processed using automatic photogrammetric techniques for point cloud and multispectral orthomosaic generation. The orthomosaic allowed for the classification of shadow and sunlit areas inside lakes, as well as to characterize water reflectance and color. The point cloud was used to derive an ultra-high resolution digital surface model (DSM), which was the basis for shadow modeling. Our analyses focused on how shadowing affects water reflectance and how this changes throughout the year (Fig. 3).

UAS data acquisition and modeling
We used a fixed-wing SenseFly eBee Standard to collect RGB and multispectral data in 2015 and 2017. The surveys differed in flight configuration, sensors and flight times ( Table 1). As a result, the RGB data of 2015 allowed to derive vegetation structure and was used for 3D and cast shadow modeling, while the multispectral data of 2017 were  The flights were done along orthogonal paths. The RGB data acquisition was done in the morning of 28/08/2015. The conditions were overcast and without cast shadows. The multispectral UAS survey was done on 01/09/2017 with clear skies and cast shadows favored by the low solar elevation (20.49° to 22.92°). In that flight, we flew a multispectral Parrot Sequoia camera with four spectral bands: Green (G: 530-570 nm), Red (R: 640-680 nm), Red Edge (REG: 730-740 nm) and Near Infrared (NIR: 770-810 nm) and an optical RGB composite. A sun irradiance sensor was used for single-image radiometric calibration and an albedo target for reflectance calculation (Lu et al. 2020).
We processed the imagery in Pix4DMapper 4.7, which automatically detects characteristic features on the ground, creating thousands of common 3D key points, which are used to generate the dense point cloud, the orthomosaics, and DSMs (McKenna et al. 2017;Pix4D, 2017). The 2015 DSM and orthomosaic show 2.6 cm/ pixel and the ones from 2017, 12.7 cm/pixel. The geolocation was determined by the Global Positioning System (GPS) of the UAS. As a post-processing procedure, we manually co-registered the models, ensuring accurate overlap among all imagery.
To generate an accurate DSM able to precisely model cast shadows onto the thermokarst lakes, we tested different UAS data processing parameters: minimum number of matches (MNM) during the point cloud generation, and the interpolation and smoothing methods during the DSM generation. For the MNM, we tested 2, 3, and 4 reprojection image matches. We then tested two interpolation methods for each MNM: inverse distance weighting (IDW) and triangulated irregular network (TIN). Finally, we tested three smoothing methods for each MNM and interpolation method results: sharp, medium, and smooth.

Validation and selection of the best DSM for shadow modeling
The more accurate DSM was selected by comparing the results of modeled and observed shadows. Ideally, data from the same flight should have been used to validate DSM performance. However, overcast conditions occurred during the acquisition of 2015, which led to the use of the data from 2017 for validation. The date and mean acquisition time were used for defining the sun position and modeling shadows using the DSMs produced for 2015.
The results were classified as true positive (TP-modeled shadows match observed shadows), false positive (FPmodeled shadows in sunlit areas), and false negative (FNmodeled shadows miss observed shadows). This allowed the calculation of quality evaluation metrics, such as precision

Shadow effects on optically diverse thermokarst lakes
To classify the main lake groups by watercolor, we performed a K-Means clustering analysis on the sunlit multispectral data. We focused on the main color classes limiting mismatches between the results and earlier studies (Laurion et al. 2010;Bouchard et al. 2011;Watanabe et al. 2011).
To evaluate the impacts of shadows depending on lake color, we compared the UAS-measured reflectance in the sunlit and shaded areas for each lake. This was obtained from sectors away from lake margins, allowing to mitigate land adjacency effects. These sectors corresponded to 10-20% of the former areas. This approach minimized signal-to-noise factors, such as vegetation scattering effects and impacts of aquatic vegetation (Sawaya et al. 2003;Watanabe et al. 2011). KWAK lakes are optically deep, as concluded by Watanabe et al (2011), which allowed neglecting bottom effects on reflectance (Zeng et al. 2017).

Monthly shadow modeling at the satellite overpass time
The best DSM was used to model the annual dynamics of cast shadows onto the lakes at the satellite overpass time. For this analysis, the solar elevation and azimuth angles for the first day of each month at 16:00 UTC were used. This time was selected as a compromise between the Sentinel-2 (ca. 16:30 UTC), Landsat 8 (ca. 16:00 UTC), and PlanetScope (ca. 15:30 UTC) overpass times. Complementarily, we analyzed the controlling factors for casting shadows by analyzing the surrounding lake environment based on 3D and 2D photointerpretation, oblique photographs from helicopter and modeling the monthly shadow footprints.

Evaluation of the quality of the DSMs for shadow modeling onto the lakes
The DSMs casted larger shadows onto the lakes when the MNM was lowest (Fig. 4). Keeping MNM lower (2) generated more TP, as well as minimizing the FN, although the FP increased slightly. IDW and TIN were inversely related with MNM when analyzing the total area of modeled shadows and TP. When the MNM was highest (4), the IDW performed slightly better than TIN, but when the MNM was lowest (2), TIN performed better than IDW. Accordingly, when MNM was 3, the performance of both methods was even. When analyzing the influence of the smoothing methods, the results for IDW always improved from the sharp to the smooth parameterization. However, for TIN, the quality decreased when the MNM was 4 and only improved for the sharp parametrization when the MNM was 3 and 2.
The best DSM was produced with a MNM value of 2, TIN, and sharp parameters. This allowed detection of the largest modeled shadow area (2,800 m 2 ), the equivalent to 72% of the observed shadows (3,912 m 2 ). This DSM had an overall accuracy of 87% according to the confusion matrix calculation ( Table 2). The precision, recall and F1 score were 0.90, 0.65, and 0.76, respectively.

Lake spectral clustering and shadow effects
The analysis of UAS multispectral data allowed for distinguishing three main lake clusters: black, brown, and light brown (Fig. 5). The black lakes cluster shows the lowest reflectance for all bands. This cluster had a mean value of 1.7% in G (SD = 0.6%), 1.3% in R (SD = 0.7%), 2.9% in REG (SD = 2.1%), and 4.6% in NIR (SD = 2.9%). The light brown lakes had the highest reflectance, with a mean value of 6.6% in G (SD = 1.5%), 7.1% in R (SD = 1.9%), 7.8% in REG (SD = 1.8%), and 9% in NIR (SD = 2.2%). Brown lakes had intermediate reflectance values, namely mean of 3.9% in G  (SD = 1%), 4% in R (SD = 1.2%), 4.7% in REG (SD = 1.3%), and 5.9% in NIR (SD = 1.8%) (Fig. 6). Shadows caused a mean absolute reflectance decrease of 2.1% in G, 2.8% in R, 2.7% in REG, and 2.5% in NIR, the equivalent of -54, -68, -54 and -40% of reflectance comparing to the sunlit sectors, respectively. The analysis of the water reflectance in shaded and sunlit areas according to lake color shows that shadow disturbances depend on the optical properties of the lakes (Fig. 7). For black lakes, the mean absolute reflectance differences between those areas varied from 0.5 to 1.1%, depending on the band. Brown or light brown lakes showed larger differences. The mean absolute reflectance differences were 1.9 to 2.6% in the first case and 3.7 to 5% in the second case.

Assessment of cast shadow impacts at satellite overpass time
The impacts of cast shadows onto the thermokarst lakes are highly dependent on solar elevation, which shows an angle of 55º on 1 July and 10º on 1 January at 16:00 UTC. As a result, the mean shaded lake area varies from 3 in July to 69% in January (Fig. 8). This shows the strong seasonality of shadows, with the smallest effects in June/July and highest effects in December/January. From 1 October to 1 March, the mean shaded lake area ranged from 16 to 69%, and from 1 April to 1 September, from 3 to 10%.
Regarding the conditioning factors, trees and erect shrubs played important and different roles in lake shadowing. Trees close to lakes are scarce in the study area but created well-defined and consistent shadows throughout the year (Fig. 9, A1-A3). In contrast, erect shrubs produced less shadowing, but had a larger overall footprint due to their widespread occurrence (Fig. 9, B1-B3). Erect shrubs were particularly important in the beginning and ending of the warm season and during the cold season. Lake morphometry also played a role, since narrow and elongated lakes were more impacted by shrub cast shadows (Fig. 9, B2). Various lakes were impacted by both, erect shrubs and trees (Fig. 9, C1-C2). Some thermokarst lakes show the importance of catchment topography, a factor especially relevant at the early stages of their formation, since the degradation of icerich permafrost mounds creates concavities in which they are deeply embedded (Fig. 9, D1).

Vegetation structure and shadow modeling
UAS image processing through Structure-from-Motion (SfM) photogrammetric principles provides an efficient and affordable approach for 3D modeling (Dandois et al 2017). However, the resulting point clouds and orthomosaics are highly dependent on the survey circumstances (McKenna et al. 2017;Aboutalebi et al. 2019). Factors such as the UAS type, sensors, conditions (atmosphere, illumination, clouds and shadows), and characteristics of the flights (configuration and overlap) need to be adjusted to ensure the quality and robustness of SfM and derived DSM results (Comba et al. 2015;Dandois et al. 2017;McKenna et al. 2017;Clark et al. 2021).
The overall accuracy of the best cast shadows model for KWAK was 87%. However, cast shadows in most lakes were   underestimated, suggesting that the SfM processing results and DSM quality varied considerably depending on the characteristics of the vegetation surrounding lakes. Modeling shadow casts from conic-shape coniferous trees was better than by erect shrubs, since these revealed irregular shapes and spaces between leaves and stems. The modeling results also suggest limitations regarding the characteristics of the image acquisition, which was always close to zenithal and non-oblique. Although our images showed well-defined shrub cast shadows, a producer's accuracy of 65% was obtained, showing that our survey was not able to fully capture this influence. As a result, shrubs may play an even more important role on casting shadows onto small lakes in the boreal forest-tundra zone than modeled here. Such shrub effect will be more important further north, where solar angles are lower and shrubification processes are known to be rapid (Stow et al. 2004;McManus et al. 2012;Beck et al. 2015;Myers-Smith et al. 2020).

Water reflectance and cast shadow influences
Analysis of the spectral properties of waterbodies is specially complex due to the bidirectional reflectance distribution function of water and the many environmental factors that affect it, such as sun glint, atmospheric scattering, wind, adjacent and bottom effects (Zeng et al. 2017;Sagan et al. 2020), as well as shadows. Finding a pure spectral signature of water using remote sensing, especially on smaller lakes, can therefore be difficult. According to Watanabe et al (2011), the different lake colors in KWAK are caused by varying concentrations of two main optically active substances, namely CDOM, responsible for spectral absorption, especially on the darker lakes (e.g., black cluster) and TSS supporting spectral scattering and absorption on the light brown colored lakes (e.g., brown and light brown cluster). Our results show shadowing changes as a function of surface characteristics in these optically diverse waterbodies. Specifically, we observed that shadow effects as measured by reflectance variations were higher for light brown lakes and almost unnoticeable for black lakes. As a result, correcting cast shadow effects in optical remote sensing imagery, will be most relevant in the light brown colored or turbid lakes, but will be dependent on the percent of lake area affected by shadows. Wauthy et al (2018) suggest an increasing dominance of erodible terrestrial-derived DOM in waters with ongoing permafrost thaw, a process known as browning with the decrease of water column transparency. In remote sensing time-series analysis, this process will act as a reflectance reductor factor, similar to shadow influences. As a consequence, correctly identifying spatiotemporal browning trends of small thaw lakes and ponds in wide regional sectors will imply that shadows are correctly considered and their impacts on reflectance are mitigated. Considering shadow dynamics will allow to avoid mismatches between at-sensor signal and what happens on the field, lowering the signalto-noise ratio and improving remote sensing time-series spectral analysis of small waterbodies in the Circum-Arctic.
The reflectance in shaded areas showed the same behavior of the sunlit areas, with higher values in light brown lakes and lower in clear black lakes, although the magnitude of these variations was comparatively much lower. In both shaded and sunlit areas, the lowest reflectance values Fig. 9 Cast shadow duration at 16:00 UTC and causal factors over KWAK lake surfaces. A. Modeled accumulated monthly footprint; B. Aerial photograph from 28/08/2015. Cast shadow causes: A1-A3-Trees, B1-B3-Erect shrubs, C1-C2-Trees and erect shrubs; D1-Topography always occurred in G and R and the highest in REG and NIR, independently of the watercolor. However, the values in REG and NIR, as true lake water properties, should be interpreted with caution, since they may be influenced by the characteristics of the surrounding environment (Watanabe et al. 2011). In addition, all bands can be affected by a combination of sun/sky glint to some extent, but this is difficult to quantify.
Modeling shadows based on ultra-high resolution DSMs produces accurate results, although such detailed height data are unavailable for the majority of the landscapes. Shadows are also known for having greater impacts when using ultra-high and very high resolution images (Movia et al. 2016;Aboutalebi et al. 2019), since in coarser resolution images they act as a complex varying mixing component of a specific pixel (Milas et al. 2017). In those circumstances, they are not explicit, but implicit, making the correction for this effect even more difficult, although some effects can be minimized using techniques such as band ratios (Zeng et al. 2017).

Conclusions
The application of UAS-derived ultra-high resolution multispectral orthomosaic and DSMs for water reflectance analysis and shadow modeling allowed the quantitative assessment of shadowing effects on thermokarst lakes and ponds. The results show that the impacts of shadow casts by vegetation and topography are greater in light browncolored (turbid) lakes, reflecting the highest absolute differences in reflectance between sunlit and shaded areas. As a result, adequately removing the spectral influence of cast shadows from light brown lakes (e.g., masking) is especially important to access spectral characteristics using remotely sensed imagery, comparatively to darker lakes, where spectral absorption prevails. The reflectance differences induced by shadows in subpolar and polar light brown or turbid lakes can have large consequences for assessing their optical dynamics, especially towards the end of the warm season when solar elevation declines. However, the shadowing effect on thaw lakes varies greatly depending on the local topographic position and its landscape context (proximity to trees, dense high shrubs, and palsa mounds). Certain lakes such as those that are in depressions as a result of thaw subsidence, or that are small, narrow and elongated, as well as closer to trees and surrounded by high and dense shrubs, will be strongly affected by cast shadows. These effects may occur even at the highest solar elevation angle and limit the usable area of unaffected water pixels for monitoring analysis. This problem will become increasingly relevant with the ongoing shrubification and accelerated tree growth that is taking place in Subarctic and Arctic regions.

Conflict of interest
The authors have no conflict of interest to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.