F17 vs. F13 SWE regression

M. J. Brodzik

Draft, 2013-10-20

Revision, 2013-10-23 (added Figure 5A)

Final, 2014-01-20 (incorporated CRREL figures, and included algorithm decisions for NISE, daily and historical data for CRREL, and possible extensions of historical SWE)

Introduction

This document contains a short recent history of the passive microwave SWE algorithm at NSIDC, and summarizes the work we have done to derive an extension to this algorithm for the current DMSP-F17 SSMIS instrument. The purpose of the extension is solely to understand F17-derived SWE in the context of the historical record. The reader should note that there are recently developed, more comprehensive, SWE algorithms, for example the one developed by the ESA GlobSnow project, that have current development funding and that Takala et al. (2011) have demonstrated to be better quantitative estimates of SWE.

There are also two efforts to reprocess the SSM/I and SSMIS Level 1b (swath data) record that include sensor intercalibrations. These swath data recalibration efforts are being conducted independently by MEaSUREs projects at Remote Sensing Systems and CSU. As part of yet another MEaSUREs grant that I have been awarded at NSIDC, I will be choosing one of these new data sets as input to a complete reprocessing of gridded passive microwave data. It is possible that the intersensor calibrations from these efforts will render this kind of algorithm regression unnecessary. However, until new new gridded data are available, this regression work is intended to be an extension of the historical NSIDC SWE algorithm.

Background

The basic SSM/I snow water equivalent (SWE) algorithm that we have generally referred to as NSIDC_S1 (Armstrong and Brodzik, 2001) is based on the observation that the temperature gradients differ by frequency as snow volume increases:

SWE(mm) = 3.0 * 1.59 * ( ( TB19H - 6 ) - ( TB37H - 1 ) )                 

for

TB19H = SSM/I descending pass 19 GHz, horizontally polarized brightness temperature (TB) (Kelvins)
TB37H = SSM/I descending pass 37 GHz, horizontally polarized TB (Kelvins)

Combining constants and rearranging terms, this algorithm can be expressed as:

SWE(mm) = A * TB19H + B * TB37H + C

where

A = 4.77
B = -4.77
C = -23.85

This algorithm is based on Chang et al. (1987), which was derived for the similar precursor sensor, SMMR, and adjusted for differences in SMMR and SSM/I frequencies (Armstrong and Brodzik, 2001). At NSIDC, we have used variants of NSIDC_S1 to produce a relatively stable global SWE record derived from SSM/I F08, F11 and F13 (nsidc-0271, Armstrong et al., 2005, updated 2007; and nsidc-0321, Brodzik et al., 2007). The snow-covered area (SCA) from this algorithm is used in our daily Near-real-time SSM/I-SSMIS EASE-Grid Daily Global Ice Concentration and Snow Extent (NISE) product (Nolin et al, 1998, updated daily). We also send daily SWE maps to CRREL that they use in operational reports to field commanders in Central Asia, who are interested in how the current year's snow pack compares to the historical record in order to predict timing and volume of runoff.

However, the F17 SSMIS, launched in late 2006, was a new generation of sensor, and initial processing that simply treated F17 TBs as if they were F13 TBs indicated that this was likely not going to be an acceptable way to extend the time series. (All of the F13 data in this report are derived from RSS V4 input data.)

In the spring of 2009, we were still using F13 data for daily NISE and CRREL SWE maps processing, but we were aware of problems with F13 that indicated its impending demise. At that time, we started comparing F13 data to the F17 data we had received, which was only available for data beginning in March 2008, nearly 1.25 years after the F17 launch date. (The earlier data had not been made publicly available by the Air Force.)

Our investigation concluded that:

  1. F17 and F13 TBs were significantly different for purposes of both sea ice and SWE algorithms, and
  2. There was a significant step change in F17 37H TBs on or around March 26, 2009. The SWE algorithm was producing far too much snow cover to be at all realistic for that spring.
  3. Although the problem we were documenting with the 37H TBs was significant enough to see in time series plots, we did not know whether it was also affecting the other channels, perhaps less noticably.

