Estimation of Surface Canopy Water in Pacific Northwest Forests by Fusing Radar, Lidar, and Meteorological Data

Surface Canopy Water (SCW) is the intercepted rain water that resides within the tree canopy and plays a significant role in the hydrological cycle. Challenges arise in measuring SCW in remote areas using traditional ground-based techniques. Remote sensing in the radio spectrum has the potential to overcome the challenges where traditional modelling approaches face difficulties. In this study, we aim at estimating the SCW by fusing information extracted from the radar imagery acquired with the Sentinel-1 constellation, aerial laser scanning, and meteorological data. To describe the change of radar backscatter with moisture, we focused on six forest stands in the H.J. Andrews experimental forest in central Oregon, as well as four clear cut areas and one golf course, over the summers of 2015–2017. We found significant relationships when we executed the analysis on radar images in which individual tree crowns were delineated from lidar, as opposed to SCW estimated from individual pixel backscatter. Significant differences occur in the mean backscatter between radar images taken during rain vs. dry periods (no rain for >1 h), but these effects only last for roughly 30 min after the end of a rain event. We developed a predictive model for SCW using the radar images acquired at dawn, and proved the capability of space-based radar systems to provide information for estimation of the canopy moisture under conditions of fresh rainfall during the dry season.


Introduction
Surface canopy water (SCW) is the water resting on the outer surfaces within the canopy layer, and plays a significant role not only on the hydrological cycle but also on the characteristics of the boundary layer of the troposphere [1]. Although SCW can refer to canopies of crops, grasses, and forests, in this paper, we refer specifically to the structurally complex forest canopies of the Pacific Northwest region of North America. To accurately describe the atmospheric processes occurring above the forest, estimates of the amount of SCW are needed. The extensive networks of weather stations make the prediction of SCW relatively easy and accurate in flat areas [2,3], but significant challenges are faced in mountainous or remote regions, as more complex extrapolations are required or crucial parameters may be unknown [4]. Nevertheless, remotely sensed data can enhance prediction of SCW in areas where the traditional modelling approaches are challenged. The critical aspect of using remote sensing data in any investigation is the selection of the appropriate wavelengths. Although the visible and near-infrared spectrum are commonly used to measure the secondary effects of SCW on the color of the canopy, radio waves can directly measure the dielectric constant from water. In many situations, radar (RAdio Direction And Ranging) is preferred over passive optical sensors, as they exhibit a series of sought-after properties, such as: (1) ability to penetrate cloud cover and vegetation (although the penetration power depends highly on the wavelength), (2) independence of local solar intensity (which allows operation irrespective of day or night), and (3) The study was carried out on six adjacent stands totaling approximately 41 ha, which are all located within the HJ Andrews forest in the Western Cascade Range of Oregon (Table 1). Additionally, five areas outside of the HJ Andrews were considered for verification. The five areas included a golf course to the south and four recently clear-cut areas to the west of the HJ Andrews forest. The HJ Andrews forest covers roughly 6400 ha of mountainous terrain and contains conifer forests typical of the Pacific Northwest of the North American continent, namely Douglas Fir and Western Hemlock [30]. The six stands are centered around the UTM zone 10 coordinates (565,600 E, 4,900,900 N) (Figure 1), and represent three stages in the life of a stand: young, mature, and old-growth ( Figure 2). Each age class is represented by two stands. The golf course will serve as control, since golf courses are generally well irrigated and have little to no geometric effects that could cause differential scattering. The clear cuts serve a similar purpose but include the scattering from geometric effects due to the abundance of slash piles and debris left from logging activities.

Study Site
The study was carried out on six adjacent stands totaling approximately 41 ha, which are all located within the HJ Andrews forest in the Western Cascade Range of Oregon (Table 1). Additionally, five areas outside of the HJ Andrews were considered for verification. The five areas included a golf course to the south and four recently clear-cut areas to the west of the HJ Andrews forest. The HJ Andrews forest covers roughly 6400 ha of mountainous terrain and contains conifer forests typical of the Pacific Northwest of the North American continent, namely Douglas Fir and Western Hemlock [30]. The six stands are centered around the UTM zone 10 coordinates (565,600 E, 4,900,900 N) ( Figure  1), and represent three stages in the life of a stand: young, mature, and old-growth ( Figure  2). Each age class is represented by two stands. The golf course will serve as control, since golf courses are generally well irrigated and have little to no geometric effects that could cause differential scattering. The clear cuts serve a similar purpose but include the scattering from geometric effects due to the abundance of slash piles and debris left from logging activities. Table 1. Summary statistics of the study areas.

Stand Area [ha] Age [Years] Number of Trees Total Height [m]
Elevation (Min-Max) [m above Sea Level] Figure 1. The six areas inside HJ Andrews for which the water canopy analysis is executed as well as the location of the weather station H15MET.

