Trends in high northern latitude soil freeze and thaw cycles from 1988 to 2002

[ 1 ] In boreal and tundra ecosystems the freeze state of soils limits rates of photosynthesis and respiration. Here we develop a technique to identify the timing of freeze and thaw transitions of high northern latitude land areas using satellite data from the Scanning Multichannel Microwave Radiometer (SMMR) and Special Sensor Microwave/Imager (SSM/I). Our results indicate that in Eurasia there was a trend toward earlier thaw dates in tundra ( (cid:1) 3.3 ± 1.8 days/decade) and larch biomes ( (cid:1) 4.5 ± 1.8 days/decade) over the period 1988–2002. In North America there was a trend toward later freeze dates in evergreen conifer forests by 3.1 ± 1.2 days/decade that led, in part, to a lengthening of the growing season by 5.1 ± 2.9 days/decade. The growing season length in North American tundra increased by 5.4 ± 3.1 days/decade. Despite the trend toward earlier thaw dates in Eurasian larch forests, the growing season length did not increase because of parallel changes in timing of the fall freeze ( (cid:1) 5.4 ± 2.1 days/decade), which led to a forward shift of the growing season. Thaw timing was negatively correlated with surface air temperatures in the spring, whereas freeze timing was positively correlated with surface air temperatures in the fall, suggesting that surface air temperature is one of several factors that determines the timing of soil thaw and freeze. The high spatial resolution, frequent temporal coverage, and duration of the SMMR and SSM/I satellite records makes them suitable for rigorous time series analysis and change detection in northern terrestrial ecosystems.


