Abstract

A procedure is described to construct time series of regional surface temperatures and is then applied to interior central California stations to test the hypothesis that century-scale trend differences between irrigated and nonirrigated regions may be identified. The procedure requires documentation of every point in time at which a discontinuity in a station record may have occurred through (a) the examination of metadata forms (e.g., station moves) and (b) simple statistical tests. From this “homogeneous segments” of temperature records for each station are defined. Biases are determined for each segment relative to all others through a method employing mathematical graph theory. The debiased segments are then merged, forming a complete regional time series. Time series of daily maximum and minimum temperatures for stations in the irrigated San Joaquin Valley (Valley) and nearby nonirrigated Sierra Nevada (Sierra) were generated for 1910–2003. Results show that twentieth-century Valley minimum temperatures are warming at a highly significant rate in all seasons, being greatest in summer and fall (> +0.25°C decade−1). The Valley trend of annual mean temperatures is +0.07° ± 0.07°C decade−1. Sierra summer and fall minimum temperatures appear to be cooling, but at a less significant rate, while the trend of annual mean Sierra temperatures is an unremarkable −0.02° ± 0.10°C decade−1. A working hypothesis is that the relative positive trends in Valley minus Sierra minima (>0.4°C decade−1 for summer and fall) are related to the altered surface environment brought about by the growth of irrigated agriculture, essentially changing a high-albedo desert into a darker, moister, vegetated plain.

1. Introduction

Long-term changes in climate response variables, such as surface temperature, are important to quantify as climate forcing parameters change. Because some of these changing forcing parameters are induced by human activity (e.g., enhanced greenhouse gas concentrations and land use changes), it is necessary to know precisely what the magnitudes of responses are so that attribution of the causes may be possible.

For surface temperature in a region a few hundred kilometers across, the long-term changes we seek to measure are small (e.g., order 0.1°C decade−1), which is of the same magnitude as the errors pervasive with the raw measurements (e.g., Christy 2002). Errors arise due to changes in location, instrument, or observational procedures, to name a few possible sources. Several studies report on the development of adjustments in order to reduce the errors such changes produce (see Folland et al. 2001). Common adjustments include the removal of the artificial effects that arise from urbanization or other forms of land use changes (Kukla et al. 1986; Karl et al. 1988; Peterson et al. 1998a).

How are adjustments determined for these and other changes that affect station data? Christy (2002) and Gallo (2005) provide evidence that generic adjustments for station moves (i.e., simple functions based on altitude or latitude) or instrument changes are inappropriate for regional time series. When a station is moved or a change occurs locally, the unique microclimate to which the thermometer responds becomes a potential source of significant new bias. Such changes can be unpredictable as shown in Christy (2002) for northern Alabama (13 stations in a region 80 km across) where magnitudes and signs of the biases associated with changes are not systematic. Indeed, Gallo (2005) found that presumed latitude and altitude relationships (i.e., that temperatures decrease as stations are more poleward and/or higher) were not sustained in paired comparisons of Climate Reference Network stations. These argue for site-specific adjustments, especially for time series representing an area less than 300 km in diameter where every station is important. [See Peterson et al. (1998a,b) for adjustments on larger spatial scales.] Adjustments for agricultural (land use) changes are relevant here and we now mention some key results.

Bonan (2001) found that warm season TMax cropland temperatures declined relative to forest lands in a comparison between the Midwest (cropland) and the Northeast (forest). A careful examination of the results of Kalnay and Cai (2003, their Figs. 2 and 3) over central California indicate that National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalyses (NNR) diagnose cooler surface trends than those of the actual station data. The implication is that since NNR diagnose the surface temperature based primarily on deeper atmospheric-layer temperatures, then changing surface effects undetected by NNR (e.g., urbanization or irrigation) have apparently acted to increase the observed surface trends beyond that expected from changes in the overlying circulation.

Small et al. (2001) documented the surface temperature changes due to the desiccation of the Aral Sea and found that increases in the warm season diurnal temperature range were generally dominated by increases in TMax. Balling et al. (1998) reported that overgrazing in the Sonoran Desert increased warm season TMax. In a general result, Gallo et al. (1999) demonstrated that as the land classification indicated more development, the greater was the negative trend in diurnal temperature range. The common theme of these results is that as water presence increases (more agriculture, more plants, etc.) TMax surface temperatures decline, especially in the warm season, resulting in a lessening of the diurnal magnitude (Karl et al. 1993). Our intent here is to examine such possible surface temperature effects of irrigation in central California, but this first requires the construction of a dataset suitable for climate analysis.

In the following sections we shall describe the study area, the selection of stations, and the collection of metadata used to determine the stations' discontinuities. We will then present the mathematical process by which generalized time series were constructed. Following this we describe several error analyses performed on the various time series. Finally we shall present the results for the San Joaquin Valley and adjacent highlands and discuss possible hypotheses to explain the results.

2. Study area

We have chosen to study the central San Joaquin Valley of California and the adjacent highlands. The 100-km-wide valley is oriented SE to NW with the Sierra Nevada rising to a 4000-m crest about 100 km NE of and parallel to the eastern edge of the valley. The elevation of the valley floor in our study region ranges from 30 to 140 m. The Coast Range assumes the western border but reaches elevations of only 1500 m along a few SE to NW trending ridges. This range is tall enough to affect a rainshadow on the western side of the valley. In this Mediterranean climate, over 90% of the precipitation falls from November to April, with annual totals in the valley between 10 cm on the western side and 30 cm on the eastern side, while May through October is essentially precipitation free. Orography enhances the precipitation in the Sierra Nevada where annual totals, much falling as snow, can exceed 125 cm of liquid equivalent.

Prior to the late nineteenth century, the valley was a vast plain, called by some the Serengeti of North America, watered in the spring by flooding rivers from the snowmelt of the Sierra Nevada. During summer and fall the valley becomes desert like, with clear skies and average daily maxima above 37°C in July.