Weather Data
The weather data used in our study were provided by the HJ Andrews Experimental Forest [31]. The closest available weather station was H15MET, which is roughly 1 km from the study areas ( Figure 1). The measurements were taken every 5 min for precipitation and every 15 min for temperature and other atmospheric variables. We computed four meteorological attributes based on the weather data: antecedent moisture condition (AMC), solar radiation, time since rain (TSR), and evapotranspiration, which was estimated with three models (i.e., Priestley-Taylor, Penman, and Brutsaert). The weather variables were summed to 15 min, 30 min, 45 min, 60 min, 90 min, and 2, 3, 4, 5, and 6 h, to match the radar imagery, according to Equation (1).
where X is one of the three statistics, xi is the value to be summed (which for AMC is the precipitation), and k is the time, which varies from 15 min to 6 h. The time since rain (TSR) estimates how many minutes have passed since rainfall has stopped (we considered that rain had stopped if <0.02 mm had fallen in the preceding 30 min).

Radar Data and Image Processing
To model SCW within the forest canopy, we have used radar data provided by the Sentinel-1 constellation, which contains two satellites in sun-synchronous orbits. Among the SAR satellites, the Sentinel-1 constellation, developed by the European Space Agency, has supplied radar images with a high spatial resolution (i.e., approximately 10 m) since

Weather Data
The weather data used in our study were provided by the HJ Andrews Experimental Forest [31]. The closest available weather station was H15MET, which is roughly 1 km from the study areas ( Figure 1). The measurements were taken every 5 min for precipitation and every 15 min for temperature and other atmospheric variables. We computed four meteorological attributes based on the weather data: antecedent moisture condition (AMC), solar radiation, time since rain (TSR), and evapotranspiration, which was estimated with three models (i.e., Priestley-Taylor, Penman, and Brutsaert). The weather variables were summed to 15 min, 30 min, 45 min, 60 min, 90 min, and 2, 3, 4, 5, and 6 h, to match the radar imagery, according to Equation (1).
where X is one of the three statistics, x i is the value to be summed (which for AMC is the precipitation), and k is the time, which varies from 15 min to 6 h. The time since rain (TSR) estimates how many minutes have passed since rainfall has stopped (we considered that rain had stopped if <0.02 mm had fallen in the preceding 30 min).

