Assessing the Impact of Land Conversion on Carbon Stocks and GHG Emissions

Abstract

With the recent thrust to convert forests in Ontario’s Clay Belt to agricultural land, a vital need arises to assess the attendant effects on carbon and greenhouse gas (GHG) emissions. This paper examines the possible effect of land conversion on soil organic carbon and GHG emissions within a study area in Northern Ontario, Canada, during the next two decades under different land management schemes. The study established a framework to conduct simulations with the DNDC model for agricultural lands and the CBM for forested areas. The methodology involves a unique change detection method for models’ land cover and disturbance inputs. The work highlights the improvement in carbon simulation accuracy from better inputs to carbon models. Furthermore, it addresses modalities to ensure fewer uncertainties are introduced while merging data from multiple geospatial data sources. The simulations demonstrated that the carbon sequestration potential in the forests was almost double the soil organic carbon accumulation in the agricultural lands. Validations done for the estimation of carbon sequestered included comparisons of the carbon model outputs from field survey data from 2018–2021. In most sites, the carbon amounts from the computer models compared to those from the field survey, within limits of error. The average uncertainties in GHG emissions ranged from ~0.5% to 12.8%.

1. Introduction

Climate change impacts have a high price tag. The Canadian Climate Institute has estimated that key climate change impacts are already costing households in Canada roughly $720 per year [1].

Over the last several years, the world has witnessed a climate that is warming. The global surface temperature from 2001 to 2020 was 0.99 °C higher than during 1850–1900, and it was 1.09 °C higher from 2011 to 2020 than during 1850–1900 [2]. A result of the warming climate has been rising Crop Heat Units (CHUs). The CHU is related to the minimum and maximum temperatures for a day. These two daily values are averaged over the corn growing season, giving the daily CHU rating. In Ontario, the CHU accumulation begins on May 1 and ends at the first occurrence of a temperature of −2 °C, though the general standard for Canada is from April 1 to October 31 [3,4]. The CHU in the Clay Belt region in Northern Ontario, Canada, has risen significantly in recent years, with the annual CHU going up 25% in some parts. The rise has come with an increase in yield per acre and allowed other crops that were not traditionally in the crop rotation in the Clay Belt to be added to the rotation, which may lead to the change in land use and land cover. The Clay Belt was chosen as the region of study because there has been a recent interest converting forests in the Clay Belt to farmland, including livestock farms, taking advantage of the more favorable farming conditions [5,6,7]. With the land transformation in the Clay Belt already underway, this work examines possible effects on the carbon stocks under different land management scenarios.

The Land Use, Land-Use Change and Forestry (LULUCF) sector is critical in calculating a country’s carbon stock, emissions, and carbon sequestration potential. It has been demonstrated that the carbon cycle on land would be considerably modified by land-use and land-cover changes [8]. In particular, terrestrial ecosystems are vital as carbon sinks. Therefore, if land-use change occurs, e.g., from natural forest to managed pastureland or to farmland, there might be a loss of sink capacity in the ecosystem [8,9]. In Canada’s drive towards net-zero emissions in 2050, this sector will play a vital role in accomplishing the targets: for instance, in 2021, the LULUCF sector removed 17 Mt from the atmosphere [10]. The most recent report to the IPCC showed that of the five most significant sources of Canada’s GHG inventory, two sources came from Agriculture and Land-Use Change and Forestry [11]. This study will assess the difference in ecosystem carbon and GHG emissions from natural forests and sites converted from forest to agricultural land. The nearby forested sites are used as the baseline for the study to estimate what the emissions would be without the land conversion.

Canada has a long history of forest carbon monitoring and accounting, which entails estimating greenhouse gas (GHG) emissions and removal. Canada has participated in various conventions such as the UN Framework Convention on Climate Change, whose aim is to stabilize GHG concentrations to prevent dangerous anthropogenic interference with the climate system [12], and the country signed on to numerous agreements, such as the 2009 Bali Roadmap and in 2015 to the Paris Agreement. These agreements are legally binding, and the federal government continues to progress toward achieving the carbon emissions and sequestration targets. More recently, Canada established the Pan-Canadian Framework on Clean Growth and Climate Change, inaugurated in 2016 [13], and the 2022 Emissions Reduction Plan [14]. This Pan-Canadian Framework was said to have been developed with the provinces and territories with input from the Indigenous Peoples. The Emissions Reduction Plan is an ambitious roadmap to have the country reduce emissions by 40–45% below the 2005 levels by the year 2030 [14]. Furthermore, Canada has a vision of being a net-zero carbon emitter by 2050 [15].

These goals will require rigorous attention in carbon accounting to measure progress against the targets and ensure that the targets are met. In addition, the international climate conventions require countries to fulfil certain reporting obligations. For example, the Kyoto Protocol requires countries to monitor any changes in carbon stocks resulting from afforestation, deforestation, or reforestation activities since 1990 [16]. Canada’s forest carbon and GHG accounting is done under the National Forest Carbon Monitoring Accounting and Reporting System (NFCMARS). The NFCMARS relies on the v.3 Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) [17] as its core ecosystem carbon model [18,19]. The model uses the best available data and scientific understanding of processes related to forest cycles. The CBM-CFS3 model adheres to the Intergovernmental Panel on Climate Change (IPCC) reporting guidelines for carbon and GHG estimates [20].

Inherent in the CBM are sources of uncertainty such as data input to the model by the user, the ecological parameters used in the simulations, processes incorrectly specified or even excluded, and the model’s algorithms [17]. Such uncertainties are not uncommon with carbon models; apart from the CBM, ecosystem carbon models have been known to have sources of errors from the uncertainty in inputs, the incorrect calibration of the models, the size and/or scope of the sample data, and even insufficient field measurements. On the other hand, these factors, such as sufficient field data and the correct calibration of the models, could contribute to more robust model estimates [21]. Models are not able to simulate all the ecosystem biogeochemical processes completely. So, approximations have to be made, although uncertainties have been said to diminish over large spatial scales as the internal variability averages out [22,23]. The international climate reports also require countries to assess the uncertainties in their carbon and GHG emissions estimates. Indeed, carbon accounting is a complex endeavor, and Canada has always reported the uncertainties associated with its estimates. However, a criticism of the reports is that the CBM-CFS3 results have been published, but the input variables are not publicly available [19].

While conducting carbon analysis, it is crucial to identify sources of error or uncertainty, quantify them, and minimize them. Rather than using default or generic parameters when estimating emissions, such as is done when using IPCC defaults [10], the site-specific or regional parameters should be used, as this increases the accuracy of the model result. This paper examines uncertainties in estimating the potential carbon and GHG emissions, as well as carbon sequestration potential, when forested land is converted to agricultural land.

The motivation for this study is the ongoing endeavor by the provincial government of Ontario, Canada, to convert some forests to agricultural land. The reasons for the land conversion are for food security, given that the population on Earth is expected to increase to 9 billion by 2050 [24], the potential for increased land productivity from the longer growing seasons brought on by warmer weather, better drainage infrastructure, new technologies and crop varieties, and the low cost of land. Furthermore, the agricultural expansion could diversify the economy in that part of Ontario for the First Nations and Metis peoples, as well as sustainably expand the livestock sector [25]. In anticipation of the potential conversion from forested to agricultural land, it is important to evaluate the impact of such conversion.

This study examines how the ecosystem carbon and GHG emissions could change over the next 20 years compared to if the forests had remained untouched. The simulations will consider land transformations from forested to cropland, pasture, or land for livestock operations. The CBM is used to model the forested land, while the DeNitrification DeComposition (DNDC) [26] simulates the carbon and GHG emissions on the agricultural land. In addition, an uncertainty analysis is performed on the results using the DNDC model’s Monte Carlo uncertainty tool.