After trying to report the problem through official channels with no interested response in the problem, my colleague Walt Meier at NSIDC reported the 37H time series issues to a personal contact, Steve Swadley, at NRL Monterey. Steve's colleague, Troy vonRentzell, at NRL in DC, did provide confirmation that the military had identified a problem with both cold and hot calibration limits for both F16 and F17 (personal communication, Swadley to Meier, 14 Aug 2009). VonRentzell confirmed that the F17 calibration values for 37H started to fall outside the defined threshold range values in early April 2009, which in turn affected the distributed 37H TBs. The Air Force eventually expanded the calibration threshold ranges in August of 2009, but had no plans at that time to reprocess earlier affected data (April-August 2009). We requested status on the missing data from launch in late 2006 to March of 2008, but delivering those data was not a high priority for the responsible facility. Swadley did, however, encourage Walt to report to NRL any other discrepancies in the the data that we noticed.

Initial F17 Regression Methods for Sea Ice and SWE

With the cessation of F13 data expected in August 2009, we needed temporary sea ice and SWE algorithms to use for F17 data in the NISE processing stream. Walt used the existing F17 data we considered to be reliable, from March 2008 to March 2009 (prior to the known calibration limits problem), to derive temporary tie points for the sea ice concentration algorithm. The 37H data are not used for the sea ice algorithm, and he did not observe any similar odd behaviors in the channels that are used for the sea ice algorithm.

Walt also derived TB regressions for the 19H and 37H data that NSIDC needed for the SWE algorithm. He used data from 28 Mar 2008 to 26 Mar 2009 (again, prior to the time when the calibration threshold range settings caused problems). Based on some earlier work we had done (Njoku et al., 2006) to identify relatively stable surface calibration locations, we selected 4 surface target locations that we expected to exhibit relatively stable TB time series:

  1. Salonga (Zaire)
  2. Summit (Greenland)
  3. Prairie location (Canada)
  4. DomeC (Antarctica)
Walt used 5x5 pixel kernels at these sites, requiring 13 or more TBs in the respective kernels from each sensor to use the data for a given date in the regression. He calculated a linear fit for daily descending pass 19H and 37H TBs, and averaged the daily slopes and offsets. The resulting values are:
Let:

F17_19H, F17_37H = actual F17 TBs
F13_19H, F13_37H = pseudo-F13 TBs

then 

F13_19H = 1.0077388 * F17_19H - 2.0615463
F13_37H = 1.0046444 * F17_37H - 2.6514931

to yield the SWE algorithm for F17 TBs as:

A = 4.807
B = -4.792
C = -21.036

In the following discussion and figures, I designate this version of the SWE algorithm as F17current. I designate the version that substitutes F17 TBs as if they were F13 TBs as F17unadjusted.

Snow-covered area (SCA) Snow water equivalent (SWE)
Time series of F13 (RSS V4), F17unadjusted and F17current, Mar08-Mar09
Time series of differences F13 (RSS V4), F17unadjusted and F17current, Mar08-Mar09

Figure 1. Time series of SCA and SWE from F13 vs. unadjusted F17 and current F17 regression, using current archive of RSS F13 (V4) and CLASS F17 TBs, Mar 2008 - Mar 2009. Coefficients for F17current were derived using the least-compromised data that we had at the time (summer, 2009).

Beginning in August 2009, we updated NISE processing to use CLASS F17 TBs, with Walt's temporary tie points for the sea ice algorithm, and F17current for the SCA. Daily CRREL SWE maps were also switched to use F17 CLASS input, using F17current. This was our best attempt at continuing the daily products with F17 CLASS data, based on the limited (and potentially low-quality) data we had at that point in time.

From Figure 1, we can say that F17current, the current algorithm we are using for F17 SWE, is a better match to F13 NSIDC_S1 SWE than using unadjusted F17 TBS. F17unadjusted undermeasures for most of the snow season with respect to the F13 algorithm results, and overmeasures near the end (May-Jun) of the season. F17current is a better match to F13 until May-Jun, when the overmeasure is slightly larger than with F17unadjusted.