Radar Data and Image Processing
To model SCW within the forest canopy, we have used radar data provided by the Sentinel-1 constellation, which contains two satellites in sun-synchronous orbits. Among the SAR satellites, the Sentinel-1 constellation, developed by the European Space Agency, has supplied radar images with a high spatial resolution (i.e., approximately 10 m) since 2016. The satellites have a revisit time of 6 days and 12 days, respectively, where the ascending (where the orbit takes the spacecraft from South to North) passes occur at local dusk, while descending (spacecraft going from North to South) passes occur at local dawn. Each satellite is equipped with a 12-m radar array that operates in the C-band (5.405 GHz) with two polarization modes (HH + HV, VV + VH). A thorough overview of the technical aspects of the satellites as well as the data collection process is provided by Torres et al. [32].
The ease and open access to Sentinel-1 data, combined with the characteristics of the C-band, makes this mission suitable for temporal investigation of canopy dynamics.
The Sentinel-1 data were obtained from the European Space Agency (ESA) portal (2014). Preprocessed imagery was downloaded in a ground-range convention for spring to fall dates (17 April through 31 October) during 2015, 2016, and 2017. We excluded all images that contain ice or snow, as frozen water behaves differently than liquid water due to differences in dielectric properties [9,26]. To exclude any residual snow or ice, we used 5 • Celsius as a minimum temperature, as we considered that, at 5 • C, all the snow and ice would have melted. In total, we used 67 dates that met our requirements. We processed the ground-range convention images with the Sentinel-1 dedicated image processing program (known as SNAP), which outputs images ready to interpret by addressing all issues, corrections, and calibrations associated with radar imagery, such as thermal noise removal, terrain correction, or de-speckling. We pre-processed the images following the recommendations of Small [9], in which the radar signal was converted to gamma-naught to account for the terrain's geometric and radiometric effects ( Figure 3). 2016. The satellites have a revisit time of 6 days and 12 days, respectively, where the ascending (where the orbit takes the spacecraft from South to North) passes occur at local dusk, while descending (spacecraft going from North to South) passes occur at local dawn. Each satellite is equipped with a 12-m radar array that operates in the C-band (5.405 GHz) with two polarization modes (HH + HV, VV + VH). A thorough overview of the technical aspects of the satellites as well as the data collection process is provided by Torres et al. [32]. The ease and open access to Sentinel-1 data, combined with the characteristics of the C-band, makes this mission suitable for temporal investigation of canopy dynamics. The Sentinel-1 data were obtained from the European Space Agency (ESA) portal (2014). Preprocessed imagery was downloaded in a ground-range convention for spring to fall dates (April 17 through October 31) during 2015, 2016, and 2017. We excluded all images that contain ice or snow, as frozen water behaves differently than liquid water due to differences in dielectric properties [9,26]. To exclude any residual snow or ice, we used 5° Celsius as a minimum temperature, as we considered that, at 5° C, all the snow and ice would have melted. In total, we used 67 dates that met our requirements. We processed the ground-range convention images with the Sentinel-1 dedicated image processing program (known as SNAP), which outputs images ready to interpret by addressing all issues, corrections, and calibrations associated with radar imagery, such as thermal noise removal, terrain correction, or de-speckling. We pre-processed the images following the recommendations of Small [9], in which the radar signal was converted to gamma-naught to account for the terrain's geometric and radiometric effects ( Figure 3). The raw images downloaded from the Copernicus web portal are typically in a slant range geometry (Figure 4a,) which corresponds to the Lagrangian coordinates in relation to the spacecraft and not in relation to a fixed datum. Pre-processing positions the images in Eulerian coordinates (i.e., a preset datum), but the raw data will be either 'flipped' along the east/west axis for ascending images or 'flopped' along the north/south axis for descending images (Figure 4b). Essentially, the raw image records backscatter as range increases from left to right, and along the azimuth track from top to bottom. Therefore, descending images appear to be inverted horizontally, as the spacecraft is travelling north to south, but looking east to west and, consequently, recording westerly locations as further away (i.e., on the right-side of the image). Similarly, on ascending images, the spacecraft is moving south to north. Therefore, the top of the image will show the southern information and the bottom will display the northern information. Since the spacecraft is looking west to east, the near-range will appear on the left-side of the image and the farrange areas will be on the right-side of the image. The raw images downloaded from the Copernicus web portal are typically in a slant range geometry (Figure 4a,) which corresponds to the Lagrangian coordinates in relation to the spacecraft and not in relation to a fixed datum. Pre-processing positions the images in Eulerian coordinates (i.e., a preset datum), but the raw data will be either 'flipped' along the east/west axis for ascending images or 'flopped' along the north/south axis for descending images (Figure 4b). Essentially, the raw image records backscatter as range increases from left to right, and along the azimuth track from top to bottom. Therefore, descending images appear to be inverted horizontally, as the spacecraft is travelling north to south, but looking east to west and, consequently, recording westerly locations as further away (i.e., on the right-side of the image). Similarly, on ascending images, the spacecraft is moving south to north. Therefore, the top of the image will show the southern information and the bottom will display the northern information. Since the spacecraft is looking west to east, the near-range will appear on the left-side of the image and the far-range areas will be on the right-side of the image.  and Processed with ESA Sentinel-1 tools (c) radar images. Note the raw images are geographically inverted due to the data acquisition process shown in (d). The columns of each raw image correspond to increasing the range distance from the spacecraft, whereas the rows of each image correspond to the spacecraft movement along its orbital track.
We de-speckled the images with a Lee Sigma filter [33]. The values used for despeckling (i.e., one look, a sigma value of 0.9, a 7 × 7 pixel window, and a 3 × 3 pixel target window) were obtained by trial and error, similarly to Foucher et al. [34]. Since terrain features influence the local incidence angle, the terrain flattening process normalizes backscatter from the ellipsoid based σ 0 to terrain flattened γ 0 [9]. The terrain impact is subsequently corrected with the assistance of an existing DEM, such as the one obtained from the Shuttle Radar Topography Mission [35]. The final output of the pre-processing workflow is an imaged terrain, with each pixel representing the γ 0 backscatter ( Figure 4c). The spatial resolution is 10 m with a horizontal accuracy of 0.2 ± 0.24 m in the slant range and 0.31 ± 1.31 m in the azimuth range [36]. We co-registered the final radar images using a bilinear geometric transformation, implemented in ArcGIS (ESRI 2018).
The data used for analysis contains the γ0 stacked by recording time (i.e., each pixel represents a point in space with multiple values of γ0 over time). We compared the radar images from the dry periods (when total rainfall accumulation <0.5 mm for two weeks preceding radar acquisition) with the images taken at varying degrees of wetness, as measured by the weather station. We expressed the wetness as total rainfall accumulation for 15, 30, 45, 60, and 90 min as well as hourly accumulations from 2 to 12 h prior to the radar image acquisition.
Considering that the radar images are not acquired continuously (i.e., more than one week but less than two weeks apart), they are suitable for models that operate with time steps of at least one week. A weekly time interval is not necessarily long in forestry, as most forestry models, such as Forest Vegetation Simulator [37] or PTaeda [38], do not operate at shorter time intervals, with the most common length being yearly. Even popular physiological models, like SORTIE-ND [39] or FORECAST [40], use a yearly timestep rather than sub-monthly. However, models like 3PG [41] or iLand [42], which and Processed with ESA Sentinel-1 tools (c) radar images. Note the raw images are geographically inverted due to the data acquisition process shown in (d). The columns of each raw image correspond to increasing the range distance from the spacecraft, whereas the rows of each image correspond to the spacecraft movement along its orbital track.
We de-speckled the images with a Lee Sigma filter [33]. The values used for despeckling (i.e., one look, a sigma value of 0.9, a 7 × 7 pixel window, and a 3 × 3 pixel target window) were obtained by trial and error, similarly to Foucher et al. [34]. Since terrain features influence the local incidence angle, the terrain flattening process normalizes backscatter from the ellipsoid based σ 0 to terrain flattened γ 0 [9]. The terrain impact is subsequently corrected with the assistance of an existing DEM, such as the one obtained from the Shuttle Radar Topography Mission [35]. The final output of the pre-processing workflow is an imaged terrain, with each pixel representing the γ 0 backscatter ( Figure 4c). The spatial resolution is 10 m with a horizontal accuracy of 0.2 ± 0.24 m in the slant range and 0.31 ± 1.31 m in the azimuth range [36]. We co-registered the final radar images using a bilinear geometric transformation, implemented in ArcGIS (ESRI 2018).
The data used for analysis contains the γ 0 stacked by recording time (i.e., each pixel represents a point in space with multiple values of γ 0 over time). We compared the radar images from the dry periods (when total rainfall accumulation <0.5 mm for two weeks preceding radar acquisition) with the images taken at varying degrees of wetness, as measured by the weather station. We expressed the wetness as total rainfall accumulation for 15, 30, 45, 60, and 90 min as well as hourly accumulations from 2 to 12 h prior to the radar image acquisition.
Considering that the radar images are not acquired continuously (i.e., more than one week but less than two weeks apart), they are suitable for models that operate with time steps of at least one week. A weekly time interval is not necessarily long in forestry, as most forestry models, such as Forest Vegetation Simulator [37] or PTaeda [38], do not operate at shorter time intervals, with the most common length being yearly. Even popular physiological models, like SORTIE-ND [39] or FORECAST [40], use a yearly time-step rather than sub-monthly. However, models like 3PG [41] or iLand [42], which operates on monthly time steps, could benefit from the information extracted from the radar images, as their environmental parameters will be more accurately estimated.

