Estimating Historically Cleared and Forested Land in Massachusetts, USA, Using Airborne LiDAR and Archival Records

: In the northeastern United States, widespread deforestation occurred during the 17–19th centuries as a result of Euro-American agricultural activity. In the late 19th and early 20th centuries, much of this agricultural landscape was reforested as the region experienced industrialization and farmland became abandoned. Many previous studies have addressed these landscape changes, but the primary method for estimating the amount and distribution of cleared and forested land during this time period has been using archival records. This study estimates areas of cleared and forested land using historical land use features extracted from airborne LiDAR data and compares these estimates to those from 19th century archival maps and agricultural census records for several towns in Massachusetts, a state in the northeastern United States. Results expand on previous studies in adjacent areas, and demonstrate that features representative of historical deforestation identiﬁed in LiDAR data can be reliably used as a proxy to estimate the spatial extents and area of cleared and forested land in Massachusetts and elsewhere in the northeastern United States. Results also demonstrate limitations to this methodology which can be mitigated through an understanding of the surﬁcial geology of the region as well as sources of error in archival materials.


Introduction
English colonization of southern New England, a region in the northeastern United States, began in the early 17th century and initiated a drastic change from previous land use activities practiced by Native American groups who had inhabited the region for thousands of years [1][2][3][4]. English-style agriculture was much different than what had previously been practiced in the region, and required widespread deforestation to create a landscape of bounded fields for tillage, pasture, and mowing, punctuated by managed woodlots [3,5,6]. This region had been glaciated up until~16,000-17,000 years ago [7], and agricultural fields were typically developed on glacial till, which was full of large stones and boulders. As forests were cleared, loose stones and boulders on the surface were organized into stone walls [6]. It was also more difficult for soils to retain moisture, and soil temperature fluctuated drastically, causing deeper freezing of the ground surface [3,6]. These temperature fluctuations combined with plowing for agriculture brought stones deeper in the soil column to the surface, where they required removal for successful agricultural pursuits by farmers [3,5,6]. Over decades, stones uncovered during plowing were moved to the edges of fields for disposal and continually added to the stone walls, demarcating property boundaries and individual fields within, and providing fencing. The material also provided a more durable, available, and affordable alternative to timber [6,[8][9][10] and was widely used throughout the region, which is now well-known for its stone walls [6,11]. During the initial stages of English settlement,  [29,30]. This census recorded the amount of improved and unimproved farmland at the town level beginning in 1850. Counties in gray did not exist in 1850 and therefore continuous data are not available. Map projection is NAD 1983 UTM Zone 18N. Abbreviations for states in the inset map are as follows: CT, Connecticut; MA, Massachusetts; ME, Maine; NH, New Hampshire; NY, New York; RI, Rhode Island; VT, Vermont.
The northeastern U.S. is now considered to be one of the most densely forested areas currently in the United States [31], and it is estimated that 65-90% is reforested [32], with only ~450 ha of old growth forest throughout Massachusetts, primarily located in the northwestern part of the state. Many of the tracts are <10 ha and located in steep, rugged terrain [33]. Research is ongoing to determine the impacts of historical widespread deforestation and subsequent reforestation since many aspects are still poorly understood. Examples include impacts to forest ecology [34,35], geomorphology [14,[36][37][38], regional climate [13], and wildlife biology [39]. For example, Hall et al. found that at local scales,  [29,30]. This census recorded the amount of improved and unimproved farmland at the town level beginning in 1850. Counties in gray did not exist in 1850 and therefore continuous data are not available. Map projection is NAD 1983 UTM Zone 18N. Abbreviations for states in the inset map are as follows: CT, Connecticut; MA, Massachusetts; ME, Maine; NH, New Hampshire; NY, New York; RI, Rhode Island; VT, Vermont. The northeastern U.S. is now considered to be one of the most densely forested areas currently in the United States [31], and it is estimated that 65-90% is reforested [32], with only~450 ha of old growth forest throughout Massachusetts, primarily located in the northwestern part of the state. Many of the tracts are <10 ha and located in steep, rugged terrain [33]. Research is ongoing to determine the impacts of historical widespread deforestation and subsequent reforestation since many aspects are still poorly understood. Examples include impacts to forest ecology [34,35], geomorphology [14,[36][37][38], regional climate [13], and wildlife biology [39]. For example, Hall et al. found that at local scales, historical land use imparted strong impacts on local vegetation, but since land use was homogenous throughout Massachusetts, averaging over broad spatial scales reduced variation [34].
Airborne Light Detection and Ranging (LiDAR) has been most commonly used to generate digital elevation models (DEMs) of the land surface and has therefore allowed for the widespread study of historical human impacts on the landscape, as it is effective at discerning extant cultural features in forested areas. Over the past~15 years, as datasets become more widespread and of higher resolution, LiDAR has revealed the extent of human impacts on the land surface in many regions of the world, including New England [40][41][42][43][44][45][46][47]. Extant historical land use features that are indicative of past deforestation, specifically stone walls and relict charcoal hearths (RCHs), are widespread throughout the landscape in the northeastern U.S., and since they have topographic signatures, they are visible using airborne LiDAR [9] ( Figure 2). Using various visualization techniques, features of interest are identified, extracted, and analyzed, and more recent studies have explored automated methods for feature extraction [20,41,[48][49][50][51][52]. Study areas vary on a global scale, as do the time periods which are being studied, especially since features from different time periods exist contemporaneously on the landscape as a palimpsest when viewed in a LiDAR DEM [53][54][55]. Availability of complementary sources that can be used in analyses of these features, such as archival data, also vary, and this can influence the types of analyses that are possible. of these features, such as archival data, also vary, and this can influence the types of analyses that are possible. Reconstructing historical forest cover and cleared land using historical land use features as a proxy has shown promise in several towns in Connecticut through comparison with improved and unimproved farmland recorded in the U.S. Census Nonpopulation Schedule for Agriculture [9]. Both stone walls and RCHs are reliable proxies in reconstructing areas of past land use and deforestation, specifically on smooth terrain that was cleared and used for agriculture, or rugged terrain which often remained forest was used specifically for harvesting timber, or used for charcoaling [9,20,34]. In this study, we demarcate and estimate areas of cleared and forested land based on extant historical land use features visible in LiDAR digital elevation models (DEMs) in several towns in Massachusetts to (1) determine the efficacy of methods published in [9] (which compared data from the 1850 U.S. Census Nonpopulation Schedule for Agriculture with LiDAR-derived estimates) on a broader scale in southern New England and the Reconstructing historical forest cover and cleared land using historical land use features as a proxy has shown promise in several towns in Connecticut through comparison with improved and unimproved farmland recorded in the U.S. Census Nonpopulation Schedule for Agriculture [9]. Both stone walls and RCHs are reliable proxies in reconstructing areas of past land use and deforestation, specifically on smooth terrain that was cleared and used for agriculture, or rugged terrain which often remained forest was used specifically for harvesting timber, or used for charcoaling [9,20,34].
In this study, we demarcate and estimate areas of cleared and forested land based on extant historical land use features visible in LiDAR digital elevation models (DEMs) in several towns in Massachusetts to (1) determine the efficacy of methods published in [9] (which compared data from the 1850 U.S. Census Nonpopulation Schedule for Agriculture with LiDAR-derived estimates) on a broader scale in southern New England and the northeastern United States in general, and (2) assess the utility of other archival sources available in the region in comparison with LiDAR-derived estimates of cleared and forested land. Specifically, we compare areas with a high density of historical land use features indicative of deforestation which have been identified using LiDAR (stone walls and RCHs) with archival map data from 1830 and improved and unimproved land areas from the U.S. Census Nonpopulation Schedule for Agriculture from 1850, 1860, and 1870. This study builds on the novel approach outlined in [9] to reconstruct areas of historically cleared and forested land using extant historical land use features identified in airborne LiDAR data, and is applicable to other regions of the world where similar historic land use features and datasets exist. While modern land cover mapping has been undertaken in this region using airborne LiDAR [56], the methods presented here provide a unique approach to estimating and mapping areas of historical land cover. High-resolution regional historical forest cover estimates can provide an important complementary dataset for use in ongoing research to understand the impacts associated with the widespread historical deforestation that this region experienced in the 17th through 19th centuries and can be used in climate models, regional mapping programs, and can assist in other uses of modeling or quantifying impacts of historical deforestation.

Study Area
The study area encompasses 18 modern towns in Massachusetts ( Figure 3, Table 1). To ensure consistency between areal estimates from the 19thcentury and present, we used an 1830 town boundary shapefile provided by Harvard Forest with their 1830 forest cover dataset [57] (see Section 2.3 for details) with which to compare modern town boundaries. For boundaries that did not match, we reviewed Massachusetts General Court records to determine which years the boundaries changed for that town [58]. We determined the state of all town boundaries in the year 1850, which is when the first agricultural census was conducted for each town, and then used a custom boundary shapefile for clipping LiDAR-derived data as well as forest data from the 1830 Harvard Forest shapefile (see Section 2.3 for details). For more information about specific town boundary changes, see Appendix A. Of note is that the modern towns of Adams and North Adams were one town (Adams) up until 1878, and for the purposes of this study, the area of both modern towns is combined in tables and calculations to represent the historical boundaries.
state of all town boundaries in the year 1850, which is when the first agricultural census was conducted for each town, and then used a custom boundary shapefile for clipping LiDAR-derived data as well as forest data from the 1830 Harvard Forest shapefile (see Section 2.3 for details). For more information about specific town boundary changes, see Appendix A. Of note is that the modern towns of Adams and North Adams were one town (Adams) up until 1878, and for the purposes of this study, the area of both modern towns is combined in tables and calculations to represent the historical boundaries.  [25,26]. Map projection is NAD 1983 UTM Zone 18N.   [25,26]. Map projection is NAD 1983 UTM Zone 18N.  4 Area in analysis reflects the modern towns of Adams and North Adams combined. 5 [60], see also Appendix C.
The majority of the towns are located in the central to western part of the state, where the bedrock geology is metamorphic schist and gneiss. Surficial deposits in these towns are primarily composed of lodgement and ablation till overlying bedrock (~75-100% of the town area; Table 1, see also Appendix C), with some fine to coarse glacial stratified (sand and gravel) deposits [61]. The exception is the town of Rehoboth, located in the southeastern corner of the state, where the bedrock geology is sedimentary rock, and is only covered by~50% glacial till and~40% glacial stratified deposits. Glacial till is the primary source of building material for stone walls [6,62]. The till varies in thickness from a few centimeters up to 60 m, and contains unsorted rock and mineral particles which are variable in distribution and range in size from clay to boulders [61]. Glacial stratified deposits primarily consist of clay to gravel sized sediment deposited by glacial meltwater [61]. Floodplain alluvium, a post-glacial deposit which generally overlies glacial stratified deposits, is located in river valleys and generally comprises less than 5-10% of the surficial deposits [61]. Some towns contain a higher percentage of post-glacial floodplain alluvium due to the presence of larger hydrological features. This includes the Hoosic River, which bisects Adams and North Adams and flows northward through Williamstown, becoming a tributary to the Hudson River in New York. The Williams River flows through West Stockbridge, eventually emptying into the Housatonic River. Rehoboth also contains a high percent of floodplain alluvium and swamp deposits due to its location within the Narragansett Bay watershed.

LiDAR Datasets, Feature Digitization, and Geospatial Analysis
While Massachusetts has full spatial coverage with LiDAR datasets currently, the datasets were collected separately over the course of several years (2002-2015) and therefore have slightly different characteristics [63]. All datasets were collected during leaf-off season in this region to obtain the most accurate measurements of the topography. The ability to accurately measure the ground surface is influenced not only by the presence of vegetation, but also by vegetation type [53,65]. The presence of coniferous vegetation, which retains its foliage year-round, can influence the number of LiDAR pulses that reach the actual ground surface in these areas [53]. Additionally, in forests with a dense understory, it is also challenging to discern the actual ground surface from low vegetation during point classification. In areas where pulses may not reach the ground surface as effectively, there may be fewer points to interpolate amongst during final DEM generation, making the morphology of fine-scale features in these areas somewhat less resolved [53,65]. The DEMs used in this study were generated by MassGIS from LAS 1.2 and 1.4 ground-classified points and have a 1 m pixel resolution [63]. All relevant tiles for each town were downloaded as DEMs from MassGIS OLIVER [64] and mosaicked using the Mosaic to New Raster tool in ArcMap.
Digitization of historical land use features from the DEMs was similar to the methods used and published in [9, 44,62], in which features were manually identified using hillshade and slope rasters derived from the DEM in ArcMap using the Hillshade and Slope tools in the Spatial Analyst toolbox. Stone walls were digitized as polylines with vertices placed at intersections with other walls, changes in wall direction, and wall endpoints. A single point was placed in the center of each identified RCH. The stone wall and RCH datasets analyzed and included in this contribution were digitized by different users over the span of several years and then edited to maintain quality. Digitized datasets may have errors associated with LiDAR dataset quality (see [66] and above discussion) or user interpretation [67]. Experienced and novice mappers may consistently over-or under-map depending on their preferences and confidence in identifying features in a LiDAR DEM [67]. For stone walls in deciduous forests and relatively simple terrain, individual digitizers agreed with the group 88-92% of the time, but frequently missed walls, and on average, mapped only 78% of walls that were verified in the field [67]. These discrepancies increase in areas of coniferous forest, rugged terrain, or if topography has been modified by roads, buildings, or other infrastructure.
To derive areas representative of cleared and forested land, we used methods consistent with those published in [9] to generate rasters with the density of historical land use features, and then reclassify areas above specific density thresholds which are representative of intensive land use areas. Area was our preferred measurement in order to compare values with those recorded in the U.S. Census Nonpopulation Schedule for Agriculture, which recorded improved and unimproved farmland area. We first generated rasters depicting the density per km 2 of the historical features, and then clipped those rasters to the town boundary. All rasters were generated using tools available from the ESRI ArcGIS suite using the arcpy library. Density rasters for stone walls were generated using the Line Statistics tool, and for RCHs using the Point Density tool, each with a circular neighborhood of 564.19 m for a resulting search area of 1 km 2 . The output rasters depicted the length of stone walls in km/km 2 and the number of RCHs per km 2 . Consistent with the intensive land area threshold determination methods described in [9], we reclassified the stone wall density raster using the Reclassify tool so that areas where the length of stone walls exceeded 2 and 3 km/km 2 were classified as improved (e.g., 2000.01-maximum = class 1) and the inverse areas were classified as unimproved (e.g., 0-2000 = class 0). For RCHs, we performed the same function for areas where RCH density exceeded 10 and 15 per km 2 in order to compare our results to the previous study, which also used those thresholds. By reclassifying the rasters using this method, we captured areas of each town where land would have been used much more intensively for these specific land use types.
To determine the percentage of glacial till and other deposits for each town, a shapefile of Massachusetts surficial geology at the 1:250,000 scale [60] was used, and this layer was clipped to the extent of each town in the custom boundary shapefile and summarized by category of interest to determine the percent of town area for each deposit type. Higher resolution surficial maps exist for the state, however, at the town scale, the 1:250,000 resolution was sufficient for estimating the percent of the town surface area contributing to construction material for stone walls.

Archival Datasets
Two archival datasets were used in this study. The first was geospatial data of all woodland in the state in 1830, which was derived from a series of maps mandated of all towns by the Massachusetts legislature in 1830 [34,68,69]. The outcome of the mandate and subsequent survey was a map of almost every town in Massachusetts during the time period, variably showing forest, fields, buildings, roads, and other features of interest. While some maps were very detailed, others were not. Harvard Forest obtained these maps from a variety of archival locations and meticulously studied and digitized them. They subsequently made the digitized geospatial products available in 2009 through the Harvard Forest Digital Archive as downloadable geospatial data [57]. Scanned copies of most of these archival maps can now be accessed online [69]. The Harvard Forest dataset provides shapefiles of woodland, buildings, town boundaries, and other features derived from these maps for almost the entirety of Massachusetts in 1830 ( Figure 4). The dataset description notes several limitations to the maps and thus geospatial data, including the fact that maps were generated by different individuals for each town and do not all conform to the same mapping standards. For example, in some towns, there are large generalized blocks of woodland, while in others, the areas are much more detailed, which could impact woodland estimates [57]. We do expect there to be some error associated with the calculated estimates we provide based on these maps (discussed in depth with results in Section 3); despite this, they do represent a valuable early resource in this region that is worth investigating and comparing to other potential historical land use estimates for this time period. form to the same mapping standards. For example, in some towns, there are large generalized blocks of woodland, while in others, the areas are much more detailed, which could impact woodland estimates [57]. We do expect there to be some error associated with the calculated estimates we provide based on these maps (discussed in depth with results in Section 3); despite this, they do represent a valuable early resource in this region that is worth investigating and comparing to other potential historical land use estimates for this time period. To calculate forested land from this dataset, the 1830lulc.shp file downloaded from the Harvard Forest Data Archive [57] was first edited to remove any polygons associated with towns denoted as missing woodland data. Next, we selected from the resulting layer any instances where the polygon attribute "WOODLAND" was 1, indicating that it had been classified as woodland. This selection was intersected with the custom town boundary shapefile described in Section 2.1 in order to calculate the amount of woodland for each town in the study using the Summary Statistics function. Using the Symmetrical Difference tool on the resulting layer allowed us to calculate areas of non-forested land in each study town, from which water that had been recorded in 1830 was removed using To calculate forested land from this dataset, the 1830lulc.shp file downloaded from the Harvard Forest Data Archive [57] was first edited to remove any polygons associated with towns denoted as missing woodland data. Next, we selected from the resulting layer any instances where the polygon attribute "WOODLAND" was 1, indicating that it had been classified as woodland. This selection was intersected with the custom town boundary shapefile described in Section 2.1 in order to calculate the amount of woodland for each town in the study using the Summary Statistics function. Using the Symmetrical Difference tool on the resulting layer allowed us to calculate areas of non-forested land in each study town, from which water that had been recorded in 1830 was removed using the Erase tool for polygons with a "LANDCOVER" class of 7 in the original 1830lulc.shp file.
The second archival dataset used was the U.S. Census Nonpopulation Schedule for Agriculture [30,70] (referred to throughout this text as 'agricultural census') for the years 1850, 1860, and 1870. The 1850 and 1860 agricultural censuses categorized farmland as "Improved" and "Unimproved" land in acres for each farm in a town, whereas the 1870 agricultural census made a distinction between "Woodland" and "Other" as subsets within the "Unimproved" category. The instructions for the agricultural census defined unimproved land as "...a wood lot or other land at some distance but owned in connection with the farm, the timber or range of which is used for farm purposes...", while improved land was defined as "...cleared and used for grazing, grass, or tillage, or which is now fallow" [71] (p. 235). We transcribed the values from the scanned ledgers of the 1850, 1860, and 1870 agricultural census available on Ancestry.com into tabular data, summed areas of "Improved" and "Unimproved" land for all recorded farms in the study towns, and then converted the acreages to km 2 (1 acre = 4046.86 m 2 ) [30].
We can expect small sources of error at the scale of this study with the use of historical census data for a variety of reasons, though studies have generally shown it to be reliable [9, [72][73][74]. The most significant expected source of error is that often the sum of "Improved" and "Unimproved" areas is not equal to the total area of the town. The average percentage of town area covered by each agricultural census is as follows: 1850, 76%; 1860, 70%; 1870, 66%. The percentage of area comprised by both "Improved" and "Unimproved" land amongst towns varies drastically; in the 1850 census, percentages for the towns in this study range from 34% to 109%, for example. These variations could be a result of many factors, including that some farms were intentionally omitted by census-takers. The 1850 schedule itself notes that small farms with production value < USD 100 should not be included [71] (p. 235), a number which changed to USD 500 by 1870 [70]. A comparison of the population census for Clarksburg in 1850 against the agricultural census shows 42 male heads of household listed with the occupation of "Farmer," yet only 32 of these men have farms listed in the agricultural census (~76%), and some men with other occupations (e.g., "Laborer", "Manufacturer") have farms noted in the agricultural census [30,75]. In addition to missing smaller farms, there were likely areas which were waterbodies, not farmland, and it is possible that some areas may have been left unenumerated or not visited. For example, the instructions for the agricultural census in 1870 specifically note "...Irreclaimable marshes and considerable bodies of water will be excluded in giving the area of a farm improved and unimproved" [71] (p. 746).

Results and Discussion
Overall, the density of stone walls varies throughout all of the study towns, with maximum densities within one individual town ranging from~5 km/km 2 in New Ashford to~13 km/km 2 in Royalston. All towns have minimum densities of 0 km/km 2 , as there are locations in all towns where no stone walls were digitized. Relict charcoal hearth (RCH) maximum densities range from 16 RCHs/km 2 in Adams to 110 RCHs/km 2 in Williamstown ( Figure 5, Table 2). For comparison, towns in Connecticut that were analyzed and published in [9] had as many as 197 RCHs/km 2 in some areas and stone wall length per km 2 exceeding 11 km/km 2 . Therefore, it appears that while stone wall density is similar to that observed in Massachusetts, the area studied in Connecticut was used much more extensively for charcoal production.   By comparing the values of improved (i.e., non-forested) land area derived from the 1830 archival maps as well as 1850, 1860, and 1870 agricultural census records (Appendix B) with area values derived from LiDAR data (Table 3), we generally found that the areas in the 1830 maps have the highest correlation with estimates derived from the distribution of historical features in LiDAR data (R 2 = 0.71), and of the agricultural census data values, those in 1850 have the highest correlation with LiDAR-derived area estimates (R 2 = 0.57) (Figure 6).
Combining the data from this study with five towns in Connecticut (see Figure 1) that were examined in a previous study [9] reveals a generally good correspondence (R 2 = 0.70) between tabulated improved land in the 1850 agricultural census and areas of intensive land clearing derived from LiDAR data using a threshold where stone wall density exceeds 2 km/km 2 (Figure 7). The trend (dotted line) also generally falls above the solid gray line that indicates a 1:1 agreement of values, suggesting that using this threshold slightly overestimates the amount of cleared land in a town in comparison to the census. In contrast with the results presented in Figure 8 of Reference [9], the correlation between 1850 agricultural data and areas derived from stone walls is not quite as strong in the subset of Massachusetts towns (R 2 = 0.57, Figure 6B) as for the Connecticut towns (R 2 = 0.97) [9]. There are several reasons this could be, including underlying surficial geology of the towns, and the timing of the historical peak of deforestation. We expect there will be more variation as more towns are included in the dataset, especially if towns are included which have lower proportions of glacial till, as there could be fewer stone walls there and higher proportions of wood historically used for fencing (see [6,8,76]).     Combining the data from this study with five towns in Connecticut (see Figure 1) that were examined in a previous study [9] reveals a generally good correspondence (R 2 = 0.70) between tabulated improved land in the 1850 agricultural census and areas of intensive land clearing derived from LiDAR data using a threshold where stone wall density exceeds 2 km/km 2 (Figure 7). The trend (dotted line) also generally falls above the solid gray line that indicates a 1:1 agreement of values, suggesting that using this threshold slightly overestimates the amount of cleared land in a town in comparison to the census. In contrast with the results presented in Figure 8 of Reference [9], the correlation between 1850 agricultural data and areas derived from stone walls is not quite as strong in the subset of Massachusetts towns (R 2 = 0.57, Figure 6B) as for the Connecticut towns (R 2 = 0.97) [9]. There are several reasons this could be, including underlying surficial geology of the towns, and the timing of the historical peak of deforestation. We expect there will be more variation as more towns are included in the dataset, especially if towns are included which have lower proportions of glacial till, as there could be fewer stone walls there and higher proportions of wood historically used for fencing (see [6,8,76]). This method of determining cleared or forested land using stone walls as a proxy appears to work more effectively in towns comprised almost solely of glacial till. This would have provided farmers a means to construct walls throughout the entire town without areas of constraint depending on the surficial materials. For example, within the plots of Figure 6, we observed several towns that fall well below a 1:1 agreement, indicating that our methods underestimate the amount of cleared land that was recorded in these towns. This may be influenced by the presence of riverine or glacial lake deposits, which would result in fewer stone walls than expected in these towns. The 1850 agricultural census values appear to be a slightly better fit with estimates from LiDAR data, as our method This method of determining cleared or forested land using stone walls as a proxy appears to work more effectively in towns comprised almost solely of glacial till. This would have provided farmers a means to construct walls throughout the entire town without areas of constraint depending on the surficial materials. For example, within the plots of Figure 6, we observed several towns that fall well below a 1:1 agreement, indicating that our methods underestimate the amount of cleared land that was recorded in these towns. This may be influenced by the presence of riverine or glacial lake deposits, which would result in fewer stone walls than expected in these towns. The 1850 agricultural census values appear to be a slightly better fit with estimates from LiDAR data, as our method seems to underestimate the amount of non-forested land in almost every town in 1830 ( Figure 6A). In the 1850 agricultural census plot ( Figure 6B), towns that fall below the 1:1 line, such as Adams, Conway, Leverett, Williamstown, and West Stockbridge, are comprised of a much higher proportion of floodplain alluvium, sand, and gravel than other towns in the study area (see Appendix C). All of these towns therefore likely contain large areas which may have been deforested or fenced, but do not have any stone walls demarcating those events given the lack of stone fencing material. It is likely that other types of fencing were used in these towns to demarcate cleared land, notably in the general area of the former glacial Lake Hitchcock and current Hoosic River watersheds. This is not unexpected, as it has been shown that surficial materials, notably glacial till, influence the presence, type, and morphology of stone walls [6,62]. It appears, therefore, that in towns where there is less glacial till and more sand or glacial lake deposits, the method does not appear to work as well, and this may be why there is such variation in Massachusetts compared with the five towns that were studied in Connecticut. In Connecticut, the towns that were studied also had extensive spatial coverage with both stone walls and RCHs, which allowed for a nearly complete sampling of space in each town. This likely contributed to the accuracy of the LiDAR-derived estimates. In Massachusetts, there are large areas of each town that are not covered by stone walls or RCHs, which leaves these areas effectively unsampled and unknown.
We hypothesize that another reason the correlation between the map or agricultural census data and the LiDAR-derived area estimates is not as strong for the Massachusetts towns could be because the height of land clearance in Massachusetts was not in 1830, or even 1850 as it was in Connecticut, but later in the 19th century, as shown in other studies [12,77,78]. For example, Foster [77] provides forest data for Petersham with a height of clearance in 1860, and for Conway/Shelburne [79], where the height of clearance was ca. 1840-1860. Despite this, comparison of LiDAR-derived areas with totals from agricultural census records from 1860 and 1870 have an even weaker correlation ( Figure 6C,D) with our estimates, suggesting that 1850 continues to provide the likely peak of deforestation for this particular group of towns on average. While Monterey and Sandisfield are not included in the 1870 plot, and those towns plus Plainfield are not included in the 1860 plot, removal of all three towns from the 1850 plot, which has the best correlation, shows it decreases to R 2 = 0.50, suggesting that this year still has the best outcome of the agricultural census values. With regard to timing and spatial distribution, it appears that the extant stone walls on the landscape are representative of all land that had ever been cleared in the town, assuming stone fencing had been used. We observe this when comparing stone wall area estimates to agricultural census estimates from time periods known to have been the height of agricultural clearance, and this is also likely why our method underestimates the 1830 values. Comparing LiDAR-derived estimates to census-derived estimates for 1850, 1860, and 1870 reveals this quite clearly.
We also found that in Massachusetts towns such as West Stockbridge, where charcoaling seems to have occurred on a fairly widespread scale, the location of RCHs corresponds well to areas demarcated as being forest during 1830, and we can infer that those areas would have also been deforested at some point in time ( Figure 8). As demonstrated in [9], rugged terrain unfit for agriculture was typically used for charcoaling during the 19th century, and these areas also correspond well to areas marked "Unimproved" in agricultural census records. However, charcoaling did not occur as frequently in some of the towns as in others, likely due to a variety of factors, including proximity to iron furnaces or ease of transport.
While the 1830 maps provide an unprecedented view of forest cover in Massachusetts during this time period, there is also error associated with them. Harvard Forest makes note of several data usage caveats in their dataset description, including that the maps were drafted by different individuals and there was no standardized cartographic approach. Additionally, it has now been almost twenty years since these maps were digitized by Harvard Forest, and it is possible that updates to GIS software could provide slightly different results. As an example of error in this dataset, see Figure 9, which depicts the eastern extent of Clarksburg, where the 19th century surveyor likely meant to include all of the mountainous area in their forested extent, but when the map is georeferenced and overlaid on a LiDAR DEM ( Figure 9B), the forested extent only accounts for a small portion of the actual mountainous terrain, and does in fact match fairly closely with the data from Harvard Forest's digitized shapefile ( Figure 9C). It is possible that a comprehensive comparison of the forest dataset with high-resolution LiDAR DEMs could potentially provide more accurate estimates. While the 1830 maps provide an unprecedented view of forest cover in Massachusetts during this time period, there is also error associated with them. Harvard Forest makes note of several data usage caveats in their dataset description, including that the maps were drafted by different individuals and there was no standardized cartographic approach. Additionally, it has now been almost twenty years since these maps were digitized by Harvard Forest, and it is possible that updates to GIS software could provide slightly different results. As an example of error in this dataset, see Figure 9, which depicts the eastern extent of Clarksburg, where the 19th century surveyor likely meant to include all of the mountainous area in their forested extent, but when the map is georeferenced and overlaid on a LiDAR DEM ( Figure 9B), the forested extent only accounts for a small portion of the actual mountainous terrain, and does in fact match fairly closely with the data from Harvard Forest's digitized shapefile ( Figure 9C). It is possible that a comprehensive comparison of the forest dataset with high-resolution LiDAR DEMs could potentially provide more accurate estimates.

Conclusions
In general, this work demonstrated that features representative of historical deforestation identified in LiDAR data can be reliably used as a proxy to estimate areas and the spatial extents of cleared and forested land in Massachusetts and elsewhere in the northeastern United States. As feature digitization continues, expanding the ability to estimate areas of deforested and forested land will broaden and become more feasible. We found that there is spatial variation amongst the towns regionally in terms of the density of extant relict land use features (i.e., stone walls and RCHs), and thus the overall ability to reconstruct areas of possible improved (i.e., deforested or cleared) land. In areas where stone might have been less available for wall building and where fencing might have been comprised partially or entirely of wood, it is also much more challenging to reconstruct areas that might have been improved, and vice versa. The same is true if stone walls have been removed or are otherwise no longer preserved or visible on the land surface. While there are some limitations with archival data, when combined with airborne LiDAR, these resources help confirm that extant historical land use features identified using LiDAR data can be used to reliably reconstruct deforested areas on regional scales.
Overall, the impacts of historical deforestation in this region were significant, and this is evident from the widespread relict land use features visible across the landscape. Continued research to understand and quantify these impacts in the northeastern U.S. is important in understanding past, present, and future aspects of forest ecology, climate, wildlife biology, geomorphology, and a range of other considerations. Recent efforts to automate digitization of these historical land use features will allow for expansion of these methods across regional scales and provide high-resolution historical land cover estimates to assist in understanding these impacts.

Conclusions
In general, this work demonstrated that features representative of historical deforestation identified in LiDAR data can be reliably used as a proxy to estimate areas and the spatial extents of cleared and forested land in Massachusetts and elsewhere in the northeastern United States. As feature digitization continues, expanding the ability to estimate areas of deforested and forested land will broaden and become more feasible. We found that there is spatial variation amongst the towns regionally in terms of the density of extant relict land use features (i.e., stone walls and RCHs), and thus the overall ability to reconstruct areas of possible improved (i.e., deforested or cleared) land. In areas where stone might have been less available for wall building and where fencing might have been comprised partially or entirely of wood, it is also much more challenging to reconstruct areas that might have been improved, and vice versa. The same is true if stone walls have been removed or are otherwise no longer preserved or visible on the land surface. While there are some limitations with archival data, when combined with airborne LiDAR, these resources help confirm that extant historical land use features identified using LiDAR data can be used to reliably reconstruct deforested areas on regional scales.
Overall, the impacts of historical deforestation in this region were significant, and this is evident from the widespread relict land use features visible across the landscape. Continued research to understand and quantify these impacts in the northeastern U.S. is important in understanding past, present, and future aspects of forest ecology, climate, wildlife biology, geomorphology, and a range of other considerations. Recent efforts to automate digitization of these historical land use features will allow for expansion of these methods across regional scales and provide high-resolution historical land cover estimates to assist in understanding these impacts.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Town Boundary Details
Adams and North Adams are today separate towns, but were the same town up until 1878; for the purposes of this study, the two modern towns were recombined and are referred to as "Adams" throughout the text. Additionally, in 1900, the boundary between North Adams and Williamstown was slightly modified to its present state [81]. Clarksburg lost its easternmost portion to the town of Florida in 1913 [82]. A small portion (~2 km 2 ) of Conway's NW corner was annexed to Buckland, an adjacent town, in 1838 [83]. Monterey was part of the town of Tyringham until 1847 when it split off, and subsequently gained more land from adjacent New Marlborough in 1851 [84] and Sandisfield in 1875 [85]. Sandisfield itself lost a different tract of land to nearby Tolland in 1853 [86]. Petersham's southern boundary was expanded in the 1930s as a result of the construction of the Quabbin Reservoir. The reservoir covers, either partially or fully, four historical towns, which were disincorporated and flooded for the reservoir's construction, and associated land was ceded to adjacent extant towns. Two of the disincorporated towns, Dana and Greenwich, are mostly encompassed by Petersham's modern boundaries today. Most of Petersham's original boundaries were not subject to flooding for the reservoir's construction, with the exception of a small piece in the west, which appears to have been a low marshy forested area in 1830 according to the maps used in this study. Petersham's land use history has been thoroughly studied and meticulously documented by researchers at Harvard Forest, which is located there. For more about impacts of historical land use on forest composition, ecology, and other factors, see the following for some examples: [80,87]. West Stockbridge's boundary was changed slightly in 1847 when part of the town was annexed by an adjacent town (Alford) [88].