At the urging of our colleagues at CRREL, the Air Force did eventually release F17 data from launch to March 2008, and reprocessed the data affected by the calibration threshold problem. NSIDC has these data from RSS, designated RSS F17 V7. It is important to understand that these data are not simply the original CLASS TBs with the threshold correction, but they also include the recent work at RSS to intercalibrate both SSMIS and the historical SSM/I record. So we have the record of F17 SSMIS from the beginning of SSMIS measurement period, but not as the NRT CLASS version.

An Improved F17 SWE Regression

My objective with this study was to ascertain whether the complete, corrected record of F17 data could be used to replace the F17current algorithm with an improved algorithm that reduces the current differences between F13 and F17 SWE. This will represent an operational improvement to the NISE data, and will allow CRREL to put the daily SWE maps into a more representative historical context, when comparing near real-time (NRT) daily data to historical SWE statistics derived from the full record of F08/F11/F13.

Ideally, I was hoping to use the same algorithm for F17 data, regardless of source (CLASS or RSS V7). However, after some discussion with Walt, I determined that there are significant differences between the CLASS data stream and the intercalibrated data stream that is RSS V7. I decided to derive two sets of algorithm coefficients, specific to each data stream. Walt took a similar approach with his adjustments to the sea ice algorithm.

Recall that the historical RSS F13 data stream we have archived at NSIDC is derived from RSS V4. Since the RSS F17 V7 data were supposed to include the Air Force corrections to the threshold range issue, I began with F13 vs. RSS F17 V7 data.

The current overlap of F13 and RSS F17 RSS TBs at NSIDC is 2006 day 348 to 2009 day 119. This spans the beginning of F17 to the end of F13. I began with the same one-year comparison time period that Walt had used, 27 Mar 2008 - 26 Mar 2009. I generated 2 sets of potential TB regression coefficients, as follows:

  1. F17stable - I used the same 4 stable target locations, for any dates during the comparison year with gridded data from both F13 and F17. I calculated a linear fit for each of descending 19H and 37H. The advantage to this approach is simplicity. It was similar to Walt's original approach. It is easy and fast to calculate. The disadvantage to this approach is that it is using data from only 4 locations, throwing out a lot of potential information for land surfaces where the SWE algorithm is usually run (always a dubious idea where regressions are concerned). Also, even if the sensors were identical, one of the targets (Canada) will necessarily exhibit differences in F13 and F17 TBs due simply to differences in atmospheric conditions at different orbital overpass times (~6:25am for F13 vs. ~7:25am for F17).
  2. 19H 37H
    Scatter plots, RSS F17 V7 vs. RSS F13 V4

    Figure 2. Scatter plots with simple regression lines, stable target TBs from F13 vs. F17, descending passes, 19H and 37H. Used to derive F17stable TB regressions.

  3. F17snowsurface - I filtered candidate TBs for the regression to only those surfaces where snow is expected to occur. I used a monthly snow climatology to filter each day of gridded TB data to "only those land surfaces where snow may reasonably be expected to occur". I derived a daily linear fit for F17 vs. F13 descending 19H and 37H TBs. The resulting daily slope/intercept values exhibited an annual variation. For each channel, I chose the slope and intercept from the date of the median slope for the year. The advantage to this approach is that it uses much more of the available TB data, on land surfaces where the SWE algorithm will actually be run. The disadvantage is that there is necessarily high autocorrelation in the input data, which violates a basic assumption of linear regression. Furthermore, since there is interannual variability in the slope/intercept functions, choosing only 1 slope/intercept pair to use for data at all times of the year will not necessarily be the correct relationship at other times of the year.
19H 37H
Daily slope/intercept, F17 vs. F13

Figure 3. Time series of daily slope/intercepts for F13 vs. F17, descending passes, daily 19H and 37H TB linear regressions for potentially snow-covered surfaces. Used to derive F17snowsurface TB regressions.

I compared the potential SWE algorithm outputs by plotting time series of SWE for the complete annual season (Figure 4).