Introduction
[2] In high northern latitude regions, excluding the North Atlantic and Greenland, mean annual surface air temperature increased at rates of up to $0.7°C per decade during the latter half of the 20th century, with the largest increase occurring during winter and spring months [Chapman and Walsh, 1993;Hansen et al., 1999Hansen et al., , 2001. Concurrently, sea ice extent and thickness decreased in the northern hemisphere [Cavalieri et al., 1984;Chapman and Walsh, 1993;Cavalieri et al., 1997;Vinnikov et al., 1999], possibly due to increases in the length of the melting season [Laxon et al., 2003]. Observational studies of snow cover in Alaska find a trend toward earlier spring snowmelt [Stone et al., 2002], and satellite measurements indicate that snow cover area decreased over the period 1972 to 1988 in both North America and Eurasia [Serreze et al., 2000].
[3] The implications of these changes for the net carbon balance of northern ecosystems remains uncertain because many of these trends may simultaneously increase rates of ecosystem respiration and photosynthesis . Soils in northern biomes are carbon rich [Dixon et al., 1994] and lie within the region expected to experience the most extreme temperature changes due to global warming over the next century [Intergovernmental Panel on Climate Change (IPCC), 2001]. Soil carbon degradation is sensitive to changes in temperature and moisture content, thus the potential for increasing temperatures to promote carbon loss from soils in high latitude ecosystems is large [Schimel and Clein, 1996;Goulden et al., 1998;Hobbie et al., 2000;Mikan et al., 2002;Dioumaeva et al., 2003].
[4] In addition to impacts on soil respiration, increasing air temperatures may also lengthen the growing season and enhance rates of net primary production (NPP) in tundra and boreal biomes. Several studies of trends in satellite-derived Normalized Difference Vegetation Index (NDVI) report both an increase in growing season length and a ''greening'' of high northern latitudes [Myneni et al., 1997[Myneni et al., , 1998Tucker et al., 2001;Zhou et al., 2001;Hicke et al., 2002;Shabanov et al., 2002], especially in Eurasia [Zhou et al., 2001;Bogaert et al., 2002]. In a recent study, the increase in growing season length (1981 -1999) was attributed to both a trend toward earlier spring leaf out (6 ± 2 days in Eurasia and 8 ± 4 days in North America) and later fall senescence (11 ± 3 days in Eurasia and 4 ± 3 days in North America) leading to an overall increase in the growing season length by 19 -24 days in Eurasia and 4 -12 days in North America [Zhou et al., 2001].
[5] Examination of trends in atmospheric CO 2 records provides additional evidence for widespread changes in northern ecosystem function. The amplitude of the CO 2 cycle increased by 40% at Point Barrow, AK, and by 20% at Mauna Loa, HI over the period 1960-1993[Keeling et al., 1996. These trends are probably caused by a combination of increases in spring and summer NPP [Keeling et al., 1996;Randerson et al., 1999], increasing rates of soil respiration during fall, winter and spring [Chapin et al., 1996;Goulden et al., 1998], and increasing levels of natural and human induced disturbance [Zimov et al., 2001]. Concurrently, the period of the CO 2 record shifted forward 7 days from 1960 to 1993 [Keeling et al., 1996] suggesting a trend toward earlier spring thaw and forward shift in the timing of CO 2 drawdown by terrestrial ecosystems.
[6] Although NDVI provides a means for evaluating regional trends in leaf area (and thus NPP), and atmospheric CO 2 records allow for analysis of trends in net ecosystem production (NEP), no comparable constraint is available for estimating trends in biome-level soil respiration. Soil respiration persists at temperatures below 0°C, but rates typically decrease by over an order of magnitude at temperatures below freezing [Goulden et al., 1998;Dioumaeva et al., 2003]. In contrast, at temperatures above 0°C, respiration is exponentially related to soil temperature [Lloyd and Taylor, 1994;Raich and Potter, 1995;Raich et al., 2002;Dioumaeva et al., 2003].
[7] As proposed by Running et al. [1999], a global satellite product that identified soil freeze and thaw boundaries would be of broad use in determining fluxes of both soil respiration and photosynthesis in northern biomes. Here we present a method to estimate the annual timing of freeze and thaw for land north of 45°N using passive microwave satellite observations. We used data from the Nimbus-7 Scanning Multichannel Microwave Radiometer (SMMR), which operated from late 1978 to mid 1987, and the Defense Meteorological Satellite Program (DMSP) Special Sensor Microwave/Imager (SSM/I), which operated from mid 1987 to present. Our algorithm uses an iterative least squares characterization of the difference between the brightness temperature at 37 and 19 (or 18) GHz to identify the transition from frozen to thawed soil. Using this algorithm we estimate the date of thaw and freeze for each 1°Â 1°grid cell north of 45°N. As part of our analysis, we assessed regional and biome level trends in soil freeze and thaw since 1988, and analyzed interannual variability since 1979. We found a trend toward earlier spring thaw in Eurasian tundra and larch biomes, a trend toward earlier fall freeze in Eurasian larch forests and a trend toward later fall freeze in North American evergreen conifer forests. The growing season length increased significantly in North American evergreen conifer and tundra biomes, but did not change in Eurasian biomes due to the parallel shift of the spring thaw and fall freeze to an earlier date.
[9] In our study we used the Equal Area Scalable Earth Grid (EASE-Grid) daily brightness temperature product that is publicly available from the National Snow and Ice Data Center (NSIDC) [Armstrong et al., 2003;Knowles et al., 2002]. The EASE-Grid is a global equal area projection that was developed at the NSIDC with a spatial resolution of 25 km Â 25 km for all channels. Data are available twice daily (ascending and descending tracks) for the SSM/I and twice every other day for the SMMR instrument. The length and temporal frequency of the data set makes it suitable for rigorous time series analysis and change detection. The radiometer channels included 6.6, 10.7, 18, 21 and 37 GHz for the SMMR, and 19, 22, 37, and 85 GHz for the SSM/I. Both horizontally (H) and vertically (V) polarized radiation were measured for all These data continue to be collected. channels except 22 GHz. In our analysis, we used the 18 and 37 GHz channels from SSMR and the 19 and 37 GHz channels from SSM/I, and considered both horizontally and vertically polarized data from ascending and descending satellite tracks. There is little atmospheric interference in this region of the spectrum, except for a water absorption band in the 22 GHz channel [Chahine et al., 1983]. Clouds are thought to have minimal impact on the channels we used, and Judge et al. [1997] determined that water vapor and oxygen corrections did not change their freeze-thaw classification results using the 18, 19 and 37 GHz channels. We therefore made no atmospheric corrections in our data processing. The EASE-Grid data were compiled for the entire period between 1978 and 2002 and were processed through several stages before development and implementation of the algorithm.

Data Resampling
[10] Our first step of data processing was to average the daily EASE-Grid data to 6-day averages. Six-day averages were necessary because the swath width of the SMMR and SSM/I instruments was not wide enough to yield complete coverage of the region north of 45°N daily. SMMR achieved complete coverage of our study area every 6 days, and SSM/I achieved complete coverage every 4 days. To maintain consistency between satellite records, we averaged both the SMMR and SSM/I data to 6-day means for our analysis. We then aggregated the global EASE-Grid data to 1°Â 1°grid cells, removing EASE-pixels with high water content and complex topography for a 'cleaner' signal prior to this step (as described below). These processing steps were performed on both H and V polarized signals from both ascending and descending satellite tracks, giving us 4 independent time series for each 1°Â 1°grid cell that spanned 1 January 1979 to 31 December 2002 (24 years).

Theory of Algorithm
[11] Brightness temperature (T b ) measured by a radiometer is the product of surface temperature (T s ) and emissivity (e) of a surface: where e depends on frequency (n), polarization (H or V), incidence angle of the satellite observation, surface structure (roughness or vegetation), and moisture. Brightness temperature measurements from microwave radiometers are sensitive to surface moisture because of the relatively high dielectric constant of water. Surface roughness and vegetation cover modulate this sensitivity. For our analysis, we used the difference between the 37 and 18 GHz brightness temperatures for SMMR and the difference between the 37 and 19 GHz brightness temperatures for SSM/I as indicators of soil freeze state.
[12] The difference between these two channels (DT b = T b(37) À T b(18/19) ) reflects the difference in emissivity (De = e 37 À e 18/19 ) that is much larger for frozen than for thawed soils ]. The SMMR and SSM/I instruments maintained fixed incidence angles, thus by using different combinations of channels (frequency) and polarization we can estimate changes in T s and e. The DT b signal is dominated by changes in surface temperature and the difference between the emissivity at 37 and 18 or 19 GHz (De).
As the soil surface thaws, De tends toward zero for pixels with low water content ]. T s remains constant between two channels measured simultaneously, leading to a plateau shaped DT b curve driven solely by changes in De at the freeze thaw boundary. Our algorithm, which is described below, was designed to detect the date of thaw and freeze by identifying the edges of this annual plateau in the DT b curve for each pixel. We define the growing season length (gsl) as the difference between the day of freeze and the day of thaw.
[13] In another study, an algorithm designed to detect snow cover from the SMMR and SSM/I radiometer data used a combination of the 19, 22, 37 and 85 GHz channels . A threshold difference between 37 and 19 GHz and a decision tree approach for the 85 -22 GHz difference were used to classify snow covered pixels ]. Although there is some correlation between snowmelt and soil thaw, they can be offset by many weeks as observed by Kimball et al. [2001]. Although we ignore snow cover in this analysis, it remains an important question that needs to be addressed further.