Since the late nineteenth century, agricultural interests sought to bring into phase the spring runoff and the abundant summer and fall sunshine to optimize the growing potential for literally hundreds of varieties of crops. Initially, small diversion projects redirected the flow of the few rivers with summer runoff onto nearby fields for on-demand irrigation. During the 1940s–1960s, the Bureau of Reclamation, Army Corps of Engineers, State of California, and local water use associations built major reservoirs, to hold back the spring runoff, and distribution systems (canals) for conveyance on demand. These simple gravity-delivery systems carry water to locations over 100 km from the impoundments.

In the six counties composing this study area (Mariposa, Merced, Madera, Fresno, Kings, and Tulare) irrigated land amounted to 242 000 ha in 1899. These farms were mostly supplied by localized, private diversion projects. By 1982, the total had risen to 1 250 000 ha (over 12 000 km2) Fig. 1, with Fresno County alone leading the nation with $4.7 billion in agricultural revenues in 2004 (Fresno Bee, 27 April 2005). The crops on these lands require different amounts of water, but most irrigated land receives 1 m of irrigation-supplied water in the course of the year, most in the dry months. Thus, the land surface and near-surface atmosphere have experienced a significant change from a sunny, half-year with dry, high-albedo surfaces, to much wetter and darker surfaces.

Fig. 1.

Land area on which irrigation was applied in five counties utilized in this study. Mariposa County had negligible land under irrigation

Fig. 1.

Land area on which irrigation was applied in five counties utilized in this study. Mariposa County had negligible land under irrigation

Does the dramatic change in San Joaquin Valley surface conditions create a response in near-surface air temperature that is measurable by the simple weather station instruments that monitored the valley throughout the twentieth century? To answer this question, we need a reasonably good time series of the valley temperature and a control case against which to test any hypothesis. For the former, we have developed a method to generate composite time series of weather station data. For the latter, we compare the valley results with those of the adjacent highlands. Our assumption here is that since the centroid of observations of the valley and mountain stations are separated by only 60 km horizontally and less than 1000 m vertically, the long-term climate trends should be very similar if no differential forcing develops. Note that we are dealing with surface temperatures only, not upper-air temperatures where trend differences may more easily occur (Folland et al. 2001). However, it is possible that some mesoscale climate change in one part of the region (valley) could impact the other so that our “control” case may not be completely independent. For example, in a very different climate regime, in which summer convection occurs, Chase et al. (1999) show evidence that Colorado high plains agriculture may have influenced summer convective events in the nearby mountains under upslope (easterly) conditions but with little impact on mountain temperatures.

3. Identifying and collecting information on station discontinuities

To determine whether this land use variation in space and time is important to temperature trends, we must first develop a dataset with sufficient precision to allow for discrimination of the effect. We accessed from the National Climatic Data Center (NCDC) daily maximum (TMax) and minimum (TMin) temperature data for all stations within the six counties of the central San Joaquin Valley (Mariposa, Merced, Madera, Fresno, Kings, and Tulare) from the valley eastward.1 We shall refer to the stations on the San Joaquin Valley floor as Valley and those in the adjacent Sierra Nevada as Sierra. Valley stations are those between the Coast Range and Sierra Nevada and generally less than 130-m elevation on the flat plain. Sierra stations are those east of the valley floor and generally above 130-m elevation. The lower-elevation foothills of the Sierra Nevada are characterized by evergreen black oak over grassland. Higher-elevation stations, above 1000 m, are generally in the yellow pine, fir, sequoia, and lodgepole pine tree range. Snow falls on these higher stations every year. We were able to utilize data from 18 Valley and 23 Sierra stations; the locations of the stations are shown in Fig. 2 and listed in Table 1.

Fig. 2.

Location of stations used in this study. Valley stations identified with crosses, Sierra stations with circles. The triangle is Darwin Glacier (Fig. 9)

Fig. 2.

Location of stations used in this study. Valley stations identified with crosses, Sierra stations with circles. The triangle is Darwin Glacier (Fig. 9)

Table 1.

Listing of stations utilized in the study

Listing of stations utilized in the study
Listing of stations utilized in the study

Through a project directed by NCDC, photographic images of metadata forms, which describe many aspects of the conditions and history of the stations, have been archived. In Fig. 3 we show a particularly useful form, 530–1, for North Fork Ranger Station that summarizes, up to the 1970s, many changes important for climate analysis. Note the comment on 18 June 1930 that “CRS [cotton region shelter] was moved 20 ft W to reduce effect of lawn sprinkling on readings.” Such events are potentially responsible for nonclimatic shifts in the temperature record.

Fig. 3.

Example of WB Form 530–1 describing the history of the station at North Fork Ranger Station beginning in March 1904. Note under “Remarks” the statement that the Cotton Region Shelter (CRS) was moved 20 ft W to reduce effect of lawn sprinkling

Fig. 3.

Example of WB Form 530–1 describing the history of the station at North Fork Ranger Station beginning in March 1904. Note under “Remarks” the statement that the Cotton Region Shelter (CRS) was moved 20 ft W to reduce effect of lawn sprinkling

We examined every form, about 1600 pages in all, for each of the stations in the region. The form identifier was recorded (e.g., 530–1, 531, 4005, 4029, B-44, etc.) and information that in some way might have bearing on the integrity of the time series was manually digitized. An example for a Valley station, Madera, is shown in Table 2. Once the information from all station forms was compiled, we determined those points in time that might be associated with a nonclimatic shift. In the case of Madera, we determined there were nine breakpoints (giving 10 homogeneous segments) based on seven station moves, a change in observation time, and an instrument change (Table 3).

Table 2.

Consolidated information for Madera (045233)

Consolidated information for Madera (045233)
Consolidated information for Madera (045233)
Table 3.

List of breakpoints determined from consolidated information of Table 2 

List of breakpoints determined from consolidated information of Table 2
List of breakpoints determined from consolidated information of Table 2