Other studies have reported on using the CBM to track carbon in agricultural lands reverting to woody land [27] or even on uncertainties calculated in carbon modeling using the CBM [18,28]. This study draws on those works while estimating the uncertainty related to forest conversions to agricultural land. In addition, the model output in the study will be compared to the soil organic carbon amounts found from field surveys in the study area. This work uses the best available knowledge and high-resolution satellite data to estimate ecosystem carbon and GHG emissions most accurately after various land management schemes. Thus, a contribution of this study is to characterize the impact of land conversion on carbon stock and GHG emissions. In previous work, a framework was generated for establishing the disturbances and land-cover changes in a study area using mostly remotely sensed data [29]. The framework was dubbed the Progressive Adaptive Method (PAM), and it was relied on to generate higher accuracy inputs for the carbon models. Thus, the framework contributed to reducing the uncertainty in the final estimates from the carbon stock simulations.

The objectives of this study were the following: (1) Identify the long-term carbon effect on the managed lands if converted to agricultural uses; (2) Determine potential emissions under various future climate conditions; and (3) Examine how the carbon and GHG flux prediction can be improved through uncertainty/sensitivity analysis.

2. Materials and Methods

2.1. Study Area

The study area is located in the Great Clay Belt area of Ontario, Canada, as shown in Figure 1. Due to rising temperatures, it is very likely that some of the forested land will be converted into pasture and/or farmland [7].

By 2007, the study region was predominantly forested, having over 70% of the area being treed and about 14% of the area being shrubland [30]. Other vegetation cover in the area consists of shrubbery, some agricultural land, and wetlands. The wetlands are both treed and wetland shrub; such ecosites can be categorized as fen, bog, and swamp. In the Clay Belt, most of the soil is heavy and rich clay loam. Within the study region, the soil is mostly of a silty clay loam texture. In addition, the tree species which were identified in the region were trembling aspen, balsam fir, black spruce, balsam poplar, eastern white cedar, eastern larch, white birch, and white spruce.

Livestock pasture sites were chosen to be simulated in this study because there will likely be an increase in livestock farming in Northern Ontario. For example, there began the 2017 initiative of a “Beef Cattle Farming Training Program for Indigenous people” [6]. Thus, the land could be used for pasture, poultry, and/or livestock rearing, as well as perhaps even growing crops since crop farming has been successfully carried out in the Clay Belt region. In fact, the soils of the Clay Belt, which are rich in glaciolacustrine and morainal calcareous clays and silts, have supported the successful growth of barley, oat, and wheat crops. More recently, the with the CHU having risen by over 25% in the region [31], more crops such as soybeans, canola, corn grain, and silage corn are being grown. The increased average temperature in the Clay Belt indicates that the region will continue to have a positive impact on crop production.

2.2. Carbon Models and Inputs for Model Simulations

Freely available and commonly used carbon models were the tools used in this study.

2.2.1. Carbon Models

Two different carbon models were used in this work. The DeNitrification DeComposition (DNDC) model v.9.5 estimated the carbon content on crop and pastureland. The DNDC model simulates carbon and nitrogen biogeochemistry in agroecosystesms. It can also predict crop growth, soil carbon dynamics, soil temperature and moisture regimes, and emissions of gases such as N₂O (nitrous oxide), NO (nitric oxide), NH₃ (ammonia), CH₄ (methane), and CO₂.

The DNDC model is a process-based model driven by four main ecological drivers: the climate, soil, vegetation, and management practices (such as crop type and rotation, tillage, fertilization, grazing, and manure application) of the area under study. Thus, the parameters required to run the model include the fertilizer scheme, the tillage applied, the climate data, and the grazing details of livestock on farmland.

The crop growth modeled by the DNDC is done on a daily time step. The soil organic matter is set in pools in DNDC, and the daily decomposition rate of the pools is dependent on such parameters as the soil clay content, nitrogen availability, and soil temperature and moisture content. The decomposed soil organic carbon (SOC) is lost as CO₂ or distributed to other SOC pools. The DNDC model also simulates the effect of rainfall, as well as NH₃ concentration in the soil, allowing for the estimation of nitrification and denitrification that occurs in the soil. Outputs from the model include the changes in water pool (from leaching, transpiration, runoff, etc.), annual change in SOC, annual flux in GHGs, and nitrogen lost due to activities such as crop uptake, weed uptake, runoff, and ammonia volatilization.

The other model used was the CBM-CFS3—the Carbon Budget Model of the Canadian Forest Sector. CBM is commonly used for assessing carbon stocks and forecasting of forest ecosystems. Several countries use this model for their carbon accounting and reporting to regulatory bodies such as IPCC. The CBM tracks carbon stocks, transfers between pools, and emissions of CO₂, CH₄, and CO (carbon monoxide).

Some of the model’s inputs include land-use classification changes, disturbances (e.g., fire, insect, deforestation, harvest impacts), and stand age. Among the numerous outputs obtainable from the tool are carbon stocks and fluxes, annual carbon transfers between pools, emissions to the atmosphere, transfers and emissions associated with disturbance types, and annual carbon stock change for the total ecosystem. CBM was primarily used to account for the carbon amounts in the forested regions of the NLP. The model simulated the carbon and GHG amounts that would have accumulated if there had been no land-use or land-cover changes. Thus, the CBM results were the baseline for the study—to be compared to the results from the DNDC models for areas that had been converted to crop, livestock, or pasture lands.

2.2.2. Inputs for Model Simulations

• Land-use type: Since each model is tailored to a particular land use, the various land uses over the study area had to be defined for both models. For most of the areas, the land use was defined by using remote sensing images from moderate- to high-resolution satellite sensors such as Landsat and Sentinel-2: First, existing land-cover/land-use maps such as Ontario’s Forest Resource Inventory [32], the Natural Resources Canada (NRCan) land cover product [33], NRCan’s Dynamic Surface Water dataset [34], and the Agriculture and Agri-Food Canada (AAFC) land-cover [33] and crop inventory [35] maps were used to derive testing and training sample areas in a reference year. Then, optical satellite data were used to predict the land use and land cover over the entire study area from a classified Random Forest classifier. The land-use type is required so that the correct carbon model is used for simulation: if the area is forested, the CBM is used, while the DNDC model is used for agricultural land. To increase the accuracy of the classification, we employed a region-based strategy in collecting training sites and testing sites: The land-cover classes produced for training the Random Forest classifier were selected from one region of the Clay Belt, while the points selected for a testing set was taken from another part of the Clay Belt [29]. Thus, the land use and land cover were derived.

The carbon models in this study require historical information such as original state of the land (or initial state from the beginning of the period to be simulated) to be defined. This historical information was obtained from the land-use data sources such as listed above, as well as Ontario’s T1 Forest Resource Inventory (FRI) for 2008–2017, which provides the inventory for all the forests in the province [32]. Furthermore, land-use changes such as deforestation (land use changing from forest to non-forested), and afforestation or reforestation (land going from non-forested to forested) need to be input in the models. Whenever there was afforestation or reforestation, the forest stand growth was also supplied, with the growth and yield curves for the trees in the area under study. The curves were obtained from various of sources, e.g., Plonski’s Normal Yield tables.

Now, the classes used for the study area were urban, water, forested, pasture, wetland, and cropland. The spatiotemporal dynamics of the forested area were obtained from classification maps, which showed if or how the forest composition changed over time, while the agricultural land changes were mainly obtained from visual observation of maps. For example, the Google Earth Maps at 30 m resolution and the PlanetScope [36] maps at ~3 m resolution were used to establish when there had been a land-use change, e.g., a deforestation for conversion to cropland. A landcover classification of the study area was performed using a Random Forest classifier. The input data consisted of a fusion of Landsat and Sentinel-2 images. For the classification, the training data were selected from one part of the study area, while the validation data were from another part of the study region. That was done to reduce spatial correlation between the training and testing data and, thus, bias in the classifier. The overall accuracy of classification using the multisource data was 90.97%. A complete description of the methodology of determining the land-use types and land-use changes is in [29].

b.