Lidar
The results of Saatchi et al. [43] suggested that pixel level radar data have limited usage in forestry applications. Therefore, to enhance the information supplied by radar and weather data, we included lidar data, which was obtained from the flights executed in 2016 [44]. The system used was a Leica ALS80 mounted on a Cessna Grand Caravan, with an average of 0.89 ground points per square meter. The vertical accuracy of the point clouds is ≤9.25 cm and the horizontal accuracy is ≤9.25 cm [45]. We incorporated lidar data in the analysis to accurately estimate the size of the tree crowns (e.g., total height or crown length), which are correlated with the amount of surface canopy water storage. Therefore, we hypothesized that tree dimensions would support radar images and weather data in estimating the water within the canopy. We did not explicitly consider the scattering associated with the trunk, or branch size as well as their arrangement because it has been shown that the needles and moisture of the canopy are responsible for the majority of scatterers in the C-band [5,17,18].
The lidar data includes the point clouds and two derived products, namely the digital surface model (DSM, the upper surface of the canopy) and the digital terrain model (DTM). We created the Canopy Surface Model (CSM) by subtracting the DTM from DSM. Field measurements using a clinometer revealed that the mean crown length was approximately 52% of the total tree height (N = 20, mean = 0.52, variance = 0.016): where L C is crown length (in meters), measured from the base of the crown to the top of the tree, and H total is the total tree height in meters.

Tree Segmentation
Individual tree crowns were extracted from the point cloud using a canopy delineation software called TrEX [46]. The centroid point of each segmented tree crown was used to extract information from the radar and topographic rasters (e.g., backscatter, aspect, and slope). The tree crowns were placed in three height classes: large (total height > 45 m), medium (total height between 30 m and 45 m), and small (total height <30 m). We considered total height in the analysis as it reflects age and size variability, which impacts radar backscatter. We assessed the accuracy of TrEx using a sample of 1.0 ha comprised from 10 circular plots of 1000 m 2 (i.e., between 10% and 21% from the stand surface), which were randomly distributed within each stand ( Figure 5). Each tree crown within the 1000 m 2 plot was visually identified from the lidar point cloud and the corresponding canopy height model. Three metrics were used to assess the TrEx performance in identifying individual trees, namely the accuracy (i.e., the number of existing trees identified by TrEx), the omission error (i.e., number of existing trees that were not found by TrEx), and the commission error (i.e., number of trees found by TrEx but not existing in reality), similar to Kaartinen et al. [47] and Ene et al. [48].

Evaporation Modelling
To estimate the amount of canopy water, the rainfall and evapotranspiration rate is needed. We evaluated evapotranspiration (the combination of evaporation and transpiration) using three different models: (Priestley-Taylor [49], Penman-Monteith [50], and Brutsaert [51], similarly to Ai et al. [52]). We considered the three evaporative models as they have been shown to approximate wet-canopy evaporation as well as dry evapotranspiration [52,53]. The evapo-transpirative loss was subtracted from the accumulated rainfall (AMC) to compute the SCW (Equation (3)).
where i represents the evapotranspiration model for a given time interval as calculated in Equation (1) where ETPT is evapotranspiration in mm, αe is a parameter of the Priestley Taylor equation (unitless), Δ is the rate of change of the saturation specific humidity (unitless), γ is the psychrometric constant [kPa], Rn is the net incoming solar radiation [W/m 2 ], and G is the surface heat flux [W/m 2 ]. We used a value of 1.05 for the αe term, as it is the recommended approximation for coniferous forests [3]. To calculate Δ, we used the equation from Murray [54].

Evaporation Modelling
To estimate the amount of canopy water, the rainfall and evapotranspiration rate is needed. We evaluated evapotranspiration (the combination of evaporation and transpiration) using three different models: (Priestley-Taylor [49], Penman-Monteith [50], and Brutsaert [51], similarly to Ai et al. [52]). We considered the three evaporative models as they have been shown to approximate wet-canopy evaporation as well as dry evapotranspiration [52,53]. The evapo-transpirative loss was subtracted from the accumulated rainfall (AMC) to compute the SCW (Equation (3)).
where i represents the evapotranspiration model for a given time interval as calculated in Equation (1) (i.e., 15 min, 30 min, 45 min, etc.).