Though we read through every available form describing each station's potential changes, it was apparent, especially before 1930, that a few climatologically important breaks were not documented. The existence of undocumented breaks is made clear when examining, for example, the unadjusted (raw) temperature data for the higher Sierra stations (above 900 m) prior to 1940. In Fig. 4 we display the unadjusted, absolute seasonal average Tmin for December–February (DJF) of each station already separated by the known breakpoints. Notice the parallel movements of the time series after 1926, indicating that anomalies during that period are evidently well characterized. However, for the segment of Huntington Lake (Hunt2) we see a sharply warmer temperature in 1925 relative to all other stations. More difficult to detect in the figure, but still significant, is the apparent error in the first segment of California Hot Springs (CaHS1) between 1916 and 1917. In this instance, all other stations indicate DJF was significantly warmer in 1917 than 1916, while CaHS1 indicates the opposite.

Fig. 4.

Unadjusted DJF TMin data for Sierra stations above 900-m elevation. Time series have been subdivided into homogeneous segments according to metadata only. Label abbreviations are of stations identified in Table 1, and the concatenated numeral is the segment number. Note in particular the unusual behavior of CaHS1 and Hunt2 as described in the text

Fig. 4.

Unadjusted DJF TMin data for Sierra stations above 900-m elevation. Time series have been subdivided into homogeneous segments according to metadata only. Label abbreviations are of stations identified in Table 1, and the concatenated numeral is the segment number. Note in particular the unusual behavior of CaHS1 and Hunt2 as described in the text

To identify these likely but undocumented breakpoints we calculate the first derivative [or first difference; Peterson et al. (1998b)] of each segment's time series (Fig. 5). Here we see that after 1926, the first differences are tightly clustered, indicating strong agreement in the progression of anomalies. Notice the characteristic signal of a single-season rogue value (Hunt2 1925), with consecutive, oppositely signed outlier values. Earlier years are not quite as tightly clustered as the post-1926 data, and likely erroneous values do appear, as in the case of CaHS1 1917 (circled, filled triangle). In this sample, it also was determined that Yosemite Valley (Yose1, 1916) experienced a break. Similar analysis was performed on the Valley stations and these additional breakpoints were then added. In the analyses to follow, we will utilize data from 1910 onward as data prior to this date in the Sierra time series were dependent on very few stations, each with several breakpoints.

Fig. 5.

As in Fig. 4, except the quantity plotted is the first difference of the time series in Fig. 4. Note the consistent character of the variations after 1926. Circled values were identified as additional segment breakpoint events

Fig. 5.

As in Fig. 4, except the quantity plotted is the first difference of the time series in Fig. 4. Note the consistent character of the variations after 1926. Circled values were identified as additional segment breakpoint events

At this point we have each station divided into several segments, which we assume to be homogeneous in time. Remaining nonclimatic trends, as opposed to sudden shifts, are assumed to be small and random. However, we are aware that spurious trends on a few relatively long segments may have undue influence on the composite time series. We shall address this issue later.

4. Computation of biases for all possible segment pairs

In our approach, each segment is treated as an independent unit relative to the other segments. Our goal is to merge these segments into one time series for each season (four cases), time of day (two cases: TMax and TMin), and region (two categories: Valley and Sierra). To create any given time series (e.g., Sierra, JJA, Tmax) we must determine the bias of each individual segment in that categorical subset relative to all others. Once the magnitude of each segment bias is determined, each segment time series may be debiased and combined with all others to generate a merged, best-guess realization of the regional time series. (In an intermediate step we can also combine the segments for specific stations into an adjusted station time series.)

We assume that stations chosen for this analysis are sufficiently similar in climate characteristics (e.g., geography of elevation bands) so that their true long-term time series would all be very similar. In other words, we assume any single station, if providing perfect observations, would be a reasonable representative of any other station in the region and for the regional average.

To solve the problem before us we appeal to mathematical graph theory. Temperature segments and their overlaps are modeled as a directed graph. A directed graph is a finite set of points, called vertices, some of which are connected by directed lines, called edges. Each edge may have a set of values associated with it describing relationships between the vertices it connects. In our case, each segment represents a vertex and an edge connects two segments that overlap. If segment Si overlaps segment Sj, the vertex corresponding to Si is connected by a directed line to the vertex corresponding to Sj. With each directed edge are associated five statistical parameters:

  • the average of the daily temperature differences δij between the overlapping portions of the segments Si and Sj (δ of TMax and TMin are separately done),

  • the number N of daily differences,

  • the standard deviation σij of the differences,

  • the autocorrelation of the difference series, r1ij so that Neff = N(1 − r1ij)/(1 + r1ij), and

  • the standard error ɛij = σij/Neff of the differences.

To compute δij, the temperatures on concurrent days of Si and Sj are subtracted and all such differences for Si and Sj averaged. The Si to Sj statistics are identical to those of Sj to Si except for the algebraic sign of the bias.

The graph model provides a convenient way to estimate these statistics for two vertices that are not connected by an edge but by a succession of edges (a path). This is similar to estimating the hypothetical outcome of a game between two sports teams who did not play each other head to head by looking at how they fared against common opponents. If Si is connected to Sj and Sj to Sk, but not Si to Sk, then we can estimate δijk, σijk, and ɛijk for the SiSk combination as follows:

 
formula

If vertex Sm provides an alternate path from Si to Sk, we take the best estimate of the bias from Si to Sk to be the one given by the path with the smaller standard error:

 
formula

The approach shown in (1) can be generalized for paths of any length:

 
formula

If there are multiple paths of possibly varying lengths from Si to Sj, we consider the “least” path to be the one with the smallest composite standard error, even if the number of its edges exceeds the number of edges in other paths. Figure 6 illustrates these concepts with a simple example. To find the biases between every possible pair of vertices, we use Dijkstra's shortest-path algorithm (Dijkstra 1959, as implemented by Standish 1995), where for our purpose “shortest” means “least” in the sense just described.

Fig. 6.

