Future (2020–2099) Carbon and Water Dynamics of Lehigh Valley Based on Land Use and Land Cover Change

: Increased urbanization has reduced the amount of green space, resulting in a reduced carbon sink potential across urban landscapes. Through the use of biogeochemical modeling, different land use scenarios have been developed and run for the future (2020–2099) to compare and quantify the potential for change in carbon and water dynamics by having more tree cover and reducing impervious surfaces or turf lawns in Lehigh Valley, PA. These results show that the effect of deforestation is larger than the effect of reforestation. Due to young-stand age trees having a lower capacity for carbon storage than mature trees, the loss of the mature trees has a more immediate impact. The conversion of lawns or impervious surfaces to forests has somewhat similar effects, although the higher nutrients of lawns allow the forest to grow better. However, replacing impervious surfaces with trees reduces runoff more. This study shows that within the city of Bethlehem, the most socially vulnerable area benefits the most from increasing the number of trees. When converting 25% of the impervious area to forest, South Bethlehem significantly increased its vegetation carbon, productivity, and carbon storage, reduced its runoff, and generally created a safer and cleaner environment for residents.


Introduction
As development becomes more extensive, planning land use and land cover (LULC) of urban areas is increasingly important, in part because urbanization has been associated with reduced carbon sequestration [1,2].On a global scale, the area available for land use and land cover change (LULCC) within urban areas is relatively small, so it has a minor effect on greenhouse gas (GHG) dynamics.However, in areas in which the LULCC exceeds 10%, there has been a noticeable change in carbon emissions [1].Realistically, there is a limited amount of land within cities that can be converted to other land use types.In this study, biogeochemical modeling is used to evaluate the effect of LULCC on carbon and water dynamics in Lehigh Valley, Pennsylvania, USA, relative to climate policy initiatives.Lehigh Valley consists of Lehigh and Northampton counties, and this study also considers a focus on the city of Bethlehem.The city of Bethlehem has developed a climate action plan [3] that outlines the specific LULCC the city plans to implement to reach zero carbon emissions by 2040.To test the effectiveness of these policies, sensitivity experiments were developed for 2020-2099 to determine the type and magnitude of urban land use change needed to meet these goals.
The terrestrial ecosystem model (TEM) is a biogeochemical model of terrestrial ecosystems that tracks the flow of carbon, nitrogen, and water among the atmosphere, vegetation, and soils.The model has been used to study the effect of human disturbances like LULCC, climate change, elevated CO 2 , ozone, and nitrogen deposition on carbon sequestration and ecosystem services [4][5][6][7][8].By inputting LULC and climate data for Lehigh Valley into the model, information about carbon fluxes, as well as moisture dynamics, can be determined.Changing the number of land use types (impervious, barren, deciduous and coniferous forests, lawn, cropland, and pastureland) and running these different regimes in TEM allows changes in water and carbon storage associated with different potential land use policies to be estimated.
When considering carbon dynamics and LULC, it is well documented that urban areas have low primary productivity and thus carbon storage [9].Ecosystems such as fertilized cropland, forests, and grazing pastures are among the highest carbon sinks [10].Urban turf lawns are fertilized and therefore can also be strong carbon sinks [11].Growing new shrubs and trees will be a net carbon sink, and over time, the sink may continue due to future disturbance and regrowth or growth enhancements from nitrogen deposition and elevated atmospheric CO 2 [12,13].However, in accounting for the entire carbon history, mowing and decomposition more than offset enhanced growth from fertilizers [10].One measure of carbon sequestration by the land that accounts for the effects of disturbance and land management is net carbon exchange (NCE), which is net ecosystem productivity (NEP) minus the fluxes associated with land conversion and product decomposition [10].A more positive NCE is indicative of a stronger carbon sink.
When considering the impact of urbanization on carbon dynamics, soils are often not considered.In areas that have a low soil organic carbon (SOC), such as the warm and dry Midwest and the western United States, urbanization results in a gain in SOC [14].For areas that have a higher SOC, such as the colder northeastern United States, urbanization decreases SOC [14].Of the cities listed in this study, Bethlehem is most similar to Baltimore, which saw a 2.2% decrease in SOC.Increasing urbanization and the replacement of forests or turf lawns with impervious surfaces also increase the amount of water runoff within these areas.
Carbon sequestration is currently not accounted for in the city's GHG inventory.However, a specific goal of the Bethlehem CAP [3] is to increase the amount of green space in the city, particularly for frontline communities (people of color, low-income residents, and those with existing health conditions), which are communities most at risk of the ramifications of climate change.These communities are meant to receive 40% of the benefits of greenhouse gas reductions.Lower socioeconomic members of the community are often disproportionately affected by the ramifications of environmental consequences, such as the urban heat island effect [15].The distribution of green spaces throughout urban areas is often not equal across the landscape and has long been associated with racist practices and redlining [16].The lack of green spaces is not only an issue of inequity but in terms of environmental function, there is a heterogeneity of ecosystem services across the landscape [16].The unequal distribution of green spaces and ecosystem function disproportionately negatively affects lower-income areas.Comparing carbon storage to social vulnerability allows for the distribution of ecosystem services across urban areas to be assessed.One outcome of this research is to provide more information on how carbon and water dynamics are affected by development decisions, how carbon storage is distributed within a city, and how it is associated with social vulnerability.
In this study, we used TEM to study the following hypotheses: (1) a 10% increase in tree cover will not significantly change carbon and water dynamics, while a 25% or more increase will significantly change them, (2) decreases in the amount of forest will have a larger impact on carbon storage than increasing the amount of forest by the same percentage due to the increased forest having a younger stand age, and (3) within the city of Bethlehem, the more socially vulnerable areas will see the largest increases in carbon storage from increases in the amount of forest due to having the highest proportion of impervious surfaces.Since previous studies have shown that at least a 10% change in LULC of the total region's area is required to achieve a significant change, it is expected that only LULC changes accounting for more than 10% of the total area of the region will produce significant changes.When comparing changes of the same magnitude (increasing or decreasing forest by the same percentage), decreasing the amount of forest will likely have a larger impact than regrowing new forests since newly forested areas are less productive initially.Within the city of Bethlehem, it is expected that the areas that are more socially vulnerable will see the greatest increase in carbon storage when having impervious areas converted to forests.Neighborhoods with higher social vulnerability tend to have a higher proportion of impervious surfaces to tree cover, which makes them of particular importance when designing LULCC scenarios.

Methods
TEM requires input variables of climate data and LULC.Remote sensing data were used to develop LULC datasets in ArcGIS Pro Version 3.0.For the climate datasets, the multivariate adaptive constructed analogs (MACA) RCP8.5 data from the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM) were used for future projections from 2019 to 2099.Sensitivity experiments were run using TEM at monthly timesteps to determine the carbon and water dynamics of a gridded representation of Lehigh Valley at the 4 km resolution.A subset of LULC focused on the city of Bethlehem was also run to compare social vulnerability and an urban heat island to the estimated changes in the carbon sink resulting from LULCC.