Freeze-Thaw Algorithm
[14] First, we fit a spline to the DT b curve. The algorithm involved a two-step process of fitting a line by ordinary least squares to the plateau of the spline-fitted DT b curve. The first step was to select a horizontal line that minimized the sum of squares of the residual between the line and the plateau of the curve. This was the first estimate (e 1 ) of the plateau level. The second step eliminated points deviating from e 1 by more than one standard deviation (s). The process of line fitting was repeated again, and a second line was chosen to represent the plateau (e 2 ). s was then recalculated and the thaw and freeze date were chosen by finding the first and last points where the spline-fitted DT b curve crossed the horizontal line 0.5 Á s below e 2 ( Figure 1). The factor 0.5 Á s desensitized the algorithm to high frequency variability in the DT b curve, which may be caused by repeated freeze-thaw events in a grid cell.
[15] To validate our algorithm and its application to a variety of surface characteristics, we compared the estimated date of thaw and freeze to that observed directly from soil temperature records at Boreas North Study Area, Manitoba (56°N, 98°W) [Dunn and Wofsy, 2003], Delta Junction, AK (64°N, 145°W) (H. P. Liu et al., manuscript in preparation, 2004), Toolik Lake, AK (68°N, 149°W) [Shaver and Laundre, 2003], and Happy Valley, AK (69°N, 148°W) [Hinkel, 2002] (Figures 2 and 3). Notice that the algorithm performs better at Happy Valley than at Toolik Lake ( Figure 3). Soil temperatures were measured at 8 cm depth in Happy Valley, and at 20 cm depth in Toolik Lake, which may explain the difference in algorithm performance. The microwave emissions measured by the SMMR and SSM/I originate in the upper few ($5) cm of the soil surface, depending on surface wetness. We would expect to see some offset between the Toolik Lake measured and satellite estimated dates of freeze and thaw because of this depth offset.
[16] For the global analysis we computed the date of freeze and thaw for each of the 4 time series (horizontally or vertically polarized and ascending or descending tracks) at each 1°Â 1°pixel, giving us 4 independent estimates of soil thaw and freeze date. Despite the fact that the data were originally aggregated to 6-day intervals, we were able to detect trends at a much finer resolution because we fit a spline to the DT b curve with a short temporal resolution. We averaged the 4 estimates of date of freeze and date of thaw for each pixel and aggregated the results to biome and regional scales for the period 1988-2002. Vegetation classifications were based on DeFries et al. [1998]. The biomes analyzed here were evergreen conifer forest and tundra in both North America and Eurasia and larch forest in Eurasia (Table 2).