A directed graph with vertices i, j, k, p, and q. There are two paths from vertex i, marked with an open circle, and vertex j, marked with an open square. One is through vertex k; the other through vertices p and q. The edges of the graph are labeled with the bias (δ) and the standard error (ɛ). Using the definitions and conventions of Eqs. (1)(3), we determine that the standard error of the path through vertex k is 0.33, whereas the standard error of the path through vertices p and q is 0.20. Hence, the least path from vertex i to vertex j is the one through vertices p and q, and we take the bias between vertices i and j to be the sum of the biases along the path through p and q, which is −0.14. For simplicity, we have not shown the edges as arrows but assume that when the direction of traversal between two vertices changes, the algebraic sign of δ is reversed

Fig. 6.

A directed graph with vertices i, j, k, p, and q. There are two paths from vertex i, marked with an open circle, and vertex j, marked with an open square. One is through vertex k; the other through vertices p and q. The edges of the graph are labeled with the bias (δ) and the standard error (ɛ). Using the definitions and conventions of Eqs. (1)(3), we determine that the standard error of the path through vertex k is 0.33, whereas the standard error of the path through vertices p and q is 0.20. Hence, the least path from vertex i to vertex j is the one through vertices p and q, and we take the bias between vertices i and j to be the sum of the biases along the path through p and q, which is −0.14. For simplicity, we have not shown the edges as arrows but assume that when the direction of traversal between two vertices changes, the algebraic sign of δ is reversed

Once the biases and standard errors have been determined for every possible pair of vertices, the results can be written as two n × n matrices, Δ, for the biases, and 𝗘 for the standard errors. Each row or column corresponds to a vertex (data segment). In Δ, the intersection of row i and column j is δij. Since δji = −δij for each possible pair, i and j, and δii = 0 for each i, Δ is antisymmetric. The bias δij is the additive adjustment for the temperatures of Sj to remove their bias relative to Si. In addition, 𝗘 is symmetric since ɛij = ɛji for each pair, i and j, while ɛii = 0 for each i. Finally, 𝗘 is irreducible since ɛij > 0, ij (see Keener 1993). We shall use this fact later.

To motivate the next step, we consider a very simple case. In the example, we will construct a single, consistent temperature time series for a single station whose temperature record can be subdivided into four homogeneous segments, S1, . . . , S4. Since these segments are derived from the same station, they are nonoverlapping. We assume we have applied the procedure described above, with the help of overlapping segments from other, nearby stations, to determine the biases, δij, i, j = 1, . . . , 4, for the possible pairings of these four segments. These biases form a square matrix Δ of order 4.

Let the adjusted time series be denoted by Sj + δij. Also, denote the operator for concatenating two nonoverlapping time series by ⊕. Then the time series S2, S3, and S4 can be brought into alignment with S1 as follows:

 
formula

This is one estimate of the total time series at the station. Three other estimates could be formed by choosing S2, S3, or S4 as the reference segment. Thus,

 
formula

It is important to realize that each of the δij were generated through different, least error paths, so that there is no constraint requiring all four estimates of the reconstructed time series to be identical. How do we know which realization is the “best” estimate?

One best estimate of the reconstructed time series for the station is simply the average of the four estimates shown in (4) and (5):

 
formula

Observe that the adjustments to each segment are just the averages of the biases in columns 1, 2, 3, and 4 of the bias matrix Δ for this station. We could generalize the adjustments by making them weighted averages, taking into account the possibility that some of the δij are more robust than others because, for instance, they are derived from sets of differences with smaller standard errors. We can also generalize the example to the case of n nonoverlapping segments at the station. In this case the adjustment aj to segment Sj is given by

 
formula

or more generally by

 
formula

where wi ≥ 0, i = 1, . . . , n, and Σni=1wi = 1. The n adjustments of Eqs. (7) or (8) form the bias adjustment vector a = (a1, . . . , an).

In Christy (2002)  a was determined by a cumulative procedure in which the columns of Δ were combined into a single column by a weighting scheme dependent upon the pooled estimate of the standard errors. In this study we use the adjustments as defined by (8) where the weights are computed in a different way.

To obtain the weights wi in Eq. (8), we rank the Si according to their ability to produce robust biases, that is, by the associated value of ɛij that represents each Si's ability to generate overlaps having differences with small errors. We can think of ɛij as being the “score” when Si “competes” against Sj. Lower scores represent less error and all possible pairs are contained in 𝗘. Let r be a ranking vector of the Si. Then we would expect the ranks of the segments to be proportional to their scores:

 
formula

where λ is the constant of proportionality. Equation (9) shows that r is an eigenvector of the matrix of standard errors associated with the eigenvalue λ. Because lower scores are more desirable, the segments with lower rankings are superior. Now we use the fact that 𝗘 has nonnegative entries. The Perron–Frobenius theorem (see Keener 1993) states that if 𝗘 is irreducible, r has strictly positive entries and λ is the largest eigenvalue of 𝗘 in absolute value.

To compute r we use the power method (Burden and Faires 1985). Since we would like for the larger weights to be applied to the segments producing the most reliable biases, we take the weights to be the reciprocals of the entries of r, normalized so that the sum of the scaled reciprocals is 1. Therefore, the ith weight is

 
formula

Once the weights are known and the bias adjustment vector a has been computed, the construction of the regional time series is straightforward. First, we debias each segment using the applicable entry from a. At this point, the data are still in daily resolution. To determine the temperature of the regional series for a particular day, we average the temperatures from every debiased segment that includes that day. From the resulting daily time series, we construct the seasonal means.

5. Error analysis

Before analyzing the results, we tested the probability distribution of each time series in a number of ways to understand the errors associated with the data and method. We will focus on the linear trend (least squares regression) because this metric is the most sensitive to changes in the procedures used and is a metric of interest for long-term changes.

a. Segment uncertainty

Do we have an adequate number of segments that are (a) consistent with each other and (b) temporally distributed in a manner that allows us to construct a robust, reproducible time series? Perhaps our set just happens to be pathologically arranged so that the final outcome contains significant error. In other words, perhaps our unique set of segments assigns great reliance on a few critical segments that, if not available or having large error, could lead to a very different solution.

