Geomorphological mapping and anthropogenic landform change in an urbanizing watershed using structure-from-motion photogrammetry and geospatial modeling techniques

ABSTRACT Increasing urbanization and suburban growth in cities globally has highlighted the importance of land planning using detailed geomorphologic maps that depict anthropogenic landform changes. Such mapping provides information crucial for land management, hazard identification, and the management of the challenges arising from urbanization. The development and use of quantitative and repeatable methods to map anthropogenic and natural processes are required to advance the science of urban geomorphological mapping. This study investigated the application of geospatial modeling, structure-from-motion (SfM) photogrammetric methods and DEM differencing as means of quantifying anthropogenic landform changes from archival aerial imagery. Anthropogenic landforms were incorporated into a detailed geomorphologic map in an urbanizing watershed located in the Washington, D.C. metropolitan suburb of Vienna, Virginia.


Introduction
Human activities have a substantial and cumulative effect on the landscape and its landforms (Marsh, 1864). Geomorphologic changes result from a range of anthropogenic activities including forest clearing, agriculture, land draining and filling, mining and quarrying, channelization, irrigation, and construction of dams or other engineering features (Nir, 1983). These activities involve both the intentional removal and deposition of material as well as the unintentional effects of hydrological changes and resulting erosion and sedimentation (Hooke, 2000).
Anthropogenic landforms are created at a variety of spatial scales. Local scale landforms result from excavation, cutting, and grading to modify slopes and drainage patterns to create level ground for development and transportation infrastructure (Nir, 1983). Broader scale landforms are created by mining and quarrying, as well as subsidence due to water or mineral resource extraction. Depositional landforms, or 'man-madeground' typically entail the use of fill material, as well as waste dumps, and may occur at a range of scales . The cumulative geomorphologic effects of anthropogenic activity are most pronounced in areas of dense human occupation (e.g. in urban spaces) and primarily occur in the early-urban to the mid-urban stages of development (Nir, 1983).
With increasing pressure on urban areas, the need for geoscience information plays a crucial role in the evaluation of resource potential and suitability for urban development, identifying hazards and informing management strategies (Cooke, 1976;Fookes et al., 2005;McGill, 1964). Situations in which development decisions were not grounded in detailed geomorphologic information often have catastrophic results including landslides (Alexander, 1989), housing and infrastructure collapse (Chirico & Epstein, 2000), and flooding and hydrological effects (Douglas, 2010). Landform changes, such as those due to road gradation and paving, have significant effects on natural processes which concentrates and alters drainage patterns (Montgomery, 1994).
There is growing recognition of the importance of urban geomorphological and geologic mapping to inform development of urban and suburban fringe areas where prior anthropogenic landform change may have occurred (Brandolini et al., 2017;Del Monte et al., 2016;Douglas, 1999;Faccini et al., 2012;Lee, 2001;Lóczy & Süto, 2011;Robinson & Speiker, 1978). There are also new methods for geomorphological mapping and landform change which are emerging in parallel with enabling technology and datasets (Otto & Smith, 2013;Paron & Smith, 2008). Stereo and monoscopic aerial photography interpretation has consistently been used as a method to map anthropogenic changes and geomorphologic units in urban areas (Del Monte et al., 2016;Maksud Kamal & Midorikawa, 2004). However, traditional stereo compilation methods are labor intensive and time consuming. New structure-from-motion (SfM) photogrammetric techniques offer a more efficient and highly accurate procedure for image matching, triangulation, and digital elevation model extraction from aerial images (Fonstad et al., 2013;Westoby et al., 2012). The development of digital elevation models (DEMs) from historical aerial photo sources for mapping geomorphologic change is a recent advance (Gomez et al., 2015) and could significantly contribute to urban geomorphological mapping applications. However, comparison of DEM data from historical sources and with different resolutions is not always a straightforward process (James et al., 2012;Smith et al., 2015).
The goal of this research is to develop and apply quantitative and repeatable geospatial mapping techniques to advance the science of urban geomorphological mapping. The research presented here specifically applies digital elevation modeling, historical aerial imagery, and SfM photogrammetry to map and integrate anthropogenic, fluvial, and regolith geomorphology for a small urbanizing watershed in the Washington, D.C. metropolitan area.