Environmental data: Climate data are a vital factor to the accurate output of the carbon models. The Dead Organic Matter (DOM) carbon and decay rate estimated by the CBM model are dependent on the specified temperature. At each time step, the decay rate is calculated based on the mean annual temperature input by the user. Although the CBM simulates the effects of temperature changes on decomposition rates, it does not consider the effects of precipitation on decomposition [17]. The mean annual precipitation is not used at all for the simulation results in CBM currently. In the DNDC model, daily weather data are used to calculate the soil climate profiles (for the soil temperature, soil moisture, oxygen concentration, etc.). The model uses climate data, as well as water and nitrogen demand and uptake, to estimate the plant growth. Daily meteorological data are entered for each year being simulated, and the user can choose whether to use minimum and maximum daily temperature or mean temperature for the temperature inputs. For this study, the daily mean temperature and mean precipitation were used. The windspeed, radiation, and daily humidity values can also be input as part of the climate profile. However, if any of these latter three parameters are used as part of the model, the daily minimum and maximum temperature values also need to be input; the mean temperature is not an option then.

c.

Disturbance events and transitions: A disturbance can be characterized as an interruption to vegetation. The disturbances that occur in the area could be brought on by anthropogenic means, as well as by natural causes, such as burns, weather events, forestry harvesting, infrastructure, and pests. In the CBM model, the disturbances can be simulated to occur over the entire study area or can be limited to certain stands. Examples of disturbances that have occurred in the Clay Belt are wildfires, logging, disease damage, and infrastructure disturbance such as from road construction. These disturbances produce corresponding changes in soil carbon [30]. Historical land disturbances (such as logging, fires, road/infrastructure creation, and insect damage) over the period of study were obtained from various governmental datasets [37,38], as well as from the farmland owners.

In the carbon models, it is possible to simulate the transitions which occur following a disturbance. For instance, if a stand consisted of mixed species such as trembling aspen, eastern larch, and black spruce, then following a harvest disturbance, the stand could now be only eastern larch and black spruce, or even eastern larch, black spruce, and balsam fir. The post-disturbance dynamics of the stand would depend on the growth and yield curves applied in the model. After a stand-replacing disturbance, the CBM model resets the age of a stand to 0. The stand will then begin to grow according to the same growth curve if the same species of tree continues growing in that place.

The yield tables input by the user are used to create the growth and yield curves in the model. The growth and yield tables permit the calculation of other parameters such as merchantable wood volume of each species, as well as the aboveground biomass components like stemwood, branches, and foliage. These tree growth curves were mainly obtained from provincial sources and some from the Plonski Yield tables, as the growth curves vary by province. For example, the growth curve for red pine in Alberta, a province in western Canada, could be different from that in Ontario, which is in central Canada. A few of the growth curves used were generic to the tree group, e.g., deciduous or coniferous. Figure 2 below shows examples of the growth and yield curves we created in CBM, having input the volume of the tree in m³/ha for each growth stage. The CBM model has volume-to-biomass conversion factors within it [17], which are used to estimate the resulting biomass carbon.

In a previous work, competition was simulated in the carbon models by adding a growth retardation factor in areas of mixed stands [30]. Thus, the effect of having different species in the stand as opposed to a uniform stand could be examined.

In this work, transitions have been modeled for changes in landcover class, e.g., forest harvest, forest management (for instance, pruning and cutting; removing dead trees), wildfires, disease damage, insect damage, and agricultural expansion. The input parameters used in our study were derived from the field study, standard practice in the study area, or the remote sensing data and geospatial databases on government websites. Hence, the input data are as close as possible to what could be obtained in the region of study.

Disturbances modeled in this study include the following:

• Land-use/land-cover changes: Land-cover change was tested annually via the spectral angle method. If a change in land-cover class was found, the simulation was updated with the new cover type. In modeling for the future, the historical trend was used to forecast the transitions until 2043. Similarly, an assumption was made about the area being allocated to cropland, pastureland, and to livestock. The agricultural expansion sequence was based on common farming practice in Ontario, as well as the information given by farmers in the Clay Belt on how they operate their farms. Urbanization was not considered in this work, as it is unlikely that there will be substantial development in the NLP within the coming years: All the NLP areas are already accessible by local and major roads and are served by existing infrastructure. Thus, an urbanization transition was not investigated in this study.
• Forest cutting and harvest: The DNDC model allows the user to input the yearly farming management information regarding grazing and grass cutting: for example, the number of hours each day the livestock grazes, the number of days in a month of grazing, and the intensity (in heads/ha) of the specified grazing animal are some parameters the user can input. Similarly, the number of times the grass is cut per month, as well as the parts of the plant cut, should be defined—root, stem, leaf, or grain. The harvest and cutting activities performed in the forested land was also simulated. Some of the events were clearcut harvesting of 60% of merchantable trees followed by salvage of the snags or clearcut harvesting of 30% of merchantable trees followed by burning organic residue (slash). Those events were stand-replacing, and the disturbance transitions were defined as such.

The Ontario government geodatasets provide information on the plans for commercial harvesting on Crown land [38]. These forest management plans are defined in two phases of five years. Following the completion of Phase I, each forest project is examined with an optional refresh of objectives or activities. During this period, a public consultation is conducted. Now, where trees that were allocated for harvest in Phase I were not logged, they can be incorporated into the Phase II plans if necessary. However, areas designated to be cut in the second phase cannot be harvested before that, i.e., in Phase I. Using the historical forest management plans in concurrence with maps produced for landcover change assessment, the deforestation activities were simulated in the CBM model. Also simulated were the different events to start farmland that was converted from forested land: e.g., clearcut with slash–burn, as well as clearcut harvesting with salvage of dead stem snags from the stand. The Progressive Adaptive Method (PAM) [29] was adopted to establish the disturbances that most likely occurred based on both the forest management plans and the remote sensing data from the Landsat and Sentinel satellites. This feedback loop between what is proposed and what actually occurs regarding forest disturbances by way of harvest is a key component of the PAM method. Thus, the maps from the provincial forestry unit were closely reviewed.

iii.

Fires: Wildfires are a significant natural disturbance in Canada’s forests. They can arise from such sources as brush burns, lightning strikes, campfires, fireworks, railway locomotives, and powerline shorts. Historically, wildfires in Ontario have burned for a variety of periods—from a few hours up to several weeks. On average, wildland fires consume 2.5 million ha. each year in Canada, roughly half the size of the province of Nova Scotia [39]. They may typically burn more of the understory and midstory—though this still depends on how severe the fire is. After the fire, the understory of the forests could be open to light, permitting subsequent growth in the area. The transitions of the forest composition following each wildfire were simulated in the CBM model.

The Ontario Wildfire Archive [40], as well as the provincial government’s database Ontario GeoHub [41], contain historical records of the fires that have occurred in the province since 1960. Thus, over the study period (from two decades ago until two decades ahead), the trend of historical fires in the region was used to simulate future wildfire patterns. Figure 3b illustrates the history of fires in the Clay Belt in the past two decades. In addition, a future effect of climate change could be an increase in extreme weather events and wildfires. Hence, the climate change impact was introduced to the model in terms of the frequency and severity of the wildfires after year 5. The different disturbances have varying impacts on carbon pools and transfers. For instance, wildfires could transfer carbon to the atmosphere (as they consume debris on the forest floor and subsequently release carbon to the atmosphere), as well as to stands and to pools like the DOM pool. However, forest harvesting transfers carbon to the forest product sector or to final wood products and stands. The annual emissions from Canada’s wildfires have been known to equal the annual emissions from burning fossil fuels [42].

Other fires that can be simulated are the regular burns done on cropland. Farmers sometimes burn leftover straw after a harvest, as fields are prepared for seeding the following spring. Controlled burning is easier and more cost-effective than tilling the straw. Furthermore, in areas where the soil has a high clay content, the soil is more likely to have compaction and draining problems. Thus, it can be difficult to till straw/residue back into clayey soil. Therefore, burning would be a choice a farmer could take to manage the residue [43]. Apart from clearing the land for planting, the burning of agricultural fields is also done in Canada to control pests, disease, and weeds [44].

iv.