Pixels With Open Water
[17] EASE 25 km Â 25 km pixels with greater than 40% open water cover were excluded from our analysis. We chose the 40% cutoff as a compromise between reducing the failure rate of our algorithm and retaining a high data volume in all of our major biomes (Table 2). With a 40% cutoff, we eliminated virtually all algorithm failures, but kept at least 86% of all pixels in each biome (Table 2).
[18] We generated a water fraction map by aggregating the pixels of a 1 km Â 1 km resolution land cover map developed by Hansen et al. [2000] to the EASE pixels and 1°Â 1°grid cells. Each EASE pixel was assigned a water percentage based on the number of 1 km land cover pixels identified as water. 1°Â 1°grid cells were constructed from EASE-Grid pixels that had less than 40% water cover. In addition, any 1°Â 1°grid cell with greater than 40% water cover was eliminated from the analysis.
[19] The physical basis of this exclusion is that De does not tend toward zero as water cover of a pixel increases ]. Instead, De becomes larger, destroying the plateau observed in pixels with low water content. Figure 4 illustrates the difference between average signals for each biome with varying water contents. At the 40% level, the shape of the DT b curve degraded to a degree that prevented accurate detection of freeze and thaw ( Figure 4) without eliminating significant portions of any major biome. Because our algorithm was based on the shape of the curve rather than the absolute value of the curve, we can safely include pixels with some water, but must remove those with high water content (Figure 4). A total of 7880 (6.7%) pixels at the EASE-Grid resolution were excluded from our analysis, and the maximum portion of any biome excluded was 14% in the tundra of North America (Table 2).

High Elevation Pixels
[20] Pixels with high elevations had a very high rate of algorithm failure compared with pixels at lower elevations. This is probably due to complex topography and in some instances, snow and ice cover that persisted over the entire year in some pixels. With a threshold of <1250 m, we excluded a maximum of 7% of any biome (Table 2). We therefore filtered out all pixels above 1250 m elevation, as identified by the United States Geological Survey's (USGS) GTOPO30 elevation map (Land Processes Distributed Active Archive Center (LP-DAAC), USGS EROS Data Center, http: //edcdaac.usgs.gov, 2003) averaged to the EASE grid resolution. A total of 2018 EASE pixels were excluded from the three biomes, or 5% of the total (Table 2).