Priestley-Taylor Model
The evapotranspiration according to Priestley-Taylor, ET PT , used in this study is: where ET PT is evapotranspiration in mm, α e is a parameter of the Priestley Taylor equation approximation for coniferous forests [3]. To calculate ∆, we used the equation from Murray [54].
where T t is the temperature in degrees Celsius at time t. The temperature was measured from the weather station at 5-min intervals over the course of three years (1 October, 2014 to 21 November, 2017). For net solar radiation, we used: where R n is the net radiation, R in is the incoming radiation, and R out is the outgoing (longwave) radiation [55]. The incoming solar radiation was measured by a nearby weather station (CENMET) and outgoing longwave radiation was calculated using the Stefan-Boltzman law.
Although an additional term for cloud cover is present in some version of the model [55], we excluded the effect of clouds on the longwave radiation, as no reliable cloud data was available. The effect of albedo was applied to the incoming solar radiation, using a value of 0.12 W m −2 , as recommended by Ogunjemiyo et al. [56].
The final term in the Priestley-Taylor equation is G, which is the ground heat flux. We have used a simplified version of Bennett et al. [57], which was concentrated on the heat flux at the surface (Equation (8)).
where I is the thermal inertia of the bulk soil material, which, according to Bennett et al. [57], is 0.66 J m −2 K −1 s −1/2 for evergreen coniferous forests.
where ET P is evapotranspiration in mm, ∆ is the slope of the saturation vapor pressure curve, γ is the psychrometric constant, R n is net radiation, and G is the surface ground heat flux, which are all as calculated above. Additionally, we have the wind function f e (u 2 ), and e * 2 and e 2 , which are the saturation vapor pressure and actual vapor pressure, respectively, at 2 m height above the surface. The wind function is the same as the simplified function used in Ai et al. [52] and Zheng et al. [58].
where u 2 is the wind speed at a height of 2 m. As described in Ai et al. [52] and Zheng et al. [58], the simplified version of the wind functions performs similarly to the more complex form, and, as such, we use the simple form in this study as well. The wind detector for this study was at an altitude of 4.5 m, and, although it was not directly located inside the stands under study, this calculation was, nonetheless, included. Since we do not have measurements of the stomatal conductivity, here, we use the original Penman equation.

Brutsaert Model
In addition to Priestley-Taylor and Penman equations that estimate evapotranspiration, we used a third nonlinear equation [51], which was explored further in different forest types [52] (Equation (11)).
where ET B is evapotranspiration according to Brutsaert, E po is the potential evaporation, and E pa is the apparent potential evaporation. Similar to Ai et al. [52], we used the Penman equation for E pa and the Priestley Taylor equation for E po . To transform the models from daily to 15-minute calculations to conform with the weather station data, the E, E po , and E pa values were divided by 96.

AMC, Time since Rain, and Time since Dry
In addition to the modelling of evapotranspiration, we also calculated the water contained in the canopy as a rolling sum of the precipitation measurements from the weather station. We used different sizes of moving windows in 15-min intervals (e.g., 15, 30, 45, and 60 min) and hourly windows (e.g., 1, 1.5, 2, 3, 4, 5, 6, and 12 h) and summed the values of precipitation within that window. This resulted in a value that reflected the accumulated rainfall, such that: where A T is the accumulated rainfall for a time window T, which is one of several (15, 30, 45, etc.) minute intervals. R i is the incrementally measured rainfall, which is measured by the weather station in 5-min increments. For example, the 1-h AMC used a moving window of width 12, meaning 12 5-min intervals were summed (i.e., i = 1 . . . 12, then 2 . . . 13, 3 . . . 14, etc.) for the length of the entire measurement schedule, which, as mentioned before, ran from 1 October, 2014 to 21 November, 2017. In addition to the evapotranspiration and rain accumulation, we calculated time since rain (TSR) and time since dry (TSD). The TSR is the elapsed time since a rain stopped whereas the TSD is the time passed since the rain started (Equation (13)). We computed TSR and TSD using a minimum rainless or rainy period of 30 min.
For TSR, we set a threshold of less than 0.02 mm in 30 min as being 'dry,' which means that the AMC for the last 30 m, AMC 30 , is ≤0.02 mm. Similarly, for TSD, values of AMC 30 greater than 0.02 mm indicate a 'wet' threshold. This gave us two additional variables that were informative of the timing of a given radar measurement in relation to rain events.