Snow-covered area (SCA) Snow water equivalent (SWE)
Time series of F13 with RSS F17unadjusted, F17current, F17stable, F17snowsurface, Mar08-Mar09
Time series of differences, F13 with RSS F17unadjusted, F17current, F17stable, F17snowsurface, Mar08-Mar09

Figure 4. Time series of SCA and SWE from F13 vs. RSS F17unadjusted, F17current, F17stable and F17snowsurface, Mar 2008 - Mar 2009. Note how much worse F17current match is when run on RSS F17 V7 TBs. Both F17stable and F17snowsurface reduce differences but reverse the bias, and are worse (greater absolute difference from F13) than F17current in the late spring.

The two new options, F17stable and F17snowsurface, reduced the average differences compared to F17current, , but were biased in the opposite direction (consistently overmeasuring with respect to F13). In the late springtime months of April, May and June, however, the positive SCA bias of both F17stable and F17snowsurface was larger than the negative SCA bias from F17current.

In the intervening time, Walt derived a newer set of tie points for the sea ice algorithm, using a bracketed iterative approach that he published (Meier et al., 20110. He defined a cost function as a metric for the differences between F13 and F17 sea ice, and assigned weight to both summer and winter sea ice extent and sea ice concentration. He iterated over various algorithm coefficients, minimizing cost, which he defined as a function of the daily differences between the algorithm results from the two sensors. His results achieved better intercalibration than those from manual tie point tuning that followed the methods of Cavalieri et al. 1999. After some discussion with Walt, I decided to adopt the iterative approach for the SWE algorithm to see if I could achieve a better fit than any of my alternatives based solely on TB regressions.

I defined my cost function with equal weight for daily SCA and SWE differences, (but I did not separate winter and summer values as Walt had):

Cost = 0.5 * f(SCA_annual) + 0.5 f(SWE_annual) 

where f(X) is RMSD (root mean squared difference):

f(X) = RMSD = ( (∑(F17Xi - F13Xi)2 ) / n )1/2

for i=1,n days in the comparison year, and F13X and F17X are
either daily SCA or SWE from the designated sensor, for
pixels where neither sensor was missing.

The cost function results for these 4 options are:

  F17unadjusted:     2025046
     F17current:     1506374
      F17stable:      272275
 F17snowsurface:      418154

As expected from Figure 4, for the period Mar 2008 - Mar 2009, the minimal Cost function was derived from F17stable, from among the options F17unadjusted, F17current, F17stable and F17snowurface. I suspected that trying Walt's iterative approach to tune the SWE algorithm could at least remove the measurement bias and possibly decrease the springtime overmeasure.

To apply the iterative approach, I used the three-coefficient version of the SWE algorithm:

SWE(mm) = A * TB19H + B * TB37H + C

adjusting each coefficient in turn, in stepwise fashion, minimizing the cost function. Several attempts yielded slightly different answers, depending on the order of the adjusted coefficients, and the step size. I took this result to indicate there is no unique solution, and that I am simply searching for a solution that is "good enough" for my purposes.

Based on F17stable showing the lowest cost so far, I began with the F17stable coefficients and minimized cost by adjusting A, B, and C by step increments of 0.001. On the first attempt at a new coefficient, if the new cost was larger than the original cost, I reversed the direction of the coefficient delta before proceeding to the next coefficient. On subsequent steps, if the new cost was larger than the last cost, I reverted to the last value fo the coefficient and proceeded to the next coefficient. The best (lowest Cost) iteration reduced the overall positive bias in F17stable to near zero, but never eliminated the positive bias in the springtime.

Given the resistant positive bias in the spring, I also calculated two further variations on Walt's sea ice summer vs. winter cost functions:

Cost = 0.25 * f(SCA_spring) + 0.25 f(SWE_spring) + 0.25 * f(SCA_nonspring) + 0.25 f(SWE_nonspring)

for spring dates in Apr-Jun 
and nonspring dates in Jul-Mar

This variation weighted the differences during April-June more heavily than during the rest of the year. And a second variation on the cost function:

Cost = 0.5 * f(SCA_spring) + 0.5 f(SWE_spring)

This variation weighted only the differences during April-June, ignoring differences during the rest of the year. Results from this variation were indistinguishable from the equal cost variation, and are not included in Figure 5.

Snow-covered area (SCA) Snow water equivalent (SWE)
Time series of F13, RSS F17unadjusted, F17current, F17stable, F17snowsurface and 2 best (lowest cost) iterative solutions, Mar08-Mar09
Time series of differences with F13 for RSS F17unadjusted, F17current, F17stable, F17snowsurface and 2 best (lowest cost) solutions, Mar08-Mar09

Figure 5. Time series of SCA and SWE from F13 vs. various RSS F17 V7 options, including best "equal" and "spring" cost iterative solutions, Mar 2008 - Mar 2009. The iterative solutions are decreasing the SCA and SWE biases during fall and winter, but are still overmeasuring compared to F17current in the late spring (May-Jun).

Snow-covered area (SCA) Snow water equivalent (SWE)
Time series of F13, RSS F17unadjusted, and 2 best (lowest cost) iterative solutions, full overlap, Dec06-Mar09
Time series of differences with F13 for RSS F17unadjusted, and 2 best (lowest cost) solutions, full overlap, Dec06-Mar09

Figure 5A. Time series of SCA and SWE from F13 vs. unadjusted RSS F17 V7 and best "equal" and "spring" cost iterative solutions, full overlap, Dec 2006 - Mar 2009.

Based on these results, for the SWE algorithm run on RSS F17 V7 inputs, I am recommending that we use the algorithm derived from the lowest spring_cost iterative solution (dark pink in Figures 5 and 5A), with coefficients:

A = 4.809
B = -4.792
C = -6.036

I also shared the full overlap time series for RSS F17 V7 TBs with T. Baldwin and E. Deeb at CRREL, who produced the following figures for the area of the Amu Darya basin.

Comparisons of "best" two RSS F17 algorithms from Figure 5A to historical data as used at CRREL.
Best "equal_cost" iterative solution (green line) compared to historical time series (red line) used by CRREL.
Best "spring_cost" iterative solution (pink line) compared to historical time series (red Line) used by CRREL.
Respective differences.

Figure 5B. Amu Darya basin comparisons. The red line data is a composite of data previously provided to CRREL:

Note the period of best agreement to F17current is the green line, from 2009 day 238 onwards.

Based on these plots, the CRREL team concluded that the "green line" on RSS F17 V7 TBs provides the best continuity of SWE with the historical (RSS F13 V4, and CLASS F17 V7 (F17current)) data that they have been using in the biweekly assessments.

An Improved CLASS F17 SWE Regression

I repeated the iterative solution approach using CLASS F17 for the same time period (Figure 6), starting all iterations from F17current algorithm coefficients. I eliminated the F17stable and F17snowsurface options. Once again, the positive bias in the spring was resistant to adjustments of the SWE algorithm coefficients.

Snow-covered area (SCA) Snow water equivalent (SWE)
Time series of F13, CLASSF17unadjusted, F17current, 3 best (lowest cost) iterative solutions, Mar08-Mar09
Time series of differences with F13 for CLASS F17unadjusted, F17current, and 3 best (lowest cost) solutions, Mar08-Mar09

Figure 6. Time series of SCA and SWE from F13 vs. iterative CLASS F17 options, including best "equal", "spring" and "spring_only" cost iterative solutions, Mar 2008 - Mar 2009.

For both SCA and SWE, the best "spring_only" solution performs marginally better than F17current in spring, but is worse than F17current during the rest of the winter. Results from "spring_cost" are almost identical to those from "equal_cost". Both of these options perform better than F17current, except for SWE from January through the spring.

If we extend the SWE products (using RSS F17 V7 TBs), I recommend that we use the coefficients corresponding to the dark pink lines in Figure 5 and 5A.

Since NISE is only concerned with SCA, I recommend that we use the coefficients corresponding to dark pink line in Figure 6 to extend the NISE processing (which uses CLASS F17 TBs).

CRREL uses the daily NRT SWE in time series plots of historical SWE data. In this case, switching from F17current to best "spring_cost" option will make F17 to historical F13 comparisons better in the early part of the winter (through January), but may be overmeasuring in the late winter/early spring. When I shared these results with our CRREL colleagues Tim Baldwin and Eli Deeb, they concluded that none of the new algorithms delivers significant improvement over the F17current algorithm when the inputs are CLASS F17.

Conclusions

Modifications of Walt's iterative approach do work to improve regression results, and in most cases produce better results than from individual regressions on component TBs.

For extending the SWE record with RSS F17 V7 TBs, I recommend using these algorithm coefficients (dark pink lines in Figures 5 and 5A):

A = 4.809
B = -4.792
C = -6.036

The NISE product only uses the SCA output from this algorithm, so for NISE processing using CLASS F17 TBs, I recommend that we use the following coefficients from now on (dark pink lines in Figure 6):

A = 4.807
B = -4.786
C = -21.036

S. J. S. Khalsa has reviewed my analysis, and concurs with this recommendation (personal communication, 30 Oct 2013). Khalsa also noted that the NISE philosophy is to apply the best current algorithm, without being unduly concerned with continuity. The NISE documentation emphasizes that NISE is not intended to be used for time series analysis.

For the CRREL daily SWE, after consultation with my CRREL colleagues T. Baldwin and E. Deeb, we agreed that F17current produces the best match for SWE time series when inputs are CLASS F17 TBs. Based on comparisons to historical data (derived from RSS V4 TBs) subsetted for the Amu Darya basin, we further concluded that the best match for RSS F17 V7 SWE is the green line algorithm from Figure 5A:

A = 4.807
B = -4.792
C = -9.036

If at some future time I produce an extension to the CRREL SWE time series using RSS F17 V7 inputs, I will use this algorithm.

References

Armstong, R. L. and M. J. Brodzik. 2001. Recent Northern Hemisphere snow extent: a comparison of data derived from visible and microwave satellite sensors. GRL, 28(19):3673-3676.

Armstrong, R. L., M. J. Brodzik, K. Knowles and M. Savoie. 2005, updated 2007. Global monthly EASE-Grid snow water equivelent climatology. Boulder, CO: National Snow and Ice Data Center. Digital Media.

Brodzik, M. J., R. L. Armstrong and M. Savoie. 2007. Global EASE-Grid 8-day blended SSM/I and MODIS snow cover. Boulder, CO: National Snow and Ice Data Center. Digital Media.

Chang, A. T. C., J. L. Foster and D. K. Hall. 1987. Nimbus-7 SMMR-derived global snow cover parameters. Ann. Glac. 9:39-44.

Meier, W. N., S. J. S. Khalsa and M. H. Savoie. 2011. Intersensor calibration between F-13 SSM/I and F-17 SSMIS Near-real-time sea ice estimates. IEEE Trans. Geo. and Rem. Sens., 49(9):3343-3349.

Njoku, E. G., S. K. Chan, R. L. Armstrong, M. J. Brodzik, M. H. Savoie and K. Knowles. 2006. Stable targets for spaceborne microwave radiometer calibration. Proc. IEEE MicroRad 2006, pp. 1-6. doi:10.1109/MICRAD.1006.1677052.

Nolin, A. W., R. Armstrong and J. Maslanik. 1998. Near-real-time SSM/I-SSMIS EASE-Grid daily global ice concentration and snow extent. Version 4. Boulder Boulder, CO:NASA DAAC at the National Snow and Ice Data Center. Digital Media.

Savoie, M. H., E. Njoku, M. J. Brodzik, K. Knowles, R. L. Armstrong. 2004. Locating Stable Calibration Targets. Eos Trans. AGU, 85(47), Fall Meet. Suppl., Abstract H13C-0444.

Takala, M., K. Luojus, J. Pulliainen, C. Derksen, J. Lemmetyinen, J.-P. Karna, J. Koskinen and B. Bojkov. 2011. Estimating northern hemisphere snow water equivalent for climate research through assimilation of space-borne radiometer data and ground-based measurements. Remote Sensing of Environment, 115:3517-3529.


M. J. Brodzik <brodzik@nsidc.org>