1) Segment uncertainty: Random removal of segments

To test this aspect of the basic structure of our method, we randomly removed 20%, 15%, 10%, and 5% of the segments 1000 times and generated time series without them. We note that at 15% elimination, some of the 1000 trials could not be completed back to 1910; hence, we stopped at 20%. We then compared the median trend of each of these four reduced-segment trials with the trend of the full-segment dataset. The magnitudes of the median's deviations from the full-segment trend ranged from 0.002° to 0.073°C decade−1 in the 16 time series with a median of the median-deviation values of 0.017°C decade−1. For the extreme case (Sierra TMin for June–August (JJA)] the median trend for 20% missing trials was −0.172° compared with −0.245°C decade−1 for the full set of segments. We conclude here that for the most part, the basic character of each time series is reproducible from random subsets of the data. There were two specific segments, however, which exerted substantial impact and we shall examine them below.

2) Segment uncertainty: Nonrandom removal of segments

To further understand the uncertainty, we repeated the calculation of composite trends for all 16 cases with the removal of each segment, one at a time. There were 137 Sierra and 112 Valley segments to eliminate individually, so we generated as many new time series for each case. The trends of the resulting time series were determined and compared with the full-segment trends. In this way we are able to determine the impact of every single segment on the trend of the entire time series.

Not surprisingly, those cases with the largest variation in trends in the randomly reduced segment trials above also were characterized by one or two segments that, when eliminated, had a significant impact on the calculation of the trend. For example, in the case of Sierra TMin JJA we discovered that one segment (Huntington Lake 1938–70), when removed, shifted the trend to be more positive by 0.25°C decade−1. Similarly, a trend increase for Sierra TMin of 0.17°C decade−1 was caused when the same segment was eliminated from September–November (SON). These two cases were by far the most extreme examples of influence by a single segment on the computed time series trend. There was only one other segment that, when eliminated, shifted the trend by more than 0.1°C decade−1: Visalia 1927–64, which affected Valley TMin JJA and Valley TMin SON. (In these four extreme cases, the recalculation of the time series without these segments caused all their associated trends to be more positive. Thus, their impact on the final results of this study regarding the difference in Valley versus Sierra TMin trends is negligible.)

With a set of trends calculated with every segment removed individually, we employ a resampling, or jackknife, method for estimating the standard deviation of the trend deviations (our test statistic) from the full-segment trend for each of the 16 time series (von Storch and Zweirs 1999). We shall focus on the time series with the largest error magnitude, Sierra TMin JJA, recognizing that the remaining time series are characterized by error magnitudes that are smaller by a factors ranging from 1.5 to 6.

We initially found that it was unreasonable to assume that every segment was an equally contributing element to the jackknife error calculation. Many segments were relatively short and therefore had little impact no matter what their error might be and therefore skewed the resulting error range to be very small by implying a large sample. We therefore reduced our sample size to those 30 segments that exhibited the largest impact on the time series, about 0.001°C decade−1 or greater. This will reduce the number of segments and therefore will increase the magnitude of the error statistic. For Sierra TMin JJA, the 95% confidence interval (CI) was ±0.127°C decade−1 and the range for the other 15 time series was ±0.018 to ±0.085, with median ±0.041°C decade−1. These error values will be included in the results later.

b. A nonparametric error analysis

To calculate the magnitudes of the trend errors in a different way, we performed another test employing only those segments whose influence on the full segment trends was at least 0.02°C decade−1 (nine Valley and eight Sierra). Because these are likely, but not necessarily, contributing error to the full-segment results, we may determine a sense of their effect by removing them in all combinations.

We therefore create all possible combinations of these critical segments (512 for Valley and 256 for Mountain) and remove these and regenerate the time series for each case. We expect that with the removal of these subsets in all combinations, an improved time series is more likely to result from the median of the cases. Thus, we include in our results presented later the median of these reduced-segment trials as a separate trend estimate.

c. Temporal sampling error

Finally, there is the issue of temporal sampling error. We are examining various seasonal time series of length 94 yr. Sampling error provides information on how much confidence one may have in the notion that the present 94-yr period is truly representative of any 94-yr period randomly selected from a large population of time series experiencing the same climate conditions. (There will be temporal sampling error even though the measurements may be perfect.) Results indicate that in only 3 of the 24 cases (16 seasonal time series and 8 difference time series) did the 95% trend sampling error exceed ±0.06°C decade−1. To determine the total confidence interval or error range of the trends, we shall combine the segment uncertainty errors and temporal sampling errors to create the error bars on the trends of the original, full-segment method. (Again, the nonparametric results will be presented as a separate trend alongside.)

In summary, our view is that the trend values of the full-segment experiment represent the best guess of the actual trends, as they were produced by the algorithm that searches for the least error path. However, we also view the difference in trend values between the full-segment method and the median of the reduced-segment trials as an indication of the likely direction of error, but perhaps not the magnitude, to which this method might be susceptible. (The magnitude of trend differences in the reduced trials is not truly representative of an unbiased value as we preselected the segments based on their large, individual impact.) The largest magnitude of this difference (median of reduced-segment minus full-segment result) occurs, as expected, for Sierra TMin JJA (+0.151°C decade−1), with the range of the others being −0.091° to +0.081°C decade−1). We conclude that the error bars displayed later for the eight difference time series are reasonable as they capture error magnitudes from all the tests performed.

6. Results

Figures 7a–e display the time series of the two elevation strata and in Fig. 8 the trend values for each time series and the trends of the difference time series are shown. It is immediately apparent that Valley TMin time series are significantly positive in all seasons and especially so in JJA and SON. Sierra trends are small (Sierra TMax trends near zero) though Sierra JJA TMin shows a tendency for significant cooling.

Fig. 7.