Disease damage: The provincial geodatabase provides the spatial extent of diseases that have afflicted the forests since the early 2000s. Thus, past forest disease damage events were simulated in CBM, and predictions for future events were made based on the trend observed. Some of the diseases that have been predominant in the Ontario forests include brown spot needle blight, ink spot of aspen, spruce needle rust, and septoria leaf spot and canker [45]. Particularly in the Clay Belt, Septoria leaf spot and ink spot of aspen were the diseases that afflicted the forests, with the extent of the damage listed as moderate to severe. (See Figure 3c).

v.

Insect damage: Insect damage within the NLP sites has come from the forest tent caterpillar, the larch casebearer, and birch skeletonizer during the study period. Again, the effects of the damage were moderate to severe. Figure 3d illustrates the insect damage experienced by vegetation in the Clay Belt in the past two decades. In the simulations conducted in CBM, salvage logging was only simulated after significant insect mortality. Significant insect infestation was taken as a stand-replacing disturbance. Partial stand mortality could also be simulated as a lesser effect of insect damage.

2.3. Simulation Framework

The DNDC model was run in the Site mode rather than Regional mode; the input parameters were entered individually for each site—five pasture sites and five cropland sites. The crops chosen for the simulation were corn, wheat, oats, and soybean, i.e., crops typically grown in the Clay Belt. The crop rotation was such that spring wheat was planted early in the year and harvested in the summer, while winter wheat was planted late in the year and harvested the following year. A fallow year was also implemented every five years. Fertilizer application was simulated, along with the tillage routine. Inputs to the model included the following:

• Latitude;
• Climate information;
• Nitrogen concentration in rainfall;
• The clay fraction in the soil;
• The soil bulk density;
• Soil pH;
• Initial nitrogen concentration at surface soil.

As per the analysis of samples from field studies, the soil nitrogen and carbon amounts were known.

The lab analysis of the soil samples gathered during the field survey provided both the pH(H₂O) and the pH(CaCl₂). However, for this study, the pH(CaCl₂) was used, as it is taken to be more accurate than the pH based on the water scale [46]. Both pH measurements were taken by mixing an air-dried soil sample with four times its weight of the fluid—water or calcium chloride—prior to measuring the pH value. The differences in pH(w) vs. pH(CaCl₂) were not always large; for instance, in site F5, the pH values were 5.865 vs. 5.593, respectively, and at the P1 site, the pH(w) was 8.070, while the pH(CaCl₂) was 7.617. Although the differences are not always large in value, the pH(CaCl₂) test is held to be more accurate, as it reflects what a plant experiences in the soil [46]. In the DNDC model, the soil textural class was required. The different classes encountered in the sites were loam, clay, sandy loam, clay loam, and sandy loam. The CBM model required the specific type of soil present in a site following conversion from forest. The Canadian soil map (REF) was used to determine soil types. The simulated areas were either of the Gleysolic, Luvisolic, or Mesiosol soil type. However, in the absence of Mesiosol as an option in the simulator, Gleysolic soil was used.

When simulating the pastureland, perennial grasses (like timothy grass) and alfalfa were selected. The study area underwent grass cutting once a year to maintain uniform pasture and promote productivity of the grasses. The pasture was also grazed from the beginning of May until the end of October each year, which is the timing of the pasture grazing in accordance with what is usually established in Ontario [47].

Another situation proposed was the forest being converted to a livestock farm. Beef farming is already widely carried out in the Clay Belt [5]. As mentioned before, an initiative for “Beef Cattle Farming Training Program for Indigenous people” was conducted in 2017 shortly after the Ontario provincial government announced the decision to convert part of the Clay Belt to agricultural land [6]. For this aspect of the simulation, it was proposed that there would be growth in the livestock, e.g., five head of cattle initially on the farm, and after five years, there would be five more head of cattle. Then, by the culmination of the twenty years of the simulation, there would be 20 head of cattle. It was simulated that the cattle would graze in the area for 10 h a day from May 1 to October 30 each year, and the grass would be cut once a year. Cattle grazing (dairy and beef cattle) is usually done in the Clay Belt [48]. Goats (simulated as sheep in the DNDC program) and pigs were simulated as the livestock in certain sites, as shown in

Table 1. Properties of the 10 sites simulated on the DNDC simulator.

Name of Site	Land Use	Plant	Planting // Harvest Dates	Texture	Clay %	pH	SOC at Surface (0–15 cm)
P1	Pasture: Dairy cattle farm	Alfalfa	Perennial	Clay	59.70	5.6	0.318811
P2	Pasture: Goat farm	Timothy grass	Perennial	Clay	47.27	6. 7	0.846974
P3	Pasture: Beef cattle farm	Alfalfa	Perennial	Clay loam	35.6	6.6	0.54909
P4	Pasture: Goat farm	Bromegrass	Perennial	Clay loam	54.6	6.2	1.147014
P5	Pasture: Pig farm	Bromegrass	Perennial	Clay	47.16	6.8	1.129642
F1	Crop Farm field	Corn	May // October to November	Clay	49.58	7.0	0.424732
F2	Crop Farm field	Spring wheat	Early April to mid-May // August	Sandy clay loam	13.93	5.9	0.110473
F3	Crop Farm field	Oat	Late April to early May // August	Sandy loam	16.91	6.5	0.457353
F4	Crop Farm field	Winter Wheat	Mid-September to mid-October // Late July or early August	Sandy clay loam	24.54	5.2	0.851421
F5	Crop Farm field	Soybean	Mid-May to early June // October	Clay loam	35.82	7. 1	0.876171

. There is provision in the model to have external manure applications simulated. However, on the pastureland, we only simulated the excreta from the grazing animals on the fields, and this manure was spread intermittently to prevent damage to the pasture. (If not done, the grass underneath the dung could deteriorate or be killed by the high nutrient amounts). Furthermore, the hooves of the animals can assist in breaking up and spreading the dung across the pasture. The speed of the dung disintegrating also depends on the climate [49]. For this scenario, there was an addition of farmyard manure, simulated consisting of about 86% carbon for each kg/ha. However, an external manure application was modeled for the cropland areas.

Another important parameter to model was the irrigation. Advances in technology, such as tiled irrigation, have played an important role in making the Clay Belt region more economically viable for crop production. An irrigation index was specified as 100% such that every water deficit is met by irrigation water [26]. The water budget in DNDC has input fluxes of precipitation, irrigation, and groundwater supply.

The crops that were simulated as being grown when the forest is converted to farmland, with their attendant parameters, are listed in Table 1 below.

Finally, the option of the forest not being converted was simulated, and the resulting carbon and GHG stocks were analyzed using the CBM tool. For these areas, the species of landcover in the NLP area was used as the input to the model. For example, if the area is not treed, the nearby forest species as taken from the historical land-cover maps are used as input to the model. The patterns of insects and disease were simulated as earlier discussed.

Table 2. Input parameters simulating carbon in sites and nearby forest for baseline comparison.

Site Name	Nearby Forest	Tree Species	Disturbance History	Conversion Year	Harvest in Past Two Decades?
P1	Coniferous	White spruce	Forest tent caterpillar infestation in 2016	No change since 2002	No
P2	Mixed	Trembling aspen, Larch, Black spruce	Forest tent caterpillar infestation in 2016	2009–2010	Yes
P3	Mixed	Larch, Black spruce	Fire in 2003, forest tent caterpillar infestation in 2015	No change since 2002	No
P4	Mixed	Balsam poplar, Eastern larch, White spruce	Forest tent caterpillar infestation in 2016	Land-cover change between 2002 and 2009	Yes
P5	Deciduous	Larch	Fire in 2003, forest tent caterpillar infestation in 2016.	Same land use since 2002; only grown smaller in area	Yes
F1	Deciduous	Trembling aspen	Fire in 2005, forest tent caterpillar infestation in 2016	2017	No
F2	Coniferous	Black spruce	Forest tent caterpillar infestation in 2015. Deforestation, then left fallow; mulched in 2017	2017	Yes
F3	Deciduous	Trembling aspen	Fire in 2003, forest tent caterpillar infestation in 2016. Burnt at deforestation in 2017. Left fallow; not mulched.	2017	Yes
F4	Mixed	Trembling aspen, White spruce, Balsam fir	Forest tent caterpillar infestation in 2016	No change since 2002	No
F5	Mixed	Trembling aspen, Balsam poplar, White birch	Spruce budworm insect infestation in 2021	No change since 2002	Yes