Cross Platform Changes
[21] The SMMR and SSM/I satellites have differences in their viewing angles, overpass time, and frequency (channels) that can introduce bias in the analysis of time series extending across the period covered by the two satellites. The viewing angles of SMMR and SSM/I are 50.2°and 53.1°, respectively. The warm and cold overpass of SMMR occur at noon (ascending) and midnight (descending), whereas for SSM/I the timing is shifted to 0600 hours and 1800 hours for the F8 satellite, and 1800 and 0600 hours for the F11 and F13 satellites. These differences are particularly important because of the variation in freezing and thawing caused by diurnal variations in temperature. Furthermore, SMMR collected data at 18 GHz, and SSM/I collected data at 19 GHz. There are small differences between atmospheric absorption and surface emissivity at these frequencies, leading to differences in measured brightness temperature.
[22] We observed a large offset between the SMMR (1979 -1986) and SSM/I (1988SSM/I ( -2002 freeze-thaw record. Derksen and Walker [2003] suggest that there is a systematic bias between the measured brightness temperatures of SMMR and SSM/I, primarily caused by differences in overpass time and instruments, but the short ($10 day) period of overlap between the two EASE-Grid records does not permit the development of correction terms that would be both spatially and seasonally valid. For this reason, we did not combine the records to study trends over the complete 24-year measurement period. Instead, we generated averages and trends based only on the 1988 -2002 period (15 years). We did, however, utilize the SMMR record to compare interannual variability of freeze/thaw transitions to that of the surface temperature record during 1979-1986.    Figure 3 from the soil temperature data (pluses) and from our satellite analysis (diamonds). For the in situ soil temperature we defined the date of thaw as the earliest date that the soil temperature is greater than 0°C for 3 consecutive days and date of freeze as the latest date that soil temperature is less than or equal to 0°C for 3 consecutive days. Each biome is divided by continent, North America (NA) and Eurasia (EA). Water fraction distribution describes the percent of pixels in each biome with a given water content. Pixels with greater than 40% water cover or greater than 1250 m elevation were excluded from our analysis. b For all areas north of 45°N excluding Greenland. calculated a correlation between seasonal surface air temperature anomalies (with a base period of 1951 -1980) and the freeze-thaw anomalies. We divided the complete 1979 to 2002 record into two periods, 1979 to 1986 for the SMMR and 1988 -2002 for the SSM/I. 1987 was excluded from the analysis because the SMMR to SSM/I shift occurred in the middle of the year. The mean for each satellite period (SMMR or SSM/I) was used separately to compute annual anomalies in thaw and freeze transitions. Thaw transitions were compared with surface air temperatures during the spring, while freeze transitions were compared with surface air temperatures during the fall. For the spring, we averaged air temperature anomalies from April and May. For the fall, we averaged air temperature anomalies from September and October for tundra biomes, and from October and November in evergreen conifer and larch biomes.

Analysis Approach
[24] We produced annual dates of freeze, thaw and length of growing season for each 1°Â 1°grid cell from 1979 to 2002. We used these data to examine the spatial and temporal patterns of spring thaw, fall freeze, and growing season length. To compare with NDVI trends and other metrics of ecosystem function, we averaged our results to biome levels and separated each biome by continent. For both the global and biome level analysis, we computed mean statistics along with trends over the period from 1988 through 2002. Interannual variability was analyzed for the entire satellite record (1979 -2002) and compared to variability in surface air temperature over the same period. Spatial patterns of the freeze-thaw anomalies were also compared to the spatial patterns of the surface air temperature record.

Mean Patterns
[25] In Eurasia, the timing of the mean thaw progressed from west to east, with a distinct transitional border at the Ural Mountains in Russia. The mean thaw date was delayed with increasing latitude. The latest mean thaw dates occurred in central Siberia along the coast of the Arctic Ocean (Figure 5a). In Eurasia, the date of thaw occurred first in evergreen conifer forests (Julian day 120 ± 25), followed by larch forests (131 ± 15) and tundra (151 ± 21) ( Table 3). In North America, the area west of the Hudson Bay in Canada had the latest mean thaw (Figure 5a). North American evergreen conifer thawed first (132 ± 20), followed by tundra (157 ± 21) ( Table 3).
[26] The timing of the freeze occurred in the reverse order of the thaw across Eurasia and North America. In Eurasia, freeze timing occurred later from east to west, again with a distinct border at the Ural Mountains. The freeze occurred earlier with increasing latitude (Figure 5b). In Eurasia the date of the freeze occurred first in tundra (272 ± 17) followed by larch (276 ± 16) and evergreen conifer forests (289 ± 23) ( Table 3). The Brooks Range in Alaska was the first part of North America to freeze, with eastern Canada following much later (Figure 5b). In North America the date of freeze occurred first in tundra (285 ± 20) then in evergreen conifer forests (293 ± 23).
[27] The mean gsl in Eurasia decreased from west to east and from south to north. It was shortest in central Siberia on the Arctic coast (Figure 5c). The gsl in Eurasia was shortest in tundra biomes (121 ± 30 days) followed by larch forest (146 ± 26) and evergreen conifer forest (170 ± 42) biomes (Table 3). In North America, the Brooks Range in Alaska had the shortest gsl. Tundra had a shorter gsl (128 ± 31) than evergreen conifer forests (160 ± 34) in North America (Table 3).

Trends
[28] In all biomes north of 45°N there was either a trend toward earlier spring thaw, or no significant change over the 15-year period from 1988 through 2002 (Figure 5d and , and 60.01-100% (dash-dot lines) water cover. The plateau shape degrades as the water content of a pixel increases, especially in the tundra ecosystem. This is due to an increase in De for water versus thawed soil. Table 3). There was a trend toward earlier spring thaw throughout much of Eurasia east of the Ural Mountains, but a trend toward later thaw in areas west of the Urals (Figure 5d). Eurasian larch forest thaw exhibited the greatest rate of change (À4.5 ± 1.8 days/decade, where the negative sign convention indicates earlier thaw) followed by tundra (À3.3 ± 1.8 days/decade) (Figure 6a and Table 3).
In North America, Alaska demonstrated a mixed response with some positive and negative trends in the date of thaw (Figure 5d).
[29] Trends in the timing of the fall freeze were substantially different for North America and Eurasia. In Eurasia there was a trend toward earlier fall freeze in the larch biome by À5.4 ± 2.1 days/decade, whereas in the North American evergreen conifer biome, there was a 3.1 ± 1.2 days/decade trend toward later fall freeze (Figure 6e and Table 3). While there were some areas of Eurasia with a trend toward later fall freeze, such as Scandinavia and the Ural Mountain region of Russia, most of eastern Eurasia experienced a trend toward earlier fall freeze similar in magnitude to the trend in spring thaw (Figure 5e). In contrast, there was a trend toward later fall freeze in many areas of North America (Figure 6e and Table 3). The area east of the Canadian Rockies between 50°N and 60°N experienced the greatest trend toward later fall freeze (Figure 5e).
[30] In Eurasian larch forests the parallel movement of spring thaw and fall freeze led to a shift in the phase of the growing season to an earlier period but resulted in no net change in gsl (Figure 6f and Table 3). In contrast, North America's growing season length increased significantly for both evergreen conifer forests (5.1 ± 2.1 days/decade) and tundra (5.4 ± 3.1 days/decade). These changes were driven by a trend toward both earlier spring thaw and later fall freeze. Some areas in Alaska showed a negative trend in growing season length, but nearly all of Canada showed a large positive trend (Figure 5f).

Interannual Variability and Comparison With Surface Air Temperature Records
[31] Mean surface air temperature patterns were similar to patterns in freeze and thaw (Table 4). There was a strong negative correlation between spring air temperature anoma-lies and thaw dates, and a strong positive correlation between fall air temperature anomalies and freeze dates (Table 5 and Figure 7). The single largest anomaly in the timing of the fall freeze occurred in 1998, where the timing of the freeze was $10 days earlier than average in all of the Eurasian biomes (Figure 7). This was matched in the surface air temperature record by a strong negative anomaly during fall (Figure 7). In addition, the spatial pattern of the 1998 freeze anomaly had a very similar structure to that of the surface air temperature anomaly during this time (Figure 8). The freeze date captured widespread cooling in central Eurasia, and more localized warming in northern Canada. Because of the high resolution of the SMMR and SSM/I data we identified an area with slightly positive freeze anomalies over the Ural Mountains in Russia, a feature that the surface air temperature record could not resolve.

Discussion
[32] Passive microwave measurements reflect changes in the emissivity of a pixel, which incorporates information from both the soil surface and vegetation. It provides a direct measurement of the freeze state of a pixel that is independent of surface air temperature measurements. Air temperature and freeze-thaw timing can, in fact, be decoupled by several factors such as ground cover (i.e., insulating mosses versus bare ground), vegetation that shades the soil surface, or water content of vegetation and soils [Baldocchi et al., 2000]. Vegetation and ground cover may be particularly important in mature evergreen conifer forests, which have a perennial canopy and thick insulating moss ground cover [Beringer et al., 2001;Chambers and Chapin, 2003].
[33] As a partial validation of our predicted trends, it was important to assess whether our results reproduced expected mean patterns of thaw, freeze, and growing season length. Generally, we expected that evergreen conifer ecosystems would thaw earlier and freeze later than tundra and that the date of thaw would increase with increasing latitude, and date of freeze would decrease with increasing latitude based on surface soil temperature measurements (Figures 2 and 3). Our biome-level results support these assumptions (Table 3),  and give confidence that our algorithm captured the largescale trends and interannual variability in the freeze-thaw record.
[34] The trends observed in this study are interesting to consider in the context of ecosystem function and the surface energy budget of northern biomes. In Eurasia, there was little net change in the growing season length, but there was a trend toward earlier spring thaw in both tundra and larch biomes and a trend toward earlier fall freeze in the larch biome. Because northern ecosystems are light and temperature limited, an earlier spring thaw may allow plants to become photosynthetically active at earlier periods, when  0.2 ± 4.2 À9.5 ± 6.5 À6 ± 2.6 0.5 ± 0.9 À0.9 ± 2.4 À0.6 ± 1.4 All evergreen conifer 2.8 ± 5.0 À3.4 ± 6.7 À1.8 ± 5.2 0.1 ± 0.8 0.6 ± 1.1 0.1 ± 0.7 All tundra À2.6 ± 4.8 À0.9 ± 4.3 À5.5 ± 6.2 0.4 ± 0.9 0.9 ± 1.2 0.4 ± 1.0 North America e 0.7 ± 5.4 1.1 ± 5.5 À1.6 ± 5.4 0.0 ± 1.0 1.3 ± .8 0.8 ± 0.7 Eurasia e 3.7 ± 6.5 À0.3 ± 6.9 À1.6 ± 7.1 0.5 ± 0.8 À0. light levels are relatively high. Earlier spring thaw may also increase soil respiration. Although there was also a trend toward earlier fall freeze, it may have impacted plant productivity less because of light limitation during this time of year. An earlier fall freeze could, however, reduce fall fluxes of CO 2 from soils. In contrast, the trend toward later fall freeze in North America could increase fall efflux of CO 2 from soils. Earlier spring thaw also allows transpiration  (1979 -1986) p-value SSM/I (1988SSM/I ( -2002 p-value SMMR (1979 -1986) p-value SSM/I (1988SSM/I ( -2002 p-value  Figure 7. Freeze-thaw anomalies (solid line) and monthly surface temperature anomalies (dotted line) for 1979 -1986 and 1988 -2002. The spring thaw is negatively correlated with spring temperatures, and the fall freeze is correlated with fall temperatures. See Table 5 for correlation coefficients.
to occur earlier in the year, increasing early season latent heat fluxes. Because much of the latent heat fluxes in evergreen conifer ecosystems are modulated by mosses and trees with shallow rooting systems [Baldocchi et al., 2000], changing freeze thaw patterns may have a disproportionately large effect on the energy budget of these biomes.
[35] Our freeze-thaw trends show substantial differences with reported NDVI time series, although the different time periods for the published NDVI analyses make it challenging to construct a direct comparison. A recent NDVI analysis by Zhou et al. [2001] over the 1981 to 1999 period concluded that there was a trend of 3.1 days/decade toward earlier spring leaf-out in Eurasia and 4.2 days/decade in North America. Zhou et al. also found a trend of 5.8 days/decade toward later fall senescence in Eurasia and 2.1 days/decade in North America. They concluded that the growing season length had increased by 9.5 days/decade in Eurasia and 6.3 days/decade in North America. Our freeze-thaw analysis over the period 1988 -2002 yields similar results for the trend toward an earlier spring, but our results differ significantly for the fall freeze. Where the NDVI analysis showed a large trend toward later fall in Eurasia (5.8 days/decade), we detected a 5.4 days/decade trend toward earlier fall freeze in the larch biome. We found no net change in the gsl of Eurasia, and a 5 days/decade lengthening trend in North American tundra and evergreen conifer biomes. A hypothesis to reconcile these results could be that in Eurasia colder temperatures were accompanied by later snow cover. In evergreen conifer and tundra ecosystems, a change in NDVI during the fall may be due in part to changes in the timing of snow cover and its impacts on moss reflectivity.
[36] High latitude atmospheric CO 2 observations suggest that there is a trend toward earlier spring uptake of CO 2 by 3.3 days/decade from 1975-1996 [Keeling et al., 1996]. There is, however, less evidence in the CO 2 record of significant changes in the zero crossing time during the fall [Keeling et al., 1996;Randerson et al., 1997Randerson et al., , 1999. This may be due to the opposite trends in the timing of soil freeze (and thus rates of soil respiration) in Eurasia and North America that would cancel out their net contribution to the seasonal cycle of atmospheric CO 2 .
[37] Limitations on the application of SMMR and SSM/I data to the study of freeze thaw transitions stems from several sources. The large pixel size prohibited including small-scale topography in our analysis, a variable that is significant for soil respiration fluxes [Goulden et al., 1998;Dioumaeva et al., 2003]. In addition, the effect of vegetation and canopy cover on our algorithm was not explicitly examined in this study. However, grouping data by vegetation type and region allowed for an internally consistent analysis of trends across years.
[38] We did not account for snow cover of pixels, although soil freeze and thaw dates chosen with this method compare well to field data in areas with seasonal snow cover ( Figure 3). The DT b curve primarily reflects changes in surface emissivity and has a strong correlation with the soil freeze state. At the four field sites we compared our algorithm to, our algorithm worked reasonably well for most years but does not capture the interannual variability perfectly at all sites ( Figure 3). Despite these limitations, this approach provides a new metric of soil freeze and thaw timing since 1988.

Conclusions
[39] The freeze-thaw product developed here represents a new and independent measurement of soil freeze state. We found that there was a trend toward earlier spring thaw in Figure 8. Maps of anomaly of 1998 (a) freeze and (b) fall temperature. Fall temperature data are the average of September, October, and November surface air temperature anomalies in 1998 from the GISS 1200 km smoothed data set [Hansen et al., 2001].
Eurasian larch and tundra biomes, and a trend toward earlier fall freeze larch forests in Eurasia but there was a trend toward later fall in evergreen conifer forests in North America. The net effect on gsl was minimal in Eurasia, but was on the order of 5 days/decade in North American tundra and evergreen conifer biomes. Although the growing season length did not change in Eurasia, the shift in the timing of freeze and thaw may have important implications for ecosystem carbon fluxes because of increased net primary productivity in spring, when light is plentiful. In contrast, the extended growing season in North America may increase plant productivity and soil respiration in spring, and increase soil respiration in fall.