Data Analysis
We initially assessed the changes in the SCW and backscatter against tree height and TSR with a simple ANOVA. Based on the findings of Dobson et al. [20], we assumed that different tree heights would have a different backscatter not only in magnitude but also in variability. Therefore, we compared the variability within the data using Levene's Test. The comparisons were executed across vegetation height classes in 1, 2, 3, 4, and 5-meter increments as well as on the orbital path (ascending or descending).
We analyzed the relationship between SCW and backscatter, TSR, and stand development with a general linear model (GLM) of the form.
where the path represents the flight direction (ascending or descending), γ is the backscatter, TSR is time since rain, and RelHum is relative humidity. The interaction between backscatter and TSR adjusts the radar reflectance in 15-min increments. The interaction term creates a set of curves that are parallel in a non-Euclidean way, which create a range of values rather than the individual curves described by backscatter and TSR. Temperature and relative humidity are the values measured from the H15MET weather station at the time of image acquisition. The 'dry' dates (i.e., no rain for more than 6 h) were dropped from analysis, as they were not pertinent to modeling wet canopy conditions. In developing the models, we used two type of inputs for backscatter: individual pixels and the average of the pixels falling within the canopy of segmented trees. In addition to Equation (14), we have constructed a multitude of linear models that assessed the relationship between SCW and some of the popular dendrometric variables, such as height, or meteorological variables, such as wind speed. However, no other variables add significant information in terms of predictive ability besides the five considered in Equation (14).
The GLM assumptions of normality, heteroskedasticity, and independence were tested using a Shapiro-Wilks test, Levene's test, and a Durbin-Watson test, respectively [59]. The differences between the amounts of SCW at different times since the rain stopped was assessed using least-square means and Tukey's multiple comparisons test.
In the eventuality that the acquisition time is significant, we will develop new models for each orbital path, to alleviate the difference in the number of rain events captured at dusk and dawn (i.e., ascending vs. descending paths). The new models would be built from Equation (14) with an emphasis on nonlinear relationships rather than linear relationships. The metrics used to assess the new models would be bias, root mean square error, and coefficient of determination.

Backscatter and Tree Height
Although the stands were comprised of different sizes of trees, they behaved in a similar way with regard to their mean backscatter (Figures 6-8 and Table 2), which decreased in the overall brightness with temperature ( Figure 7) and during rain events ( Figure 9). The decrease in brightness showed a substantial spread in backscatter values (~5 dB). Stands containing larger trees displayed a smaller variation in backscatter overall than stands with smaller trees (Figure 8).

Backscatter and Tree Height
Although the stands were comprised of different sizes of trees, they behaved in a similar way with regard to their mean backscatter (Figures 6-8 and Table 2), which decreased in the overall brightness with temperature ( Figure 7) and during rain events (Figure 9). The decrease in brightness showed a substantial spread in backscatter values (~5 dB). Stands containing larger trees displayed a smaller variation in backscatter overall than stands with smaller trees (Figure 8).         Figure 9. Comparison of backscatter during or not during a rain event. Note that there were only six images captured during a rain event. Although we found significant differences (p < 0.001) in backscatter for all tree height increment levels (1-5 m intervals), the differences in R 2 values between each class were minute (R 2 < 0.01) when comparing models that used tree height against models without a height variable. Similar results were obtained when the vegetation was grouped into the broad classes (i.e., small, medium, and large), as the Tukey-Kramer test revealed significant differences between classes but with a dismal coefficient of determination. The model as a whole performs only slightly better when tree height characteristics (such as total height, crown depth, or crown base) are included in the model with an increase in overall R 2 of only about 0.01 with the inclusion of these measures.
The canopy-level analysis, which complemented the pixel-level analysis, depends on the tree canopy segmentation algorithm. We found that the automatic procedure identified more than 75% of the tree canopies, irrespective of the stand (Table 3). TrEx was documented to have reduced commission errors but larger omission errors [60], which Although we found significant differences (p < 0.001) in backscatter for all tree height increment levels (1-5 m intervals), the differences in R 2 values between each class were minute (R 2 < 0.01) when comparing models that used tree height against models without a height variable. Similar results were obtained when the vegetation was grouped into the broad classes (i.e., small, medium, and large), as the Tukey-Kramer test revealed significant differences between classes but with a dismal coefficient of determination. The model as a whole performs only slightly better when tree height characteristics (such as total height, crown depth, or crown base) are included in the model with an increase in overall R 2 of only about 0.01 with the inclusion of these measures.
The canopy-level analysis, which complemented the pixel-level analysis, depends on the tree canopy segmentation algorithm. We found that the automatic procedure identified more than 75% of the tree canopies, irrespective of the stand (Table 3). TrEx was documented to have reduced commission errors but larger omission errors [60], which was confirmed by assessment plots, as only one out of the six stands had accuracy less than 80%, supporting the usage of the segmented crowns in the subsequent analyses.