contains the properties of the forested sites, which were simulated in CBM.

2.4. Field Survey Data for the Analysis of Simulation Results

We embarked on fieldwork in 2018, 2019, and 2021 within the Clay Belt to collect soil and vegetation data. The field survey was conducted to support the research of digital soil mapping and to assess the impact of land conversion on soil properties and GHG emission. This data collected were input data for some of the parameters the carbon models required. Furthermore, as discussed in the next section, the field data were compared with the carbon model results. Sample sites were selected at various locations within the Clay Belt in order to have a mix of landcover types—e.g., farmlands, woodlots, forested, abandoned pasture, etc. The understory and overstory observed in the forest areas were always recorded for the plots where sampling was conducted. Figure 4 is an example of some of the sites visited for field study.

The coordinates of the locations sampled were recorded using a high accuracy Leica GPS receiver with a Total Station. In order to determine the locations of the plots accurately, the base station was set up (wherever and whenever possible) in an open area near the access points and allowed to run for about 1.5 to 2 h, during which time the soil sampling was completed.

At almost every sample site, very dense clay was encountered. All the soil samples collected were sent to a lab for analysis. The lab results derived the soil carbon content and bulk density data.

Given the bulk density (in g/cm³) and percentage of carbon in each transect from the soil analysis, the carbon amount was calculated from Equation (1). Within a specific volume of soil collected, we have the following calculation:

$$\begin{array}{c}\mathrm{Carbon} \mathrm{amount} = \mathrm{Percentage} \mathrm{of} \mathrm{carbon} \times \mathrm{Mass} \mathrm{of} \mathrm{soil}\hfill \\ \mathrm{or}\hfill \\ \mathrm{Mass} \mathrm{of} \mathrm{carbon}/\mathrm{area} = \mathrm{Organic} \mathrm{carbon} \% \times \mathrm{Bulk} \mathrm{density} \times \mathrm{horizon} \mathrm{thickness}\hfill \end{array}$$

(1)

Similarly, the nitrogen amount was estimated, with the percentage of nitrogen in each sample provided by the lab analysis results. Thus, the field data served as validation to the simulated carbon and nitrogen content where possible.

2.5. Uncertainty Analysis Framework

One of the methods to calculate uncertainty in carbon models is to select several parameters related to the factor one is investigating, such as soil organic carbon. The parameters are perturbed within a range of variation, and their effect on the factor can be estimated. White et al. [28] used the Gaussian Emulation Machine for Sensitivity Analysis (GEM-SA) software v.1.1 [50] to calculate the uncertainty of dead organic matter. GEM-SA permits the user to build an emulator based on Bayesian statistics. GEM-SA proposes to perform uncertainty and sensitivity analyses with fewer code runs than Monte Carlo-based methods [50]. Metsaranta et al. calculated uncertainty by ratioing individual uncertainties to the total uncertainty: running the model over the entire study area of 2.3 × 10⁶ sq. km five times while varying one factor each time. Five factors were assessed for uncertainty in that study [18].

As a way to combat climate change due to atmospheric carbon, the Paris Agreement recommends maintaining or increasing carbon sequestration and storage in forests. Yet, under a changing climate, there could be challenges to maintaining forests as carbon sinks. Minor temperature increases (for instance, between 1 °C and 2 °C) would increase the rate of forest growth and productivity. Hence, there would be more carbon sequestration in the forests. However, temperature increases of more than 2 °C can trigger a decline in growth due to higher mortality [51]. Therefore, rising or reducing temperature could be used to examine the effect of climate change on carbon sequestration.

In that vein, for the DNDC model, the temperature and precipitation were the climate parameters modified in the Monte Carlo uncertainty run. The Monte Carlo test can be performed on different categories of variables like Crop parameters (such as maximum yield, water requirements, plant carbon/nitrogen ratio), Soil parameters (e.g., pH, porosity, bulk density), and Grazing management (such as grazing days and grass cutting amount). With the climate variables chosen to quantify the uncertainty of the analysis, a range was specified for the variables—for air temperature, it was 1.1, and for atmospheric pressure, it was 0.9. Thus, for a 5-year cycle, the Monte Carlo test was run 4000 times for each time step while varying both parameters. For example, the pressure was run for values ranging from −0.9 to 0.9, and the temperatures ranged from −1.1 to 1.1.

In the forested areas, the quantification of sensitivity was related to the quality of input data: we ran the model with disturbance data and landcover information that was not correct, and we compared the simulation results to the original results obtained from correct data. The difference between the two results could be indicative of the importance of such input data when monitoring carbon dynamics. It could also be used to measure uncertainty in the modeling from the input data.

3. Results

This section presents the ecosystem carbon and GHG emissions simulated using the different carbon models. Whereas the DNDC model tracks the flux of GHGs such as CO₂, CH₄, NH₃, NO, and N₂O, the CBM model captures the flux of only CO, CO₂, and CH₄ emissions. Thus, the comparisons and analyses of GHG emissions being performed in this work will only focus on the GHGs that overlap between the two models.

3.1. Comparison of Field Survey to Simulated Ecosystem Carbon—To Present

As mentioned earlier, field survey data were acquired, which served as validation for the simulated carbon from both models. Our simulations, which utilized historical data to predict the current carbon stock, have proven to be a reliable method for generating realistic results. Thus, we can also predict the future carbon stocks with confidence.

We highlight the comparisons in the carbon estimated from some of the sites in Figure 5, taking the pasturelands and farmlands as examples. The farm labelled #3 had its conversion in 2017 according to the high-resolution satellite imagery examined. That corresponded with what we observed when we visited the site: it was a clearcut field, ready to seed. The pasture site #4 was an oat field when we visited for our field campaign. Obtaining some local knowledge to combine with our remote sensing techniques—having the area’s history—allowed us to better simulate the processes that occurred, leading up to the amount of carbon stock that we estimated at the current time.

During our field visit, the pasture site #2 was a mix of different grasses. It was a livestock pasture, and we observed about 30 heads of livestock grazing in the area. Performing the comparison of the calculated amount of soil carbon from our field studies with the simulations done, in all the cases, the results were comparable within standard deviation except at pasture site #2. The discrepancy might have arisen because, at that location, we did not ask the farmer how often he did his grass cutting nor confirm the livestock information. The difference in soil carbon estimated at that site from the simulation vs. the amount determined by soil analysis underscores the importance of obtaining sufficient and correct data when performing simulations. Furthermore, having some local knowledge to combine with our remote sensing techniques (to learn about the history of the study area) allowed us to better simulate the processes leading up to the current carbon and GHG stock estimated. This, again, is a principle of our PAM method—to use the best available data in generating a depiction of the past and a prediction for future processes.

3.2. Simulations Results of Ecosystem Carbon—Forested Sites

First, the simulation result of the historical ecosystem carbon stocks and GHG emissions for the sites is presented. All the sites had undergone conversion by 2017. Figure 6 shows the future prediction for about two decades, with the disturbances expected based on historical trends. (Table 1 gives the descriptions of the sites labeled F1–F5 and P1–P5).

If no land conversion had taken place, the pattern of ecosystem carbon and dead organic matter accumulated over the four decades of simulation would be as shown in Figure 7 and Figure 8 below. For instance, the total ecosystem carbon came out to 21,769 t/ha in F4, 24,264 t/ha in P4, and 7656 t/ha in the site adjacent to P5.

The dead organic matter always shows a change following a disturbance. For instance, in the forest annexing site P2, when there was a wildfire in 2003, there was an increase in the dead organic matter carbon from 315 to 400 tC/ha. In addition, after insect damage, there was an increase in DOM carbon from 342 to 352 tC/ha. Also noted is the impact of different disturbances on the forested areas. For example, forest sites near P2, P4, and F2 experienced harvest and insect damage, whereas sites by F4 and P1 experienced only insect damage. The slopes of DOM carbon change for P2 and P4 were −11.5 and −9.3, respectively, while the latter group was −5.4 and −1.8 following the initial disturbances. Similarly, when there was fire and insect damage as in sites F1 and P3, the slope of carbon change was between −4.0 and −6.1.