Remote Sensing
To run TEM, LULC datasets are required to provide the amount of each cohort (described in Section 2.2) in each grid cell.To create the dataset, aerial imagery was classified to separate the pixels into different cohorts and then overlain with 4 km grids.The aerial imagery used for this project was provided by USGS Earth Explorer (Table S1) [17].A Landsat 8 Operational Land Imager (OLI) C1 L1TP imagery (Collection 1 Level 1 Terrain Precision Correction) taken on 30 July 2017 was used to perform a classification to determine the fractional amount of each cohort within each grid in Lehigh Valley, as well as the city of Bethlehem.An additional classification was conducted on an image from 21 December 2017 to determine the proportion of coniferous and deciduous trees within Lehigh Valley (Table S1, Figure S1).L1TP data were radiometrically calibrated and orthorectified through the use of ground control points and digital elevation models (USGS).The image is in tier 1 (T1), which is the highest available quality and processing level for Landsat data [17].
For both of the images, the remote sensing data were used to categorize and quantify the amount of land within Lehigh Valley that is developed (impervious), barren, forest, turf lawn, cropland, pastureland, and water.For the image from December, the forest cohort was replaced by coniferous and deciduous cohorts, using a forest area that was bare of leaves to classify the deciduous areas.In both images, water was removed from the LULC dataset because it is not a cohort used within TEM.Training areas were created to match pixels of similar reflectance into the land use type groups.The training areas were used in the pixel-based, random tree classification to group the pixels.Multiple iterations of the classification were performed to improve accuracy by visually checking that the pixels were being classified correctly.Using a fishnet tool, 4 km grids were created to determine the amount of each land use type in each of the grids.The classified dataset was overlain with the fishnet, outlining Lehigh Valley, and the outline of the city of Bethlehem was used for the experiments involving the subset of Bethlehem (Figure S1).
As a result of the winter image, the forest cohort was split into 15% coniferous and 85% deciduous forests.The forest cohort division was determined by comparing an image of Lehigh Valley in the winter (21 December 2017) to an image from the summer (30 July 2017).The total area of the forest in the winter image was compared to the total area of the forest in the summer image to ensure accuracy.To verify the ratio of coniferous to deciduous forest, additional classifications were conducted and produced similar results.The ratio of forest in the winter to the summer was then used to determine what percentage of the total forest is coniferous and deciduous, as the winter image would only show the area represented by coniferous trees.Seventy percent of the lawn cohort was split into 65% cropland and 35% pastureland in grids deemed non-urban (<35% impervious surfaces) as per the Lehigh Valley Planning Commission [18].The remaining 30% was left as lawn.The classified data were split into 4 km grid cells.After gathering data for the categories, the data were input into TEM and run under current climate conditions to serve as a control run.The most dominant land use types were forest area in most cohorts, urban in the areas with the most impervious surfaces, and cropland in more rural areas of Lehigh Valley (Figures S1 and S2).
Using the same classified image from Lehigh Valley, four grids in the shape of the city of Bethlehem were extracted to represent the north, east, south, and west sections of the city (Figure S1).A forested area in the southern section of the city was removed because it is part of the Lehigh University campus, no one lives in the area, and our primary goal within the city was to focus on how the potential LULCC could impact socially vulnerable populations.

