Hydrological Mapping in the Luquillo Experimental Forest: New Local Datum Improves Watershed Ecological Knowledge

: Streams and rivers of the Luquillo Experimental Forest, Puerto Rico, have been the subject of extensive watershed and aquatic research since the 1980s. This research includes understanding stream export of nutrients and coarse particulate organic matter, physicochemical constituents, aquatic fauna populations and community structure. However, many of the streams and watersheds studied do not appear in standard scale maps. We document recent collaborative and multi-institutional work to improve hydrological network information and identify knowledge gaps. The methods used to delimit and densify stream networks include establishment and incorporation of an updated new vertical datum for Puerto Rico, LIDAR derived elevation, and a combination of visual-manual and automated digitalization processes. The outcomes of this collaborative effort have resulted in improved watershed delineation, densiﬁcation of hydrologic networks to reﬂect the scale of on-going studies, and the identiﬁcation of constraining factors such as unmapped roadways, culverts, and other features of the built environment that interrupt water ﬂow and alter runoff pathways. This work contributes to enhanced knowledge for watershed management, including attributes of riparian areas, effects of road and channel intersections and ridge to reef initiatives with broad application to other watersheds.


Introduction
Hydrological maps are used for decision-making, and there is an implicit expectation that they are accurate representations of the landscape conditions they illustrate. For those who work along forested montane streams, hydrological maps have until recently only partially conveyed information about perennial and non-perennial channels in the landscape. Mapping of freshwater resources and their associated landscape features, including potential catchment areas, small-scale watersheds, and floodplains is of vital importance to understanding the Caribbean region [1,2] ridge to reef ecosystem connections and more globally for developing research on disturbance along hydrologic networks [3]. Freshwater in the Caribbean is a complex resource that can be scarce and limited for consumption [4], affect forest ecosystems and aquatic fauna [5,6], or become a management challenge during storms and floods [7]. In mountainous areas, headwater streams are prone to flash floods with a rapid onset of high flows and then an equally rapid return to baseflow conditions [8]. It is precisely these smaller order streams in headwater systems that are often not documented in standard or continental scale (1:20,000) maps, despite their global importance to ecosystem processes such as greenhouse gas exchange [9]. The maps at 1:20,000 were meant for a specific designation or fitness for use because of the scale and the associated level of detail. These maps were valid for the municipal and regional levels, and not meant for specific sites such as research watersheds, and their associated hydrological networks. As the LEF has been under forest cover consistently and under cloud cover often, mapping this steep and complex terrain has always been a challenge for land surveyors to conduct field work. With the availability of improved and accessible remote sensing technologies (e.g., Light Detection and Ranging, LiDAR) in combination with long-term stream ecology research, we are now able to rethink our understanding of these aquatic systems and move forward to identify gaps in the research needed to improve watershed ecological knowledge.

Study Site
The Luquillo Experimental Forest (LEF) is located in the Luquillo mountain range, northeastern Puerto Rico (Figure 1). The LEF is home to the Luquillo Long-Term Ecological Research (LUQ-LTER) program and the Luquillo Critical Zone Observatory (LCZO), both of which have been funded by the US National Science Foundation (NSF) [10,11]. The 11,330 hectares (28,000 acres) area is the only tropical forest in the USDA Forest Service National Forest System and bears two official names: El Yunque National Forest and LEF, which refers to the forest's additional designation as a research experimental forest [10], encompassing the entire National Forest and nearby forested and developed areas as a research area.
Hydrology 2021, 8, x FOR PEER REVIEW 2 of 12 maps at 1:20,000 were meant for a specific designation or fitness for use because of the scale and the associated level of detail. These maps were valid for the municipal and regional levels, and not meant for specific sites such as research watersheds, and their associated hydrological networks. As the LEF has been under forest cover consistently and under cloud cover often, mapping this steep and complex terrain has always been a challenge for land surveyors to conduct field work. With the availability of improved and accessible remote sensing technologies (e.g., Light Detection and Ranging, LiDAR) in combination with long-term stream ecology research, we are now able to rethink our understanding of these aquatic systems and move forward to identify gaps in the research needed to improve watershed ecological knowledge.

Study Site
The Luquillo Experimental Forest (LEF) is located in the Luquillo mountain range, northeastern Puerto Rico (Figure 1). The LEF is home to the Luquillo Long-Term Ecological Research (LUQ-LTER) program and the Luquillo Critical Zone Observatory (LCZO), both of which have been funded by the US National Science Foundation (NSF) [10,11]. The 11,330 hectares (28,000 acres) area is the only tropical forest in the USDA Forest Service National Forest System and bears two official names: El Yunque National Forest and LEF, which refers to the forest's additional designation as a research experimental forest [10], encompassing the entire National Forest and nearby forested and developed areas as a research area.