Comparing the baseline ecosystem (the nearby natural forest) to the total ecosystem carbon in areas that are being used for agricultural purposes, Figure 9 below illustrates the difference in amounts. It shows that, for instance, on site P3, there is a −46.5% difference between the baseline and the pastureland, and on F1, the difference is −7.6% between the baseline forest and the farmland; at site F4, it is −64.9%. Similarly, the difference in GHG emissions for these scenarios came out to −74.8%, −44.9%, and −97.7%, respectively in CO₂ amounts, as well as −99.9%, +0.05%, and +100% in CH₄ amounts. The GHG values on field P3 are very high because there was a fire in that area in 2003. The CH₄ amount simulated in all the years after was negligible. Site F1 was converted to a farm field in 2017, and the total CH₄ emissions thereafter were 0.02 tC/ha and 41.42 tC/ha of CO₂ in the next 27 years under investigation. However, F4 was a wheat field for the entire duration of the simulation compared to a mixedwood forest that never underwent fire nor harvest. Thus, the CO₂ emitted from the forest was ~98% higher than what was produced on the farm field, while the CH₄ emissions were negligible in both scenarios (farm or forest), but the farm had a total of 0.03 tC/ha over the four decades to the forest’s roughly nil CH₄ emissions. The carbon sequestered by the forest is also ~65% more than that on the farm.

In Figure 10 below, we illustrate the carbon sources in the forest stands. (The figure shows the amount of carbon simulated when the forest remains unconverted). Prior to disturbances, the carbon stocks are more in the biomass than in the DOM. Following disturbances, the biggest contributors to the carbon stocks are the litter and soil carbon in the mixed and deciduous forests. For the coniferous forest, the dead organic matter and soil carbon were the largest sources of ecosystem carbon following the disturbances. Three sites are shown in Figure 10 below, with a conifer forest, a deciduous forest, and a mixed forest. Table 2 gives the descriptions of the forested sites depicted in Figure 10a–c.

3.3. Simulation Results of Carbon Stock and GHG Emissions—Agricultural Sites

Considering the agricultural lands, simulations were performed to determine the sequestered soil organic carbon and the GHG emissions if the land usage was persistently agricultural for four decades. Figure 11 below illustrates how the soil carbon and CO₂ amounts varied over the years. It is observed that there was always a drop in soil carbon and GHG emission values during the sixth year when the croplands were left fallow. Furthermore, sites P2 and P4 had the same animals grazing over the simulation period, with the only difference being the type of grass planted on the sites. Their difference in soil organic carbon was, on average, a difference of 3.85%/year, i.e., a total of 27,098 kgC/ha over the four decades simulated. Similarly, P1 and P3 both had cattle being grazed each year, with identical growth patterns of the herds. The grass type was alfalfa in both areas, and the difference in the sites was that the animals being raised were dairy cattle on site P1 and beef cattle on site P3. The total GHG emissions in the former site was 633,652 kgC/ha of CO₂ and 0.00985 kgC/ha of CH₄, while the beef cattle’s total emissions over the simulation period was 620,348 kgC/ha of CO₂ and 0.01298 kgC/ha of CH₄. (It should also be noted that the soil properties were different at the sites, e.g., different pH, clay percentage, and SOC amount at the surface).

Assessing the baseline forest simulations, the mixedwood forests have considerably more ecosystem carbon than the forested areas consisting of single species: whereas the uniform-species areas have between 170 and 431 tons of carbon, the mixedwood forests range in ecosystem carbon from 291 to 881 t of carbon. In fact, over the duration of the period simulated, apart from the forest by site P3, the lowest amount of ecosystem carbon is a little less than 450 t of carbon, where the highest amounts of ecosystem carbon in the single-species sites was 431 t during the period under study.

The impact of the changing climate was also examined. This was conducted by introducing a multiplication factor on some climate variables in the DNDC model. The carbon sequestered on the corn and soybean farm fields is shown in Figure 12a,b (where the stocks simulated under the changed climate are represented as the climate-impacted part of the figures). The soybean soil carbon stock was mostly lower under the changed climate (average of −1.0%), while the soil carbon in the corn fields was −2.0% different on average. However, the CO₂ content on the soybean farm was mostly greater when the climate was changed (an average of +10.65% increase), while the corn fields differed in that respect (decrease of −0.55% CO₂ on average). Figure 12 illustrates the impact of the changed climate compared with the current climate scenarios.

The pasturelands also saw differences in soil organic carbon and GHG emissions under climate change. Site P5, which has livestock and where farmyard manure application was simulated, showed significantly less soil carbon under climate change during the simulation period. Yet for site P3, the SOC would be higher under the new climate regimes. Overall, the difference on average came out to +0.33%. The CO₂ emissions, however, were on average more for P5 under the changed climate conditions and less for P3, as illustrated in Figure 11c,d. For the 41-year study period, the difference in the GHG emissions was found to be less on site P3 under climate change by an average of −25.6%.

3.4. Sensitivity and Uncertainty Analysis

The average uncertainties of carbon stock estimated on the agricultural lands were not very high, as seen in

Table 3. Average uncertainty of carbon and GHG (CO₂) stock estimates for the different agricultural lands.

Crop	Soil Organic Carbon Uncertainty	GHG Uncertainty	Livestock Grazed	Soil Organic Carbon Uncertainty	GHG Uncertainty
Corn	0.022%	0.553%	Beef cattle	0.332%	25.559%
Spring wheat	0.0073%	5.357%	Goats	1.549%	16.157%
Oat	0.129%	0.489%	Goats	0.860%	11.011%
Winter wheat	0.0169%	0.489%	Pigs	1.073%	15.506%
Soybean	0.102%	10.651%	Dairy cattle	0.478%	42.481%

, while those of the CO₂ emissions were much higher. Additionally, the ranges of CO₂ emissions uncertainties reveal great fluctuations that were not seen in the carbon stock uncertainties. The disparity in average uncertainty might arise because the GHG amounts are smaller than the carbon stock values, so there is a greater impact of differences when the comparisons are made between the original and the result from the uncertainty analysis. That might also play a role in the ranges of uncertainty; for instance, the spring wheat field’s CO₂ emission’s uncertainty ranged from −17.57 to +22.52, and the pig field’s CO₂ emission’s uncertainty ranged from −7.99 to +28.51. Table 3 lists the outcome of the uncertainty analysis in estimating SOC stock and CO₂ emissions.

Another test was done for the forested areas using the CBM model. We sought to determine whether there would be significant differences when approximated data inputs were used in the simulation as compared to the more accurate data obtained from remote sensing tools. Figure 13 below shows some of the results obtained.

Table 4. Comparing carbon and GHG amounts simulated under approximate and more accurate inputs.

Stock	Approximate (t)	Accurate (t)
Total Ecosystem Carbon	360,068	304,919
Soil Carbon	156,413	118,696
GHG (CO, CO2, CH4)	6228	4245

is a summary of the simulations when approximate values using generic data were input to the carbon models instead of the more accurate data from classification maps obtained from remotely sensed data. Comparing the carbon stock and GHG amounts estimated by the carbon models when approximate values from generic data were used as input to the carbon model as opposed to the more accurate data from classification maps obtained from remotely sensed data, a pattern was revealed: the more accurate data input yielded lower amounts. That is, the total ecosystem carbon was overestimated by 15.32%, the SOC by 24.11%, and the GHG emissions by 31.84%.

More specifically, we also used four sites to demonstrate the effects of accurate input data by comparing the SOC and total ecosystem carbon amounts when some of the disturbance data were omitted:

• On site by F1—simulation performed without disturbances;
• On site by P5—simulation performed without disturbances;
• On site by F4—only one species (white spruce) was simulated as being present;
• On site by P4—simulation performed without land use conversion (whereas it should have been simulated as occurring in year 80).

