The Influence of Interannual Carbon Variability on Long-Term Sequestration in Proximate Northern Forests and Wetlands

Carbon dioxide (CO 2 ) levels are rising dramatically as a result of increased anthropogenic activity. One way of countering excessive CO 2 emissions is by restoring natural ecosystems that have historically been found to be efficient carbon (C) sinks. To be economically viable, these efforts must consider biomes with long-term sustained C sequestration capacities. Low interannual variation in this sink capacity minimizes risk of sequestration reversal. The goal of this study was to compare the interannual variability (IAV) of C at eight proximate Ameriflux eddy covariance sites across northern Wisconsin, Michigan’s Upper Peninsula, Saskatchewan, Alberta, and the Northwest Territories with up to two decades of observations per site. Two wetlands (Allequash Creek (US-ALQ) and Lost Creek (US-Los)) and an unmanaged and managed forest (Sylvania Wilderness Area (US-Syv) and Willow Creek (US-WCr), respectively) were considered in the temperate region while boreal sites consisted of a bog (CA-SCB), peatland (CA-WP1), evergreen needleleaf forest (CA-SCC), and deciduous broadleaf forest (CA-Oas). To consider the fuller C budget, stream discharge data from the United States Geological Survey was also incorporated for some sites. In most of the measured years, on average, net ecosystem carbon dioxide exchange (NEE) in all ecosystems was negative, indicating C uptake by the ecosystem. The standard deviation of the yearly NEE cumulative sums for US-Los was 63 g C m -2 yr -1 while for US-Syv and US-WCr it was 111 g C m -2 yr -1 and 154 g C m -2 yr -1 respectively, implying greater variability for the deciduous forests than the wetlands. A similar result was found for the boreal sites. Mutual information analysis was used to determine influences of carbon components (gross primary productivity (GPP) and ecosystem respiration (RECO)) and drivers (photosynthetic photon flux density, air temperature, latent heat, and streamflow) on NEE. A larger influence on boreal NEE was found relative to temperate NEE on seasonal and yearly scales. NEE was also more impacted by GPP on hourly and diel scales and somewhat equally influenced by both GPP and RECO on mulitday, seasonal, and yearly scales. Our results demonstrate that for these regions, wetlands are a more reliable biome for C storage on decadal scales than forests.


Introduction and Background
One of the most urgent problems facing humanity today is global climate change. Carbon dioxide (CO 2 ) levels are rising dramatically as a result of increased anthropogenic activity. In its Fifth Assessment Report, the Intergovernmental Panel on Climate Change showed that atmospheric CO 2 concentrations had risen from 280 ppm in 1950 to 400 ppm in 2000 (Were et al., 2019). Additionally, computer models from the U.S. Climate Change Science Project have projected that in order to stabilize atmospheric CO 2 at 550 ppm (around two times the industrial level), annual global emissions must be reduced by 75% in the next century (Sundquist et al., 2008). If not combated, climate change will remain a major threat not only to the survival of species and ecosystems, but also to civilization as we know it (Erwin, 2009). Flooding and drought-induced damages have been estimated to cost billions of US dollars annually, with the latter being responsible for the 2007 California wildfires which displaced nearly one million people (Trenberth, 2005).
One way of countering excessive CO 2 emissions is by restoring natural ecosystems that have historically been found to be efficient C sinks. In order to make informed decisions on restoration measures, an analysis of C flux interannual variability (IAV) must be conducted. This is a useful approach for determining the long-term C sequestration potential of biomes. In addition, water and energy fluxes must also be considered since these drivers govern the magnitude of variability over time. It is important to note, however, that these efforts, known as natural climate solutions (NCS), must be complemented by drastic mitigation efforts in industry and energy if we are to meet the Paris goal (Anderson et al., 2019).

The Global Carbon Cycle
C is a vital element that helps maintain life on Earth. It is responsible for creating organic molecules that are crucial for cellular reproduction and metabolism, provides a "heat blanket" in the atmosphere that regulates global temperatures, and is dominantly used as a source of energy via the burning of C-based fossil fuels. In the past few centuries, human-induced emissions due to rapid economic development have caused anthropogenic radiative forcing to increase by 11% from 1750 to 2017, resulting in unprecedented increases in global temperatures (Bruhwiler et al., 2018). The effects of these changes may be seen in the form of more severe floods and droughts (threatening global food production) and intense acidification of the world's oceans (threatening marine life and biodiversity) (Bruhwiler et al., 2018). C is stored in many different reservoirs across the Earth system ( Table 1) and is constantly being transferred mechanically, chemically, or biologically from one stock to another (Figure 1). One of the relatively rapid exchanges occurs between the terrestrial biosphere and atmosphere, where vegetation uptakes and emits CO 2 via photosynthesis and respiration respectively. Anthropogenic disturbances (degradation and deforestation), biological interruptions (insect outbreaks), and natural hazards (fires) may cause a rapid increase in respiration, therefore increasing C emission rates and exacerbating global warming.
The direction of exchange between the ocean surface and atmosphere depends on 1) the relative concentrations of C between each reservoir and 2) surface water temperature (Bruhwiler et al., 2018). Due to spontaneous gas flow from higher to lower concentrations, regional surface waters and atmospheres supersaturated with CO 2 will experience transfer towards the atmosphere and surface water respectively. Additionally, regions with upwelling warm water are known to experience C outgassing while those with sinking cold water are C absorbers (Bruhwiler et al., 2018).
Terrestrial C is a highly dependent variable, owing to natural or anthropogenic changes in water availability (section 1.2) and energy (section 1.3). As a result, a thorough evaluation of drivers will aid in better understanding C's spatial and temporal variability across biomes.

Carbon and Water
The water budget of a particular ecosystem depends on inflows via precipitation (P) and outflows via evapotranspiration (ET) and discharge. Water storage is usually ignored due to its negligible contribution on long temporal scales (Gedney et al., 2006).
Water has a significant influence on the local ecosystem C budget. An increase in atmospheric CO 2 concentrations, for instance, has been shown to close stomatal apertures which decreases transpiration rates and causes water outflow to be dominated by discharge (Gedney et al., 2006). C can also be physically transported through runoff. In fact, a study has shown that lateral hydrologic C inflows and outflows through discharge have a significant impact on wetland C budgets .
C in streams and rivers can take the forms of organic and inorganic. Organic C includes methane (CH 4 ) or is synthesized from photosynthetically-assimilated CO 2 . It can enter inland waters directly as plant detritus (particulate organic C (POC)) or leached material (dissolved organic C (DOC)). Inorganic C can take the form of dissolved CO 2 , carbonic acid, bicarbonate, or carbonate and is the product of chemical and mechanical weathering or rock erosion. These processes release particulate inorganic C (PIC) as calcium carbonate, calcite minerals, and C alkalinity (Drake et al., 2017).
Due to the tightly bound relationships between C and water, it follows that water availability plays a crucial role in determining an ecosystem's C source/sink status. Wetlands such as small lakes, floodplains, and marshes, contain a large portion of the world's C despite their small global land cover (~2-6%). Wetland C has been described as originating from five main reservoirs: plant biomass, POC, DOC, microbial biomass, and gaseous CO 2 and CH 4 . The latter four of which reside in the water or soil (Kayranli et al., 2009).
Another important link can be found between CH 4 and water table level. A study by Harriss et al. (1982) on the Great Dismal Swamp in Virginia has revealed low water table levels to correspond to a greater CH 4 sink. On the other hand, during high water table levels, the swamp was observed to become a significant CH 4 source. This was later confirmed via CH 4 monitoring at a German minerotrophic fen by Augustin et al.
Along with wetlands, another common ecosystem to consider is the deciduous forest. Water cycling at these biomes mainly depends on daily P patterns and soil-water storage capacity. Due to a variety of factors, the interplay between C and water can be very different than that of a wetland.
Temperate forests are known to hold ~20% of the global plant biomass and ~10% of the terrestrial C (Bonan, 2008). Due to deforestation, many of these forests have become C sources. There have been efforts, however, that aim to restore these ecosystems and bring them back to C sink status (Bonan, 2008). Temperate deciduous forests respond to high atmospheric CO 2 concentrations through reducing stomatal conductance, leading to a decrease in transpiration which ultimately augments discharge (Leuzinger and Körner, 2010).