Model Description
The terrestrial ecosystems model version hydro (TEM-Hydro [4][5][6] consists of multiple pools for vegetation carbon and nitrogen and a single pool for soil organic carbon and nitrogen, as well as a pool for inorganic nitrogen.The vegetation pools include leaves, active and inactive stem tissues (e.g., sapwood and heartwood), fine roots, and a labile pool for storage.Fluxes into and out of the pools include photosynthesis, nitrogen uptake, respiration, litterfall, allocation, and dissolved organic carbon (DOC).The more recent version of the model (TEM-Hydro2) includes improvements to the modeling of land use and a reduced-form open nitrogen cycle [4], which allows for N deposition, N fixation, denitrification, and leaching of dissolved organic and inorganic nitrogen (DON and DIN).Felzer et al. [4,5] provide a complete description of the model, illustrated in a summary figure (Figure S3), along with how human disturbance is treated (Figure 1 from Felzer [4]), which is relevant to this paper.
The cohort approach is designed to keep track of land use transitions so that when converting impervious surfaces or lawns to forests, for example, the new forest is tracked separately from the existing forest as it is younger and contains soil nutrients from the original disturbed land, not the existing forest.When a forest is converted to an impervious surface or lawn, new cohorts are created to distinguish the new impervious surface or lawn from the original.A full description of the cohort approach is described in Felzer [19] and Felzer and Jiang [10], but its application here is relatively simple as there is only a single year of disturbance in each experiment.Since forest consists of both deciduous and evergreen broadleaf types, adding additional forests or creating new cohorts from the forests creates two new cohorts.In these experiments, there are therefore 7 cohorts prior to 2030 (2 forests, impervious surface, barren land, lawn, crop, and pasture) and 9 cohorts after the disturbance in 2030.Following a disturbance, the area of the original cohort that has been depleted is adjusted accordingly.The output is then area-weighted for each of the cohorts.
The initial model runs started with the 7 land cover cohorts, as follows: impervious surface, barren field, deciduous forest, coniferous forest, lawn, cropland, and pasture.The area of each cohort was used in TEM as postprocessing to area-weight the output.To examine land use effects, a change to the land use cohort areas was implemented in 2030.The decision to implement LULCC in a single year is a simplification of a more realistic gradual change, but an advance of the assumption of two different states of LULC, which would not capture the carbon effects of that land use conversion and decomposition of its products.The changes were 10%, 25%, and 50% of the area of the cohort that was being converted.All the experiments involved conversion to or from forests in the year 2030, and the forest cohorts were split into two cohorts to represent coniferous and deciduous forests.When forest area is being converted into either impervious surface or lawn, we assume the equivalent of agricultural clearance in the year 2030.During the year of the disturbance, 60% of the vegetation is burned.The remaining 30% goes into the 10-year pool, and 10% goes into the 100-year pool, where they decay at a rate of 10% and 1% per year, respectively [10].Felzer and Jiang [10] tested several assumptions about the distribution between the amount of timber going into the fuelwood and the 10-and 100-year product pools for both agricultural clearance and timber harvesting and concluded that the net effect on cumulative NCE was relatively minor, while there is a redistribution between the pools.Harvested crops go into the 1-year pool, where they are assumed to be consumed within the same grid where they are grown within the year.The seed is harvested, so the seed carbon goes into the 1-year product pool.The residual crop from the harvest returns to the soil as stubble.TEM-Hydro has undergone extensive validation.Felzer et al. [5] validated the model for several watersheds in the eastern U.S. for evapotranspiration (ET) and runoff.ET showed a strong correlation, and most basins had a positive Nash Sutcliff value for runoff.Felzer et al. [6] validated the model at several grassland and forest eddy covariance sites in the western U.S. for ET and NEP and discussed the particularities at each site that may lead to disagreements.Generally, the model, like other models, underestimates peak summer NEP in forested sites but is closer to biometric estimates.Validation for western watersheds [6] generally shows a good fit for runoff, especially for forested watersheds.
Felzer et al. [4] show how including disturbance in forests increases the NEP to more closely resemble observed values, reasonable crop yields for Pennsylvania, and reasonable dissolved inorganic nitrogen (DIN) leaching for watersheds in Lehigh Valley, which are within the range of observed values.Felzer and Sahagian [20] compared gridded TEM-Hydro trends from 1982 to 2000 of ET, GPP, and NEP to trends based on an amalgamated remote sensing and an eddy covariance-based product for the conterminous U.S. and showed an agreement in trend direction in the majority of grids.Felzer and Jiang [10] assessed TEM-Hydro NEP and net carbon exchange (NCE) estimates relative to other modeling (forward and inverse) and inventory studies and discussed how to appropriately compare the different data and concluded that TEM-Hydro compares favorably when accounting for the full range of processes and is in the range of other models.

Datasets Input for the Terrestrial Ecosystem Model (TEM)
The input climate variables for TEM are temperature, diurnal temperature range, vapor pressure, incoming shortwave radiation, precipitation, and surface wind speed.MACA-downscaled and bias-corrected data were used to provide future climate data.The area sizes in the LULC dataset provided from the imagery classification were used to run the model based on several different landscape regimes.CO 2 levels were taken from historical data (1979-2019) [21] and future data from the RCP8.5 (2020-2099) scenario [22].TEM also requires soil texture and elevation, ozone, and nitrogen deposition, which were taken from standard datasets [10], with ozone and nitrogen deposition levels held constant for the future.The products of TEM include carbon, nitrogen, and water fluxes.
Crops were fertilized but not irrigated, because Lehigh Valley is in a mesic climate.The lawn cohort was fertilized and irrigated.Pasture was neither fertilized nor irrigated but received natural fertilization from the livestock, as explained in Felzer [4].The pasture, lawn, and crop cohorts were irrigated.For irrigated cohorts, if grid precipitation in a month during the growing season is less than 200 mm, the cohort was irrigated to the equivalence of 200 mm of precipitation.Cropland and lawn were fertilized starting in 1979, at around 15 gN/m 2 yr −1 , while the pasture cohort was fertilized with 5 gNm −2 yr −1 .Management processes are further described in Felzer [4].Impervious surfaces were treated as 100% clay, as water cannot penetrate as easily through clay; additionally, the size of the bucket was set to zero to prevent any water from penetrating.MACA uses a statistical method that produces downscaled and bias-corrected climate data from global climate models from the Coupled Model Intercomparison Project 5 [23] to capture patterns of daily near-surface meteorology (surface temperature, daily temperature range, vapor pressure, solar radiation, precipitation, wind) [24] for the future.The MACA statistical downscaling method utilizes a training dataset (GridMET) to remove historical biases and match spatial patterns in climate model output [24].Data provided by MACA were converted to monthly climate data and provided future climate data ranging from 2020 to 2099 under the RCP8.5 scenario for the NCAR CCSM4 r6i1p1 ensemble (Figure 1a-f).The data provided from MACA is ~4 km grid resolution (0.047 • latitude by 0.036 • longitude).

GridMET
GridMET is a historical dataset of daily high-spatial resolution climate variables ranging from 1979 to the present [23].GridMET was chosen for the model spin-up because MACA was originally bias-corrected using GridMET, so no further bias correction was needed to connect the historical and future datasets.The GridMET data were converted to a monthly timescale.GridMET data are also available at ~4 km resolution.The climate variables used from GridMET were the same as for MACA and included surface temperature, daily temperature range, vapor pressure, solar radiation, precipitation, and wind.This dataset was used solely for the purpose of providing model initial conditions in 2020, rather than re-equilibrating at that year.The model was initially equilibrated based on the 1979-2019 GridMET climate and then run in the transient version of the model, with the final conditions in 2019 serving as the initial conditions for the future runs, starting in 2020.All results are reported solely from the starting and ending period of the future runs.

Experiment Design
For the future (2020-2099) period, a control run without LULCC and experiments with LULCC were run to examine the impact of LULCC on carbon and water dynamics.The control experiment consisted of the future model run in which the LULC remains at 2017 values.Sensitivity experiments were then performed using monthly transient climate from 2020 to 2099, to experiment with different percentages of LULCC for the city of Bethlehem based on the types of LULCC outlined in the CAP [3].Different percentages of land use conversion (10%, 25%, and 50%) were used to determine what threshold of change results in a significant difference from the control.The change to LULC will be implemented in the year 2030 when new cohorts are created from the source cohort.The new cohorts are land use types created from the preexisting land use types following the disturbance.While the Bethlehem CAP is designed to increase sustainability (experiments 1 and 2), two experiments were added to see what the impact of decreased forest (experiments 3 and 4) would be and simulate further development (Table 1).When a forest is converted to an impervious surface or turflawn, a disturbance is created in the forest cohorts, as described above; the area is reduced, and new cohorts of impervious surface or turflawn are created.When impervious surfaces or turflawns are converted to forest, the area of the impervious surface or turflawn is reduced, and new forest cohorts are created.The newly created cohorts have different properties from the original cohorts because the forest stand age is reset and, in the case of forest to impervious surface, the soil texture is changed to 100% clay, and water penetration into the surface is allowed.Additionally, no evapotranspiration was allowed for sitting water on impermeable surfaces, assuming the water runs off fast enough to prevent evapotranspiration at the source location.From the original classification image used for Lehigh Valley, a subset of Bethlehem was extracted.The Bethlehem run used the same inputs as Lehigh Valley after being cropped to only retain the grids of Bethlehem.The city was divided into four regions, north, south, east, and west, to study how the same percentage of LULCC could affect the regions with their differing proportions of each land use type.The northern section of Bethlehem has more green space than the other sections, and South Bethlehem (locally referred to as the South Side) has the highest proportion of impervious surfaces.The values for the city of Bethlehem were recalculated to account for accurate estimates of LULC for the four grids representing the four quadrants of the city.

Statistical Analysis
To analyze the significance of the results, a two-tailed Student's t-test was used, as the difference of means between the control run and the run with LULCC change for each time period (2030-2049, 2050-2069, and 2070-2099).The results were analyzed at an α = 0.05 level of significance.For the carbon storage fluxes for Lehigh Valley (NEP and NCE), the t-test was performed on the cumulative values to show the compounding effect of the change in carbon storage throughout each time period.Additionally, analysis of variance (ANOVA) using a post hoc Tukey approach was used to determine the significance of the difference among the 10, 25, and 50% cases.ANOVA analysis is shown by the different letters either above or below the bars on the graphs.For each subset of experiments, the letters that are the same do not have a significant difference.Differing letters have statistically different means, e.g., a bar with "a" and a bar with "b" have statistically different means from each other.The statistical significance of the climate trends and the base run were also assessed at the α = 0.05 level using the approach described in Felzer and Sahagian [20].The statistical analysis was performed using IBM SPSS 29.0.1.0statistics software.

Climate Conditions
For the future (2020-2099), climate conditions indicate a warmer (Figure 1a) climate.Air temperature and vapor pressure are both modeled to increase significantly (Figure 1e).The daily temperature range is predicted to increase significantly, indicating higher daytime air temperatures (Figure 1b).The reduced wind speeds may be a result of overall warming as the equator-to-pole temperature gradient is reduced (Figure 1f).The small increasing trend in precipitation is not statistically significant.Along with increasing CO 2 in the atmosphere (resulting in 927 ppm CO 2 by the year 2099), warmer air temperatures can be expected to increase productivity and biomass in a region that is not moisture-limited, like Lehigh Valley.However, the same climate conditions can also lead to a decrease in soil carbon; as the temperature increases, the conditions become more favorable for heterotrophic respiration in the soils.The future climate for Lehigh Valley with this scenario includes higher air temperatures and more humid conditions, which are consistent with expectations from a warming scenario.These climate trends would suggest increases in productivity and biomass, as long as ecosystems do not become moisture-limited.

Lehigh Valley Control Run: Changing Climate and CO 2
The control run shows the state of current carbon and water dynamics, with the LULC of 2017 remaining constant, but with changing climate and CO 2 .For the future period, solely as a result of the changing climate and CO 2 , the vegetation carbon (VegC), and vegetation carbon fluxes (net primary productivity (NPP) and gross primary productivity (GPP)) show a statistically significant increase during the 21st century (Table 2).There is also a statistically significant decrease in soil carbon.While the ecosystem carbon fluxes (NEP and NCE) have a slight increase, the increase is not significant.Overall, the water variables (evapotranspiration, runoff, and soil moisture) do not change by a large amount, with only ET showing any significant increase.These results are consistent with a warmer climate and higher CO 2 levels increasing vegetation productivity, warming causing more decomposition, and no large change in precipitation.Experiment 1 (L2F) involved the creation of new forests from areas that were previously lawns.Results are shown as percent differences from the control (Table 3), as displayed by bars (Figure 2), time series (Figure 3), values for different PFTs (Table 4), and maps (Figures S2 and S4-S7).By the end of the century, results show a significant increase in vegetation carbon by 5.7% and a decrease in soil carbon by 8.1%.The soil carbon also decreases in all cases (Figure 2, Table 4a).There are nonsignificant increases in the carbon fluxes and decreases in all the water variables.Through time, only the vegetation carbon increases for the 25 and 50% cases and all the soil carbon decreases are generally outside the range of interannual variability (IAV) of the control, while there is not a significant increase in cumulated NCE (Figure 3).The increasing vegetation carbon throughout the future shows the growth of the forested area over time from the newly created forest cohorts since the climate conditions are the same.Soil carbon decreases with the conversion to forest, however, which is due to the calibrated parameters used to describe the lawn cohort (see Section 4), resulting in higher heterotrophic respiration in those areas when they become newly forested.The largest increases in productivity and carbon storage occur in urban and surrounding areas, as they tend to have more lawns (Figures S2, S4 and 5a).

Experiment 2 Lehigh Valley: Conversion of Impervious Areas to Forest
Experiment 2 (I2F) involved the conversion of impervious areas to forested areas by 10, 25, and 50%.The only significant changes are increasing soil carbon in the 25% and 50% cases and reduced ET and increasing soil moisture in the 25 and 50% cases (Figure 2).Other nonsignificant changes are an overall increase in carbon stocks and fluxes and less runoff (Table 4b).The increase in soil moisture is consistent with no soil-holding capacity for impervious surfaces.When converting from impervious to forested areas, urban grids that are located near the center (the majority impervious) have the largest increase in carbon storage (Figure S4b).These are also the grids that generally see the largest increase in NPP and vegetation carbon (Figures S5b and S6b).The time series (Figure 4) shows that all the changes are within the control IAV.
It would be expected that converting impervious to forest would have a larger magnitude of change for a similar amount of land change than converting lawn to forest.However, due to the fertilization of lawns, forest regrowing in areas that were impervious is more nitrogen-limited than those regrowing in previously lawn areas.The resulting nitrogen limitation in previously impervious areas inhibits the growth and vegetation in these regions.

Experiment 3 Lehigh Valley: Conversion of Forested Area to Lawn
The runs in experiment 3 (F2L) converted forested areas to lawns.The LULCC decreases all the carbon variables but increases the water variables (Figure 2, Table 4c).The vegetation carbon and soil carbon decreases are significant in all cases, while the carbon fluxes are all significant for 50%.There is a strong decrease in cumulative NCE in all cases.The ET increases are significant for 25 and 50%, and the soil moisture increases for all cases.In the initial 2030-2049 period, there is a large decrease in NCE due to the conversion of forest into lawn (Figure 3f).When forests are cut down, it is assumed that 60% of forests are burned, resulting in a large negative NCE.The NCE is still negative in the other periods following 2030-2049; however, it is not as negative since the change occurred in 2030 and a majority of the products have burned or decayed in that period.The cumulative NCE is significantly negative in all cases (Figure 3f).The increase in soil moisture is due to the irrigation of lawns, which also produces increases in runoff, although the increase in runoff is not significant.The northern and southern sections of Lehigh Valley are more heavily forested and see the largest decreases in NCE, NPP, and vegetation carbon when converting forested areas to lawns (Figures S4c, S5c and S6c).

Experiment 4 Lehigh Valley: Conversion of Forest to Impervious Area
The runs in experiment 4 (F2I) converted forested areas to impervious surfaces.Unsurprisingly, there are large decreases in vegetation carbon and the vegetation and ecosystem carbon fluxes, with all decreases significantly negative for all cases except for 10% NPP and NEP (Figure 2, Table 4d).The ET and soil moisture decrease in all cases, and runoff increases in all cases.This increase in runoff is expected as water cannot percolate the impervious surface.
There is an increase in soil carbon, however, which is due to the assumption that 85% of the material that is burned remains in the soil, and there is no heterotrophic respiration in impervious surfaces.The burning of wood for the conversion, as well as the carbon entering the soil, can be seen immediately following the conversion in 2030 (Figure 3a,f).When the forested area is burned during the conversion and a large portion of the burned material is added to the soil, it appears that there is a large increase in soil carbon, even though the impervious surface is calibrated to assume no carbon.Although soil carbon increases for the cases in which the forest area is burned, the trend for soil carbon throughout the future period is still negative (Figure 3b).The 50% change case produces results outside the IAV of the control for most cases, with the cumulative NCE changes being negative and significant for all cases (Figure 4f).The forested areas along the northern and southern sections of Lehigh Valley see a decrease in vegetation carbon (Figure S5d) and productivity decreases that are larger in magnitude than the experiment in which the forest is converted to lawn (Figure S4c,d).The decrease in NCE is somewhat similar between the two cases (Figure S3c,d), with imperviousto-forest having more grids showing a larger decrease in NCE than when converting from lawn.Since lawns still have the capacity for carbon storage and productivity, they see less of a decrease in carbon than impervious surfaces.However, turflawn still does not have the capacity for carbon storage and productivity that forests do.

Experiment 1 Bethlehem: Conversion of Lawn to Forest
The results are shown as percent differences from the control (Table 5), as displayed by bars for the 50% case (Figure 5) and values for different PFTs (Table 6), but divided between the four regions of the city.For the first experiment (L2F), in the city of Bethlehem, vegetation carbon increased while soil carbon decreased, and the carbon fluxes all increased (Figure 5, Table 6a) in all the regions.In all regions, vegetation and soil carbon saw significant changes for the 10%, 25%, and 50% cases (Table 6a), while GPP and NPP generally saw significant increases for the 25% and 50% cases.Cumulative NEP and NCE are significantly positive in all regions.The region that saw the largest percent change increases in biomass, carbon storage, and productivity is the South Side of Bethlehem, which is usually significantly different from the other regions.While the South Side sees the largest percent increase, the region still has the lowest absolute values of carbon storage and productivity of all the regions (Table 5) due to having the highest proportion of impervious surfaces to forest area.The change from lawn to forest results in a decrease in soil moisture and runoff due to areas defined as lawns being irrigated, so the lawn cohorts have higher soil moisture than the forest cohorts per unit area.

Experiment 2 Bethlehem: Conversion of Impervious Areas to Forest
Experiment 2 (I2F) for the city of Bethlehem converted impervious surfaces to forest by 10, 25, and 50%.Across the regions of Bethlehem, all the carbon variables increase (Figure 5, Table 6b).This change resulted in increases in soil moisture and ET, with reduced runoff in all regions.For the same percentage of the area changed during the same period, the South Side has the largest increases in carbon stocks and fluxes, with increases of 77, 33, and 510%, respectively, for vegetation carbon, NPP, and NCE.However, this region still has the lowest biomass and carbon when compared to the other regions due to the higher proportion of impervious surfaces (Table 5).Cumulative NEP and NCE are significantly positive for all regions, and especially pronounced for the South Side (Figure 5).

Experiments 3 and 4 Bethlehem: Conversion of Forest to Lawn and Impervious
In the third (F2L) experiment for the city of Bethlehem, forest area was converted to lawn area by 10, 25, and 50% to simulate deforestation.The LULCC resulted in decreases in the carbon stocks and vegetation and ecosystem carbon fluxes (Figure 5, Table 6c).The percent decrease (Table 6c) for all the regions in vegetation carbon is similar for each case.However, the absolute values show that North Bethlehem decreased the most in vegetation carbon, NPP, and carbon storage due to the region having the highest proportion of forest.Experiment four (F2I) produced similar results in carbon to experiment three, except the decreases in GPP and NPP are larger.(Figure 5, Table 6).For F2L, soil moisture and ET increase, while in F2I, they decrease, and runoff increases in both cases, but much more so when converting to an impervious surface, while lawn irrigation leads to an increase when converting to lawn.Cumulative NEP and NCE are significantly positive or negative for all regions for both F2L and F2I.

Biomass of Bethlehem
A goal of the Bethlehem CAP is to gain an estimation of canopy cover for the entire city of Bethlehem, including both street and park trees (measured by the city), as well as private property and large forested areas, such as Lehigh University.After extracting the forested areas in Bethlehem, the average vegetation carbon for the year 2020 was calculated to provide an estimation of the biomass for the city of Bethlehem.The section of Bethlehem that contains the most biomass is North Bethlehem, which has the highest proportion of forested area to impervious surfaces (Table 7).The estimated total biomass of forests in Bethlehem is 240 TgC.Future work is to calculate the biomass for Bethlehem using the tree inventory collected by the city and compare it to the calculated biomass from this study.This process would require applying species-specific allometric equations to the diameter at breast height (dbh) values measured by the city's tree inventory.The biomass calculated from this study is expected to be higher due to the inclusion of trees in private properties and other areas not accessible in the ground survey conducted.

Assumptions
This study explored the effects of land cover change among impervious surfaces, turflawns, and forests in Lehigh Valley and within the city of Bethlehem.The results were based on model runs with and without LULCC using the MACA-downscaled and bias-corrected NCAR CESM model RCP8.5 scenario from 2019 to 2099.The results were all presented as differences from a control run without any LULCC change.Lehigh Valley is the equivalent of one grid of the CESM climate model, which means that for future data, MACA only captures the trend of a single grid.The bias correction is important to ensure that the range of climate variation throughout the time period being modeled is reasonable and therefore responsible for the correct range of ecosystem responses.The bias correction in MACA is performed with respect to the GridMET data, which provided the historical climate data [23].We examined the historical and future climate input variables to ensure they were properly bias-corrected.MACA has been shown to perform extremely well in temperature, humidity, wind, and precipitation through the methods of calculation, allowing for the joint downscaling of temperature and dew point temperature [24].This study uses the same climate conditions for all the experiments; thus the climate serves as a control, and the only thing changing between the experiments is the LULC datasets.
Implementing the change to LULC in a single year may seem unrealistic, but the total area change, if spread out over 80 years, is consistent with historical rates of change.The area change generally falls within Pennsylvania's average annual land cover change of 115,981 km 2 yr −1 [25].While the 50% conversions of forests to lawn and impervious cohorts in experiments 3 and 4 fall within that range when accumulated over 80 years, the 10% experiments are less than that average rate.When breaking down the land use change provided by NRCS into land use types, the average rate of land use change in developed areas of 15,302 km 2 yr −1 is less than the amount of area converted in this study in a single year, with even the 10% conversion of impervious to forest changing 18,823 km 2 in a single year (Table S2).However, the rates of change we imposed are consistent with historical rates but just implemented in a single year rather than spread out over 80 years.
To explore the effects of applying the LULCC in a single year, we repeated experiments 6 and 12 with the same total change applied on a decadal basis.When converting from impervious to forest (Figure S8a), there is virtually no difference in carbon or water variables (Figure S8b,c) because the new forest is still very young and therefore does not substantially increase the vegetation carbon.While young, growing forests are stronger carbon sinks, the relative area each decade is still fairly small.When converting from forest to impervious (Figure S9a), however, the change in vegetation carbon and cumulative NCE (Figure S9b) takes a different path, with a more gradual change than the original results, but ends up in the same place by the end of the century, although cumulative NEP is virtually the same.The runoff (Figure S9c) is less at first since not as much impervious surface is being created immediately but ends up at the same values by the end of the century.
Studies have shown that at least 10% LULCC is required to see a significant change in carbon storage [1].In this study, the only experiments in which over a 10% change in the total region occurred were in the 25 and 50% conversion of forest to either lawn or impervious (Table S3).The 25% and 50% experiments for F2L and F2I were outside the IAV range for the stocks (vegetation C, soil C, and soil moisture), but the 10% case was also outside the IAV for the two carbon stocks, and water variables for the F2I case (Figures 3 and 4).The F2I cases for GPP and NPP are also outside the IAV for the 25% and 50% cases.Generally, stocks have less IAV than fluxes do, resulting in the fluxes remaining within the IAV, but not the stocks.Changes in NEP or NCE are generally smaller as they are the difference between two fluxes (GPP and ecosystem respiration), which often change in tandem.The IAV of the vegetation and ecosystem carbon fluxes and the water variable fluxes tend to be quite large.The resulting variability in climate from year to year impacts the ecosystem's carbon fluxes.Due to this large IAV, the change in LULC would need to be large in order to exceed the naturally large IAV in the ecosystem carbon fluxes.Both temperature and moisture have large impacts on productivity [26][27][28], which can create large IAV dependent on the climate conditions of that year [29,30].For instance, if a year is warmer and wetter than the previous year, productivity will also increase, resulting in increased IAV.

Young vs. Old Forests
Young, fast-growing forests do not possess the same capacity for carbon storage as old-growth forests [31,32].The newly forested areas do not have the same carbon storage that mature forests do, resulting in lesser immediate impacts from reforestation than deforestation.While young to middle-aged forests have higher productivity (NPP) than mature forests [33,34], mature forests have higher carbon storage than younger forests [35].Reforesting efforts have been found to only make up for 10% of the carbon storage lost in deforestation events [32].The impact of carbon loss from deforestation is much more profound than the sequestration of carbon from reforestation in the period right after the LULCC.When forested area is reduced, there is a larger negative impact on carbon storage than when forest area is created, due to both the lower rate of carbon storage by the young forest and the immediate release of carbon upon timber harvest disturbances, especially if some of the wood is burned.The assumption that 60% of the forest area is burned in the same grid in which it is harvested results in a larger decrease in NCE initially than what would be expected.In a real-world scenario, it is possible that more of the harvested wood would likely be exported outside of the grid as timber harvest.
However, younger forests have important positive consequences for ecosystem services.In terms of carbon, they sequester more carbon than mature forests, even if they store less, although mature forests continue as small carbon sinks [36].They will also eventually grow into mature forests and store the carbon they have sequestered during their lifetime both in their wood and soils [37].Mature forests provide a source of oxygen, fix nitrogen, improve air quality, and serve an important ecological niche, as they provide unique microclimates and habitats [38].The carbon stored in larger trees will be sequestered even longer if the products are used as timber and replaced by younger trees to take up even more carbon, which is included in the TEM model.Conversely, the destruction of mature forests without sequestering the wood releases a large amount of CO 2 into the atmosphere.There is also important cultural, religious, and medicinal value to older mature forests [38].

Conversion to Forest
When comparing L2F to I2F, a primary reason that carbon storage and productivity increase more for the same percentage change in L2F is due to the fertilization of lawns.Since the lawn cohort is fertilized, the new cohorts created in those areas continue to receive nitrogen, allowing them to grow better than I2F.In an additional model run to examine this effect, deciduous trees grew better with nitrogen fertilization and reached a higher mature biomass.While coniferous trees are also not nitrogen-saturated, the fertilization effect was not as impactful in these trees and will require further examination.In a realworld scenario, the soils in formerly impervious areas have fewer nutrients than the soils in greenspaces [39].In order to maximize carbon storage and productivity in formerly impervious areas, fertilization would be required.For both cases in which new trees are grown, the future period of 80 years is not enough to reach maximum vegetation carbon.To enter the old-growth stage, forests often need to exceed the 100-year stand age [40].Based on the definition of old-growth forest, the newly created forest from both lawns and impervious areas will need more time to reach full maturity.However, existing trees in Lehigh Valley are not old-growth, so the biomass of the trees at initialization may be larger than reality, which would further enhance the loss of biomass due to deforestation.
The relation between vegetation and soil carbon is sometimes counterintuitive.In the L2F case, there is increasing vegetation carbon but decreasing soil carbon.While the percentage differences are great for the reduced soil carbon, the absolute differences are actually less for the changes in soil carbon due to the larger absolute value of vegetation carbon.However, turflawn is calibrated to a much larger value of soil organic carbon (16,260 gCm −2 from Jo and McPherson [41], 1995 vs. 8290 gCm −2 for temperate forests from Gaudinski et al. [42]).Therefore, when lawn is converted to forest, the heterotrophic respiration rates are much higher for the new forest cohort, leading to a reduction in soil carbon in spite of the increase in forest litter.On the other hand, when forest is converted to impervious surface, heterotrophic respiration rates are set to 0, so there is no loss of any existing soil below the surface and the amount of soil carbon is larger than a run without the conversion.Studies have shown that urban soils, especially in temperate and colder climates, accumulate 3-5 times more carbon than natural soils [43].The conversion of impervious surface to forest would release carbon due to years of sealed and cultural layers, which is not considered in this study.
Altering the amount of area covered by lawns or impervious surfaces can cause profound changes to carbon and water dynamics.Overall biomass, productivity, and carbon storage do not show a strong increase when converting impervious surfaces to forests but show a strong decrease when converting the forests to impervious surfaces (Figure 2, Table 4b,d).The experiments in which forests are being converted to impervious surfaces show a larger impact than when impervious surfaces are converted to forests.However, it is important to note that more total area is changed when forests are converted due to Lehigh Valley having more forested area than impervious area.The newly created forests are also less productive due to their young stand age, while the conversion of forest to impervious produces an immediate result.Similarly, converting forest to lawn has a larger effect than converting lawn to forest.Due to the treatment of lawn being irrigated in the model, the experiment with a higher proportion of lawn has higher soil moisture, evapotranspiration, and runoff than the newly created forests.
Through the experiments, the conversion of lawns and impervious areas to forests has only a small and often insignificant increase in carbon stocks and fluxes within Lehigh Valley.Other studies suggest that increasing the amount of forested area can be effective in increasing carbon uptake from the atmosphere [1,9].The results of this study, while agreeing with an increase in vegetation carbon, offset this increase with reduced soil carbon due to the higher heterotrophic respiration rates.The increase in atmospheric CO 2 contributes to the increased biomass as forests grow back under elevated CO 2 levels [44].While the CESM scenario does not display a significant increase in precipitation for Lehigh Valley, the increased CO 2 concentrations in an environment that is not moisture-limited (since soil moisture remains steady) still results in increased biomass.

Water Dynamics
The use of a bucket model for moisture in TEM allows the model to account for the effect of soil texture on porosity and saturation [45].Since clay still has the ability to hold soil moisture, we have chosen to entirely eliminate the water-holding capacity of impervious surfaces.However, without allowing the water to run off without evaporation, most of it evaporates rapidly in the model, given enough net irradiance and low humidity, so we have had to eliminate ET for standing water on impervious surfaces in order to allow the water to run off.
When converting impervious surfaces to forests, the runoff would be expected to decrease [46].The benefits in reducing runoff, however, are dependent on the location of the reforestation and can be diminished if there are impervious surfaces upland from where the reforestation is occurring [46,47].When increasing forests from impervious surfaces, the bucket becomes larger, so soil moisture also increases and runoff decreases.In both experiments involving forest and impervious surfaces, evapotranspiration decreases.In the case of converting impervious to forest, soil evaporation decreases due to shading while transpiration increases due to tree cover, but the effect of the soil evaporation decrease dominates, especially since transpiration is small while the forests are still young with a less developed canopy and smaller leaf area index.The conversion from forest to impervious eliminates both transpiration and soil evaporation, although there can be consideration of how much sitting water evaporates from the surface, which we have eliminated in this study.

Bethlehem
The results within the city of Bethlehem show that the area with the largest potential for increasing carbon storage and productivity is the South Side of Bethlehem.Compared to the map of social vulnerability, this area aligns with South and East Bethlehem, which are the most vulnerable regions within the city based on the map of social vulnerability and have the highest heat island index [3].The high proportion of impervious surfaces to greenspaces, which is typical in lower-income neighborhoods, results in lower productivity and carbon storage.In terms of the Bethlehem CAP, the city should prioritize increasing greenspace in the South Side and east Bethlehem, as those are areas that will see the largest improvement in carbon storage.However, unlike east Bethlehem, which is a heavily industrial area, the South Side is a densely populated residential area.Thus, even though both the South Side and east Bethlehem would largely benefit from increasing greenspace, the South Side should be prioritized.Increasing the number of street trees has been found to reduce rates of violent crimes, with underserved neighborhoods benefiting the most [48,49].The neighborhoods that would gain the most in carbon storage from the increase in tree cover are also the most socially vulnerable [3].In addition, with these areas being the most vulnerable, increasing green space will help bridge the gap in environmental equity.The high proportion of impervious surfaces also results in a higher urban heat island effect, which creates dangerous health conditions for the citizens within those neighborhoods [50].By increasing greenspace within South and East Bethlehem, there is also a potential for reducing the impact of the urban heat island effect [51].
However, there are some negative consequences to increasing green space too much (for example, in the 50% case) in lower-income neighborhoods, in terms of affordable housing.Bethlehem currently has only a 2% vacancy rate, far below the 6-8% considered sustainable to provide affordable rents and prevent homelessness.Home sale prices and rents have risen dramatically from 2019 to 2023, by 50% and 40%, respectively, creating unaffordable rents for much of the population.The city has developed a comprehensive plan [52] to address this housing crisis.However, environmental gentrification often comes at the expense of affordable housing by driving real estate prices up and forcing low-income residents to leave [53].On the other hand, affordable housing can be built with significant areas preserved as green and recreational space, as recently proposed by Bethlehem's First Presbyterian Church [54].There needs to be a balance between the two needs that best provides environmental justice for lower-income residents.

Magnitude of Carbon Sink
The average person has a carbon emission of 4.73 Tonnes Cyr −1 , or 4.725 × 10 −6 TgCyr −1 [55].As of 2019, the population of Lehigh Valley is 840,000 citizens, resulting in 3.97 TgC being produced annually.For the 2070-2099 period, the 50% change from impervious to forest results in an 8.2 TgCyr −1 increase in NCE, which is equivalent to twice the annual carbon emissions from the current population.When comparing the 50% change of lawn to forest for the same time period, there is a 30.0 TgCyr −1 increase in NCE, which is equivalent to 7.5 times the annual carbon emissions from the current population.The L2F change produces a larger increase than I2F, likely due to the higher nutrient content in the L2F soils providing better growing conditions for the new trees.Certainly, adding fertilizer when converting impervious surfaces to forests would increase forest growth even more.Note that these increases are the result of the 50% L2F and I2F LULCC changing 4 and 5% of the total area of the region, respectively (Table S3).To further develop greener cities and mitigate the carbon footprint of cities, increasing greenspace even more will be necessary.
For Bethlehem, PA as of 2021, the population is 75,000, equating to 0.35 TgC being produced annually.For the 2070-2099 period, combining the regions of Bethlehem, there is an increase in NCE of 0.8 TgCyr −1 for the 50% change in L2F, which is equivalent to 2.3 times the annual carbon emissions of the current population.For the same period and percentage of land use change, the change in NCE for I2F is 0.75 TgCyr −1 , which is equivalent to 2.1 times the carbon emissions of the current population.The values are much closer for Bethlehem because of the extremely high proportion of impervious surfaces (Figure 6, Table S4) compared to that of lawns within the city.While the growing conditions are better in L2F, the amount of land converted in I2F for the city of Bethlehem outweighs the lack of nutrients and therefore has a larger increase in carbon storage as a result.

Hypotheses Confirmation
The original hypotheses were partially supported.However, there were some significant changes when less than 10% of the total LULCC of the region was changed (Table S3), for example, particularly for the stocks of vegetation carbon, soil carbon, and soil moisture, as well as ET in several of the L2F and I2F cases (Figure 3).The results showed that the decrease in forests has a larger impact than increasing them by the same amount.The effect of the burning of the products from deforestation makes the immediate impact of the experiments in which forests are converted to lawns or impervious surfaces more profound.The low carbon storage of young trees reduces the immediate impact of reforestation (Figures 3 and 4).Reforestation had the largest effect on the South Side in terms of carbon stocks, vegetation carbon fluxes, and ecosystem carbon fluxes (Figure 5).Since this region is the most socially vulnerable section of the city, increased carbon sequestration is an ancillary benefit to greater tree cover for the residents.

Hypotheses Confirmation
The original hypotheses were partially supported.However, there were some significant changes when less than 10% of the total LULCC of the region was changed (Table S3), for example, particularly for the stocks of vegetation carbon, soil carbon, and soil moisture, as well as ET in several of the L2F and I2F cases (Figure 3).The results showed that the decrease in forests has a larger impact than increasing them by the same amount.The effect of the burning of the products from deforestation makes the immediate impact of the experiments in which forests are converted to lawns or impervious surfaces more profound.The low carbon storage of young trees reduces the immediate impact of reforestation (Figures 3 and 4).Reforestation had the largest effect on the South Side in terms of carbon stocks, vegetation carbon fluxes, and ecosystem carbon fluxes (Figure 5).Since this region is the most socially vulnerable section of the city, increased carbon sequestration is an ancillary benefit to greater tree cover for the residents.

Conclusions
Replacing lawns or impervious surfaces with tree cover can be an important part of sustainable development in urban areas.However, a relatively large percentage of 10% of the area will have to be converted to see a noticeable change in carbon storage, although changes in some variables are evident with as little as 3% land use conversion.Due to the maximum amount of each land use type that could realistically be converted, to achieve more sustainable goals, multiple land use changes would need to occur in tandem, such as replacing large areas of both impervious surfaces and turf lawns with forests.Converting lawn to forest has a larger effect in this study than converting impervious surface to forest due to the higher nutrients in the soils of the lawn, but fertilizing converted impervious surfaces should increase their carbon storage capacity.While these transitions may be difficult, the replacement of parking lots or the inclusion of more trees in residential areas and lawns or even city parks will be crucial for reaching these goals.The South Side also has many parks that are solely abandoned fields, making these areas ideal candidates for reforestation efforts, although conversations with residents also reveal the importance of open space for children to play.
Based on the modeling experiments, higher proportions of forest to impervious surface and lawn result in higher productivity and carbon sequestration, which agrees with other studies [1,2,9].With urbanization expected to increase in the future due to increasing populations, there is a potential for a reduction in the amount of carbon storage due to the removal of trees.Even after increasing tree cover, the forests will not have the same carbon storage capacity as mature forests, which reduces the immediate impact, although younger forests sequester carbon at a faster rate than mature forests.Therefore, further deforestation will have a larger impact than planting new trees in the same size area.While the younger trees will not store as much carbon in the short term, the negative impacts of deforestation are felt immediately.
With an overarching goal of producing more sustainable cities capable of sequestering more carbon and decreasing the effects of climate change, it is clear that increasing the number of trees within cities will be necessary.The sensitivity experiments show that even moderate changes to land use (10%) can have a positive impact on carbon sequestration within cities.In several scenarios with less than 10% change in total area, the results end up being significant, suggesting that even smaller changes can be impactful.Urban areas with high proportions of impervious surfaces or abandoned fields have the largest potential for change and should be prioritized when considering implementing land use changes with the goal of sustainability.In addition to enhanced carbon sequestration, added tree cover reduces crime rates and reduces the urban heat island, making cities more livable and reducing environmental inequities.

Figure 1 .
Figure 1.(a-f) Time series showing the future climate data under RCP8.5 using the CCSM4 model as downscaled and bias-corrected by MACA for the period 2020-2099.Bolded equations denote a statistically significant trend at the α = 0.05 level.(a) Temperature, (b) daily temperature range, (c) net irradiance, (d) precipitation, (e) vapor pressure, (f) wind speed.

Figure 2 .
Figure 2. (a-i) Bar charts showing the percent change for 10, 25, and 50% land use changes in Lehigh Valley.Within each subgroup, the changes in the area are 10, 25, and 50%, respectively.From left to right, the categories are lawn to forest (L2F), impervious to forest (I2F), forest to lawn (F2L), and forest to impervious (F2I).Post hoc ANOVA testing is displayed, indicating which experiments are statistically different from each other in each category.Bold letters indicate if the run was statistically significant compared to its control value, based on a two-tailed Student's t-test at the 95% confidence level.(a) Vegetation carbon, (b) soil carbon, (c) GPP, (d) NPP, (e) cumulative NEP, (f) cumulative NCE, (g) ET, (h) runoff, (i) soil moisture.

Figure 2 .
Figure 2. (a-i) Bar charts showing the percent change for 10, 25, and 50% land use changes in Lehigh Valley.Within each subgroup, the changes in the area are 10, 25, and 50%, respectively.From left to right, the categories are lawn to forest (L2F), impervious to forest (I2F), forest to lawn (F2L), and forest to impervious (F2I).Post hoc ANOVA testing is displayed, indicating which experiments are statistically different from each other in each category.Bold letters indicate if the run was statistically significant compared to its control value, based on a two-tailed Student's t-test at the 95% confidence level.(a) Vegetation carbon, (b) soil carbon, (c) GPP, (d) NPP, (e) cumulative NEP, (f) cumulative NCE, (g) ET, (h) runoff, (i) soil moisture.

Figure 5 .
Figure 5. Bar charts showing the percent change for 50% land use changes in Bethlehem.Within each subgroup, the bars represent West, South, North, and East Bethlehem, respectively.From left to right, the groups are lawn to forest (L2F), impervious to forest (I2F), forest to lawn (F2L), and forest to impervious (F2I).Post hoc ANOVA testing is displayed adjacent to the error bars, indicating which experiments are statistically different from each other in each subgroup.Bold letters indicate if the run was statistically significant compared to its control value, as described in Figure 2. (a) Vegetation carbon, (b) soil carbon, (c) GPP, (d) NPP, (e) cumulative NEP, (f) cumulative NCE, (g) ET, (h) runoff, (i) soil moisture.

Figure 6 .
Figure 6.Percentage of each cohort for each of the four regions of Bethlehem, PA.

Figure 6 .
Figure 6.Percentage of each cohort for each of the four regions of Bethlehem, PA.

the
Lehigh Valley shown on 1/24th degree grids.Results reflect the difference between the control run and the 50% land use conversions for each of the experiments for the 2070-2099 period.(a) Lawn to Forest, (b) Impervious to Forest, (c) Forest to Lawn, and (d) Forest to Impervious. Figure S5.a-d.Mapped results of Net Primary Productivity for the Lehigh Valley shown on 1/24th degree grids.Results reflect the difference between the control and the 50% land use conversions for each of the experiments for the 2070-2099 period.(a) Lawn to Forest, (b) Impervious to Forest, (c) Forest to Lawn, and (d) Forest to Impervious. Figure S6.a-d.Mapped results of Vegetation Carbon for the Lehigh Valley shown on 1/24th degree grids.Results reflect the difference between the control and the 50% land use conversions for each of the experiments for the 2070-2099 period.(a) Lawn to Forest, (b) Impervious to Forest, (c) Forest to Lawn, and (d) Forest to Impervious. Figure S7.a-d.Mapped results of runoff for the Lehigh Valley shown on 1/24th degree grids.Results reflect the difference between the control and the 50% land use conversions for each of the experiments for the 2070-2099 period.(a) Lawn to Forest, (b) Impervious to Forest, (c) Forest to Lawn, and (d) Forest to Impervious. Figure S8.Repeat of experiment 6 (50% Impervious to Forest) with same areal change implemented each decade instead of in a single year.(a) aerial change of forest and impervious surface, (b) nep, nce, and vegetation carbon for the run with single disturbance and the run with decadal disturbances, and (c) runoff for the single and decadal disturbance runs.Figure S9.Repeat of experiment 12 (50% Forest to Impervious) with same areal change implemented each decade instead of in a single year.(a) aerial change of forest and impervious surface, (b) nep, nce, and vegetation carbon for the run with single disturbance and the run with decadal disturbances, and (c) runoff for the single and decadal disturbance runs.Author Contributions: Conceptualization, B.S.F. and C.A.; methodology, B.S.F. and C.A.; software, B.S.F. and C.A.; validation, B.S.F.; formal analysis, B.S.F. and C.A.; investigation, B.S.F. and C.A.; resources, B.S.F.; data curation, B.S.F.; writing-original draft preparation, B.S.F. and C.A.; writing-review and editing, B.S.F.; visualization B.S.F. and C.A.; supervision, B.S.F.; project administration, B.S.F.; funding acquisition, B.S.F.All authors have read and agreed to the published version of the manuscript.Funding: This research was partially funded by a 2019 Lehigh Faculty Innovation Grant (FIG).

Table 2 .
Trends over all PFTs.Carbon units in TgCyr −1 ; Water units in 10 12 mm H 2 O yr −1 .Bold values indicate a trend statistically significant at α = 0.05 level.

Table 4 .
Percent differences from control run for 10, 25, and 50% land use changes from (a) lawn to forest, (b) impervious to forest, (c) forest to lawn, and (d) forest to impervious for Lehigh Valley covering the future period (2070-2099).VEGC and SOLC are TgC, GPP, NPP, NEP, and NCE are TgCyr −1 , soil moisture is 10 12 mm H 2 O, and evapotranspiration and runoff are 10 12 mm H 2 Oyr −1 .Bolded values denote a statistically significant difference at the α = 0.05 level.Note that NEP and NCE are not cumulative.

Table 6 .
Percent differences from control run for 10, 25, and 50% land use changes from (a) lawn to forest (b) impervious to forest, (c) forest to lawn, and (d) forest to impervious for the city of Bethlehem covering the future period (2070-2099).(a)Lawn to forest, (b) impervious to forest, (c) forest to lawn, and (d) forest to impervious for the city of Bethlehem, PA covering the future period (2070-2099).Results are split into the north, south, east, and west regions of the city.VEGC and SOLC are TgC, GPP, NPP, NEP, and NE are in TgCyr −1 , soil moisture is 10 12 mm H 2 O, and evapotranspiration and runoff are 10 12 mm H 2 Oyr −1 .Bolded values denote statistically significant differences at the α = 0.05 level.Note that NEP and NCE are not cumulative.

Table 7 .
Biomass for the four regions of Bethlehem shown in gC/m 2 and TgC as of 2017.