In Figure 14, the total ecosystem carbon from the forests by the farm and pasture sites F1, F4, P4, and P5 are plotted vs. the estimated ecosystem carbon without the disturbances or correct landcover information. (The simulations done with incorrect information are labeled as the ‘modified’ sites in Figure 14). When the disturbances were not included, the estimated ecosystem carbon amounts were higher than the actual values—by 12.5% by site P5 and by 14.7% by site F1. This is expected, since the disturbances such as insect damage and wildfires create more dead organic matter, thus influencing the total ecosystem carbon amounts. On the other hand, having the forest by site F4 labeled as coniferous with only one tree species, rather than a mixed forest, there was a 58.6% reduction in the total ecosystem carbon over the period simulated. Similarly, when the deforestation by site F4 was not included in the simulation, the estimated carbon amount was incorrect by about 26%. The SOC comparisons showed a similar trend, with sites by P5 and F1 overestimating the SOC by 2.7% and 2.5%, respectively. Such results from the uncertainty/sensitivity analysis underscore the importance of having the right input for the models. It is vital to provide realistic data, from the best available knowledge sources, as input for simulation to obtain accurate ecosystem carbon stock and emissions estimates.

4. Discussion

Ontario’s boreal forest is said to be a fire-dependent ecosystem: after a wildfire, new vegetation grows [52]. Wildfires can produce enormous amounts of carbon and greenhouse gases. According to the principle of uniformitarianism, future wildfires are predicted based on their past occurrence, as tracked and recorded by governmental databases. In some years, there could be anomalies, such as in 2023 when Ontario experienced a severe wildfire season.

This study took the baseline case as a situation where the nearby natural forests are not converted to agricultural lands. Comparing the sequestered carbon simulated in those scenarios shows the forests to have resulted in 1.96 times more sequestered carbon in total. Furthermore, the forests continually acted as a carbon sink, not a source, over the four decades simulated, as demonstrated in Figure 10.

Depending on their age, older forests are usually better at sequestering carbon—old forests have more stored carbon. However, they sequester at a slower rate than young forests. Most of the deforestation in Canada occurs when forests are converted to other land uses such as agriculture, resource extraction, industrial development, and urban expansion [10,16]. Thus, emissions from forest conversion between 2005 and 2021 have fluctuated around 16 Mt [10].

It was observed that when the forests were converted to cropland, the average GHG emissions were 60.14 tC/ha CO₂ and 0.781 kgC/ha CH₄. Furthermore, the average ecosystem carbon sequestered in these scenarios was 9252.9 tC/ha. (The average ecosystem carbon simulated in the forested areas was 16,616 tC/ha). Amongst the crops examined in this study, soybean showed the most ecosystem carbon sequestered, while winter wheat had the least. Other studies showed that conversion to cultivated land from forests similarly reduced the soil organic carbon [53,54] and in pasturelands used for grazing [55].

In addressing the GHG emissions that the country will experience in the coming decades, it is essential to examine the input from methane: CH₄ has 28 times the heat-trapping ability of each molecule of CO₂. Interestingly, roughly 25% of the global CH₄ emissions are from cattle [56]. Goats and cattle have high methane emissions because they are ruminants—they produce a lot of methane while digesting food. Another result in this work was that the beef cattle produced 1.32 times more CH₄ emissions than the dairy cattle. Other studies have shown that beef cattle produced 6.9 times the direct CH₄ emissions of dairy cattle [56]; that beef cattle produced 2.9 times the dairy cattle CH₄ emissions [57]; and in 2012 in the U.S. that beef cattle produced 71% of the total enteric emissions (mainly CH₄), while dairy cattle produced about 25% of the country’s enteric emissions [58]. Now, dairy cattle emit less methane than beef cattle because lactating cows intake primarily for milk production, while beef cattle consume feed entirely for meat [59].

Pastureland was also simulated, with it being grazed for 10 h a day from May to October. Under these conditions, the average GHG emission was 303.71 tC/ha CO₂ and 0.334 kgC/ha CH4, with the sequestered carbon being, on average, 7708.7 tC/ha. The pasture field with the highest emissions was alfalfa pasture with dairy cattle grazing.

Although continuous grazing needs less management and incurs minimal costs, rotational grazing requires more management and costs for watering and fencing certain parts of the fields. However, the rotational grazing practice is more beneficial in terms of the quality of grasses in the pasture. Disadvantages of rotational grazing include the loss of an area for production each cycle and more work in managing the rotation system; however, the benefits of allowing the plants to rest, not losing forage from trampling, and having more efficient animal distribution and forage use would outweigh the costs. Intensively managed pasture grazing is also known to produce more meat or milk per unit of land [49], as well as control disease and pest disturbances to the livestock, since they are constantly being moved away from fresh waste [60].

There are also environmental benefits to rotational grazing: environmentally sensitive areas such as wetlands, or even pastures that are too wet or too dry, can be fenced off so as not to be disturbed by the livestock. In addition, areas that are near riverbanks should be grazed in the late summer (dry period) rather than in spring (rainy periods) to protect the banks from erosion [60]. Stable banks help maintain water quality and improve fish habitats. Finally, healthy pasturelands not only reduce the risk of soil erosion from snow, rain, and wind, but they are also known to be more efficient at storing soil carbon than annual crops [49,60].

The simulations indicate that there were higher CO₂ emissions when forests were converted to pasture than when converted to cropland. However, on average, croplands had higher amounts of SOC than pasturelands following their conversion from forests. In addition, different crops exhibited varying amounts of sequestered ecosystem carbon in the croplands, and the pasture sites also had different amounts of SOC and GHG emissions depending on the grass being grazed, the livestock present, and the management practices employed on the sites.

Although some of the results obtained in this study demonstrated higher SOC or GHG emissions on farm sites compared to pasture sites, it should not be assumed that the conversion to cropland will necessarily produce more carbon and GHG than a livestock farm. It is important to note that the crop rotation chosen for this study, as well as its management, have affected the outcome. In addition, periods of fallow were incorporated into the simulations. Not surprisingly, the emissions were reduced in the fallow years. Furthermore, there was less sequestration during the fallow years. If portions of the Clay Belt are converted to croplands, farmers might choose a more diversified agrifood production. That, in turn, will affect the soil carbon and GHG emissions that arise. Further research will examine different land management scenarios following initial conversion from forested land, e.g., a reversion to forest cover if agricultural land is abandoned. In addition, future work will focus on minimizing uncertainties in the land-cover classification and the change detection results of data fused from multiple satellite sensors. Efforts to generate sufficient and more accurate input data to carbon models will in turn reduce uncertainties and produce more accurate simulation results.

The uncertainty analysis in this study also served to examine the effects of climate change on the soil organic carbon or GHG emissions in the various land uses. In this study, the climate parameters used for the uncertainty analysis were the atmospheric temperature and pressure. The Monte Carlo test of sensitivity was conducted by modifying the pressure values by a multiplication factor of −0.9 to 0.9, and the temperature values had a multiplication factor ranging from −1.1 to 1.1 of the original values. For instance, when the temperature was 0.56 times the current air temperature specified in the input file, the pressure was 0.39 times the original. Or, when the temperature was at the −1.01 factor, the pressure was at −0.83 of the original pressure value. The outputs of soil organic carbon and CO₂ emissions were then compared to the original values.

To continue demonstrating that the methodology is easily accessible to an everyday user, the computation times of the simulations were recorded. Processing times were compared for the activity that took the longest time—the Monte Carlo uncertainty runs. The simulations for uncertainty analysis were performed on three different desktop computers. PCs 1 and 2 have an i7-6700 Intel processor, with CPU @3.40 GHz; the only difference between them that is PC1 has an installed RAM of 16 GB, while PC2 has 8 GB of RAM installed. PC3 has an i7-8700K Intel Core processor, whose frequency is 3.70 GHz. The installed RAM is also 16 GB. All the systems have a 64-bit operating system. We conducted several Monte Carlo runs on each computer and noted their processing times. The average time for completing the run on a 5-year cycle was 5.88 h on PC1, 5.90 h on PC2, and 5.05 h on PC3. For the 1-year fallow runs, the runtimes were similar on the three computers, with processing times ranging from 57 min to 1 h and 33 min to complete. (A combined 23 fallow runs were tested across the three machines). Thus, we observe that a standard computer can perform the simulations in this work, though higher-powered computers give a slight difference in processing time.