GLM Model
To develop the SCW model, we examined TSR of up to three hours but did not extend TSR further in the model, as the backscatter signal plateaued after this time. Our findings agree with Keim and Link [61], which found that mean residence time of intercepted rainfall does not exceed 1 h. Furthermore, our results indicated that, after 30 min of the end of a rain event, the predictive capability of our model dropped substantially. Additionally, the models that were based on individual pixels of the radar data performed poorly when compared to those using a tree-level segmented from lidar.
The GLM model using the simplified measure of AMC (Equation (14)) has an R 2 value of 0.686, which was the highest among all the models tested in this study, and all terms are deemed significant (p < 0.001). All the stands behave similarly with regard to how long it has been since the last rainfall, despite having different canopy heights and depth ( Figure 6). We also find that the change in backscatter from rainfall disappears between 15 and 30 min, and the backscatter returns to normal 'dry' levels.
Since the interaction term was found to be significant, we executed two separate GLM analyses, including one for the ascending orbital path and one for the descending path. However, there were insufficient wet ascending passes, so only descending passes were used in the modelling. Unfortunately, the linear model that pointed toward the difference in the path estimates (Equation (14)) did not provide encouraging results. Therefore, we developed a new GLM (Equation (15)) that can be used for predicting the surface canopy water, as it exhibits a significant improvement in precision compared with Equation (14).
where wind is the accumulated windspeed of the previous 45 min. All the variables were significant (p-value < 0.01), irrespective of the model used to account for evapotranspiration (i.e., Brutsaert, Penman, and Priestley-Taylor), with an R 2 of 0.64 (i.e., Penman adjustment had the largest R 2 among the three model). The significance of the variables was also true for the unadjusted SCW, which had the largest overall R 2 of 0.69. The situation was mirrored by the root means square error, which ranges between 0.155, for the unadjusted SCW model, to 0.689, for the Priestley-Taylor adjustment. The bias for all models was 0, as they were obtained as a linear combination of nonlinear functions. Considering the short amount of time after the rain for which the SCW model has been developed, it is expected that the adjustments should be minimal. Nevertheless, all three adjustments seem to increase data variability, which is reflected in the reduced ability of the models to capture the relationship between SCW and radar backscattering, TSR, Temperature, relative humidity, and windspeed. Therefore, we considered the unadjusted SCW model is more suitable for prediction, than the Brutsaert, Penman, and Priestley-Taylor rectifications. The SCW computed considering the evapotranspiration were not significantly different among themselves (p > 0.1).

Orbit Path
The Levene's Test for homogeneity of variance and Welch's test found a significant difference (p < 0.0254) in backscatter of a different stand size between the ascending and descending orbital path. The mean backscatter for ascending orbits was −7.848 dB (standard deviation 0.957 dB) and for descending orbits was −7.865 dB (standard deviation 1.115 dB).
The clearing area revealed a significant difference (p < 0.0001) between the ascending and descending orbit, showing that dawn and dusk have different backscatter behavior (i.e., for ascending passes, it is −8.404 dB, whereas, for descending passes, it is −7.246 dB).
The golf course, as expected, displayed a narrower variation in backscatter, even though the mean showed a larger degree of variability than the stands in the HJA. The ascending and descending passes exhibited mean backscatter differences that were significantly different statistically speaking, but still relatively close (i.e., −15.918 dB and −15.683 dB for the ascending and descending passes, respectively, p < 0.0002).
ANOVA on temperature effects on backscatter showed a significant difference between orbital paths (p < 0.001), but the R 2 , with a value of 0.228, was not deemed to be useful.
Considering that radio waves are differentially reflected by wet surfaces [27,29], we tested whether or not the backscatter changes from dry to rainy conditions. We found the backscatter to be significantly different during the rainy periods from dry periods (p < 0.0001), but the R 2 was a mere 0.078 ( Figure 9).

Discussion
Backscatter appears to be correlated with the presence of rain as expected (i.e., when the image is taken during or immediately after a rain event), but a stronger correlation was exhibited between the time of day when the image was acquired (i.e., ascending or descending orbital pass). The backscatter from the clear-cut areas displays an even greater degree of variability than the vegetated stands. The increase in variability associated with clear cuts is most likely due to the lack of scattering by vegetation cover and abundance of woody debris, which create a surface dominated by complex scatterers and negligible attenuation.

Orbital Path and Local Time of Acquisition
The orbital path at acquisition (either ascending or descending) clearly plays a critical part in the overall backscatter of the image. Although all images are corrected for terrain effects, there is not a generalized calibration for diurnal environmental effects, such as dew and atmospheric conditions. The ascending passes, taken during local dusk, have higher backscatter overall. Although this effect is statistically significant, it is not large enough to produce a practically significant change ( Figure 8). Compared to local dusk during the summertime, local dawn in the Pacific Northwest typically has a lower temperature and more dew present, which could account for the decrease in backscatter. After correcting for dates with no prior rainfall, the difference between ascending and descending images is reduced even further.
Although the temperature did not heavily influence the backscatter directly, its effects on the physical processes governing water storage and uptake may have led to the increase in radar brightness at lower temperatures. Although sub-freezing and near-freezing temperatures were excluded, the effects on water droplet shape and canopy storage capacity with temperature [62] are likely to be a factor in the reflectivity of water to C-band microwaves.