(a)–(d) Time series of seasonal anomalies of Valley and Sierra TMin and TMax through 2003. (e) Time series of annual anomalies of mean temperature for Valley and Sierra stations

Fig. 7.

(a)–(d) Time series of seasonal anomalies of Valley and Sierra TMin and TMax through 2003. (e) Time series of annual anomalies of mean temperature for Valley and Sierra stations

Fig. 8.

Seasonal trends for each time series of Figs. 7a–d and trends of the difference time series. Solid gray bars represent trends from the original, full-segment calculations with error bars deduced from the combination of the 30-segment jackknife method and temporal sampling error. The lightly hashed bars represent the median trends of the nonparametric experiments

Fig. 8.

Seasonal trends for each time series of Figs. 7a–d and trends of the difference time series. Solid gray bars represent trends from the original, full-segment calculations with error bars deduced from the combination of the 30-segment jackknife method and temporal sampling error. The lightly hashed bars represent the median trends of the nonparametric experiments

The trends of the differences between the Valley and the Sierra time series are also significant and provide the key results for this research (Fig. 8, right). The TMax differences are greatest in JJA with the Valley trend being more negative than that of the Sierra stations. The most striking result is the highly significant relative positive Valley TMin trends, peaking in JJA, at over +0.5°C decade−1 in the trend of the differences. This amounts to a relative warming of 5°C in Valley JJA TMin versus the Sierra stations over the 94-yr period.

The correlations of seasonal anomalies between Valley and Sierra time series for the 1910–2003 period are given in Table 4. The TMax anomalies are highly correlated in March–May (MAM) and SON. The TMax correlation is lower in DJF when multiday periods of inversion events occur that are characterized by valley fog and low cloudiness (“high” fog or tule fog), decoupling the Valley and Sierra temperatures. The JJA TMax correlation is low due to smaller variance magnitudes in this season and the relatively greater magnitude of the variance contained in the trend differences (negative for Valley, neutral for Sierra), which begin to overwhelm the interannual variabilty. The TMin anomalies are poorly correlated in JJA and SON as, again, the magnitude of the variance carried by the differing long-term trends approaches that of the small interannual fluctuations. However, when detrended the time series are highly correlated (Table 4).

Table 4.

Correlation (Pearson product moment) of seasonal anomalies between Valley and Sierra time series, 1910–2003 (detrended in parentheses)

Correlation (Pearson product moment) of seasonal anomalies between Valley and Sierra time series, 1910–2003 (detrended in parentheses)
Correlation (Pearson product moment) of seasonal anomalies between Valley and Sierra time series, 1910–2003 (detrended in parentheses)

a. Hypotheses to explain results

The purpose of this paper is to present a methodology that generates a regional temperature record for climate applications using a technique that adjusts for the numerous discontinuities in the individual station records. In viewing the results we see the robust significance of the Valley versus Sierra trend differences, which beg explanations. The results are consistent with the hypothesis that the massive growth of irrigation in the San Joaquin Valley has impacted long-term trends and an inspection of Fig. 1 and Figs. 7a,b shows the irrigated acreage versus temperature correspondence, especially in the warmer seasons. One would expect the seasonal cycle of trend differences to coincide with irrigation deliveries, which are largest in JJA; furthermore, irrigation would be expected to have its largest impact in JJA due to the phase of solar forcing. Likewise, our largest trend differences between Valley and Sierra for both TMax and TMin were in JJA.

If our results accurately reflect the near-surface air temperature changes over the past century, we may hypothesize here about the causes for the dramatic warming observed in Valley TMin relative to nearby Sierra. (Though our hypothesis focuses on irrigation as an obvious cause, we have not ruled out the effects of subtle circulation changes that might differentially affect the Valley and Sierra stations.) Agricultural development with irrigation is the one most likely to do so, since it is prevalent around the Valley stations but not the Sierra stations, in the following ways:

  1. enhanced greenhouse warming due to increased water vapor concentrations in the atmosphere from evaporation and evapotranspiration into the valley boundary layer,

  2. enhanced nighttime downward infrared flux due to swelling of aerosols as humidity increases toward morning, and

  3. enhanced nighttime sensible heat flux from the surface due to the increased heat capacity of the vegetation and moist soil, both of which more readily absorb and store solar energy due to lower albedo, relative to the original desert surface, and a larger heat storage capacity due to existing water mass.

We have calculated trends of the 3-h, synoptic dewpoint observations for the period 1950–2001 at Fresno/Yosemite International Airport (Table 5). In all synoptic times, the dewpoint trends are positive during this period, with general maxima in the warm season afternoons (> +0.4°C decade−1) while Valley JJA TMax fell –0.26°C decade−1. Daytime moistening and dry-bulb cooling are thus observed and support the results of studies cited earlier.

Table 5.

Decadal trends (°C decade−1) of hourly dewpoint temperatures for 1950–2001 at Fresno/Yosemite International Airport. Columns are local hour [Pacific standard time (PST)]

Decadal trends (°C decade−1) of hourly dewpoint temperatures for 1950–2001 at Fresno/Yosemite International Airport. Columns are local hour [Pacific standard time (PST)]
Decadal trends (°C decade−1) of hourly dewpoint temperatures for 1950–2001 at Fresno/Yosemite International Airport. Columns are local hour [Pacific standard time (PST)]

It is important to note that the time series of surface moisture for Fresno contains some uncertainties. For instance, the instruments that measure moisture content have changed from manual psychrometers to analog-to-digital hygrometers, and observations of moisture were taken at different times in different periods. In other words, the time series is not truly homogeneous. In addition, the increase in dewpoint temperatures (if any) may be near the absolute precision of a single instrumental time series as deployed here (a shortcoming of our homogeneous segment technique is that it may not be applied to a single station). Furthermore, it is difficult to conclude from observations at one station, which is situated in a metropolitan area that has grown 10-fold over this period to >500 000 population, that large-scale agricultural irrigation is the direct cause of any moisture increase. In short, moisture trends presented here should be viewed with caution. (As a side note, it is evident that with increasing background dewpoints, the efficiency of household evaporative cooling systems, widespread when the lead author grew up in Fresno, will have declined.)

