Like all reanalysis efforts, the Modern Era Retrospective-Analysis for Research and Applications (MERRA) must contend with an inhomogeneous observing network. Here the effects of the two most obvious observing system epoch changes, the Advanced Microwave Sounding Unit-A (AMSU-A) series in late 1998 and, to a lesser extent, the earlier advent of the Special Sensor Microwave Imager (SSM/I) in late 1987 are examined. These sensor changes affect model moisture and enthalpy increments and thus water and energy fluxes, since the latter result from model physics processes that respond sensitively to state variable forcing. Inclusion of the analysis increments in the MERRA dataset is a unique feature among reanalyses that facilitates understanding the relationships between analysis forcing and flux response.
In stepwise fashion in time, the vertically integrated global-mean moisture increments change sign from drying to moistening and heating increments drop nearly 15 W m−2 over the 30 plus years of the assimilated products. Regression of flux quantities on an El Niño–Southern Oscillation sea surface temperature (SST) index analysis reveals that this mode of climate variability dominates interannual signals and its leading expression is minimally affected by satellite observing system changes. Conversely, precipitation patterns and other fluxes influenced by SST changes associated with Pacific decadal variability (PDV) are significantly distorted. Observing system changes also induce a nonstationary component to the annual cycle signals.
Principal component regression is found useful for identifying artifacts produced by changes of satellite sensors and defining appropriate adjustments. After the adjustments are applied, the spurious flux trend components are greatly diminished. Time series of the adjusted precipitation and the Global Precipitation Climatology Project (GPCP) data compare favorably on a global basis. The adjustments also provide a much better depiction of precipitation spatial trends associated with PDV-like forcing. The utility as well as associated drawbacks of this statistical adjustment and the prospects for future improvements of the methodology are discussed.
Steady progress in data assimilation efforts has produced a legacy of important contributions regarding atmospheric dynamics (Dole 1989; Simmonds and Keay 2000; Thompson and Wallace 1998, 2000; Thompson et al. 2000; Hoskins and Hodges 2005), tropical water and energy fluxes, and climate variability on interannual time scales (Trenberth et al. 2000, 2001, 2008; Bove et al. 1998; Wang 2002). Yet, as a tool for capturing decadal-scale variability and detection of trends, the reliability of reanalysis products has yet to be realized (Bengtsson et al. 2004; Sterl 2004). In part, this is because water and energy (and other) fluxes produced by the assimilation procedure are subject to uncertainties in the physical parameterizations. Perhaps more importantly, because our satellite observing system has been characterized by an evolving capacity to correct model state variables of heat, moisture, and momentum on which the physical parameterizations critically depend, the fluxes determined from reanalyses contain variability that arises as a result of these discreet, aperiodic observing system improvements. In this paper we will examine how these observing system changes have affected the depiction of climate variability within the 30-plus year period covered by the Modern Era Retrospective-Analysis for Research and Applications (MERRA).
Bias correction methods developed since initial reanalysis efforts have significantly improved the consistency of more recent products. Dee and Todling (2000) used rawinsonde humidity data to correct slowly evolving background guess systematic error in the forecast model component of the Goddard Earth Observing System (GEOS) data assimilation system. Andersson et al. (2005) were able to diagnose a posteriori the incompatibility of European Centre for Medium-Range Weather Forecasts (ECMWF) model background dry humidities in cloud-free areas with satellite moisture retrievals that ultimately resulted in increased rainout of assimilated water vapor and unrealistically high precipitation in the 40-yr ECMWF Re-Analysis (ERA-40). More exacting demands to accommodate the monitoring of decadal climate signals and the growing but ephemeral mix of remotely sensed data have spurred great progress in so-called bias aware methodologies (Dee, 2005). Variationally constrained bias estimation strategies have been developed (Derber and Wu 1998; Dee 2004) that embed the correction within the cost function minimization process (i.e., reconciling observations with the background guess). This allows the detection of time-dependent biases from the entire suite of observing systems (e.g., Dee and Uppala 2009). Like the Japanese 25-year Reanalysis (JRA-25; Onogi et al. 2007) and the ECMWF Interim Reanalysis (EC-Interim; Berrisford et al. 2009), the MERRA employs this bias correction approach.
Despite these advances, additional improvement in bias correction of datasets and models is still required as evidenced by the results of Trenberth et al. (2009) and Bosilovich et al. (2011) who highlight uncertainties in water and energy fluxes among the most recent reanalysis efforts. For example, net downward energy flux at the surface over oceans is approximately 14 W m−2 for MERRA (but −18 W m−2 for JRA) over the period 2000–04 when modern remote sensing measurements should have the greatest impact. Global precipitation trends in MERRA exceed 0.5 mm day−1 or nearly 20% of the climatological global amount over this period, similar to the older ERA-40 reanalysis. In contrast, the EC-Interim reanalysis exhibits a decrease in rainfall of about 0.3 mm day−1 since 1990. Clearly, there are remaining issues, not only with model physics inadequacies but also their interactions with constraining datasets whose bias properties are varied and time dependent in nature.
Bosilovich et al. (2011) have presented a comparative summary of MERRA fluxes with those of other recent reanalyses and also highlighted the relationship between the MERRA trends and the initial availability of Advanced Microwave Sounding Unit-A (AMSU-A) data. The goal of this paper is to provide a more detailed characterization of these time-dependent biases and explore the efficacy of a statistical correction strategy to reduce these artificial trends in water and energy fluxes. After reviewing the data used in the study (section 2), we begin by analyzing the temporal variability of the water vapor and heat budget increments and relate these forcings to the stepwise evolution in passive microwave sensing associated with the AMSU-A and Special Sensor Microwave Imager (SSM/I) (section 3). Section 4 provides some diagnostics on what aspects of climate variability within the 30-yr period survive or are distorted by these artifacts. We then use principal component analysis to characterize space–time variability of the increment terms. In section 5 we show how applying a statistical regression methodology to identify and eliminate trend artifacts associated with satellite sensor changes provides an effective correction. An assessment of the adjusted fluxes is presented in section 6 along with concluding remarks.
2. MERRA data and diagnostics
All data used in this study are monthly-mean assimilated products taken from the MERRA collection archived at the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC) (http://disc.sci.gsfc.nasa.gov/). These include monthly mean vertical integrals of heat and water budget components (tavgM_2d_int_Nx and tavgM_2d_int_Nx) and individual radiative flux components (tavgM_2d_rad_Nx) for the period 1979–2009. These quantities at their native resolution (0.50° latitude × 0.67° longitude) were regridded via box averaging to a 2.5° latitude × 2.5° longitude grid to facilitate computing covariance statistics. More detailed documentation on the quantities in these datasets can be found in the MERRA File Specification Document (Suarez et al. 2010). For a more detailed discussion of data assimilated into MERRA, see Rienecker et al. (2011) and the Global Modeling and Assimilation Office Web site (http://gmao.gsfc.nasa.gov/research/).
From these data we can reconstruct the MERRA budget for enthalpy, H ≡ cpTυ, written as
where cp is the heat capacity of dry air at constant pressure and Tυ is the virtual temperature. Overbars denote a vertical integration over the mass of the atmosphere and the subscript H denotes reference to the virtual enthalpy budget. The first two terms on the right-hand side of (1) represent the convergence of enthalpy and the release of total potential energy, the global integral of the latter being equivalent to the generation of kinetic energy. The other tendencies involve moist processes (MSTH, which is essentially precipitation), radiative heating (RAD), and turbulent diffusion, that is, sensible heat flux (SENS). The analysis increment (ANAH) is the tendency that is added to the prognostic budgets due to the observational analysis: resH is a small value that includes the gravity wave drag and a residual that results from maintaining energy balance in the presence of numerical dissipation, each of which is also included in the output diagnostics (Suarez et al. 2010).
The vertically integrated radiation heating term, RAD, can be expanded further into the difference of the top-of-atmosphere (TOA) and surface shortwave fluxes (SWT and SWS, respectively) minus the sum of outgoing longwave radiation (OLR) at the TOA and the net longwave flux to the surface (LWS):
The MERRA vertically integrated total water vapor budget can be written as
where the change of water vapor is related to the dynamical convergence of water vapor; the net moist physical processes (MSTq), essentially the negative of the sum of convective, large-scale and frozen precipitation, since precipitation is a sink of water vapor; and evaporation (EVAP). Here the subscript q denotes reference to the water vapor budget. Note that −(L/cp)MSTq = MSTH. Two nonphysical terms affect the moisture budget. ANAq represents the analysis increment of water vapor, and resq represents a very small amount of negative filling, ensuring positive water vapor content.
3. MERRA T, q increments and their evolution
In this section, we examine the statistical properties of the MERRA vertically integrated moisture and temperature increments, focusing in particular on signatures of satellite impacts. Climatological mean maps of the moisture and heating increments are presented in Fig. 1. The degree of systematic spatial structure found in these fields is consistent with our understanding that model physics deficiencies can vary as a function of climate regime (e.g., regions of predominantly upward or downward motion). In general, positive moisture increment values maximize over the tropics, most notably in the Northern Hemisphere summer warm pool areas (western Pacific, northern Indian Ocean, and inter-American seas). These are preferred areas of deep moisture and frequent convection (Adler et al. 2003), and the added moisture from the increments acts to systematically increase precipitation (not shown). There are also positive values in oceanic subtropical ridge locations in eastern ocean basins where precipitation is quite small but the transition from stratocumulus to trade wind cumulus regimes are found. In the eastern Pacific a region of negative increments (drying) just south of the equator forms a pronounced dipole structure with the positive values lying immediately to the south. In the region of negative moisture increments, boundary layer air over the cold upwelling water off the coast of South America is systematically too wet in the background field. Presumably these features indicate that the delicate interplay between subsidence drying and shallow convective moistening of the equatorward moving air masses is not handled as well in the model as it should be.
Heating increments are predominantly positive over the continents and adjacent waters with an equatorially symmetric cooling pattern over the central subtropical Pacific. That the moisture and heating increments have such different spatial patterns indicates that dynamical transports of heat and moisture significantly affect the thermodynamic adjustments. The off-equatorial cooling pattern is reminiscent of Rossby modes emanating from forced heating over the western Pacific (Gill 1980; Hoskins and Karoly 1981) where ANAq has a relative maximum. One possible explanation for this structure is that the precipitation and associated Gill-like dynamical heating response, poleward and downstream of the heat source, is somewhat excessive compared to reality. This is consistent with the comparisons to Global Precipitation Climatology Project (GPCP) precipitation shown by Bosilovich et al. (2011). The extra heating would then require a systematic cooling increment for the analysis to match satellite radiances.
Given a uniform observing system, ANAH and ANAq would be essentially stationary in a statistical sense with some constant bias owing to model physics deficiencies.1 Unfortunately, this is not the case in reality due to the significant stepwise evolution of the satellite observing system. Figure 2 presents Hovmoeller plots of zonally averaged ANAq and ANAH anomalies for the 1979–2009 period (henceforth, all anomaly quantities are understood to be departures from their monthly varying seasonal climatology defined by the period 1979–2009). From this perspective discontinuities in the increments resulting from the discreet changes in the observing system are now evident, particularly for ANAq. The most obvious of these is the abrupt change to positive moisture forcing in November 1998, when the assimilation of NOAA-15 AMSU-A radiances begins. This moistening is prominent at tropical latitudes and also over the southern oceans. Also apparent is the effect of the SSM/I sensors, whose data begin in late 1987 and tend to decrease the drying from other observations in the subtropics but dry the high latitudes, particularly in the SH. The increasing amplitude of the apparent SSM/I-induced changes over the Southern Ocean likely results from the increasing number of those sensors deployed with time during the 1990s. Presumably, the SSM/I and AMSU-A effects in the NH are muted because of the far denser conventional observing system there compared to the SH and tropics.
The ANAH change associated with the onset of AMSU-A data ingest is prominent, particularly in the equatorial region and in the Southern Hemisphere. It is evident that, as the moisture forcing increases with time in these regions, heating increments generally trend in the opposite direction; less heating is needed to keep the analyzed temperature consistent with observations. Beyond this first-order change in ANAH with time series, other features are also noted such as the presence of “blockiness” or discontinuities of multiyear duration in the Northern Hemisphere, particularly poleward of 30°N. It is not clear why these shorter duration features should occur in a relatively data rich portion of the globe. It is likely that significant surface emissivity uncertainties affecting the assimilation of window channel data, especially the Microwave Sounding Unit (MSU) channel 1 and AMSU-A channels 1–3 and 15, exist, varying from sensor to sensor in such a way as to affect the inference of atmospheric temperature, but the explanation of these patterns will require more in-depth analysis.
The presence of an annual cycle power in both increments despite the removal of the respective mean values over the reanalysis is also striking and indicates variations both in amplitude and in phase. After the build up of SSM/I sensors beginning in 1987, the increments near 50°S have a maximum drying in SH winter [June–August (JJA)]; yet after the AMSU-A ingest begins, a dramatic SH winter moistening is evident. A summer heating–winter cooling cycle in ANAH centered south of 45°S also reverses phase after the start of AMSU-A data. Another shift in the annual variation of ANAH begins in the Northern Hemisphere in 1998 and strengthens substantially in 2002. Clearly, the changes in sensors are associated with changes in the required annual cycle of forcing. For purposes of climate analysis this is a complicating factor.
All of these structures noted here emphasize the important interactions of a changing observing system interacting with a model having its own specific hydroclimatic biases. To explore the origin of these discontinuities further we turn first to monthly mean statistics of observation-minus-forecast (OMF) and observation-minus-analysis (OMA) residuals for the various AMSU-A channels of the NOAA polar orbiting suite. Simulated radiances based on the model 6-h forecast and analyzed fields are generated using the forward observation operator of the Community Radiative Transfer Model (CRTM) (Han et al. 2006). The difference between the monthly mean values of OMF and OMA residuals, in turn, provides an estimate of the monthly mean contribution of the radiance increments [analysis minus forecast (AMF)] in observation space. Figure 3 shows time-mean maps of AMF from 1999 to 2002 for NOAA-15 AMSU-A channel 15, one of the AMSU-A window channels (1, 2, 3, 15) sensitive to water vapor and liquid water emission, and for channel 5, the oxygen absorption channel sensitive to mid- and lower-tropospheric temperature and surface temperature and emissivity. It is primarily through the window channels that the AMSU-A instrument can affect the water vapor increments. Over the ocean the channel 15 AMF pattern shows a strong resemblance to various features of the climatological moisture increment pattern in Fig. 1 with significant amplitude over the northern Indian Ocean/warm pool region, along the South Pacific convergence zone, and around the Gulf of Mexico and surrounding regions. The AMF values in channel 5 have a pattern of cooling over the NH summer warm pool position—consistent with the heating increment Hovmoeller plot in Fig. 2.
Over land, where the window channels are sensitive to surface emissivity variations, the signals of both channels are quite large. In MERRA, the window channels are assimilated over land also, but only for footprints that can be classified as having a single surface type (i.e., mixed land/water or ice covered/ice free points are eliminated). The challenges of forward modeling for window channels include the complexities of vegetation and soil attributes and their interplay. A related issue is that model rainfall driving of soil moisture likely departs sufficiently from reality to further degrade first guess emissivity and complicate the subsequent retrieval of atmospheric behavior from these channels. It should also be noted that biases in the model diurnal cycle over land may generate significant differences between the first guess and observation. To the extent that AMSU-A multiple channel data differ from the single 50.3-GHz channel on the MSU in their ability to correct the first guess emissivity, systematic changes may propagate to the atmospheric temperature. Conversely, we suspect that because of the far denser conventional observation network over land, the increment contributions from these channels end up having much less impact on the increments ANAq and ANAH. However, as noted in the discussion of Fig. 2, there is some evidence to suggest that their effects on temporal variations over land are not negligible.
Confirmation of the central role that AMSU-A window channels play in producing a change in the moisture increments is demonstrated in an assimilation experiment in which the window channel data are withheld from the assimilation (Fig. 4). Within two weeks of the onset of AMSU-A data assimilation (2 November 1998), the moisture increments in MERRA have led to systematically larger precipitation (~0.1 mm day−1) compared to the experiment in which these moisture-sensitive channels are withheld. This result is comparable to that from an experiment in which all AMSU-A channels are withheld (not shown). Accounting for the radiometric warming of the window channels by liquid water (after removal of likely precipitating footprints) is necessary to correctly interpret the observed radiances in terms of water vapor present. A second experiment in which the window channels are retained but the standard bias correction for liquid water is eliminated results in an even larger forcing of precipitation than seen in MERRA. This illustrates the importance of the liquid water correction and raises the question of whether the present correction in MERRA is sufficiently accurate. Other possible causes of these artifacts might include differences in the way the cross-track scanning AMSU-A and the conical scanning SSM/I instruments detect vertical moisture structure, or perhaps its ensuing bias correction. In particular, the radiometric effects of footprint filling by liquid are nonlinear and may differ significantly between SSM/I and AMSU-A scenes. Differences in bias correction efficacy between the sensor systems may result. A full analysis of these issues is beyond the scope of the present paper but is the subject of ongoing investigation. Regardless of the ultimate cause, there is strong evidence here that the problem begins with nonphysical moisture forcing changes leading to condensation and heating increment changes.
4. Modes of variability in MERRA
a. Globally averaged time series
The presence of these discontinuities in the ANAq and ANAH forcing, and their demonstrated relationship to the microwave components of the satellite observing system, naturally raises concerns regarding the ability to detect real variability in the climate system from MERRA. In this section, we present an analysis of the major modes of variability in the vertically integrated MERRA water and energy budget variations in terms of global means over separate ocean and land areas (Figs. 5 and 6). The water vapor increment and precipitation2 in Fig. 5 are strongly correlated over ocean and are the largest terms in the moisture balance equation. Note that the time-mean moisture increment is small over ocean but the increment changes sign dramatically beginning in November 1998 when NOAA-15 AMSU-A assimilation begins. Moisture convergence and evaporation also respond in time with the increment and decrease over the oceans after AMSU-A enters the datastream in November 1998. Moisture transport from ocean to land thus increases at this time. Over land the predominant balance is between increasing moisture convergence and precipitation anomalies with the increment playing a much smaller role. The moisture increment correlates negatively with moisture convergence, suggesting that over land the influence of the moisture transport from oceanic regions runs counter to the conventional observations over land.
Especially prominent in the over-ocean budget is an annual cycle component. Despite the removal of a mean annual cycle over the 30-yr record, Fig. 5 provides evidence of nonstationarity at this frequency. These signals change phase over the 30-yr period.
The most prominent aspect of the enthalpy budget in Fig. 6 is the negative correlation between the heating increment and precipitation over ocean. This is consistent with the interpretation that as the moisture increment forcing increases, it triggers more precipitation, so less net atmospheric heating is needed to force the model temperature to agree with the observations. From this global perspective the stepwise increases in moisture increments appear to be the driving mechanism. However, the enthalpy increment has a more complicated steplike structure than does the moisture increment, and the upward precipitation steps in late 1987 and 1998 are less obvious in the heating increment. It is also possible that other factors, such as intercalibration uncertainties of the MSU and AMSU temperature sounding channels, may be inducing additional increment changes that complicate the signal.
Like the precipitation heating, dynamical heat transport has strong variations with power at the 1-yr time scale that results from heating increments that have changing phase and amplitude. These are particularly dominant over land. Note that for both moisture and heat the increments over land are strongly anticorrelated with the dynamical transports, indicating that increment-induced precipitation heating changes over the ocean lead to circulations and transports that must be compensated over land to keep Tυ and qυ near observed values there. Sensible heat flux variations are small and unremarkable when averaged over global domains, but in the following section we will note significant regional variability.
For averages taken over ocean and land areas at the global scale the trend and nonstationary annual cycle signals appear to completely dominate any physical variations that might be present on interannual and longer time scales.
b. EOF analysis
For convenience as a diagnostic, we have also used principal component analysis (PCA) to extract major patterns of coherent variability in exploring the analysis increments, ANAq and ANAH. The leading five EOF patterns and PCs for ANAq and ANAH are presented in Figs. 7 and 8 and explain roughly 30% and 40%, respectively, of the monthly variance after the climatological mean plus annual cycle are removed3 (note that the product of the EOF and principal component time series of a given mode produces the contribution of that mode to the total anomaly signal). The leading EOF of the moisture increment has as its primary signal the increase in water vapor forcing over the northern Indian Ocean, tropical Pacific, and inter-American seas. Variations in the phase and amplitude of the annual cycle dominate the second mode and relate to the pre-SSM/I and AMSU-A epochs. The third mode appears to capture the tropical moistening and Southern Ocean drying as the number of SSM/I sensors grows after August 1987, and then the near reversal of this forcing with the beginning of AMSU data. For the enthalpy budget, the leading trend mode carries the signal of decreased (increased) heating over the warm pool and western Pacific (eastern Pacific) with time. Other modes modify the linearity of this trend.
In both ANAq and ANAH, each of the modes has superimposed on it signals of evolving departures from a mean annual cycle. These features change both phase and amplitude between the pre-SSM/I and the Advanced Television and Infrared Observation Satellite Operational Vertical Sounder (ATOVS) epochs, as discussed in connection with Fig. 2. These signals of nonstationarity are also present, to a smaller extent, in the remaining modes, which we have not shown. Although most of the dominant signals in both increments are found over the global oceans, there are also significant centers of action over land. The moisture increments over land are largely confined to the Southern Hemisphere, while heating increments are present over Northern Hemisphere land areas as well. Central Africa is an area of particular interest where Bosilovich et al. (2011) have shown that strong decreases here in MERRA precipitation, evaporation, and moisture convergence occur after 1995 in conjunction with a change in equipment at radiosonde station Bangui (station ID 64650). They also show that GPCP precipitation gives no indication of the drop seen in MERRA precipitation. The leading mode of ANAq (Fig. 7) indicates a change in sign from negative to positive in moisture forcing in central Africa as the AMSU sensors become available after 1998. We suspect that the onset of more moisture data with the window channels may enable the analysis to partially combat the erroneous drying produced by the radiosonde instrumentation change.
While PCA is able to isolate a large fraction of the variability associated with sensor changes, one cannot say that a given mode reflects a specific sensor change. It is important to remember that, by construction, the orthogonality of the successive EOFs is purely a mathematical property. Nevertheless, the PCA diagnostics presented here highlight the discrete temporal discontinuities in the MERRA budget terms and, by virtue of their coincidence with sensor changes, illustrate the dominance of the sensor-related artifacts in explaining the total variance of the increments.
c. Interannual and near-decadal signals
Given that a large fraction of the interannual-to-decadal variability in monthly mean flux anomalies can be explained by the effects of the AMSU-A and SSM/I observations entering the datastream, what are the implications for the utility of MERRA in terms of climate variability within the period analyzed? The El Niño–Southern Oscillation phenomenon is the largest globally coherent variability signal in the climate system. To determine whether this signal can be extracted despite the trends in the data, we have regressed the monthly flux anomaly fields against an ENSO SST index. Our base SST time series is constructed using the extreme least lag methodology (Chen et al. 2008a), which varies similarly in time with Niño-3.4 SST but also accommodates both positively and negatively lagged SST anomalies over the global oceans. The results of regressing the moisture and heat flux anomalies against the time series of this globally averaged index, normalized by its standard deviation of 0.061 K, are shown in Fig. 9. Clearly, the precipitation patterns (−MSTq or MSTH) agree closely with well-known indices of ENSO-induced rainfall anomalies (Ropelewski and Halpert 1987; Curtis and Adler 2000). Increased precipitation over the equatorial central and eastern Pacific with suppressed precipitation over the Maritime Continent and South Pacific convergence zone is well captured in the reanalysis. Teleconnected precipitation anomalies over the Amazon and Indian Ocean regions are also quite realistic. Radiative heating anomalies associated with the precipitation patterns capture the anomalous heating in the precipitation regions and the increased energy loss via OLR in the surrounding subtropics and tropical region west of the date line. The former is due principally to increased greenhouse effects of water vapor and clouds, while the latter is due to a warmer radiating atmosphere with somewhat enhanced aridity combining to enhance OLR in anomalous subsiding regions. Both moisture and heating dynamics terms indicate the strong role in bringing moisture to regions of anomalously large precipitation and supporting the release of potential energy through anomalous vertical overturning. The latent heat flux exhibits a complex pattern with double maxima in the eastern tropical Pacific and broadly negative values west of the date line and over the subtropical North Atlantic Ocean. Both changes in wind speed and near-surface moisture anomalies act to produce this complex pattern. Sensible heat flux anomalies are largest over land and show elevated values in the Australasia region and Amazon basin; these are areas well known to have excessive aridity during warm ENSO events. Reduced sensible heat fluxes are found in eastern equatorial Africa, associated with elevated cloudiness, precipitation, and reduced solar forcing during warm ENSO events. Additional future studies can examine these gross patterns in more detail—for example, aspects of the lag relationships between SSTs and fluxes, the asymmetry of warm and cold events, and details of energy transformations and transports. Our point here is just to show that, despite the pronounced trend in many of the fluxes, coherent relationships between SST variations and the regional flux patterns on the interannual scale appear to remain intact.
The relatively high frequency of ENSO (2–7 yr) allows effective separation of the SSTs and flux signals from the trend signal; furthermore, in the case of ENSO the signal is well known, a priori. For lower amplitude, lower frequency variations in fluxes, the problem becomes more difficult. For example, Pacific decadal variability (PDV) (Zhang et al. 1997; Deser et al. 2004; Chen et al. 2008b; Burgman et al. 2008) is a well-known low-frequency mode of behavior that produces systematic variations in precipitation, radiation, and circulation in conjunction with broad changes in SST structure over the Pacific basin. Figure 10 shows the leading SST EOF and PC after the ENSO signals have been removed from the data during the 1979–2009 period via the method of Chen et al. (2008a). After performing an EOF analysis on this filtered SST data, the first mode contains the bulk of the PDV signal. There is the characteristic trend of cooler water over the eastern Pacific and the horseshoelike increase in SST in the western subtropics of both the Northern and Southern Hemispheres. Note that some signature of the Atlantic meridional overturning (Schlesinger and Ramankutty 1994; Enfield et al. 2001) and the cooling induced by the Mt. Pinatubo eruption are also present. Some of the trend in this PC is also related to the century-scale rise in SST, but overall the PDV signal dominates regional trends. Also shown in Fig. 11 are GPCP and MERRA precipitation anomalies regressed onto this PC. The GPCP precipitation pattern is very similar to that reported by Burgman et al. (2008) in an observational analysis of PDV. MERRA bears some similarity in the large-scale patterns but notable differences, such as those in the central Pacific and Atlantic Oceans, are prominent. Clearly, the nonphysical long-term trend signal in MERRA gets picked up in the regression procedure. While no 30-yr reanalysis is long enough to truly quantify decadal-scale variability, low-frequency variability patterns within the satellite era are distorted in MERRA by the satellite-related artifacts.
5. A statistical adjustment methodology
Is it possible to adjust the fluxes in any defensible way so as to minimize these artifacts? One might argue that simply removing the trend modes by means of the EOF analysis would be useful, but these modes explain only a portion of the detectable sensor change signal—it is not linear. Furthermore, additional analysis of the EOF decomposition of the fluxes (not shown) indicated that for many terms the trends are mixed with real variability (e.g., the ENSO signals). Attempting to remove these modes would damage real signals of variability.
A more viable approach is to use the leading modes of ANAq and ANAH discussed above to remove the artifacts in the fluxes by linear regression, that is, projecting the flux data onto the PCs of the individual increment modes and then removing these components from the fluxes. This technique is commonly referred to as principal component regression (PCR). Conceptually, the vertically integrated equation for moisture or enthalpy can be written as
where X is moisture or enthalpy, Fi is any one of n physical terms (e.g., surface and TOA fluxes and dynamical transports), and Ix is the moisture or enthalpy increment. All quantities are understood to be vertically integrated, mass weighted departures from their respective monthly resolved climatologies.
We can formulate an adjusted physical budget term as
where the Fi,j are the spatial vectors of regression coefficients of the ith budget term on the PCs of the m leading increment modes. The PCs have only temporal variability.
We also formulate a modified increment,
composed of the original increment minus the sum of m modes, each of which is composed of the EOF spatial vectors Ij multiplied by their respective principal components, PCj.
The adjusted budget equation now reads
where the modified increment and the physical terms are now free of the leading modes that relate to the satellite artifacts.
Of course, one potential drawback to this procedure is that we are subjectively declaring which subset of increment modes is affected by artifacts. In general, we cannot expect all of the artifact-related variability to be collected into just a few modes. It is also likely that some physical variability will project to some extent onto the modes that we select, mixing with the nonphysical signal. As noted earlier, removing sensor change effects from the fluxes is complicated by the fact that we also expect the analysis increments to have some sensitivity to climate variability (e.g., ENSO or other events that systematically perturb the distribution and intensity of precipitation and radiation).
The methodology outlined here is closely related to redundancy analysis (RA) (von Storch and Zwiers 1999; Wang and Zwiers 2001; Bakalian et al. 2010), a variant of several techniques widely used to find covariability between two datasets. Like canonical correlation analysis (CCA), RA finds a hierarchy of paired modes, predictors and predictands that maximizes some metric of variability. But, rather than maximizing correlation, RA maximizes the explained predictand variance and associates that with a predictor mode. We performed the analysis with both RA and PCR methods, with virtually identical results.
Based on examination of the increment EOFs and PCs, we selected five PCA modes from each of the increments to define the artifact signals. These were the first four and the seventh moisture increment modes and the first five enthalpy increment modes (i.e., Figs. 7 and 8). Experimenting with variations on this selection showed that the first three modes of each increment are crucial and that using more than 10 modes began to weaken the retained variance noticeably.
In centering the data for the regression analysis, we must define an annual cycle and mean values. We have opted to use the entire record from 1979 onward to define the climatology because choosing any shorter period might not yield adequate sampling of the annual cycle. Obviously this is subjective and one could argue that the most recent period with the most robust measurements from satellite, say 2000 onward, would be a more defensible choice. Given the uncertainties associated with the AMSU-A window channel effects on the assimilation of moisture, we believe this choice will remain ambiguous until specific satellite algorithm and bias correction issues are studied more thoroughly with numerical experimentation.
6. Results and discussion
a. Adjusted time series for heat and water budget terms
Figures 11 and 12 show area-averaged anomaly time series of the corrected budget quantities over separate global ocean and land domains. By construction the large trends in the original moisture and heating increments, as well as precipitation and evaporation, are now absent after applying the PCR procedure. These reductions are substantial over ocean but much smaller over land (cf. to Figs. 5 and 6, noting differences in the scales). Perhaps the most significant result here is that, when averaged to global domains, the remaining increments are of the same size as the leading physical terms. Furthermore, whereas the unadjusted heating and moistening increments were also strongly correlated, they now have a much weaker relationship to each other. These remaining increments should be due largely to model physics inaccuracies interacting with time-varying flows, as well as remaining observational error. The latter may have both random and flow-dependent components. The moisture increments and evaporation tend to be anticorrelated, as do the heating increments and dynamical heat transport. This is not surprising since we have removed the variability in these fluxes that correlates with the part of the increment signal induced by satellite sensor changes. But, it also suggests that, when moisture is added (or subtracted) by the increments to follow the observations, it acts to decrease (or increase) the near-surface humidity deficit and alter the evaporation accordingly. Likewise, if enthalpy increments constraining the analysis to temperature observations are positive, the dynamical response to this heating should be that heat is exported to land areas and potential energy lost to kinetic energy via vertical overturning. The opposite argument holds for negative increments.
Bosilovich et al. (2011) noted significant changes in the increments, transports, and fluxes over tropical continental regions between the decades before and since the availability of AMSU-A data. Because the changes in radiosonde moisture measurements in central Africa also result in a local steplike change in a datastream, they also project onto the low-frequency variance of the leading increment and constitute a trend that gets removed.
Ideally, the next step would be to “explain,” as much as possible, the remaining increments in terms of statistical corrections to the model physics terms. This is the objective that Schubert and Chang (1996) sought in their analysis of a much more limited data sample. We leave this type of adjustment for future work but point out here the likely direction of such adjustments. Note that the corrected moisture increment has some anticorrelation with ENSO signals—that is, it is subtracting (adding) moisture when SST is warm (cold), as does precipitation. If the moisture increment is used entirely to adjust the precipitation rate, that is, , the resulting global mean value agrees better with GPCP precipitation in terms of interannual variability (Fig. 13). There are some remaining departures in agreement between and GPCP, such as in the mid-1990s and post-2005, but here the agreement with the Hilburn and Wentz (2008) SSM/I precipitation over ocean (not shown) is better, underscoring the fact that significant uncertainties remain with precipitation retrievals.
b. Component radiative fluxes
Another interesting aspect of the PCR adjustment is that, although the radiative heating term in the atmospheric heat budget, RAD, has no large global trend in the raw MERRA data, this is decidedly not so for the individual SW and LW fluxes at the TOA and the surface. Significant corrections through the PCR procedure are made to these component fluxes. Figure 14 shows time series of these quantities before and after the correction. Except for OLR, significant stepwise trends are present in each flux component. The strong correlation between TOA and surface SW variations is present, as expected (Ramanathan 1986; Cess et al. 1991).
Selected non-ENSO EOFs of the unadjusted MERRA radiative flux components that dominate their trend and annual cycle behavior are shown in Fig. 15. The leading SW mode corresponds to broadly increased precipitation-driven cloudiness over the warm pool region, the eastern Pacific, and the Southern Ocean. Mode 3 of the TOA net SW captures a change in annual cycle phase that maximizes at higher latitudes. The second OLR mode shows a decrease with time, consistent with increasing cloudiness over the Indo-Pacific warm pool. Surface LW increases are ubiquitous over the global oceans, and the PC of the leading mode correlates well with column water vapor increases (not shown), suggesting that the onset of SSM/I and eventually NOAA-15 and NOAA-16 AMSU-A moisture data have a dominant impact on surface downward LW. It is notable that the largest surface LW signal is present over central Africa where the decreasing moisture convergence associated with the evolving rawinsonde signal in the unmodified MERRA data projects strongly and indicates the effects of lower-tropospheric drying.
The opposition of SW and LW fluxes and the significant correlation between SW at the surface and the TOA in Figs. 14 and 15 have consequences for net atmospheric radiative cooling and warming. The increase in water vapor increments with time sustains higher water vapor in the lower troposphere which, in addition to enhancing downward LW at the surface, supports more condensation and cloudiness. The subsequent SW reduction at the surface from the increasing cloud fraction more than offsets the increase in LW there and, in fact, results in the decreasing net radiation to the surface of about 5 W m−2. High cold cloudiness associated with increased tropical precipitation induced by increasing water vapor increments slightly reduces OLR and offsets the increasing longwave loss to the surface. So, the net adjustments in SW and LW brought about by the changing increments of water vapor (and temperature) are near canceling in terms of atmospheric radiative loss. But, as in the case of surface changes, the increased TOA SW cloud reflection overcomes the OLR signal to produce a net decrease of about 5 W m−2 in net radiation absorbed by the earth–atmosphere system, which would be felt by the ocean in a coupled system.
Adjustments derived from the PCR procedure are effective in identifying the step functions in the fluxes related to the SSM/I and AMSU-A datastream changes (Figs. 11 and 12). Net surface and TOA fluxes now have essentially no trend. Although the net effects of SW on the atmosphere are small even before the correction, the spurious decrease in net SW at the surface of ~8 W m−2 is removed by the PCR procedure. Of course the opposing net LW trend of about half that magnitude is also removed. In light of the mean 13.8 W m−2 net energy flux bias to the ocean in the 2000–04 period noted by Bosilovich et al. (2011), the O(5 W m−2) net ocean surface radiative flux decrease over the 30-yr MERRA period could mean significant biased forcing for ocean and land models forced offline using these unadjusted fluxes. The signals associated with ENSO events are now readily identifiable in the individual LW components with elevated (suppressed) atmospheric LW cooling to space and to the surface during warm (cold) SST events.
c. Regional decadal signals
With the removal of global mean trends from the fluxes, there is concern that any long-term climate signals may also have been removed. The statistical nature of the PCR procedure is incapable of distinguishing a physical trend in the fluxes from that induced by one or more steplike artifacts in the increments. However, this does not necessarily mean that regional trends or variability that do not project on a near-global domain are removed. Here, we reexamine the PDV signal in MERRA to see how the PCR adjustment has affected the realism of this feature. Figure 16 shows the regression of original and corrected precipitation, evaporation, and TOA net SW flux onto the PC. There are striking differences between the adjusted and original quantities. In the adjusted data there is a distinct decrease in precipitation over the equatorial western Pacific and increases in the eastern equatorial Pacific that agree well with the result using GPCP precipitation (Fig. 10) and, thus, with the results of Burgman et al. (2008). These features are distorted in the unadjusted MERRA data. The net radiative heating, RAD, responds to cloudiness changes and mirrors the effects of the PCR adjustments over the western Pacific, but it also indicates the substantially different response in the cloudiness in the eastern Pacific. The unadjusted evaporation pattern is essentially the same as that for its leading EOF (not shown), whereas after the correction much of the same structure remains but with smaller amplitude.
Because the PDV SST signal used in the regression is to first order a trend, any other signal having a trend will be picked up to some extent, whether or not it is physically related to the PDV mechanism. We see that in central Africa the downward trend in precipitation and evaporation, which is likely erroneous (Bosilovich et al. 2011), has also been greatly reduced. This has occurred because of the temporal coherence between the analysis increment, ANAq, and the physical budget terms that oppose the increment, as noted in section 4b. It is perhaps fortuitous that both the AMSU-A timing of data availability and the change in rawinsonde instrumentation induce low frequency trend components of opposite sign.
Clearly the PCR procedure has resulted in significant differences. Based on the better precipitation comparisons, it appears that the PCR adjustments enable a much clearer picture of the low frequency behavior of fluxes during the MERRA period. More analysis is needed to confirm the degree to which this is so.
Like all reanalysis efforts, MERRA has to contend with an inhomogeneous observing network. From a climate perspective, the satellite era (~1979 onward) arguably presents the largest challenges to assimilation efforts because of the significant, but discreet, advances in remote sensing of water vapor and clouds. Here we have addressed the two most obvious observing systems epoch changes: the AMSU-A series in late 1998 and, to a lesser extent, the earlier advent of SSM/I in late 1987. More specifically, the window channels of AMSU-A have been found to produce a strong influence on MERRA water vapor increments that changes their global mean sign from negative to positive. By virtue of the moisture increment influence on net condensation (precipitation), the heating increment anomalies turn negative. The total increment of heating (climatological mean plus anomaly) thus falls to its smallest positive magnitude during the post-1979 period. The dominant pattern of variability for the moisture increment is the stepwise increase in moisture over the tropical oceans, particularly the warm pool and inter-American seas during NH summer. Drying of the Southern Ocean region during the tenure of SSM/I is reversed with the onset of AMSU-A data. As tropical ocean precipitation increases in response to this additional moisture forcing, the heating increment anomalies change from warming to cooling over the western Pacific and much of the Indian Ocean while remaining positive over land areas.
We also find that aspects of the annual cycle are affected by the sensor changes. Removing a mean annual cycle defined over the entire period leaves signals in the fluxes that still have power at the annual time scale, indicating that a nonstationary annual cycle has been induced by the satellite changes. This leads to the question of which period should be regarded as the one over which the climatological mean is defined. One would logically expect the most recent era from 2000 onward as having the most accurate data with which to constrain the GEOS-5 model. But this choice carries the implicit assumption that the degrading effects of possible biases in the algorithm, instrument, or forward radiative transfer modeling are outweighed by the increased ability of the augmented observations to correct model physics-induced biases. We do not believe this issue is settled and are currently analyzing diagnostics of the bias corrections to understand more fully the origin of the increased moisture increments from SSM/I to AMSU-A.
EOF analysis reveals that the increment artifacts are largely captured by approximately five modes (the first four plus the seventh in the case of the moisture budget, and the leading five for the heat budget). Principal component regression was found useful in isolating and greatly reducing artifacts produced by changes of satellite sensors. The method uses the patterns of variability contained within a subset of increment modes that capture the discontinuities in increment forcing to predict a corresponding sequence of modes in each of the flux terms. We have shown that, after the adjustments are applied, the spurious trend components related to the assimilation of the SSM/I and AMSU data are largely eliminated. Although the net radiative heating term in the enthalpy budget does not have nearly the trend amplitude compared to the precipitation term (or the forcing ANAH), the TOA and surface net SW and the surface net LW components do. The adjustment process is effective in removing these spurious signals. Perhaps the most encouraging evidence for the utility of this approach is present in Fig. 13, where the corrected global mean precipitation anomalies were found to relate well to the variability in GPCP precipitation. This also means that, to the extent of our confidence in global mean GPCP values, the global area-averaged time series of evaporation estimates recovered from the corrected MERRA are also more accurate. It was also shown that simply using the raw assimilated fields from MERRA to diagnose pacific decadal variability leads to distorted signals. The patterns of precipitation covariance with the PDV SST mode obtained from the adjusted data were in much better agreement with recent diagnostics (Burgman et al. 2008).
The PCR technique (and ultimately any linear regression method) is not without problems, as it is not possible to distinguish between trends associated with physical processes and those arising as a consequence of the step functions. The selection of those increment modes deemed to be associated with the satellite-induced artifacts is ultimately a subjective process. As is frequently the case, individual EOFs and their PCs cannot usually be equated purely with physical modes. Here, too, we cannot expect individual events where sensors change to be captured by a single mode of the increments. Instead, it takes several modes to effectively describe the collective effects of the sensor changes. One might consider alternative strategies such as a regression of the fluxes on SST to recover the physical variability associated with climate variations. In the case of ENSO, this is possible since the leading mode has very little trend. In contrast, signals such as that of the PDV, which dominates the global SST trend, would end up retaining also the nonphysical trends induced by stepwise increment changes. While we have not attempted to rigorously quantify the error associated with the PCR procedure, a posteriori comparisons to independent observational data may provide the best assessment of how sensitive the results are to the selection of predictor modes.
It is anticipated that the analysis of the MERRA fluxes and the statistical adjustments presented here will add to the utility of the MERRA dataset in climate diagnostics studies. Many types of studies using MERRA data should be largely immune to the problems of time-dependent bias. These would include, for example, compositing studies, various diagnostics on processes at the synoptic-to-intraseasonal scale, or any study in which the characteristic time scale is short compared to time intervals between significant sensor changes. For studies of variability on interannual and longer time scale, the issues presented here should be considered when interpreting the results. The degree to which our statistical adjustments can be taken a step further to partition the remaining systematic increment variability among the physical budget terms is currently under study. The approach of Schubert and Chang (1996) provides one potential methodology; however, whether the assumptions made in that study, which used 6-h data over a limited region for several months, can be extended to global monthly datasets is unknown.
An important complement to the adjustment strategy analysis presented here involves the use of multiple data withholding experiments to determine the contribution of each observing system change to the evolving time series of fluxes. This avenue is actively being pursued. Ultimately the answer to the problem of uncertainties induced via observing system changes lies both in reducing model errors, in more robust bias correction strategies and more mature understanding of sensor retrieval algorithms that will be employed in future reanalysis efforts. These are ongoing challenges being addressed at all institutions involved in data assimilation activities.
This study was supported by the NASA Energy and Water Cycle Study (NEWS) and by the NASA Modeling, Analysis, and Prediction (MAP) Program. The authors thank Dr. Emily Liu, GSFC, for producing the MERRA withholding experiments and for discussions of these results. Also, Drs. Max Suarez and Siegfried Schubert provided many thoughtful comments and discussions concerning the results in this paper. We are especially grateful to Dr. Ron Gelaro, whose comments substantially improved an earlier version of this paper.
This article is included in the Modern Era Retrospective-Analysis for Research and Applications (MERRA) special collection.
MST in the MERRA archived data is the net moisture source due to moist physics, which is essentially precipitation and is negative in the moisture budget and positive in the heat budget. Thus, the budget terms in Fig. 5 are all additive in the sense of contributing to positive qυ and Tυ tendencies.
For the moisture increment, the EOF7 is shown in place of EOF 5 since the latter contain a significant projection of the physical structure associated with ENSO.