Carbon and Energy
The interplay between C and energy can be described in terms of radiation absorption and chemical reactions. One example is the role of latent heat (LE) in causing turbulent energy exchange in the form of ET. As a result, larger LE contributes to a stronger linkage between heat fluxes and water transport.
Solar energy and vegetation also play a crucial role in ecosystem C production via photosynthesis. C uptake typically peaks during the warm season and when sunlight is available.
The dependence of C accumulation on sunlight can be quantified in terms of net ecosystem exchange (NEE) and photosynthetic photon flux density (PPFD). One finding indicated that NEE increases and then levels off at a positive value as PPFD increases. The same study also found these trends to be similar for different ecosystems within the same category (rich and poor fens) (Frolking et al., 1998). Furthermore, daytime CO 2 uptake and PPFD were found to be positively correlated for a mixed conifer forest, mixed short-grass prairie, and sagebrush shrubland (Kelly et al., 2002).

Water and Energy
The interplay between water and energy has not been as well documented as those of C/water and C/energy. However, one notable relationship between these two components is encompassed in the LE term. In this context, by definition, LE represents the amount of heat added to a substance to change its phase while keeping its temperature constant. In ecosystem dynamics, LE is a mechanism that cools the surface through water evaporation. This form of energy is mainly responsible for the vertical transfer of water vapor from the surface to the atmosphere.
Due to the strict tie between LE and water, energy partitioning is also governed by water availability. For example, a study conducted on two ecosystems in Florida, scrub oak and pine flatwoods, found that the partitioning of net radiation into its sensible and latent heat flux components was mainly driven by fluctuations in soil moisture and leaf area. Specifically, a decrease in Bowen ratio was found with increasing soil moisture and leaf area (Bracho et al., 2008). The rate of ET has also been shown to vary with the development of leaf area index, generally increasing and decreasing as leaves develop and senesce respectively (Zhou and Zhou, 2009). Moreover, increases in canopy interception have been shown to increase ET. This reduces discharge and eventually leads to a decline in water availability (Liu et al., 2016).

Wetlands
Due to their high potential to sequester atmospheric CO 2 , wetlands are regarded as one of the most important ecosystems when it comes to global warming mitigation (Were et al., 2019). Atmospheric CO 2 is transferred, accumulated, and trapped into wetland soils as soil organic matter. Other major inputs of C include organic matter from senesced vegetation and lateral transport of dissolved C via inflowing waters. C may also exit the wetland through outflowing water or, most commonly, be released directly as CO 2 and CH 4 through decomposition (Villa and Bernal, 2018).
Other factors such as water table depth (WTD), temperature, and oxygen availability also have an influence. WTD is the most significant driver of greenhouse gas fluxes in peatlands. Agricultural activity has caused peatlands to be drained for cultivation, converting the land from a net C sink to a net C source.
These ecosystems are also extremely sensitive to small changes in WTD (Evans et al., 2021). In addition, during dry and wet periods, the land has been found to become a net CH 4 sink and source, respectively (Kayranli et al., 2009).
Temperature is also an important controller of C sequestration due to its role in accelerating organic matter decomposition. Higher temperatures have been associated with higher CH 4 emissions but do not override the overall C sink capacity of most wetlands (Mitsch et al., 2012;Olsson et al., 2015;). Aerobic conditions only allow CO 2 formation during decomposition while both CO 2 and CH 4 are formed in anaerobic conditions (Kayranli et al., 2009).
Although these dynamics have been found to be generally true for the typical midlatitude wetland, it is important to acknowledge that these findings may not hold for other regions. Temperature and water table level in tropical wetlands and ultra-wet wetlands, for example, have a weaker influence on C sequestration (Christensen et al., 1998;Sjögersten et al., 2014;Villa and Bernal, 2018).

Terrestrial Forests
Terrestrial forests (hereon referred to as "forests") have also been known to play a key role in climate warming mitigation. Like wetlands, many variables play a role in determining their net C sink/source status. Over the past few decades, deforestation and forest degradation have been responsible for a fifth of global greenhouse gas emissions (Schrope, 2009), while intact forests are important global C sinks (Pan et al., 2011). Specifically, these ecosystems store 40% of the belowground C and 80% of the aboveground C while also carrying 90% of the C exchanges between the land and atmosphere (Wei et al., 2014).
A majority of the C storage in forests takes place in soil and vegetation. Boreal forests store a majority in soil while temperate forests sequester mostly in woody biomass. Furthermore, influences on plant and soil C sequestration potential are driven by forest management and soil moisture respectively (Ma et al., 2015).
The ability of a forest to efficiently sequester C is sensitive to factors such as temperature, precipitation, topology, human activity, and natural disturbances (Ma et al., 2015). Soil organic C, for instance, is known for its high sequestration capacity but is believed to have been reduced due to anthropogenic land cover change (Murty et al., 2002). Climatic changes over the past few decades are responsible for elongating the growing season and have a positive impact on forest productivity (Boisvenue and Running, 2006). In some cases, this has even resulted in higher growth rates (Cole et al., 2009).

Interannual Variability
Researchers have been studying the interactions between the land and atmosphere for multiple decades. Throughout this endeavor, new technologies have been implemented to facilitate the quantification of meteorological, hydrological, and gas flux variables. Eddy Covariance (EC) (see section 2.4), which debuted in the early 1990s, has proved to be the most ubiquitous method for obtaining flux measurements. The community's confidence in this approach eventually led to the integration of international/continental flux networks across the globe (Baldocchi et al., 2017). Since these stations have been sustained for decades, it is now possible to conduct analysis on C interannual variability (IAV) between sites on multi-decadal scales.
Due to the climatic impacts on carbon components, NEE may experience significant year-to-year fluctuations. For example, the spatial variability of ecosystem respiration (RECO) and gross primary productivity (GPP) has been shown to be strongly linked to temperature and precipitation while radiation availability only influences GPP (Yi et al., 2010).
Based on location and ecosystem type, C flux IAV can vary greatly. For datasets that cover temperate deciduous forests, the standard deviation of interannual NEE is approximately ± 100 g C m −2 yr −1 (Baldocchi et al., 2017). Furthermore, a mixed forest in southern Ontario, Canada was found to switch relatively quickly between being a C source or sink. Over the 17 years of the study, net ecosystem productivity increased by approximately 15.7 g C m −2 yr −1 due to decreases in RECO and increases in GPP. It was also found that PPFD and TA were the main drivers of C fluxes at the site (Froelich et al., 2015). On a seasonal scale, earlier start to leaf emergence caused by earlier spring TA has also been deemed responsible for these interannual variations (Saigusa et al., 2005).
Wetlands are also prone to NEE IAV. On average, these ecosystems mainly act as C sinks and their dynamics are driven mostly by plant phenology and WTL. Lower WTL reduces CO 2 uptake. Therefore, future drier conditions due to climate change may lead to NEE that is less negative during the growing season (McVeigh et al., 2014). In addition, warmer TAs during the winter have been associated with increases in CO 2 uptake (Helfter et al., 2015). It is important to note, however, that one wetland has been found to change its C sink/source status in a period of as little as two years, with its NEE variability being mostly attributed to the transition period between vegetation growth and senescence (Serrano-Ortiz et al., 2020). Other drivers of IAV include vapor pressure deficit, annual maximum leaf area index, and growing season mean stomatal conductance.
A study by Zscheischler et al. has also found that IAV is caused by the most active (high percentile) fluxes at temperate forests on an hourly and daily scale. Therefore, it was concluded that IAV is governed by only a small handful of short-term fluxes (less than 20% of the total dataset).
Lastly, various models have been introduced to predict future trends in biome IAV. While these have been shown to reproduce IAV magnitudes, they fail to agree with the timing of observations made in mid-latitude forests across North America. GPP and RECO have been found to be drastically underestimated in deciduous and evergreen forests respectively, indicating model errors due to processes related to vegetation type (Keenan et al., 2012).

