Check my arithmetic

In the last post, I forecast next year’s September sea ice extent using not only extent data, but latitude of the ice edge. To compute that from extent, I applied a formula given by Eisenman in the supplemental info of his publication.

A reader pointed out that there seemed to be a discrepancy between his numerical values and mine — and he’s right. I’ve double-checked the approximation formula and my program, but I can’t resolve the discrepancy.

Therefore, I’d be grateful if some of you would reproduce the calculation of latitude based on extent data, using Eisenman’s approximation formula. The formula as given is this:

The data I used for sea ice extent is here (it’s monthly averages from NSIDC):


If you could report the range of latitude values you get from the extent data in that file, using the formula as given by Eisenman, that would confirm or deny that I’ve implemented his formula as stated.


P.S. The figure I used for the radius of the earth is 6371 km.

About these ads

15 responses to “Check my arithmetic

  1. I get max 79.855, min 60.269.

    [Response: Thanks -- that's exactly what I got. That means that either we both implemented his formula correctly or we both made a mistake which ended up giving us identical results.

    But those numbers don't agree with the data used in Eisenman's paper, which range from 65.382 to 80.689. There are some differences, namely that he used daily data whereas these are monthly, and his data only go through Jan. 2010. Also, he states that the ice edge latitude was *not* computed by this formula, that it's just a convenient approximation.

    But -- I used his formula to compute ice edge latitude for daily data and compared that to his stated values. Not only do they not match, there's a consistent bias between the results. The biggest discrepancy is on the low end -- his latitude range bottoms out at a little over 65N (winter max) whereas we both got down to near 60N. Also, if I plot his values for extent vs. latitude, the scatterplot does *not* follow the approximate formula -- it's really not that close.

    I'm puzzled. Perhaps I should email Eisenman and ask his opinion. In the meantime, if others want to check -- there is still a small chance that we both made the same mistake (but at this point I doubt it).]

    • Same numbers here.

      I’m not quite happy with the 0.49 value chosen by Eisenman from a numerical perspective, I think he cut off the number too early. I am assuming that Eisenmann intended the formulas to smoothly blend into each other. However, if you evaluate the polynomial in the denominator at 16 * 10^6 km² it has a value of 0.4969420… . Cutting it off as 0.49 for the simplified formula introduces a jump of about 0.2° at this threshold value between the latitudes of the ice edge as calculated by the long and the short formula, respectively. The resolution of the sea ice extent data implies a resolution of the latitudes of the ice edge of 0.01°, a factor of 20 less than the original jump. Using the rounded value of 0.497 instead reduces the jump to insignificant 0.0016°.

      So, assuming the complex formula is the better approximation and the simplified formula is supposed to be a smooth continuation I would recommend to use 0.497 (or go the extra mile and use 0.496942) instead of 0.49 in the simplified formula.

  2. Yeah, after coding this I looked at the previous post and found the same discrepancies. I’d say it’s worth an email to Eisenman.

    [Response: Sent.]

  3. I get max 79.855, min 60.269 as well.
    Had different min at first, then realised I forgot to calc the lat differently for large extents.

  4. Looking into details, the formula seems to be pretty rough. o1..o5 are the five terms in the parentheses of the given formula, o1 being the order 1 term of 1.56*eps. The contributions of these terms are given below for extends in the range of 4.0 to 8.0e6 km^2. lat4 uses the terms up to o4 to estimate the latitude. The order 6 term would be a negative one, meaning the formula underestimates the latitude in the range of interest. The most optimistic guess on the latitude uncertainty would be at least 1 deg.

    ext: 4.0
    o1: 0.1537, o2: -0.5866, o3: 0.4383, o4: -0.1289, o5: 0.0129
    lat5: 79.1, lat4: 79.2

    ext: 5.0
    o1: 0.1922, o2: -0.9165, o3: 0.8561, o4: -0.3148, o5: 0.0394
    lat5: 77.4, lat4: 77.7

    ext: 6.0
    o1: 0.2306, o2: -1.3198, o3: 1.4794, o4: -0.6527, o5: 0.0980
    lat5: 75.5, lat4: 76.4

    ext: 7.0
    o1: 0.2690, o2: -1.7964, o3: 2.3492, o4: -1.2092, o5: 0.2119
    lat5: 72.8, lat4: 75.2

    ext: 8.0
    o1: 0.3075, o2: -2.3463, o3: 3.5066, o4: -2.0628, o5: 0.4131
    lat5: 67.3, lat4: 74.1

    The formula is probably not really useful to get a realistic latitude of the ice edge with the accuracy you look for. It’s hard to tell if it is a useful statistical measure to forecast sea ice extend. It’s biased, but the bias vanishes with the decreasing extend.

    Further, the formula has problems for high ice extends:

    ext: 14.0, lat5: 64.9
    ext: 15.0, lat5: 62.5
    ext: 16.0, lat5: 60.9
    ext: 17.0, lat5: 61.3
    ext: 18.0, lat5: 64.2

    I’m afraid it’s not really useful for these extends. The zero order alternative offered in the paragraph might help a little bit.

    ext: 14.0, lat5: 64.9, lat0: 62.6
    ext: 15.0, lat5: 62.5, lat0: 61.6
    ext: 16.0, lat5: 60.9, lat0: 60.7
    ext: 17.0, lat5: 61.3, lat0: 59.8
    ext: 18.0, lat5: 64.2, lat0: 58.9

    Switching the formula for some extends effects the bias of the latitude anomaly.

  5. “Also, he states that the ice edge latitude was *not* computed by this formula, that it’s just a convenient approximation.” Perhaps the ice edge latitude values in his paper are determined by some observational criteria, whereas that formula is simply an approximation of geography? [Response: Yes, I think that's the case.] (OTOH, if the formula assumes solid ice above some latitude and no ice below it, I’d expect the formula to give higher latitudes than using some observational criteria, which is the opposite of what you found. Curiouser…)

  6. David B. Benson

    60N is Nunivak in the Bering Sea and Cape Farewell, the southern tip of Greenland, in the North Atlantic.

    I doubt that in the 21st century.

  7. FWIW, the “equivalent extent” column in Eisenman’s dataset
    seems to have been computed using the following equation:
    lat = arcsin(1-ε/2π)

  8. Also FWIW, a reasonably good fit to Eisenman’s data can be found with:
    lat = arcsin(1-ε/(2π(2.2-26.91ε+204ε^2-607ε^3+608ε^4)))
    … for all extents. A fifth order equation does little to improve the fit.

  9. Doesn’t the radius of the earth depend on where you are?

    • Yep, but for almost all calculations you can throw in your favorite average radius and be accurate to within a few %. The difference between the polar and equatorial radii for earth is only 20 km or .3 %….

  10. We are starting to get into the nitty gritty here, so has anyone one taken the shading effect of Greenland into account? It casts a very long shadow and so significantly affects insolation for a period when the sun is low in the sky.

  11. Final word (maybe):
    Eyeballing Eisenman´s data 1978-2010 shows that it does seem to be a two-domain function, but the breakpoint occurs at about 10 million km², not 16 million. (typo?).

    Using gnuplot´s fitting, for extents 9.8
    lat = arcsin(1-ε/(2π(1.642-2.076ε)))

  12. … trying again …

    Using gnuplot´s fitting, for extents greater than or equal to 9.8e6, I get
    lat = arcsin(1-ε/(2π(4.546-61.94ε+361.2ε²-772ε³+401ε^4)))
    … and for extents less than 9.8
    lat = arcsin(1-ε/(2π(1.642-2.076ε)))

  13. OffTopic… continuing on gallows humor on … Looking at US Drought Monitor by state, it is clear that Mississippi state is an outlier in the correlation between the amount of republican representatives in the state senate and the drought area of the respective state. What are Mississippians doing wrong?