Our initial look at the three hypotheses indicates the most likely explanation is option 3 above, though all may contribute to some extent. Regarding the enhanced water vapor greenhouse effect (option 1), the additional moisture seems to be of a small enough amount in a relatively shallow layer that there would be little impact, and certainly not as much as 5°C. For option 2, the valley generally does not experience relative humidities greater than 80% in the warmer seasons, which is the general threshold at which aerosols begin to swell. At present, therefore, we hypothesize that the significant increases in Valley TMin are related to the darkening and moistening of the formerly dry, high-albedo desert surface (option 3). The darker surface allows for more absorption of solar energy while the additional water mass in plant material and wet ground increases the heat capacity, providing a daytime repository of energy to be lost via sensible heat flux at night. In future work, it is our intent to test these hypotheses with a high-resolution, boundary layer model to quantify the possible impacts of these irrigation-related perturbations.

b. Related findings

We note that Cayan et al. (2001) examined hydrologic data beginning in 1950 and discovered a trend toward earlier spring snowmelt, or peak discharge dates, in the Sierra Nevada. The Sierra trends for MAM TMax (+0.18°C decade−1) and TMin (+0.14°C decade−1) for 1950–2003 in our dataset are highly consistent with Cayan et al.'s result. The MAM Sierra trend in TMean (+0.16°C decade−1) is the most positive of all seasons since 1950. However, as implied in Fig. 8, once the entire century is considered, the Sierra MAM mean trend (+0.01°C decade−1) is not significantly different from zero.

Even though our century-scale Sierra trends are fairly unremarkable, there is clear indication of change in this region. Figure 9 displays photographs taken in 1908 and 2003 from the same view of Darwin Glacier (37.1702°N, 118.6771°W) near the crest of the Sierra Nevada in Fresno County. The elevation of the upper ice line is approximately 3960 m. The reduction in extent is obvious and indicates that the conditions that support this glacier have changed during the twentieth century. As there is no evidence of significant long-term temperature changes in our Sierra time series, though trends may be different at 4000 m for unknown reasons, other factors are likely involved, for example, decreases in cloudiness or precipitation. However, some proxy indicators suggest the twentieth century was wetter than previous centuries (Graumlich 1993) while measurements show a general increase since 1900 (USGCRP 2000). In any case, the causes of glacial mass balance changes in this region are evidently more complex than can be inferred from simple temperature records.

Fig. 9.

Southward looking views of Darwin Glacier (37.1702°N, 118.6771°W) near the crest of the Sierra Nevada, Fresno County, CA, taken in (top) 1908 and (bottom) 2003. The elevation of the upper ice line is approximately 3960 m. The 1908 photo was merged from two USGS file photographs by H. Basagic, Portland State University. The 2003 photo was taken by N. L. Stephenson, Research Ecologist, USGS Western Ecological Research Center, Three Rivers, CA

Fig. 9.

Southward looking views of Darwin Glacier (37.1702°N, 118.6771°W) near the crest of the Sierra Nevada, Fresno County, CA, taken in (top) 1908 and (bottom) 2003. The elevation of the upper ice line is approximately 3960 m. The 1908 photo was merged from two USGS file photographs by H. Basagic, Portland State University. The 2003 photo was taken by N. L. Stephenson, Research Ecologist, USGS Western Ecological Research Center, Three Rivers, CA

Finally, we note that our TMean trends of both Valley and Sierra composites are less positive than implied by assessments based on larger-scale analyses for this region (e.g., USGCRP 2000; Folland et al. 2001). Indeed, our trends are in closer agreement with unforced model hindcasts of twentieth-century climate than with human-enhanced forcing (e.g., Tett et al. 2002). A comparison of the time series of six stations common to this study and version 1 of the United States Historical Climatology Network (USHCNv1) dataset [Fresno, Hanford, Merced, Visalia, Yosemite Valley, and Lemon Cove; Karl et al. (1990)] indicate the composite USHCNv1 TMean trend is 0.10°C decade−1 more positive than calculated here. [A similar comparison for the JJA TMax trend in North Alabama also indicates a 0.10°C decade−1 more positive trend in USHCNv1 than Christy (2002).] This suggests that utilizing as much data as possible, and applying site-specific adjustments, may yield lower rates of surface temperature increases, though our sample here is small. In any case, the highly significant warming in Valley TMin does suggest that land use changes have had a substantial impact on the local climate.

7. Conclusions

We have demonstrated a technique to create regionally consistent time series of temperature data based on the assumptions that we are able to identify all significant discontinuities in station records and that the stations are situated in a climatologically homogeneous region. We composited the temperature records of 18 stations in the San Joaquin Valley of central California and 23 stations in the adjacent Sierra Nevada into, respectively, two regional time series for each season. Our analysis of trends begins in 1910 though records are available in earlier years from fewer stations. Our results indicate that the central San Joaquin Valley has experienced a significant rise of minimum temperatures (∼3°C in JJA and SON), a rise that is not detectable in the adjacent Sierra Nevada. Our working hypothesis is that the rapid valley warming is caused by the massive growth in irrigated agriculture. Such human engineering of the environment has changed a high-albedo desert into a darker, moister, vegetated plain, thus altering the surface energy balance in a way we suggest has created the results found in this study. Additionally, if these results are confirmed, the lack of long-term warming in the generally undeveloped Sierra Nevada (annual mean trend, 1910–2003, −0.02° ± 0.1°C decade–1) coupled with significant, nighttime-only warming in the valley, suggests a regional inconsistency compared with twentieth-century simulations of climate forced by human influences other than land use changes.

Acknowledgments