Since climate variables were chosen for the uncertainty analysis, the runs were also used to project how the land would respond to an increase in temperature (at 1.1 times the historical temperature) or a reduction in pressure (a fraction of 0.9 the historical or projected future values). Thus, these runs simulated the impact of climate change with rising air temperature.

Several factors affected the simulation result. Typically, carbon stock predictions generate uncertainty from the carbon models and from the scenarios being simulated. Uncertainty from the carbon models has to do with their responses and the model inputs. One of the ways to reduce model uncertainty is to select models that have been well validated and shown to be robust under a variety of situations, such as the models used in this study. It is crucial that the model correctly simulates the processes in the ecosystem. So, in this study, different models were explored regarding suitability for the required task. The carbon models rely greatly on the inputs given by the user. Thus, as this study demonstrates, the simulations should be performed with the best or most reliable input data. For instance, fire is known to be a key driver of carbon dynamics and GHG fluxes in Canada’s managed forests [61]. Hence, it behooves the carbon model user to seek accurate input data for a model’s parameters and algorithms. In this study, we used remote sensing data to obtain realistic model inputs. The estimation of uncertainty within and between carbon models is still a complex issue; for instance, there is no consistent method for calculating confidence intervals in ecosystem models [18]. Yet, it is essential to make uncertainty estimates, as they will enable policymakers to know where to direct resources to reduce uncertainties in carbon estimates. Furthermore, the results of carbon models can affect the adaptations made by the policymakers under the changing climate. Hence, there is potentially a sizable economic value in reducing the models’ uncertainty when predicting future carbon and emissions.

Scenario uncertainty arises because the climate events of the future are unknown. Usually, the future is predicted from historical records. However, the simulations in this body of work were conducted knowing that unprecedented events could occur and cannot be properly modeled. For instance, in 2023, Canada’s eastern province of Nova Scotia experienced three months’ worth of rain in a day [1]. Future climate conditions were simulated to determine the potential emissions and ecosystem carbon stock under various conditions. (The possible disturbances that could occur, based on the disturbance history in the area, were also modeled). Using the SSP1.9 pathway from IPCC, which predicts the temperatures to be elevated above the most recent decade [2], future temperatures were increased by a factor to account for the effects of the changing climate. Similarly, the precipitation was modeled to be different in the future based on historical trends. Climate change is projected to cause changes in precipitation [52].

In this study, the estimated average uncertainty in emissions was 0.46% in carbon and 12.83% in CO₂ emissions from the agricultural lands. Other studies also vary in their uncertainties, e.g., ±15% for carbon [18] and ±19% for the uncertainty of emissions on cropland [10]. The carbon and GHG flux estimates typically have large uncertainties, especially for large areas, since modeling all the physical, chemical, and biological processes entailed could be difficult. For example, the national uncertainty range for GHG fluxes from non-converted forest land was ±39 MT of carbon dioxide (CO₂) between 1990 and 2021 (where non-CO₂ emissions do not contribute much to the total uncertainty) [10]; or, the uncertainty of emissions/removals on cropland was reported as ±19%; and the uncertainty in estimates of CH₄ emissions from enteric fermentation in 2020 was between −3 Mt and +5 Mt CO₂ eq. [10]. It is not unusual for the results to have uncertainties up to 20% of the mean value for carbon modeling [18,28]. Due to the various inputs and the complexities of the carbon models, uncertainties are bound to result.

Some limitations were observed in this study. Although we used reliable data sources for disturbances in the simulation, more comprehensive information could have been used in the models. For instance, when providing the occurrence of fires, the extent of the fire was not considered, even though this would impact the simulation results. Future work will consider not only the presence or absence of a disturbance but also the extent or severity of each disturbance that will be incorporated in simulations. Another area that could be expanded on from this study is the spatial response of the models: the carbon analysis conducted on the plots was broad, in a sense, because of the forest carbon model used—the version we used of the CBM model is aspatial. As such, the model can represent a stand of trees (or a group of stands), but it does not recognize the spatial relationships of stands within the spatial unit (a subdivision of the study area). Therefore, we were constrained to only looking at plots individually. In future work, we will use the new version of the CBM, which is spatially explicit. Then, a map of the emissions and carbon stock can be created.

Another limitation encountered was in the land practices explored. In the current study, only one practice of agriculture and forestry was considered. By specifying the forest being managed differently, the resulting carbon stock would also be different. Furthermore, the practice on the agricultural land after conversion is important in predicting the emissions. Additional studies will be conducted on various land management schemes to provide more insight into the ensuing soil carbon stock and GHG emissions.

5. Conclusions

This work has demonstrated how the accuracy of carbon estimates could be increased—by using remote sensing tools to generate input data for carbon models. For instance, for the carbon and GHG emission estimates in the study region, the total ecosystem carbon was overestimated by 15.32% when generic data (rather than more accurate input data) were used for simulations. Similarly, the generic data overestimated the SOC and the GHG emissions by 24.11% and 31.84%, respectively.

This study utilized multisensor image classification data for detecting land-cover and land-use changes. The multisource land-cover classification technique we developed [29] sought to use the best available data for the area’s training and subsequent land classification. In addition, the geospatial data sources provided high-quality inputs of disturbance information to the carbon models. Hence, the estimation of ecosystem and soil carbon flux was found to be more accurate than the model’s estimate using generic data.

The conclusions drawn from this study are specific to the type of land use and properties of the agricultural rotations simulated. In addition, the disturbance-related processes introduced in the simulations were based on the historical trends in the region. Thus, if different parameters are introduced, the uncertainties calculated would likewise differ. The carbon models used in this study were used in accordance with the land-cover type they were designed for (forested or agricultural). And the CBM model in particular has been used for Canada’s carbon accounting [18,19] to fulfil the reporting requirements of the IPCC.

Carbon accounting is crucial for proper resource management. However, uncertainties arise in the carbon estimates from sources such as inconsistent accounting methods and inadequate data inputs. As in this study, a limitation for accurate carbon simulations is accessing accurate and realistic data for model parameters in order to break the cycle of ‘garbage-in’ producing ‘garbage-out’. Geospatial accounting methods, such as those used in the approach in this study, would give a boost to natural resource monitoring. Furthermore, consistent carbon counting methods across regions and provinces would also instill more confidence in the estimates calculated from carbon models. This work adjusted two climate parameters, omitted disturbances, and entered an inaccurate land cover type to observe their effects on the simulations. Yet, further work should be undertaken to investigate the sensitivity of the carbon models to other parameters.

In exploring the long-term effects of different land management practices on the carbon cycle and GHG emissions, our next task will be to simulate emissions up to the year 2100 under the Representative Concentration Pathway 8.5, where the carbon cycle–climate feedback is largest. Under this scenario, no extra climate policies or limits will be instituted, yet land-cover change will be present with agricultural expansion.

An area that has not been considered in this work is the contribution from understory vegetation. We intend to extend the modeling to first identifying the understory using remote sensing tools and then incorporating the understory vegetation in the carbon models. In addition, we aim to increase our understanding of how the use of mechanized farm equipment in the farms will affect the GHG emissions estimated. Although past infrastructure disturbances (such as road developments) were observed in the Clay Belt area, this study did not include urbanization as a possible disturbance in the upcoming years. Future research will make room for urbanization as a factor in predicting its effect on ecosystem carbon and GHG emissions in the area. Future studies will explore these subjects in order to reduce the uncertainty surrounding the calculations of ecosystem carbon and GHG on agricultural lands.

Such findings are significant, as they support a comprehensive understanding of the changes that could occur as Ontario’s forests are converted into agricultural land over the coming years.