Data Needs for Densification of the Hydrographic Network
The combined investment of various agencies through time has made possible the current and ongoing work of hydrological network densification taking place in the LEF. The cumulative process among State and Federal agencies and integration of products to establish the baseline for producing the densification of the hydrographic network started two decades ago ( Figure 2). In the past, to generate detailed topographic maps involved methods such as conducting plane table surveying and taking elevations in the field. Although these resulted in maps with an appropriate level of detail, this was cumbersome as it required moving people and equipment through forested areas to take x, y, z readings plus bearing angles. In a steep and forested landscape such as the LEF this work could take years from the initial field campaign to data compilation, and the final actual placement on a map at the scale of 1 in 20,000. Starting in 2002, the National Geodetic Survey (NGS-NOAA) began to carry out the geodetic leveling process for Puerto Rico [12]. This geodetic leveling process was financed by three agencies in Puerto Rico; PR Office of Management and Budget (agency initials in Spanish/English, OGP/PROMB); Electric Power Authority (AEE/PREPA); Department of Transportation and Public Works (DTOP/PRDOTPW), and the Engineers and Land Surveyors Corps [13].

Data Needs for Densification of the Hydrographic Network
The combined investment of various agencies through time has made possible the current and ongoing work of hydrological network densification taking place in the LEF. The cumulative process among State and Federal agencies and integration of products to establish the baseline for producing the densification of the hydrographic network started two decades ago ( Figure 2). In the past, to generate detailed topographic maps involved methods such as conducting plane table surveying and taking elevations in the field. Although these resulted in maps with an appropriate level of detail, this was cumbersome as it required moving people and equipment through forested areas to take x, y, z readings plus bearing angles. In a steep and forested landscape such as the LEF this work could take years from the initial field campaign to data compilation, and the final actual placement on a map at the scale of 1 in 20,000. Starting in 2002, the National Geodetic Survey (NGS-NOAA) began to carry out the geodetic leveling process for Puerto Rico [12]. This geodetic leveling process was financed by three agencies in Puerto Rico; PR Office of Management and Budget (agency initials in Spanish/English, OGP/PROMB); Electric Power Authority (AEE/PREPA); Department of Transportation and Public Works (DTOP/PRDOTPW), and the Engineers and Land Surveyors Corps [13]. The geodetic survey of Puerto Rico was the fundamental and critical component to generating a reliable mathematical model (geoid) that would allow the production of elevation data ( Figure 1). The development of the local geoid/geodetic model (Puerto Rico Vertical Datum 2002, PRVD02) was a necessary precursor, that established sufficient elevational benchmarks, to enable collection of and integration with other spatial data resources (e.g., USGS, LiDAR) and create a finely detailed 1 m 2 DEM. After the new vertical datum for Puerto Rico was published in 2010 (PRVD02), the United States Geological Survey (USGS), through the 3D Elevation Program (3DEP [14], financed the generation of elevation data via aerial LiDAR surveys. This survey was able to cover greater than 98% of the Puerto Rican archipelago with only a small portion of the southeast of the main island not surveyed due to constant cloud cover. Once the LiDAR data set was produced, then it could be downloaded from the USGS site, and used to generate the digital bareearth elevation model which has a spatial resolution of 1 square meter per cell. For the first time in Puerto Rico we were able to obtain a digital model of high resolution and accurate elevations, from which multiple levels of data can be extracted ( Figure 3). The geodetic survey of Puerto Rico was the fundamental and critical component to generating a reliable mathematical model (geoid) that would allow the production of elevation data ( Figure 1). The development of the local geoid/geodetic model (Puerto Rico Vertical Datum 2002, PRVD02) was a necessary precursor, that established sufficient elevational benchmarks, to enable collection of and integration with other spatial data resources (e.g., USGS, LiDAR) and create a finely detailed 1 m 2 DEM. After the new vertical datum for Puerto Rico was published in 2010 (PRVD02), the United States Geological Survey (USGS), through the 3D Elevation Program (3DEP [14], financed the generation of elevation data via aerial LiDAR surveys. This survey was able to cover greater than 98% of the Puerto Rican archipelago with only a small portion of the southeast of the main island not surveyed due to constant cloud cover. Once the LiDAR data set was produced, then it could be downloaded from the USGS site, and used to generate the digital bare-earth elevation model which has a spatial resolution of 1 square meter per cell. For the first time in Puerto Rico we were able to obtain a digital model of high resolution and accurate elevations, from which multiple levels of data can be extracted ( Figure 3).

Processing Decisions and Trade-Offs in Generating a Hydrographic Network
Since the elevation model generated based on LiDAR and the new datum provides a high degree of detail, and in the absence of a densified hydrographic product for the LEF area, we decided to experiment by generating a much denser hydrographic network than currently available. This step allowed us to identify and illustrate hydrologic networks where long-term ecological research has been taking place in the LEF [15,16].

Processing Decisions and Trade-Offs in Generating a Hydrographic Network
Since the elevation model generated based on LiDAR and the new datum provides a high degree of detail, and in the absence of a densified hydrographic product for the LEF area, we decided to experiment by generating a much denser hydrographic network than currently available. This step allowed us to identify and illustrate hydrologic networks where long-term ecological research has been taking place in the LEF [15,16].

The Sonadora Sub-Watershed and Stream Systems, a Long-Term Research Area Densification Case Study
Part of the Espíritu Santo Watershed, the Quebrada Sonadora (QS) stream flows from El Yunque Peak (1010 m elevation) through the El Verde Research Area (EVRA). Improved mapping of this sub-watershed and its hydrological network will allow researchers to better apply findings from long term observations to other headwater streams and watersheds in the LEF and elsewhere in the Caribbean. Initial field work in the 1980s on the Quebrada Sonadora, for example, revealed that the 2 stream channels shown on USGS topographic quadrangles [17] greatly underestimate perennial channels in this basin, limiting the ability to characterize the entire drainage network (Figure 4).
Use of a flow accumulation of ≥2000 cells results in a drainage network of just over 1200 elements, identifying both perennial and non-perennial runoff. Using the automated method as a guide, on a flow accumulation of ≥10,000 cells (one cell = 1 m 2 ), the network consists of 123 elements within the sub-watershed. Generating and validating items with this level of detail using a flow accumulation threshold ≥ 2000 cells, and an average scale of 1:200, can be cumbersome and time consuming within each square kilometer.

The Sonadora Sub-Watershed and Stream Systems, a Long-Term Research Area Densification Case Study
Part of the Espíritu Santo Watershed, the Quebrada Sonadora (QS) stream flows from El Yunque Peak (1010 m elevation) through the El Verde Research Area (EVRA). Improved mapping of this sub-watershed and its hydrological network will allow researchers to better apply findings from long term observations to other headwater streams and watersheds in the LEF and elsewhere in the Caribbean. Initial field work in the 1980s on the Quebrada Sonadora, for example, revealed that the 2 stream channels shown on USGS topographic quadrangles [17] greatly underestimate perennial channels in this basin, limiting the ability to characterize the entire drainage network (Figure 4).
Use of a flow accumulation of ≥2000 cells results in a drainage network of just over 1200 elements, identifying both perennial and non-perennial runoff. Using the automated method as a guide, on a flow accumulation of ≥10,000 cells (one cell = 1 m 2 ), the network consists of 123 elements within the sub-watershed. Generating and validating items with this level of detail using a flow accumulation threshold ≥ 2000 cells, and an average scale of 1:200, can be cumbersome and time consuming within each square kilometer.
LiDAR was flown in 2010-2011 as part of the Luquillo Critical Zone Observatory [18] and watershed delineations were generated based on a bare earth Digital Elevation Model (DEM) created from that LiDAR data. This dataset originally had a number of issues such as lack of coverage at high elevations due to cloud cover when the flights were conducted, low LiDAR point cloud densities (also due to cloud cover) and other artifacts. As a result, older USGS quadrangle 10 square meter DEMs were used to fill in the gaps. With the availability of the new USGS LiDAR based on a 1 square meter DEM, these watersheds were delineated again. This resulted in better defined watersheds for many of the sampling sites and corrected flow paths. In some cases, due to the rugged terrain in the LEF, narrow ridges that are the boundaries between watersheds were sometimes missed and flow paths were misrouted when using the 10 m 2 DEM. Using the flow accumulation from the 10 m 2 DEM results in the flow path for Quebrada Prieta A (QPA) being misrouted into the Quebrada Toronja (QT) ( Figure 5). Use of the 1 m 2 DEM successfully identified the ridge that forms the boundary between Quebrada Prieta (QP) and QT, and the watershed flow paths are corrected for the first time in a map product that that agrees with field observations ( Figure 5).  LiDAR was flown in 2010-2011 as part of the Luquillo Critical Zone Observatory [1 and watershed delineations were generated based on a bare earth Digital Elevation Mod (DEM) created from that LiDAR data. This dataset originally had a number of issues su as lack of coverage at high elevations due to cloud cover when the flights were conducte low LiDAR point cloud densities (also due to cloud cover) and other artifacts. As a resu older USGS quadrangle 10 square meter DEMs were used to fill in the gaps. With t narrow ridges that are the boundaries between watersheds were sometimes missed and flow paths were misrouted when using the 10 m 2 DEM. Using the flow accumulation from the 10 m 2 DEM results in the flow path for Quebrada Prieta A (QPA) being misrouted into the Quebrada Toronja (QT) ( Figure 5). Use of the 1 m 2 DEM successfully identified the ridge that forms the boundary between Quebrada Prieta (QP) and QT, and the watershed flow paths are corrected for the first time in a map product that that agrees with field observations ( Figure 5).

Other or Ongoing Challenges of Work within the LEF and Surrounding Landscape
The work of densifying the hydrographic network will not be limited to the research example watersheds. The goal is to cover the entire area of the LEF. After consulting LUQ LTER and USDA Forest Service scientists who conduct research in the study area, a flow accumulation threshold of ≥5000 cells or 0.5 hectare (ha) was established. Flow paths below 0.5 ha could be derived but a flow accumulation of 0.5 ha is likely well below the threshold for perennial streams and channel formation; overland flow still occurs along flow paths with areas below 0.5 ha during high intensity rainfall events [19].
We generated the flow accumulation model process using the ArcGIS Pro 2.4 program: DEM > Depression-less DEM > Flow Direction > Flow Accumulation > Linework derivation (stream to feature).

Other or Ongoing Challenges of Work within the LEF and Surrounding Landscape
The work of densifying the hydrographic network will not be limited to the research example watersheds. The goal is to cover the entire area of the LEF. After consulting LUQ LTER and USDA Forest Service scientists who conduct research in the study area, a flow accumulation threshold of ≥5000 cells or 0.5 hectare (ha) was established. Flow paths below 0.5 ha could be derived but a flow accumulation of 0.5 ha is likely well below the threshold for perennial streams and channel formation; overland flow still occurs along flow paths with areas below 0.5 ha during high intensity rainfall events [19].
We generated the flow accumulation model process using the ArcGIS Pro 2.4 program: DEM > Depression-less DEM > Flow Direction > Flow Accumulation > Linework derivation (stream to feature). This last step, the stream to feature geoprocessing tool converts the flow accumulation raster to vector format from a specified accumulation threshold. Limitations to this method include bridges or road culverts that are not registered in the base map nor in the DEM. However, with the high-resolution DEM we can find grooves that indicate a path or route for flow and manual adjustments can be made.
Raster-based flow accumulation of <1000 square meters (1 m 2 cells) becomes too dense ( Figure 6) and impractical to field validate, and thus the threshold >5000 square meter provides the optimal details needed for current applications for watershed and aquatic research needs. method include bridges or road culverts that are not registered in the base map nor in the DEM. However, with the high-resolution DEM we can find grooves that indicate a path or route for flow and manual adjustments can be made.
Raster-based flow accumulation of <1000 square meters (1 m 2 cells) becomes too dense ( Figure 6) and impractical to field validate, and thus the threshold >5000 square meter provides the optimal details needed for current applications for watershed and aquatic research needs.

Figure 6.
Here we show the raster based flow accumulation down to 300 m 2 , with the manually adjusted stream flow path features in brown. This level of detail is not practical at the landscape scale due to the amount of effort it would require. However, was done here for an intensively studied Quebrada Sonadora watershed.  Figure 7 shows a manually validated hydrological network with a threshold ≥ 5000 accumulated cells (5000 square meters). Previous work in the BEW used a 6-ha threshold [20] from field knowledge. However, it is evident that this threshold is too high in many parts of the LEF and use of a lower threshold such as 5000 cells is more appropriate as shown in Figure 7. Here we can see clear flow paths in the hillshade up to near the watershed boundary. Figure 7 shows a manually validated hydrological network with a threshold ≥ 5000 accumulated cells (5000 square meters). Previous work in the BEW used a 6-ha threshold [20] from field knowledge. However, it is evident that this threshold is too high in many parts of the LEF and use of a lower threshold such as 5000 cells is more appropriate as shown in Figure 7. Here we can see clear flow paths in the hillshade up to near the watershed boundary. To organize and record daily work over the LEF, a grid of one square kilometer was arranged for the study area. Considering the threshold ≥5000 accumulated cells, each 1 km frame can take several hours of computing time and researcher review to process and validate the hydrological network, depending on the number of generated flow paths and flow diversions caused by road/trail surfaces and roadside channels/culverts. The presence of anthropogenic channels and surfaces requires a significant investment of time and researcher validation to interpret and correct flow paths through the built environment. Sometimes it is not possible to find the groove or depression that connects flow networks, making this work susceptible to errors. Follow up field work can help validate these flow path corrections.

Ongoing Work with Densified Hydrological Networks
In the case of the LEF, the densified hydrological networks can be combined with other products like the high-resolution Goddard's Lidar, Hyperspectral and Thermal, G-LiHT, [21] aerial photography from flights flown by NASA in 2017 and 2018 and rendered together using tools such as ArcGIS story maps or journal maps to show more detailed views of study areas [15]. Such story maps are now available on the Luquillo LTER website, illustrating the immediate integration of densified hydrological networks with other ecological site information [https://luq.lternet.edu/ (accessed on 3 January 2021)]. More To organize and record daily work over the LEF, a grid of one square kilometer was arranged for the study area. Considering the threshold ≥5000 accumulated cells, each 1 km frame can take several hours of computing time and researcher review to process and validate the hydrological network, depending on the number of generated flow paths and flow diversions caused by road/trail surfaces and roadside channels/culverts. The presence of anthropogenic channels and surfaces requires a significant investment of time and researcher validation to interpret and correct flow paths through the built environment. Sometimes it is not possible to find the groove or depression that connects flow networks, making this work susceptible to errors. Follow up field work can help validate these flow path corrections.

Ongoing Work with Densified Hydrological Networks
In the case of the LEF, the densified hydrological networks can be combined with other products like the high-resolution Goddard's Lidar, Hyperspectral and Thermal, G-LiHT, [21] aerial photography from flights flown by NASA in 2017 and 2018 and rendered together using tools such as ArcGIS story maps or journal maps to show more detailed views of study areas [15]. Such story maps are now available on the Luquillo LTER website, illustrating the immediate integration of densified hydrological networks with other ecological site information [https://luq.lternet.edu/ (accessed on 3 January 2021)]. More accurate hydrological networks and watershed delineations will aid researchers in identifying study sites, correctly identifying a location within a reach or watershed, and integrating complementary data on aquatic ecosystems and their ecological communities.
Initial work with the densified hydrological networks included analyzing watersheds to determine the percent area of different lithologies [22], mean elevation, mean precipitation [23], mean annual runoff based on estimated annual rainfall [24] as shown in Table 1. Many of these datasets are available on Hydroshare [25,26] and presented as a story map https://arcg.is/1XfiO4 (accessed on 3 January 2021) [15]. Watershed delineations have been completed for more than 36 watersheds above long-term water quality research locations with the aim to better understand the variability in water quality due to factors such as variability in lithology, landcover, forest cover type, slope and other factors.

Opportunities and Challenges to Improve Watershed Ecological Knowledge
A suggested next step is to create a more appropriate definition of flow accumulation areas to establish perennial streams based on concurrent variation in mean annual rainfall across the landscape [24,27]. Establishing a threshold beyond the previously established 6 ha value developed for specific conditions in the watersheds of the BEW [20] is needed to enhance the value of this newly available 1 m DEM. This advances our understanding of headwater streams, which have ephemeral responses to flow and are tightly coupled to catchment hillslope processes, including the accumulation of organic matter and sediments [28,29]. Drought conditions in the Caribbean region and specifically the LEF are predicted to increase in frequency [30,31]. One outstanding issue is how changing environmental conditions may alter the dynamics of riparian forest organic matter inputs to stream ecosystems [32]. As a basal resource for stream fauna communities the effects of accumulation of leaflitter in non-perennial network streams and the pulsed input of nutrients and organic matter when either seasonal or flash-flood event flows are established can now be quantified. For example, in temperate systems, non-perennial streams have been estimated to contribute up to 10% of the daily CO 2 emission from perennial rivers and streams [33]. The newly densified hydrology in the LEF long-term research catchments/watersheds will allow the quantification of ecosystem process such as CO 2 emissions from tropical headwater perennial and non-perennial streams, which are important elements in carbon cycling assessments [9].
Future research in hydrological and ecosystem science will now be able to investigate differences in the distribution and conditions of intermittent and ephemeral streams, stratified by geological substrates, elevation, and forest types. Flash-flood hazard potential and risks assessment identification in headwater areas could also be improved [34]. Additionally, work quantifying connectivity and drought conditions of stream reaches in relation to the ecology of migratory aquatic fauna [4,5,35] can elucidate the role of headwater streams previously not identified in hydrological maps. The re-delineated and densified hydrological networks presented here provide more accurate depictions of discharge and stream habitats, opening up new lines of knowledge development and practical applications.