Study Area
The area selected for this research was the Piney Branch watershed in Vienna, Fairfax County, Virginia, a suburb of Washington, D.C. (Figure 1). Piney Branch is a tributary of Difficult Run which feeds into the Potomac River. The watershed is approximately 10 km 2 in size with 70 m of total relief.

Geology
The watershed is underlain by Early Cambrian and/ or Late Proterozoic metasedimentary, metavolcanic, and transported metaigneous rocks of the Piedmont geologic province, including the Piney Branch Complex and the Peters Creek Schist (Figure 2). The Piney Branch Complex (CZpb) is characterized by an intermixed complex of peridotite, pyroxenite, and gabbro now represented by serpentinite, soapstone, actinolite schist, and metagabbro which constitute an autoclastic mélange (Drake & Morgan, 1981). The unit also contains small sheets and dikes of quartz-albite plagiogranite. The Peters Creek Schist is characterized by fineto coarse-grained, quartzrich schist, and lesser mica gneiss. Interbedded with these rocks are fineto medium-grained, well-bedded metagraywacke and semi-pelitic schist (CZpg), as well as a few layers of calc-silicate rock. Also present is the Old Mill Branch metasiltstone Member of the Popes Head Formation (Cpo), which is characterized by thin-to medium-bedded, medium-to very finegrained, very mature, micaceous metasiltstone. These units are intruded by the Bear Island Granodiorite, a fine-grained, muscovite-biotite, granodiorite and related pegmatite composed largely of quartz. Bear Island Granodiorite holds up irregular elongated ridges roughly perpendicular to the direction of streamflow (see Drake & Lee, 1989). The bedrock units have a significant regolith weathering profile (described in Pavich et al., 1989) and are overlain by quaternary alluvium in the floodplain of Piney Branch and its tributaries.