AMC and Evapotranspiration Adjusted AMC
The unadjusted AMC model performed slightly better than the models, which were adjusted using evapotranspiration methods. Although the difference in model predictability is surprising, such disparity in the two approaches to moisture modelling suggests that evaporative loss is not as crucial in estimating the effect of precipitation on the C-band radar backscatter as the precipitation itself. There are two other possible interpretations of the weaker ability to predict SCW when evapotranspiration adjustments are considered. First, the parametrization is not appropriate, which likely is not the case as all the values were taken from similar studies, and, second, the adjustments are not warranted in any climate or problem, which is likely the situation of our study. The results of our research mirror the finding of Paun et al. and Seppelt and Richter [63,64], which pointed towards algorithms' impact on the estimates, and Strimbu et al. [65], which found that modeling details depends on the studied problem.
As shown by the difference between the two approaches to precipitation modelling, the timing of rain events is important to the effect on radar imagery. We found that approximately 30 min after a rainfall ends, the change in backscatter from moisture is negligible. Our findings suggest that the six-day temporal resolution of Sentinel-1 does not produce radar images sensitive enough to adequately detect rainfall for the dry summers of the Pacific Northwest. We should point out that the results are based on a relatively reduced dataset, as for the three summers (April to October 2016-2018), only six dates out of a total of 67 satellite passes occurred during rain events including three in the morning and three in the evening. Our study provides evidence that SCW can be modeled using radar imagery, but, to develop a robust and reliable model, more data are needed. The GLM indicates that identification of a formal relationship between SCW, radar imagery, and other relatively easy-to-measure attributes of freely available data is possible.
Additionally, the factor of wind on evapotranspiration right before image acquisition was crucial to modelling evapotranspiration and, thus, the effect of moisture on backscatter [66]. We found that windspeed between 30 and 45 min prior to imaging had a pronounced effect on the model. However, beyond 45 min, the effect is lost. Since we were using timeframes of 15-min intervals in measuring the wind, we could not pinpoint the optimum time to investigate. Further study should not rule out windspeed, and we recommend narrowing down the time window that is most crucial.
The main finding of our study is that foliar SCW should be modelled using values relevant to this investigation, namely the tree and not the radar raster pixel as the elementary unit. We were able to provide evidence of a relationship when we changed the perspective of the investigation from pixel to tree. The change in modeling SCW was possible because of the inclusion of canopy lidar models, which proved to be crucial in our modeling of SCW.

Conclusions
Our study provides evidence that C-band radar is sensitive to summertime tree canopy moisture in the hilly Pacific Northwest forests. Therefore, it can be used to estimate the surface canopy water. Furthermore, we found that values during or within, at most, a halfhour of rain events are the most helpful, since the leaves and branches have a similar radar signature to dry leaves or branches after 30 min because the water evaporated, dropped, or was absorbed. Local time is important, as summers in the Pacific Northwest, particularly in the Cascades, are very hot, which suggests that dusk acquired data (i.e., ascending path of the satellite) requires short intervals after the rain stopped because of increase evaporation compared with dawn satellite passes (i.e., descending path).
The limited summertime rain events corresponding to radar acquisition led to the development of a predictive model for the amount of the surface canopy water only for the descending path of the satellite (i.e., local dawn). The model is nonlinear and predicts surface canopy water from radar backscatter, time since rain, relative humidity, and wind speed. We were able to reach operational results (i.e., R 2 = 0.69) because we changed the unit of analysis from a radar pixel to an individual tree canopy. Therefore, we suggest that further forest studies using radar imagery should use a combination of radar with lidar or high resolution multispectral derived products (such as a projected tree crown).
This study is the first of its kind in the Pacific Northwest and can serve as the basis of an array of projects in radar remote sensing. To begin with, this study explicitly avoided wintertime observations due to the previously mentioned issues with water-ice and wetsnow backscatter. Further investigation in the Pacific Northwest forests should expand the space-time components of this study, beginning with an inquiry into the wintertime behavior of canopy backscatter as well as examining summertime backscatter beyond the HJ Andrews forest. The climatological duality of the Pacific Northwest, both temporally (wet winters and dry summers), and spatially (coast, valley, and the Cascade Mountains) suggests limited generalization abilities, such as findings from one area or season may not be applicable to another area or season. Nevertheless, our study provides a framework for development of models among regions and seasons.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
AMC Antecedent Moisture Condition, in this paper, is taken to mean the cumulative amount of water over a given period of time (e.g., AMC 30 is the accumulated precipitation over 30 min).

ET
EvapoTranspiration, the combined effect of evaporation and transpiration that accounts for the loss of AMC over time.

ESA
European Space Agency LUT Look Up Table, each Sentinel-1 Level-1 product includes in the metadata a look up table with calibration operators specific to that image.

PNW
Pacific NorthWest, the region of North America comprising Northern California, Oregon, Washington, and British Colombia. This region is typified by thick evergreen forests and high annual rainfall.

SAR
Synthetic Aperture Radar, a form of radar that uses either (1) the motion of a sensor over time, or (2) radar images taken at different times to produce a radar image of finer spatial resolution than traditional beam-scanning radar. SAR is commonly used in most modern radar sensors.

SCW
Surface Canopy Water, how much water is contained within the canopy on the surface (as opposed to water contained within the biological structure of the canopy components) SNAP SeNtinel Application Platform, a computer program for processing ESA Sentinel mission images TSD Time Since Dry, a variable that indicates how long it has been since a rain event began (defined as <0.02 mm rainfall in 30 min) TSR Time Since Rain, a variable that indicates how long it has been since a rain event ended (defined as >0.02 mm rainfall in 30 min) VH Vertical Horizontal polarization, one of several polarization modes used in radar imaging. This means vertically polarized microwaves are sent out through the antenna and only horizontally polarized backscatter is measured.
VV Vertical-Vertical polarization means vertically polarized signals are sent by the antennae, and vertically polarized backscatter is measured.