Support for Christy and Norris was provided by the National Science Foundation. We thank Pam Waisanen, USGS National Center for Earth Resources Observation and Science, for assistance with irrigation data displayed in Fig. 1. We thank Nathaniel Stephenson (USGS, Ash Mountain, California) for contributions to the discussion on glacier variations. The staff of NCDC, Tom Ross and Joe Elms in particular, was very helpful in deciphering some confusing station information and in accessing some of the more obscure data files.

REFERENCES

REFERENCES
Balling
Jr.,
R. C.
,
J. M.
Klopatek
,
M. L.
Hildebrandt
,
C. K.
Moritz
, and
C. J.
Watts
,
1998
:
Impacts of land degradation on historical temperature records from the Sonoran Desert.
Climatic Change
,
40
,
669
681
.
Bonan
,
G. B.
,
2001
:
Observational evidence for reduction of daily maximum temperature by croplands in the midwest United States.
J. Climate
,
14
,
2430
2442
.
Burden
,
R. L.
, and
J. D.
Faires
,
1985
:
Numerical Analysis.
3d ed. Prindle, Weber, and Schmidt, 676 pp
.
Cayan
,
D. R.
,
S. A.
Kammerdiener
,
M. D.
Dettinger
,
J. M.
Caprio
, and
D. H.
Peterson
,
2001
:
Changes in the onset of spring in the western United States.
Bull. Amer. Meteor. Soc
,
82
,
399
415
.
Chase
,
T. N.
,
R. A.
Pielke
Sr.
,
T. G. F.
Kittel
,
J. S.
Barron
, and
T. J.
Stohlgren
,
1999
:
Potential impacts on Colorado Rocky Mountain weather due to land use changes on the adjacent Great Plains.
J. Geophys. Res
,
104
,
16673
16690
.
Christy
,
J. R.
,
2002
:
When was the hottest summer? A state climatologist struggles for an answer.
Bull. Amer. Meteor. Soc
,
83
,
723
734
.
Dijkstra
,
E. W.
,
1959
:
A note on two problems in connection with graphs.
Numer. Math
,
1
,
269
271
.
Folland
,
C. K.
, and
Coauthors
,
2001
:
Observed climate variability and change.
Climate Change 2001: The Scientific Basis, J. T. Houghton et al. Eds., Cambridge University Press, 99–181
.
Gallo
,
K. P.
,
2005
:
Evaluation of temperature differences for paired stations of the Climate Reference Network.
J. Climate
,
18
,
1631
1638
.
Gallo
,
K. P.
,
T. W.
Owen
,
D. R.
Easterling
, and
P. F.
Jamason
,
1999
:
Temperature trends of the U.S. Historical Climatology Network based on satellite-designated land-use/land cover.
J. Climate
,
12
,
1344
1348
.
Graumlich
,
L. J.
,
1993
:
A 1000-year record of temperature and precipitation in the Sierra Nevada.
Quat. Res
,
39
,
249
255
.
Kalnay
,
E.
, and
M.
Cai
,
2003
:
Impact of urbanization and land-use change on climate.
Nature
,
423
,
528
531
.
Karl
,
T. R.
,
H. F.
Diaz
, and
G.
Kukla
,
1988
:
Urbanization: Its detection and effect in the United States climate record.
J. Climate
,
1
,
1099
1123
.
Karl
,
T. R.
,
C. N.
Williams
Jr.
,
F. T.
Quinlan
, and
T. A.
Boden
,
1990
:
United States Historical Climatology Network (HCN) serial temperature and precipitation data. Environmental Science Division Publ. 3404, Carbon Dioxide Information and Analysis Center, Oak Ridge National Laboratory, Oak Ridge, TN, 389 pp
.
Karl
,
T. R.
, and
Coauthors
,
1993
:
A new perspective on recent global warming: Asymmetric trends of daily maximum and minimum temperature.
Bull. Amer. Meteor. Soc
,
74
,
1007
1023
.
Keener
,
J. P.
,
1993
:
The Perron–Frobenius theorem and the ranking of football teams.
SIAM Rev
,
35
,
80
93
.
Kukla
,
G.
,
J.
Gavin
, and
T. R.
Karl
,
1986
:
Urban warming.
J. Climate Appl. Meteor
,
25
,
1265
1270
.
Peterson
,
T. C.
, and
Coauthors
,
1998a
:
Homogeneity adjustments of in situ atmospheric climate data: A review.
Int. J. Climatol
,
18
,
1493
1517
.
Peterson
,
T. C.
,
T. R.
Karl
,
P. F.
Jamason
,
R.
Knight
, and
D. R.
Easterling
,
1998b
:
First difference method: Maximizing station density for the calculation of long-term global temperature change.
J. Geophys. Res
,
103
,
25967
25974
.
Small
,
E. E.
,
L. C.
Sloan
, and
R.
Nychka
,
2001
:
Changes in surface air temperature caused by desiccation of the Aral Sea.
J. Climate
,
14
,
284
299
.
Standish
,
T. A.
,
1995
:
Data Structures, Algorithms, and Software Principles in C.
Addison-Wesley, 748 pp
.
Tett
,
S. F. B.
, and
Coauthors
,
2002
:
Estimation of natural and anthropogenic contributions to twentieth century temperature change.
J. Geophys. Res
,
107
.
4306, doi:10.1029/2000JD000028
.
USGCRP
,
2000
:
Climate Change Impacts on the United States: The Potential Consequences of Climate Variability and Change (Overview).
National Assessment Synthesis Team, U.S. Global Change Research Program, Washington, DC, 154 pp
.
von Storch
,
H.
, and
F. W.
Zweirs
,
1999
:
Statistical Analysis in Climate Research.
Cambridge University Press, 484 pp
.

Footnotes

* Visiting scientist at National Center for Earth Resources Observation and Science, U.S. Geological Survey, Sioux Falls, South Dakota

Corresponding author address: John R. Christy, ESSC, University of Alabama in Huntsville, Cramer Hall, Huntsville, AL 35899. Email: christy@nsstc.uah.edu

1

One station utilized in the study is just outside these counties. Kern River Power House resides in Kern County, 4 km from the southern edge of Tulare County.