Natural Climate Solutions
NCS involve cost-effective restoration and conservation practices that increase C sequestration and consequently mitigate impacts of climate warming. Different biomes require different strategies. For example, forests may see increased fire management and reforestation while wetlands would need peat and coastal restoration.
NCS will play a major role in limiting increases in global temperatures to 2°C due to its high C sequestration potential in relatively short amounts of time. When taking into account the constraints of food security, fiber security, and biodiversity conservation, it has been found that the maximum potential for NCS is ~23.8 Pg CO 2 equivalent (CO 2 e) yr -1 . In addition, these practices are expected to deliver ~37% of CO 2 e mitigation from now until 2030 and ~20% from now until 2050. These estimates are more likely if fossil fuel emissions are held constant for ten years, brought down to 7% of current levels by 2050, and completely abolished by 2095 (Griscom et al., 2017).
Forests are one of the most important biomes for C sequestration. In 2018, 11.6% of the total annual greenhouse gas emissions were offset by forests in the contiguous United States (EPA, 2020). The most effective mitigation practices have been found to be reforestation (307 Tg CO 2 e yr -1 ), tree management (267 Tg CO 2 e yr -1 ), and fire management (267 Tg CO 2 e yr -1 ) (Fargione et al., 2018).
When stands reach the climax of mean annual increment, they are usually logged for economic purposes. If the harvesting interval is increased (i.e. less harvests during the same period of time), it would allow for potential increases in forest C stocks (Kaarakka et al., 2021).
Wetlands are also known to have great sequestration potential and are therefore a prime candidate for NCS. Inland wetland restoration is known to have a net cooling effect on decadal and century time-scales and therefore would not prove to be effective for rapid mitigation. These biomes also require a major investment (~$4200-$49200 per ton C) rendering them economically inferior compared to other restoration efforts (Taillardat et al., 2020). On the other hand, conservation and restoration of coastal wetlands (such as mangroves) is more cost effective (~$1800 per ton C). The downside, however, is that these ecosystems cover a small percentage of the global land, making them appropriate for mitigation efforts on a national rather than global scale (Taillardat et al., 2020).
Despite evidence of their effectiveness as a mitigation strategy, NCS will undoubtedly face challenges in its implementation. Some potential roadblocks include complexities in social-ecological systems and doubt regarding the costs of ecosystem services. Furthermore, its ongoing uncertainty as a long-term C sink may lead to hesitation of action by policy makers (Seddon et al., 2019).
In Wisconsin, NCS has not been explored as much as in other regions despite the relatively high density of forests and wetlands in the area. Our work on the sites in northern Wisconsin and Michigan's Upper Peninsula (UP) will attempt to fill in these knowledge gaps and paint a clearer picture regarding C sequestration potential in the upper Midwest.