Land-use and land cover history
Colonial and post-colonial crop and livestock-based small farms began in the area in 1754 and expanded through the 1850s. Colonial-era deforestation and the introduction of European-style agricultural practices resulted in significant erosion of topsoil from deforested hillslopes resulting in deposition of sedimentary material inherited from previous anthropogenic activities or 'legacy sediments' (James, 2013(James, , 2018. Subsequent transport and deposition of legacy sediments into valley bottoms has elevated the floodplain 1-2 m above the previous natural flood plain (Hupp et al., 2013). The establishment of several known mill dams located along Piney Branch may have also compounded development of legacy sediments as noted in other regional research (Hupp et al., 2013;Merritts et al., 2011).
The Washington and Old Dominion (W&OD) railroad was constructed between 1855 and 1860 for regional passenger and freight travel. Running roughly parallel to Piney Branch and bisecting the watershed, the development of this transportation corridor required significant land modification through excavation and fill and resulted in continued land clearing for development along its route. The W&OD rail line operated until the 1960s, when it was purchased by the Virginia Department of Highways and Virginia Electric and Power Company for use as an electric power line transmission right of way. The railroad line itself was eventually converted to a paved recreational rail trail beginning in the late 1970s (Virginia.gov, 2019).
Concurrently, from the late 1940s onward, small farms have gradually been replaced by expanding urbanization and development which is a consistent pattern displayed in other similar urban areas in the US (Hooke, 1999). The watershed is primarily suburban, with housing communities and major transportation corridors. The transition from rural agricultural land use to suburban developed land use began in the downtown core area of Vienna and radiated outward through subsequent decades. The development of medium-or high-density housing was accompanied by significant land modification to excavate basement foundations for houses, to build sewer and stormwater infrastructure, to grade existing slopes, and to grade and pave new roadways. A map of the progression of development on a parcel-by-parcel temporal scale is shown in Figure 3.

Materials and Methods
Due to the varied land-use and land cover transitions that have occurred, multiple data sources were incorporated to create a detailed urban geomorphologic map highlighting anthropogenic landform changes (Main Map).

Data
Stereo aerial imagery and ortho-imagery ranging from 1937 to 2017 were acquired. Examples of this imagery, which documents the 'farm to housing' transition, are shown in Figure 4 for 1950, 1973, and 2018 respectively. Figure 4(B) shows the grading, filling, and excavation of building foundations and basements typical of post-1950 housing construction methods utilized in this area.
Prior to processing, all digital aerial images were cropped using Adobe Photoshop to remove the photo frame and fiducial information. SfM processing was conducted using Agisoft PhotoScan [Metashape v. 1.6] to produce a DSM and orthophoto from each set of images. Table 1 shows the dates and scales of imagery and DEMs used in this study. Eleven ground control points (GCPs) were added to improve photo alignment and to control the final orthophoto and DSM models. Table 2 shows the data and estimated error for the GCPs as modeled in Agisoft Metashape.
The DSM was then filtered to remove vegetation and manmade features above the land surface. This was done in PCI Geomatica v.2018 using the DSM2DTM algorithm. DSM2DTM functions by applying a series of filters to derive a local minimum based on a user-defined window size and compares DSM values to the local minimum or ground surface elevation. Elevation values exceeding the defined ground surface elevation are removed and replaced with the minimum for the patches identified (Barbarella et al., 2017). In this way, the DSM was filtered to produce a DTM and was prepared to be directly compared with the 2012 LiDAR data.
The geologic map of the Vienna Quadrangle in Fairfax County, Virginia (Drake & Lee, 1989) was digitized and used as a data input to modeling regolith thickness in combination with analysis of borehole data, fieldwork, and literature on regolith properties in the larger region (Pavich et al., 1989). Specifically, a database of 20 boreholes at 7 locations around the Piney Branch watershed was developed from previously acquired hollow-stem auger borings. These borehole logs record sediment characteristics, material strength properties, degree of weathering, and depth to un-weathered bedrock. The boreholes range in depth from 6 m to 34 m. Field mapping also included shallow hand-augered sediment samples (locations not shown) which generally penetrated less than 2 m and were confined to floodplain areas. Figure 5 shows the distribution of GCPs, independent check points (ICPs), and borehole data used in this study.
Datasets depicting topographic and other anthropogenic development features such as roads, building footprints, and hydrologic features were acquired from various sources including municipal government agencies.

DEM of Difference (DoD) modeling for anthropogenic landforms
Quantification of anthropogenic topographic change through direct comparison of elevation surfaces is done by differencing two elevation surfaces. In this study, the 1950 SfM DSM was subtracted from the 2012 DEM, resulting in a DEM of Difference (DoD) where positive values represent land accretion and negative values indicate land removal occurring between those time periods (Etzelmuller, 2000). DEMs have unique error characteristics based on the method used to develop them which are compounded when multiple DEMs are directly compared. Therefore, it is necessary to account for the error in both of the input DEMs to ensure that landform changes can be distinguished from the combined error in the DoD dataset (Lane et al., 2003;Burrough & McDonnell 1998). In this study, the accuracy of the filtered SfM DTM created from 1950 aerial imagery was calculated using the Root Mean Squared Error (RMSE) statistic from a set of 32 ICPs at locations that exhibited little evidence of landform change in sequential aerial image analysis. These 32 ICPs were attributed with elevation values from the  Drake & Lee, 1989) 2008 USGS LiDAR DEM. The 2008 DEM was used for the elevations of these points so that the check point elevations were independent of the 1950 and 2012 DEM elevation datasets used in the DoD modeling. Table 3 shows the RMSE accuracy results of the 1950 DSM for the ICP locations.
After RMSE error modeling, a change threshold value was developed using the Root of Sum in Quadrature (RSiQ) method (Chandler, 1999;Lane et al., 2000;Taylor, 1997;Westoby et al., 2012;Wheaton et al., 2010). The RSiQ method characterizes the combined vertical inaccuracies in the DoD, where values above or below the threshold constitute significant topographic change and values within the threshold cannot be distinguished from noise or possible error in the respective DEMs.

Geomorphological modeling of fluvial geomorphology features
Using the 2012 USGS LiDAR DEM, hydrologic characteristics such as flow direction and flow accumulation were calculated and used to model streamflow and fluvial processes. A relative elevation model and path-distance model were completed using the methods presented by Chirico (2010) for surficial geologic and regolith modeling. It assumes that fluvial erosion and slope processes are the main contemporary natural geomorphic processes at work in the study area. The relative elevation modeling develops a base DEM dataset using the elevation values of points along the hydrologic network. It then subtracts the base elevation data from the 2012 DEM which results in a dataset showing all elevation values of the terrain above the fluvial hydrographic network (erosion surface). The combination of the relative elevation model and the path distance model enable the classification of alluvial and colluvial deposits along the hydrologic network (Chirico, 2010).

Geomorphological modeling of regolith units
The geologic units in the study area have developed regolith profiles as a result of chemical and physical weathering over time. Thus, regolith was mapped by developing a model of regolith thickness based on underlying geologic rock unit properties and supported by analysis of available borehole data (for example see Pavich et al., 1989) and field studies. The relative elevation dataset was used in modeling regolith depth as a function of weathering, groundwater level, and elevation of local stream base level. This analysis revealed that the regolith profiles of the Peters Creek Schist and the Piney Branch Complex can be over 30 m thick and progressively thin from the interfluve downslope to the stream channels eroding them. The intrusive rock units have a thinner (generally less than 3 m thick) weathered bedrock surface.   Figure 6 shows a small area of the 1950 and 1973 SfM DEMs in comparison to the 2012 LiDAR data for reference.
Calculation of a change detection threshold using the RSiQ method indicated that values of the 2012-1950 DoD between −1.46 m and +1.46 m cannot be distinguished from error or noise. Elevation change exceeding −1.46 m or +1.46 m can be considered as areas of substantive change. Substantive geomorphologic change was classed into the following categories.
Areas of negative change are represented by two categories on the map. One showing −1.46 m to −3.5 m (one standard deviation from the mean) and a second showing −3.5 to −10 m (two or more standard deviations) of land removal such for road   through integration of other data sources and detailed aerial image analysis. An overlay of abandoned and buried pre-1950 gravel surfaced roads that were mapped from aerial photo interpretation and verified during fieldwork are shown in light gray. Contemporary paved surfaces are indicated in medium gray, and buildings and structures are indicated in dark gray. Retaining walls present along some transportation infrastructure are mapped using a hatched line.

Geomorphological mapping of fluvial geomorphology features
The results of the fluvial geomorphologic model were used to develop three fluvial geomorphologic map units including alluvial floodplain deposits (in yellow), alluvial levee deposits (in light tan-brown), and colluvial hillslope deposits (in peach). A combination of the relative elevation and path distance models of the DEM data were used to segment floodplains and colluvial deposit zones and then were refined based on field observations along Piney Branch Creek and its tributaries. Alluvial floodplain deposits occupy the valley bottoms where slope values range from 0.0°to 2.5°. Alluvial levee deposits occur in direct proximity to the stream channel and were mapped as higher relative elevation to that of the floodplain deposits and low path-distance values (slope-distance from the stream channel). These deposits were verified in the field through selected hand auguring and field observations from recent overbank flood events. Colluvial units occur at the base of the interfluves on a wider range of slopes (2.5-5.0°) with some local colluvial accumulation on slopes higher than 5.0°. The colluvial deposits occur mainly in the upper portions of the watershed's tributaries but there are also some deposits along Piney Branch's downstream floodplain margins where slope deposits accumulate as a result of surface wash near steep interfluves. Field mapping revealed the presence of legacy sediments in each of the fluvial geomorphologic map units but predominantly in the alluvial and colluvial units. Legacy sediments in alluvial floodplain deposits were observed in streambank profiles along Piney Branch as an approximately 1-2 m layer of silt/clay, fine sand/sand atop a historical soil layer. The bedload of Piney Branch indicates many fragments of building materials (bricks, shingles, etc.), pebble and gravel sized concrete and asphalt nodules, as well as large anthropogenic debris including iron rail ties and spikes from the former railroad (Figures 7 and 8).

Geomorphological mapping of regolith units
Correlation of relative elevation modeling, previous regolith mapping (Pavich et al., 1989), and field observations resulted in the development of 4 map units depicting the presence and thickness of regolith. The first regolith unit (ri) is generally less than 3 m thick and is underlain by units composed of the Lower Cambrian, Cambrian, and Ordovician Piedmont intrusive rock units including amphibolite, quartz bodies, and the Bear Island Granodiorite. These units are more resistant to weathering and thus occur as weathering resistant ridges running perpendicular to Piney Branch and are composed of large boulder to cobble size corestones with abundant quartz fragments mixed with thin sandy to clayey sands at the surface ( Figure 9).
Modeling of the terrain overlying Early Cambrian and/or late Proterozoic and Cambrian Piedmont rocks resulted in three map units (rm1, rm2, rm3) corresponding to the depth of regolith profiles. The regolith profiles are composed of a soil layer which is generally less than 0.6 m thick of sandy-clays to clayey-sands with abundant quartz fragments. The soil is underlain by a mixture of micaceous saprolite  and large corestones of the host rock with abundant quartz fragments throughout. Unit rm1 (in light green) is less than 5 m thick and occurs in the toeslope position. Unit rm2 (in medium green) ranges in thickness from 5-10 m along mid-slope areas, and unit rm3 (in dark green) exceeds 10 m in thickness up to a maximum of 35 m and occupies the highest terrain position on upslope interfluves.

Discussion/Conclusion
This study focused on geomorphological mapping in an urbanizing watershed in the metropolitan Washington, D.C. region of northern Virginia. The study developed a quantitative and repeatable process for mapping anthropogenic and natural landforms integrating several geospatial techniques to map the presence and thickness of anthropogenic cutting and filling activities as well as fluvial geomorphologic and regolith deposits.
In the Piney Branch watershed, colonial-and postcolonial agricultural practices and associated forest clearing resulted in significant erosion and deposition of legacy sediments above alluvial materials in the floodplain, which has affected contemporary fluvial processes. Post-1950s urbanization and development led to significant topographic change in the form of excavation and filling to construct roads, utility infrastructure, basement foundations, and slope grading for housing developments. These changes are highly evident and quantifiable from the image analysis conducted through SfM photogrammetry and topographic change detection. The map presented in this study clearly shows that the spatial distribution of anthropogenic landform change is associated with the most densely urbanized areas and the areas which have undergone the longest periods of urbanization.
Further, this study demonstrated the value of using SfM photogrammetry, historical aerial photography, and topographic change techniques as efficient and repeatable methods for quantifying elevation changes resulting from anthropogenic activities in urban areas. As such, the study contributes towards the goals of advancing the science of urban geomorphological mapping by expanding the technological tools and methods available to scientists, engineers, and planners for the preparation of these mapping products.

Software
Historic aerial images were cropped as necessary using Adobe Photoshop. A digital surface model (DSM) and orthophoto were developed from historic stereo imagery using the structure-from-motion (SfM) algorithms in Agisoft PhotoScan [Metashape v. 1.6]. Fluvial geomorphologic landforms were modeled using hydrologic and raster analysis tools in ESRI ArcMap. The final geomorphologic map was also created using ESRI ArcMap.

Disclosure statement
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.