This Study
To better understand how local biomes can contribute to global warming mitigation, accurate descriptions of the C cycle must be made on the regional scale. Previous studies on C dynamics in wetlands have considered lateral transport (Gao et  Taking into account everything that is already known about the C sink capacity of forests and wetlands, this study will aim to conduct a direct comparison between four proximate forests and wetlands in northern Wisconsin and Michigan's UP. The water and energy budgets of each ecosystem will be analyzed with an attempt to determine their broader role in influencing long-term C sink capacity. Comparisons will be carried out on multiple time scales (hourly, diel, multi-day, seasonal) to determine C IAV and its drivers. These temperate sites will also be compared to boreal forest and wetland sites to determine the influences of location on carbon-driver dynamics. This study will aim to answer the following questions: • How does the NEE IAV of the forest sites compare to that of the wetland sites? Is this more a function of GPP or RECO? • Are there any significant differences in C drivers between or within biome types on multiple time scales? • What do these findings imply for NCS?
Our hypotheses are the following: • Stream C export will be more pronounced at the wetland sites, therefore leading to lower IAV than the forest sites • Both biome types will have PPFD serve as the dominant driver on short time scales (hourly, diel) while water availability and temperature will be responsible for variability on longer time scales (multi-day, seasonal) for the wetlands and forests respectively.

Allequash Creek (US-ALQ)
This site is located in the Allequash Wetland near Allequash Creek (Latitude: 46.0308, Longitude: -89.6067) (Figure 2). The wetland is dominated by peat and covers ~32 hectares of the Trout Lake basin. It is one of the regions within the National Science Foundation's North Temperate Lakes Long-Term Ecological Research and is also included in the United States Geological Survey's (USGS) Trout Lake Water, Energy and Biogeochemical Budgets site (WEBB). The basin is monitored with a network of 60 observation wells and 4 stream gaging stations.
Since the soil consists of highly conductive outwash sand and gravel on top of crystalline bedrock, groundwater discharge to Allequash Creek is highly promoted (Anderson and Lowry, 2007). The creek flows downstream through the wetland and drains into Allequash Lake. The vegetation comprises a mix of broad-leaved evergreen/deciduous shrub/trees and narrow-leaved persistent emergent and wet meadow (Turner et al., 2019).

Lost Creek (US-Los)
This site is situated in a shrub wetland near Lost Creek (Latitude: 46.0827, Longitude: -89.9792; Figure 2) at an elevation of ~480 m and was established in September 2000. Since the creek and associated floodplain provide a consistent source of water, the wetland experiences a large amount of peat accumulation. Vegetation near the flux tower mostly consists of alder, willow, and sedges. Since it is located near a stream and has a long, narrow shape, this ecosystem shares many of the characteristics of typical minerotrophic wetlands in the Great Lakes region (Sulman and Desai et al., 2009). Its climate is characterized by warm, wet growing seasons and cold, dry winters. The growing season in the region starts at around June and ends in August (Pugh et al., 2017).

Sylvania Wilderness Area (US-Syv)
The Sylvania Wilderness Area site (hereafter referred to as Sylvania) (Latitude: 46.2420, Longitude: -89.3477; Figure 2) is an old-growth forest that was established in late 2001. Tree age ranges from 0 to 350 years old with the dominant vegetation consisting of sugar maple and eastern hemlock. Other types of vegetation include yellow birch, basswood, and ironwood. The climate can be considered northern continental (short growing seasons and cold winters) and the average elevation is ~540 m. Most of Sylvania's area is sheltered by hemlock-northern hardwood forest (66%), along with lakes and forested wetlands/marshes covering 21% and 13% of the area respectively (Desai et al., 2005).

Scotty Creek Watershed
CA-SCC and CA-SCB are characterized as a boreal forest and wetland respectively. These sites will be used to compare C dynamics and drivers with the four main temperate sites. Both are located in Scotty Creek (Latitude: 61.3, Longitude: -121.29; Figure 2) which is a 152 km 2 watershed near Fort Simpson, NT in the southern Taiga plains in northwestern Canada. The region experiences an average TA of -2.8 °C and mean total precipitation of 388 mm (149 mm of which is snow). The wetlands in this landscape are composed of bryophytes, ericaceous shrubs, and to a lesser extent black spruce and tamarack. The forests consist of a dense cover of black spruce along with ericaceous shrubs, bryophytes, and lichens dominating the understory (Helbig et al., 2017).

Alberta Western Peatland (CA-WP1)
CA-WP1 (lat: 54.9538; long:-112.4670; Figure 2) is one of many ecosystems located in the southern boreal forest of Canada. The regional climate has been characterized as having long, cold winters and short, cool summers with a mean annual temperature of 2.1 °C. Average yearly precipitation is composed of both rain (382 mm) and snow (122 mm). The terrain is relatively flat and is composed of nearly homogenous vegetation. Dominant tree and shrub species include   . No data was available for US-WCr or the boreal sites since no gages were found within a reasonable distance to be considered as a proxy. The distance from the gage at Cisco Lake Outlet and US-Syv is 8 km and the distance from the Bear River gage to US-Los is 4 km. The Allequash Creek gages are relatively closer to US-ALQ being a distance of 3.7 km, 0.6 km, and 0.24 km away for County Trunk Highway M, Sayner, and Site No. 3 respectively (Figure 3).

The Eddy Covariance Method
The Eddy Covariance (EC) method has been used widely by researchers due to its relatively few theoretical assumptions and its potential for measuring micro-meteorological variables. Gas, energy, and momentum fluxes may be directly measured, giving direct observations of land-atmosphere interactions within the flux tower footprint. EC has also benefited from technological advances in the past two decades, making it one of the most popular methods for data acquisition among the micro-meteorological community (Liang and Wang, 2020).
Measurements of ecosystem variables are determined through the covariance between vertical wind velocity and the quantity being measured. An anemometer is used to measure wind speed and direction, an infrared gas analyzer measures gas concentrations in the air, and radiation sensors measure solar radiation (section 2.5).

Theory
This overview will follow the discussion outlined in (Baldocchi, 2003). Vertical turbulent motions in the atmosphere are responsible for transporting CO 2 and other constituents between the atmosphere and biosphere. The EC technique involves keeping track of these turbulent motions in order to determine net inflows and outflows of these constituents with respect to the land and atmosphere. To achieve this, the Reynolds averaging method is applied on the following flux equation: where is the instantaneous mass flux density, is the vertical wind velocity, and is the CO 2 density. This will yield the following expression: where represents the average flux density of the constituent over some specified period of time (i.e. hourly, monthly, yearly, etc.), is the air density, and is the ratio between and ( / ). The overbars represent the time averages and primes indicate deviations from the mean (i.e. ). Positive and negative values of covariance denote a loss and gain of constituents by the ecosystem respectively.
Interpretations of EC measurements can be done through assumptions and manipulation of the mass conservation equation: where is defined above, is the three dimensional flux vector, is the source/sink term, and , , and are the wind velocities in the , , and directions respectively. Term I represents the change in CO 2 mixing ratio with time, II is the advection of CO 2 by the zonal, meridional, and vertical winds, III is the flux divergence, and IV takes into account any biological sources or sinks and is dependent on location.
With two key assumptions, equation 3 can be simplified such that only one term remains on each side. The assumptions are the following: the average CO 2 mixing ratio is constant with time and the terrain being measured is flat and horizontally homogeneous. This implies that term I and the horizontal parts of terms II and III may be set equal to zero. As a result, only the vertically-dependent source/sink term and the vertical flux divergence terms survive, yielding the following: Multiplying by and then taking an integral from zero to an arbitrary height above the canopy gives an expression for the mean vertical CO 2 flux density at that location: where is the mean vertical CO 2 flux density at the surface (i.e. fluxes from the soil or underlying vegetation). One of the main purposes for the EC method is to evaluate in equation 5 for numerous ecosystems around the globe.

Filtering
High quality data acquisition via the EC method requires high turbulence in the measurement area. During stable conditions, CO 2 flux measurements have been found to be underestimated using this technique (Goulden et al., 1996). Therefore, filtering methods must be applied to discard data taken during periods of low turbulence.
Goulden et al. proposed that friction velocity ( ) may be used as a means of determining the validity of flux measurements. Once a threshold ( ) is determined, any fluxes below that threshold would be discarded. One shortcoming of this technique is that the same cannot be used universally for all sites due to its dependence on local topography, surface roughness, and heterogeneity. Despite this, its effectiveness has made it the most-used flux-filtering method.
To obtain a lower threshold, one must know the point at which NEE becomes insensitive to changes in . Therefore, nighttime NEE is plotted as a function of several different classes. Afterwards, statistical comparisons are made to determine when the difference between the averaged nighttime NEE of a class is not statistically significant with the NEE at higher values of . The value where this occurs is deemed to be the . If it is found that is co-varying with other drivers of respiration, the data is normalized to remove any covariation between respiration and the variable in question. Due to the likeliness of covariance with temperature in temperate regions, an NEE-temperature function is typically used for normalization (Aubinet et al., 1999;Aubinet et al., 2012). (with one having zero slope). The point at which these lines cross is noted as the "change point" and is used as a guide for estimating appropriate thresholds. On the other hand, the MP method compares the average NEE of the twenty classes with the average NEE of the ten classes with the greatest magnitude. A threshold is then determined if the class contains a mean nighttime NEE greater than 99% of the mean NEE at the ten highest classes. An improvement of this method has also been implemented that reduces the effects of noise. More details may be found in Pastorello et al., 2020.

NEE Uncertainties
To determine NEE uncertainties for this study, the REddyProc Software was used from the Max Planck Institute for Biogeochemistry. The software was utilized as a means to create reasonable bounds for NEE. The percent error between the exact (measured) and approximate (gap-filled) cumulative annual NEE was calculated to determine uncertainty. Most uncertainties were within 10% of the measured values for all sites and years except for CA-Oas. This discrepancy is believed to be due to the disparate flux processing methods used at the site.

Gap Filling
Due to the great temporal resolution of EC flux measurements, gaps in data are not uncommon. These may be due to various reasons (i.e. winter power failures, dusty solar panels, etc.). Since one of the many purposes of EC is to compile averages of C fluxes over long time scales, it is essential to have sophisticated techniques for filling in missing data. According to the discussion in Papale, 2012, there are multiple methods for gap filling, including mean diurnal variation (MDV), look-up tables (LUT), artificial neural networks (ANNs), and nonlinear regressions.
The MDV approach fills gaps by taking the average of valid values on adjacent days during the same period. It is recommended that the length of the averaging period doesn't exceed two weeks since after this time nonlinear dependencies may introduce significant uncertainties and errors. The position of the averaging periods may be fixed or may vary depending on the data. Preferably, the latter should be used since periods will be defined based on each gap, resulting in more accurate results.
LUT looks at meteorological conditions for missing fluxes and then gap-fills with the average of valid flux measurements during similar conditions. Several classes of meteorological drivers (e.g. radiation, air or soil temperature, soil water content, vapor pressure deficit, etc.) are considered depending on site characteristics. These are then tabulated and "looked-up" in order to gap-fill fluxes accordingly. Linear interpolation is used if no flux data is available for a specific combination of the drivers.
ANNs consist of empirical nonlinear regression models that predict missing fluxes through various algorithms. These operate in a set of nodes and are connected to weights that represent different regression parameters. To utilize an ANN, it must first be "trained". Training involves feeding sets of data to the ANN (usually meteorological drivers) and receiving the associated outputs. The dependency of the output on the meteorological drivers is mapped onto the weights which ANN then uses to gap-fill the data. This method is used as a standard technique for gap-filling fluxes in the FLUXNET global network (along with marginal distribution sampling which is an extension of LUT; see Reichstein et al.).
Similar to the LUT method, the nonlinear regressions approach bases its gap-filling on relationships between CO 2 fluxes and its associated drivers. Instead of tabulating data, however, parameterized nonlinear equations are used to estimate missing data. Daytime gap-filling typically takes PPFD to be the main driver and uses a rectangular hyperbola or exponential function for estimation (Falge et al., 2001). Due to lack of PPFD during the nighttime, equations where temperature is the main driver are used instead (Lloyd and Taylor, 1994;Falge et al., 2001;Moffat et al., 2007). To obtain estimations for gap-filling, the inputs for these equations must all be valid fluxes.
It's important to note that while there are a variety of gap-filling techniques, methodological comparison studies have shown most approaches to agree with synthetic data on several timescales (Moffat et al., 2007;Desai et al., 2008). However, other factors, such as ease of use, may lead researchers to employ one method over another.

Partitioning Fluxes
Due to technological limitations not allowing for direct measurements of GPP and RECO, these dependent variables must be inferred from measured NEE via equation 6: To obtain an estimation of NEE, the following equation is used: where is the estimated turbulent CO 2 flux via EC through the horizontal plane at an arbitrary height above the canopy and is the change in CO 2 storage below the horizontal plane (see sections 1.4.2 and 2.5 in Eddy Covariance (Foken et al., 2012;Munger et al., 2012)). Many methods have been proposed on how to partition NEE into its constituent fluxes. One of the most popular is the night-time data based method.
The night-time based approach exploits GPP's heavy dependence on radiation. After data has been quality filtered, all night-time data points for NEE are assumed to represent RECO (i.e. GPP is zero). A simple model of RECO's dependence on TA is then formulated. One of the most common methods is using the equation: • where is the base respiration at 10°C, is the TA, and is the temperature sensitivity parameter which represents the change in RECO for every 10°C change in temperature. After calculating RECO, equation 6 is used to calculate GPP.
Unfortunately, RECO is driven by a multitude of other variables (e.g. nutrient levels and soil moisture) causing this method to introduce possible errors and biases. The lack of measurements of these variables at FLUXNET sites is one of the main reasons why nutrient level and soil moisture-dependent models are lacking (Reichstein et al., 2012).
Other models include using variants of the Arrhenius equation (Lloyd and Taylor, 1994) or using soil temperature-dependent models such as the Eyring function (

Flux Instrumentation
The instrumentation at all flux tower sites was virtually the same.

USGS Gage Instrumentation/Methods
All streamflow measurements made by the USGS follow a similar protocol. Streamgages are unable to measure discharge directly. Instead, these instruments directly measure stream velocity, stream depth, and stream width and then mathematically compute streamflow using the following equation: Where is discharge and , , and are stream depth, width, and velocity respectively. The units for discharge are typically given in . While depth and width can be calculated relatively easily across many sites, stream velocity measurements are less trivial. The most common methods used include The Mechanical Current-Meter Method and The Acoustic Doppler Current Profiler Method (ADCP). The former involves partitioning depth and width measurements in a way that resembles the Riemann Sum approximation. The stream velocity is then measured by placing a wheel of metal cups in each segment. As water flows, the cups rotate around a vertically oriented axis which then prompts an electric signal transmission. This signal counts the number of revolutions per unit time which is then converted to linear water velocity. Each partitioned width and height is multiplied by its respective water flow quantity and are then summed together to determine the stream velocity.
The latter method improves measurement accuracy by taking advantage of the doppler effect. ADCPs are hydroacoustic instruments that send pulses into the water and then measure the amount of time it takes for signal retrieval. The relative change in frequency is then translated into water velocity.

Data Acquisition and Analysis
The Ameriflux data repository was used to acquire all C, water, energy, and micro-meteorological data for all sites. Therefore, all data measurement and processing methods adhere to standard Ameriflux protocols (https://ameriflux.lbl.gov/). All temperate sites are For variables where three-index position qualifiers were used (i.e."VariableName_HorizontalPosition_VerticalPosition_Replicate"), the position qualifier "VariableName_1_1_1" was chosen. More information on positional qualifiers and their use in data labeling can be found on the Ameriflux Website (https://ameriflux.lbl.gov/) → Data → About Data → Data Variables.
The Streamflow data was obtained from proximate USGS gage locations at each site (excluding US-WCr and boreal sites; see section 2.3; https://www.usgs.gov/). Relatively minimal data gaps were present in the gage data.
Regression analysis between PPFD and NEE was plotted via numpy's polyfit command in Python. We then obtained numerical values for the Pearson correlation coefficient and p-value using the "pearsonr" function in the Scipy Python library. Due to the non-linear relationship between C components and TA (Fei et al., 2017), we utilized a quadratic fit for GPP and NEE and an exponential fit for RECO (see section 3.2.4). The numpy polyfit function was then used to implement fits across years. To determine the r 2 and p values for these fits, the "Real Statistics Data Analysis Tool" external Microsoft Excel package was utilized along with the scatter plot function.

Wavelet-Based Time Scale Decomposition
The maximal overlap discrete wavelet transform (MODWT) was used to decompose the time scales of variability for NEE, GPP, RECO, and various C drivers. These half-hourly fluxes (hourly for discharge) were reconstructed for scales 1 (2 1 measurements equal to 1 hour) to 14 (2 14 measurements equal to 341 days). Scales 1-2, 3-6, 7-10, and 11-14 represent hourly (small perturbations such as clouds passing overhead), diel (day-night cycles of radiation and temperature), multiday (variabilities in water table and synoptic weather), and seasonal (vegetation phenology, solar cycle, and hydrologic cycle) scales respectively (Sturtevant et al., 2016). The Wavelet Methods for Time Series Analysis (WMTSA) as part of the MATLAB wavelet toolkit was utilized to wavelet decompose data into these scales (Cornish et al., 2003). Due to this toolkit's incompatibility with data gaps, all gaps were dropped during analysis.

Mutual Information
In information theory (Shannon 1948), mutual information describes the mutual dependence between two random variables. Its purpose is to measure how much information may be received from one variable through observation of the other variable. In other words, it is defining the mean tendency for a single paired state of two variables to co-exist (Fraser & Swinney, 1986). Since this method calculates the statistical dependence of Y on X, it lends itself to be an appropriate approach for resolving C flux-driver relationships (Knox et al., 2021).
In this study, we input one of the variables to be NEE and the other to either be GPP (no half-hourly data for CA-Oas), RECO, PPFD (incoming shortwave radiation was used for CA-SCB), TA, LE, or discharge (streamgage data only available for US-ALQ, US-Los, and US-Syv; no yearly scale analysis was conducted due to limited resolution). A mutual information score (MIS) was then calculated between NEE and the other variable to determine the extent of dependence. High and low magnitudes of MIS indicate stronger and weaker links between NEE and the other chosen variable, respectively.
One weakness with utilizing this method is the MIS's high dependency on the number of sampled data points. To account for this, we used a histogram approach. A set number of bins were determined depending on each site's number of data points via the equation On average, equation 10 provides five points for each cell on the histogram for two uniformly distributed random variables. This technique mirrors the adaptive partitioning approach discussed in Cellucci et al. Since all MISs are now characterized based on , direct comparisons between sites is now possible.

Interannual Variation of NEE
Among the temperate sites, the deciduous broadleaf forests had the largest range of annual cumulative NEE relative to the wetland sites. US-ALQ (2 years of data) and US-Los (15 years of data) had ranges of 17 and 193 g C m -2 yr -1 respectively. For US-Syv (12 years of data) this was 362 g C m -2 yr -1 and for US-WCr (16 years of data) it was 551 g C m -2 yr -1 . Similar results were found for the boreal sites. CA-SCB (3 years of data), CA-WP1 (6 years of data), and CA-Oas (10 years of data) had ranges of 85, 265, and 382 g C m -2 yr -1 respectively. CA-SCC (2 years of data), however, saw a range of only 10 g C m -2 yr -1 . This site is composed of mostly evergreen needleleaf forest and therefore has different carbon dynamics than that of a deciduous broadleaf. An alternative and more likely explanation for this discrepancy is a lack of sufficient data. All sites transitioned from C source to sink in a similar window period (between DOỸ 150-200). For the temperate region, yearly NEE variability was found to be substantially smaller at the wetland site (US-Los) than at the forest sites (US-Syv and US-WCr). On the last DOY, NEE at US-Los varied between -209 and -9 g C m -2 yr -1 across years, while at US-Syv and US-WCr, it was between -280 and 97 & -575 and -4 g C m -2 yr -1 respectively (Figure 3).
Annual sums of GPP at the temperate sites were greater at the deciduous broadleaf/mixed forest sites relative to the wetland sites. On the other hand, RECO was similar in magnitude across the sites with the exception of US-Syv. None of the sites were dominators in terms of NEE annual sums across all years. Boreal sites were also found to be carbon sinks across all years and were comparable to their temperate counterparts (Figure 4). The deciduous broadleaf forest site (CA-Oas) was recorded to have a larger carbon sink potential relative to the wetland site (CA-WP1). CA-SCC and CA-SCB had limited data but the former saw higher rates of GPP and RECO relative to the latter for adjacent years (2014-2016). Due to the near equal partitioning of NEE, both biomes became either carbon neutral or modest carbon sinks during the measured years (Figure 4).

NEE Zero-Crossing Day
Both types of temperate ecosystems had an NEE zero-crossing day (ZCD) within the range of the 150th (May 30th) and 210th (July 29th) day of the year (DOY). US-Los saw an earlier NEE sign change for each year from 2001-2020, while US-Syv and US-WCr were more variable (Figure 5). From 2002From -2005, the minimum and maximum DOY when US-Syv became a net C sink was 138 (May 18th) and 189 (July 8th) respectively. The site then saw a steady earlier NEE sign change from 2014-2018 and became a C sink relatively later in 2020. US-WCr's DOY sign change was also very variable with no specific trend throughout the measured years (not shown).
For the boreal sites, CA-Oas experienced moderate variability with its average ZCD being around day 170 (June 19th). On the other hand, CA-WP1 became a net C sink during later parts of the year after 2005. We also found a steady earlier ZCD for CA-SCB and CA-SCC, but more data is needed to determine long-term trends (not shown).

Component Fluxes
GPP annual sums were also found to be larger than most of the RECO annual sums at the sites. In general, US-WCr and CA-Oas experienced higher GPP values than US-ALQ, US-Los, and CA-WP1 (Figure 6). US-Syv had approximately the same yearly sum of GPP and RECO over the studied years. Cumulative GPP showed a similar trend across all sites and all studied years. Each site initiated GPP at ~150th DOY. There was year to year variability, however, for GPP annual sums, with US-Syv showing more IAV on average than US-WCr and US-Los. In comparison, average year-to-year variability for RECO annual sums were almost identical across all sites. RECO began earlier at US-WCr for most years (between DOY 50 and 100) and had a smoother curve. US-Los and US-Syv, on the other hand, initiated RECO at around DOY 150 and saw steeper slopes as they entered summer months (Figure 4).

Radiation
One of the main drivers of C accumulation across all sites was PPFD. As PPFD increased, NEE decreased most significantly during the months of June through September for the measured years. GPP's sensitivity to solar radiation caused a dramatic increase in production during these months, causing the second term on the right side of equation 6 to dominate. This is also reflected in (Figure A.3) and (Table A.3), with the relationship between the two variables exhibiting a steeper positive slope (increased r-value) in the summer months for both wetlands and forests. Ultimately, this led to a general decrease in NEE during the summer months across all sites and years (Figure A.2

; Table A.2).
Even though PPFD was found to be a driver for C accumulation on a multi-year scale, the correlation coefficient experienced substantial interannual variation. For the forest sites, the correlation between NEE and PPFD was somewhat constant over the years, with the strongest dependency occurring in August (not shown). The wetland site experienced the most variability across years, with a general increase in PPFD reliance for C accumulation over the two decade period. On average, Pearson correlation coefficients for NEE vs PPFD and GPP vs PPFD were -0.799 and 0.804 respectively from 2001-2010.
Of the three summer months, US-Syv and US-WCr had consistent periods of higher correlation during July and August throughout the two decades. US-Los, however, was more variable in this regard, showing increases and decreases in radiation dependence across all months and years (not shown).

Precipitation
Average cumulative precipitation across all years was greatest at US-WCr (mean ± standard error: 685 ± 38 mm) followed by US-Los (677 ± 69 mm) and US-Syv (378 ± 83 mm).  P and carbon component data were divided into daytime and nighttime segments to determine relationships with and without solar radiation. Daytime and nighttime were defined as hours from 9:00:00 to 17:00:00 and from 17:30:00 to 24:00:00/00:30:00 to 8:30:00 respectively. CA-SCC, CA-SCB, and CA-WP1 utilized rainfall data instead of precipitation data. All P at 0 mm was dropped. GPP did not have any significant relationship with P during the daytime or nighttime across sites (not shown). RECO was generally positively correlated with increases in P. Daytime correlation with RECO was higher among the boreal sites while during the nighttime, RECO showed higher correlation among the temperate sites. The RECO-P relationship was not significantly affected by the type of ecosystem (Figure 7).

Air Temperature
Monthly averaged TA was found to be nearly identical across all sites, with temperatures generally residing between -16 and 22°C. Yearly averaged TA, however, was more variable between sites. US-WCr was the warmest amongst the sites from 2001 through 2008. From 2015 through 2020, US-WCr starts out warmer but is then surpassed by US-Los beginning in 2017. US-ALQ had an abnormally warm year during 2015, but this is likely due to relatively high gaps in TA data during that year (Table A.5). Excluding 2015, US-ALQ was cooler relative to US-Los and US-WCr during 2016, 2019, and 2020. Due to high amounts of data gaps during 2017 and 2018 for US-ALQ, TA was not considered during those periods ( Table A.

5).
To quantify the responses of GPP, RECO, and NEE to TA, quadratic (equation 10) and exponential (equation 11) regression models were used based on Fei et al., 2017: where , , , , and are fitted parameters, is ecosystem respiration, and is air temperature. The quadratic model was applied for GPP/NEE versus TA and the exponential model was utilized to evaluate the relationship between RECO and TA. Overall, the regression models between the C components and TA agreed well across all sites and years (Figure 8). RECO was observed to increase exponentially with TA, while a quadratic relationship was found between NEE/GPP and TA. The regression equations were also variable from year to year. For equations representing RECO vs TA, the y-intercept term showed small variation (0.55 -1.75 across all sites and years).
Due to its limited impact on the steepness of the regression curve, only the exponential growth constant was considered when determining TA's influence on RECO. When considering sites that had the least amount of data gaps (US-Los and US-WCr with n=13 and n=18 respectively), we observed the greatest impact of TA on RECO at US-Los relative to US-WCr (mean value for exponential growth constant was .09 and .06 for US-Los and US-WCr respectively). Furthermore, the variability of impact across years proved to be higher at US-WCr (SD = .012) than at US-Los (SD = .006).
To determine the extent of TA's impact on NEE/GPP, the year to year variabilities in the quadratic coefficient (QC) were considered. The average QC for GPP was .006, .008, .009, and .01 for US-ALQ, US-Los, US-Syv, and US-WCr respectively. The variability in QC (standard deviation) was similar across all sites. For NEE, all QC was found to be negative. The average QC was largest at US-WCr (-.006), followed by US-Los (-.003), US-Syv (-.0026), and US-ALQ (-.0018). The greatest variabilities in QC occurred at the forest sites (US-Syv: 0.001; US-WCr: 0.004) relative to the wetland sites (US-ALQ and US-Los: .0005).

Evapotranspiration
GPP was more impacted by increases in LE than RECO across all sites and years. RECO experienced a positive correlation with LE, with US-WCr having the lowest correlation relative to other sites (Figure A.5). In particular, increases in LE promoted GPP. Most Pearson correlation coefficients were calculated to be greater than 0.65 across all years and sites (Table A.6). Consequently, this resulted in an anti-correlation between NEE and LE across sites and years with most coefficients being less than -0.6. Some years experienced changes in RECO when LE was at 0 W m -2 , suggesting a more complex relationship between water balance and RECO at the sites. On hourly (1-2 hours) and diel (4 hours to 1.3 days) scales, GPP was the dominant predictor of NEE regardless of region or ecosystem type. On multiday (2.7-21.3 days), seasonal (42.7-341 days), and yearly scales, however, each variable showed indistinguishable magnitudes of predictability for NEE. The boreal sites saw more NEE predictability by biophysical variables relative to the temperate sites on long timescales (seasonal and yearly). Discharge did not have a drastically different NEE impact between sites and timescales (Figure 9).

Interannual Variability and Ecosystem Type
One of the objectives of this study was to determine differences in C flux IAV between biome types in temperate and boreal regions. When partitioning NEE fluxes into its components, GPP was greater than RECO for both ecosystem types, resulting in negative NEE (Figure 6). GPP was more of an NEE predictor than RECO on short timescales (hourly and diel) for both the temperate and boreal wetland and forest sites (Figure 9).
Overall, the forest sites were found to have higher GPP. This is likely due to larger photosynthetic capacity and greater leaf area index, promoting high photosynthetic uptake relative to the wetland sites. These results agree with a comprehensive study that found deciduous forests to be one of the largest C sinks across the conterminous United States (Xiao et al., 2011).
In addition to higher uptake, the deciduous broadleaf forest sites experienced larger GPP variability with respect to the wetland sites. Due to productivity's high reliance on meteorological variables, year-to-year fluctuations in local weather patterns play a vital role in GPP IAV. This finding suggests an ecosystem type dependence as opposed to a regional dependence for year-to-year carbon uptake potential.
Ecosystem respiration was also found to have high IAV across sites, with similar variabilities at US-ALQ, US-Los, and US-WCr and the highest occurring at US-Syv (Figure 3). The boreal sites experienced similar RECO IAV with respect to one another.
A study observing a temperate mixed forest in northeastern China concluded that only leaf respiration increased among other modes of respiration (i.e. stem and soil) when TA increased (Guan et al., 2006). Even though RECO is not dominant when partitioning NEE fluxes, leaf respiration's high temperature sensitivity has the potential to influence IAV and must be taken into account. Research on 15 European forests from 1996 to 1998 has shown variabilities in RECO to be the main determiner of NEE, indicating its important role in influencing IAV magnitudes (Valentini et al., 2000).
The C sink/source status of an ecosystem depends on the relative magnitudes of GPP and RECO, with the biome becoming a source or sink if RECO or GPP dominates respectively. In general, cumulative annual NEE was more variable from year-to-year at the forest sites than at the wetland sites regardless of region. GPP and RECO were also subject to variabilities (to a lesser extent) based on ecosystem characteristics. This result agrees with the finding that NEE IAV is larger than that of GPP and RECO (Xie et al., 2014).
Due to higher GPP potential at the forest sites relative to the wetland sites, US-Syv and US-WCr were found to have more years where the magnitude of negative NEE was greater at the end of the year than that of US-Los. The cumulative NEE at the forest sites also did not follow a specific temporal trend on decadal scales. At US-Los, however, C uptake generally increased as years passed, with the greatest C uptake occurring during 2018 and 2019 and the least during 2000 and 2005 (Figure 3). Similar results were found for the boreal sites.
An increase in cumulative yearly NEE at the wetland site may be indicative of the fertilization effect as a result of increased atmospheric CO 2 concentrations. More years of C fluxes need to be measured, however, to be certain.
In addition to uptake on yearly scales, the forest sites generally had greater C uptake on seasonal scales. During the growing season, US-Syv and US-WCr had negative NEE that was 1.5 and ~1.7 times more than US-Los on average. Year-to-year variabilities in growing season NEE were also larger at the forests versus the wetlands.
While NEE magnitude is helpful in determining the extent of C uptake, temporal variations in C fluxes must also be accounted for. One such measure is determining when NEE shifts from source to sink, indicating the first time when production rates surpass respiration rates. Due to significant increases in global temperatures in the past century, some ecosystems have been shown to respond with earlier leaf-out times and later senescence, therefore extending the growing season (Helfter et al., 2015). Only US-Los saw an overall earlier date of NEE sign change as years progressed. US-Syv and US-WCr were more variable in this regard. This was consistent with the boreal wetland and deciduous forest sites as well. CA-SCC (boreal evergreen needleleaf site), however, experienced similar yearly variabilities to that of the boreal wetland sites (CA-WP1 and CA-SCB).
Lengthier growing seasons have implications for higher rates of productivity. Therefore, consistently experiencing longer growing seasons as years progress may cause the wetland sites to reach GPP levels that are comparable to that of the deciduous forests. On the other hand, high IAV in NEE zero crossing day may prove inefficient for sustained C uptake on long timescales.

Drivers of Interannual Variability and Implications
Ecosystems owe their year-to-year C variability to external factors such as PPFD, TA, and water availability. Our results agree well with the notion that increases in PPFD and TA promote higher ecosystem productivity in both wetlands and forests (Froelich et al., 2015;Helfter et al., 2015). Therefore, longer periods of warmer temperatures and larger areas of foliage are amongst the most ideal conditions for enhanced GPP. It is important to note that linear regression has obvious shortcomings as an analysis tool, in that it is blind to non-linear relationships. Therefore, a lack of linear correlation with water parameters (i.e. precipitation and discharge) indicates a more sophisticated relationship between water availability and NEE (i.e. through changes in soil moisture (Xu et al., 2004)).
Precipitation's influence on C fluxes at US-Los, US-Syv, and US-WCr were observed to be modest when compared to those of other variables. When partitioned into daytime and nighttime sections, precipitation had more explanatory power for changes in NEE across all temperate sites, while GPP was not affected (Figure 7). RECO at US-Los was found to have a higher reliance on precipitation relative to US-Syv and US-WCr during the daytime, while nighttime influence on respiration was similar across all sites. This suggests that daytime C accumulated via photosynthesis (GPP) is being offset by enhanced RECO in the presence of precipitation. Alternatively, nighttime precipitation was observed to lead to more positive NEE values whereas no correlation was found during the daytime. Out of the temperate sites, NEE at US-Los was influenced most by this variable.
One study concluded that variability in net C balance depends more on PPFD than precipitation (Xie et al., 2014), indicating how different biomes respond differently to these drivers. As a result, for the sites considered, increased precipitation can be regarded as a way of advancing these biomes towards becoming a net source of C.
Due to precipitation's small influence on GPP, this leads NEE to be strongly characterized through precipitation (affecting RECO) and monthly PPFD (affecting GPP). Therefore, we would expect wet periods to lead to higher RECO and gradually lead these sites to becoming net sources of C. On the other hand, more sunlight exposure would promote higher rates of GPP and lead to a higher likelihood of the sites becoming a net C sink. These results were found at sites where zero precipitation was recorded 80%, 50%, and 94% out of all recorded days at US-Los, US-Syv, and US-WCr respectively at half-hourly scale. Hence, more studies need to be conducted at sites where precipitation rates are higher to determine the significance of this relationship across biomes. The boreal sites also saw increases in respiration rates during higher precipitation events (except for CA-Oas during the daytime). This may indicate that the RECO-P relationship is not regionally dependent, therefore pointing to a similarity between the temperate and boreal biomes considered in this study.
It is important to note that linear regressions have limitations when determining these relationships. A study has found that nonlinear fits may better capture year-to-year variabilities, although this was with regards to precipitation and aboveground net primary production (Knapp et al., 2016).
LE also had a significant impact on C fluxes. Increasing GPP with LE suggests higher rates of photosynthetic production with higher evapotranspiration magnitudes. Relating to what was found with precipitation, this creates a "tug-of-war" between drivers. Precipitation enhances RECO while LE promotes GPP. GPP was found to be less correlated with LE at US-ALQ whereas most other sites experienced correlation coefficients within 0.6 and 0.9. A relatively low correlation of LE with RECO across all temperate sites and years led NEE to decrease dramatically with increases in LE, with NEE at US-ALQ having the smallest correlation coefficient. This suggests a high dependence of C IAV on the water balance of each ecosystem. When the ecosystems lose water to the atmosphere, they are more prone to become net C sinks. Precipitation, however, promotes these biomes to become a net C source. No significant differences were found in evapotranspiration as a driver for C IAV when comparing temperate forests and wetlands. A similar relationship was found for the boreal sites regardless of ecosystem type. Furthermore, a study in China comparing C and water fluxes also found similar results, indicating that the C and water cycles are tightly coupled (Xiao et al., 2013).
In general, TA acted as a way of driving up rates of ecosystem production and respiration. Therefore, magnitudes of NEE remained largely unchanged with increasing temperature. RECO's exponential increase and GPP's quadratic correlation with increasing TA was found not to have significant IAV (Figure 8). This indicates that to maximize productivity, these biomes must be exposed to large amounts of solar radiation and also reside in warm climates. On the other hand, a warm climate also promotes respiration, and is maximized when combined with high amounts of precipitation. A study on multiple ecosystem types in China also found an exponential increase in respiration with increases in TA (Yu et al., 2012). Furthermore, a study found GPP at multiple deciduous forests to be negatively correlated with annual water balance (Law et al., 2002). Lastly, a synthesis report including 49 sites concluded NEE to be correlated linearly with mean annual temperature and had a logarithmic correlation with precipitation. It also determined PPFD and temperature to be the most significant predictors of annual NEE which is in agreement with our results (Kato and Tang, 2008).
Precipitation data was sporadic during the course of the study period and therefore we would expect each biome, on average, to be dominant in GPP versus RECO due to the greater abundance of sunlight relative to precipitation on yearly scales. During most of the years, on average, NEE was indeed found to be negative for all sites, giving each biome a net C sink status. Our findings agree with a study measuring terrestrial ecosystem respiration sensitivity to TA across 60 diverse FLUXNET sites, where no difference was found among biomes (Mahecha et al., 2010).

Study Limitations
As mentioned previously, most sites did not have stream gages located directly at the tower locations. Therefore, remote discharge data was used as a proxy for on-site discharge, which may introduce uncertainties when compared to NEE responses. The largest uncertainties may come from the Bear River (US-Los), Cisco Lake Outlet (US-Syv), and County Trunk Highway M (US-ALQ) stream gages since all were positioned 3 km or more from their respective Ameriflux towers.
Data inavailability was also a limitation in this study. For precipitation, US-Syv had the highest data gaps (2008-2013; 2019-2020) ( Table A. (Table A.

5).
As a result, when comparing C flux and micro-meteorological variables, this data gap problem was exacerbated. In some cases, this has inevitably caused an inability to conduct analysis on a continuous set of data on a yearly scale, resulting in some years not being considered entirely. This was also the case at the boreal sites. We believe this to be one of the main sources of uncertainty in this study.

Implications for Natural Climate Solutions
As mentioned in section 1.8, NCS is one of the most cost effective ways to increase C sequestration to limit global temperature increases due to greenhouse gas emissions. To fully take advantage of NCS's mitigation potential, ecosystem conservation and restoration practices must be implemented such that C uptake is maximized on decadal scales. To achieve this, each practice must be tailored with both ecosystem type and its corresponding C drivers in mind.
Increasing C storage requires increases in GPP while also reducing RECO. According to this study, the most effective way to promote increases in GPP across ecosystem types is through an increase of PPFD. Therefore, NCS efforts must implement afforestation and reforestation to increase leaf radiation interception such that PPFD-use is maximized. However, it is important to note that photosynthetic uptake has temporal limitations, peaking only during the summer months. On the other hand, reports of a lengthening growing season due to warmer temperatures may allow for significant C uptake as soon as early spring that lasts as late as mid-fall (Helfter et al., 2015).
To maximize C uptake and reduce risks of sequestration reversal, restoration efforts must also limit ecosystem RECO. This study shows RECO to be somewhat modulated by precipitation. NEE also became more positive with increased precipitation, indicating these biomes to likely become net C sources during high precipitation events. The dependency of the hydrologic cycle on RECO agrees well with a previous study on these sites (Sulman et al., 2009).
This correlation also goes beyond forests and wetlands. For instance, a mixed grass prairie in Wyoming was found to experience lower and higher respiration rates during decreases and increases in summer precipitation respectively (Chimner and Welker, 2005). This renders precipitation to potentially be a universal driver of RECO regardless of ecosystem type. However, other studies have shown that the magnitude of precipitation is also relevant, even causing net C uptake at a shortgrass steppe when total precipitation is above a certain threshold (Parton et al., 2012). This opens the possibility of the wetlands and forests to experience net uptake during periods of high rainfall (>10 mm day -1 ), but this cannot be determined with certainty due to limited precipitation data at the study sites.
To effectively implement restoration strategies to minimize RECO, one must understand the underlying mechanisms behind the RECO-precipitation relationship. High amounts of soil water has been found to promote microbial decomposition of organic matter and consequently, the release of CO 2 . Furthermore, high soil water content may also facilitate plant root/shoot growth which is associated with higher respiration rates (Chimner and Welker, 2005). It has also been documented that heterotrophic respiration is enhanced during small precipitation events (<10 mm day -1 ) (Parton et al., 2012).
Therefore, restoration and conservation efforts as part of NCS must prioritize controlled drainage after minor precipitation events to reduce RECO. This reduction will cause higher GPP, resulting in more negative NEE. Methods such as tile drainage may be used to reduce soil water content. We believe these efforts will create optimal conditions for C sequestration at these sites.

Conclusion
Here we have observed that two wetland and forest ecosystems spread across northern Wisconsin and Michigan's UP differ drastically in terms of C IAV. Results show higher C IAV magnitude in the forests relative to the wetlands which has implications for ecosystem restoration and conservation prioritization for NCS. A mix of boreal wetland and forest sites showed similar results, revealing small dependency of these processes on latitudinal location.
GPP and RECO were found to have a quadratic and positive exponential relationship with increases in TA, respectively across all sites. In addition, PPFD was found to be highly correlated with productivity (mostly during the summer months) with little to no correlation with RECO. On the other hand, respiration rates increased during relatively high precipitation events with little to no correlation with GPP. This indicates both PPFD and precipitation to be major contributors to C IAV. Year-to-year variability in GPP was also found to contribute more to NEE IAV than yearly variabilities in RECO.
Mutual information analysis revealed NEE at the boreal sites to have more dependency on biophysical variables relative to the temperate sites on seasonal and yearly scales. This suggests that the extent of influence on C fluxes may be more regionally-dependent on long time scales. Moreover, NEE was more dependent on GPP on short timescales (hourly and diel) and dependent on GPP and RECO somewhat equally on multiday, seasonal, and yearly scales.
Due to large NEE IAV at the forest sites, we believe wetlands to be the most reliant biome for long-term C sequestration. Although the potential for large uptake is higher at the forest sites, the high variability is a cause for high uncertainty. On the other hand, while the wetlands were characterized by lower C uptake, their tendency for more constrained C flux magnitudes from year to year leads to a more reliable C sink in the long-term.
Due to the spatially limited nature of this study, more work needs to be conducted at other biomes, especially those in the tropics, to determine whether similar results hold in other regions.  Table 2). Figures 8,9,10,11 from top to bottom display US-ALQ, US-Los, US-Syv, and US-WCr respectively.  Table 3.