From: Phil Jones To: Gabi Hegerl , "Michael E. Mann" Subject: Re: Mann and Jones (2003) Date: Tue Aug 10 15:47:04 2004 Cc: Tom Crowley Gabi, No second attempt - don't know what the first was? We'll be doing a new instrumental data set (surprisingly called HadCRUT3), but that's it at the moment. Attached is a good review of corals - just out. Cheers Phil At 10:36 10/08/2004 -0400, Gabi Hegerl wrote: Hi Mike and Phil, Thanks! Yes, factor 1.29 will get me closer to my best guess scaling (factor 1.6 to same-size signals). The scaling is a tough issue, and I think there are lots of possibilities to do it depending on what one wants to do. For comparing underlying forced signals, I think tls is best. To get a conservative size paleo reconstruction (like what part of instrumental do we reconstruct with paleo), the traditional scaling is best. I'll write up what Myles and I have been thinking and send it. Phil, if there is a second attempt at that with the Hadley Centre, let me know, I don't like racing anybody! Gabi Michael E. Mann wrote: Dear Phil and Gabi, I've attached a cleaned-up and commented version of the matlab code that I wrote for doing the Mann and Jones (2003) composites. I did this knowing that Phil and I are likely to have to respond to more crap criticisms from the idiots in the near future, so best to clean up the code and provide to some of my close colleagues in case they want to test it, etc. Please feel free to use this code for your own internal purposes, but don't pass it along where it may get into the hands of the wrong people. In the process of trying to clean it up, I realized I had something a bit odd, not necessarily wrong, but it makes a small difference. It seems that I used the 'long' NH instrumental series back to 1753 that we calculated in the following paper: * Mann, M.E., Rutherford, S., Bradley, R.S., Hughes, M.K., Keimig, F.T., [1]Optimal Surface Temperature Reconstructions using Terrestrial Borehole Data, Journal of Geophysical Research, 108 (D7), 4203, doi: 10.1029/2002JD002532, 2003. (based on the sparse available long instrumental records) to set the scale for the decadal standard deviation of the proxy composite. Not sure why I used this, rather than using the CRU NH record back to 1856 for this purpose. It looks like I had two similarly named series floating around in the code, and used perhaps the less preferable one for setting the scale. Turns it, this has the net effect of decreasing the amplitude of the NH reconstruction by a factor of 0.11/0.14 = 1.29. This may explain part of what perplexed Gabi when she was comparing w/ the instrumental series. I've attached the version of the reconstruction where the NH is scaled by the CRU NH record instead, as well as the Matlab code which you're welcome to try to use yourself and play around with. Basically, this increases the amplitude of the reconstruction everywhere by the factor 1.29. Perhaps this is more in line w/ what Gabi was estimating (Gabi?) Anyway, doesn't make a major difference, but you might want to take this into account in any further use of the Mann and Jones series... Phil: is this worth a followup note to GRL, w/ a link to the Matlab code? Mike p.s. Gabi: when do you and Tom plan to publish your NH reconstruction that now goes back about 1500 years or so? It would be nice to have more independent reconstructions published in the near future! Maybe I missed this? Thanks... ______________________________________________________________ Professor Michael E. Mann Department of Environmental Sciences, Clark Hall University of Virginia Charlottesville, VA 22903 _______________________________________________________________________ e-mail: [2]mann@virginia.edu Phone: (434) 924-7770 FAX: (434) 982-2137 [3]http://www.evsc.virginia.edu/faculty/people/mann.shtml % COMPOSITENH" % % (c) 2003, M.E. Mann % % THIS ROUTINE PERFORMS A RECONSTRUCTION OF NORTHERN HEMISPHERE % MEAN ANNUAL TEMPERATURE BASED ON A WEIGHTED COMPOSITE OF LONG-TERM TEMPERATURE % PROXY RECORDS SCALED AGAINST THE INSTRUMENTAL HEMISPHERIC MEAN TEMPERATURE % SERIES, AS USED IN THE FOLLOWING TWO PUBLICATIONS: % % % Jones, P.D., Mann, M.E., Climate Over Past Millennia, Reviews of Geophysics, % 42, RG2002, doi:10.1029/2003RG000143, 2004 % % Mann, M.E., Jones, P.D., Global Surface Temperatures over the Past two Millennia, % Geophysical Research Letters, % 30 (15), 1820, doi: 10.1029/2003GL017814, 2003 % % % 1. READ IN INSTRUMENTAL RECORD % % Read in CRU instrumental NH mean temeperature record (1856-2003) load nh.dat; yearinstr=nh(:,1); % calculate both warm-season and annual means warmseason=(nh(:,5)+nh(:,6)+nh(:,7)+nh(:,8)+nh(:,9)+nh(:,10))/6; annualmean=nh(:,14); % use annual mean record in this analysis nhmean=annualmean; % % 2. READ IN PREVIOUSLY PUBLISHED PROXY-RECONSTRUCTIONS OF NH ANNUAL MEAN % RECONSTRUCTIONS AND FORM APPROPRIATELY SCALED COMPOSITE % % Read in Mann et al (1998), Crowley and Lowery (2000), and Jones et al (1998) % NH temperature reconstructions load nhem-millennium.dat; load crowleylowery.dat; load joneshemisrecons.dat; nhmbh=nhem_millennium(1:981,2); nhjones=joneshemisrecons(1:981,2); nhcl=crowleylowery(1:981,2); yearmillen=nhem_millennium(1:981,1); % since some reconstructions are only decadally resolved, smooth each on % decadal timescales through use of a lowpass filter with cutoff at % f=0.1 cycle/year. Based on use of the filtering routine described in: % % Mann, M.E., On Smoothing Potentially Non-Stationary Climate Time Series, % Geophysical Research Letters, 31, L07214, doi: 10.1029/2004GL019569, 2004. % % using 'minimum norm' constraint at both boundaries for all time series nhsmooth=lowpass(nhmean,0.10,0,0); nhmbhsmooth=lowpass(nhmbh,0.10,0,0); nhjonessmooth=lowpass(nhjones,0.10,0,0); nhclsmooth=lowpass(nhcl,0.10,0,0); % Mann et al (1998) already calibrated in terms of hemispheric annual mean temperature, but % reference mean has to be adjusted to equal that of the instrumental series % over the 1856-1980 overlap period (which uses a 1961-1990 reference period) admbh=mean(nhsmooth(1:125))-mean(nhmbhsmooth(857:981)); newmbh=nhmbhsmooth+admbh; % need to adjust and scale Jones et al (1998) and Crowley and Lowery (2000) % reconstructions to match mean and trend of smoothed instrumental series % over 1856-1980 t1=1856; t2=1980; x=(t1:t2)'; nhlong=nhmean(1:125); smoothlong=lowpass(nhlong,0.10,0,0); amean0=mean(smoothlong); y=smoothlong; [yc,t,trend0,detrend0,xm,ym] = lintrend(x, y); % y=nhclsmooth(t1-999:t2-999); [yc,t,trendcl,detrendcl,xm,ym] = lintrend(x, y); % y=nhjonessmooth(t1-999:t2-999); [yc,t,trendjones,detrendjones,xm,ym] = lintrend(x, y); % multjones=norm(trend0)/norm(trendjones); adjustedjones=nhjonessmooth*multjones; offsetjones=amean0-mean(adjustedjones(t1-999:t2-999)); newjones=adjustedjones+offsetjones; newjones=newjones'; % multcl=norm(trend0)/norm(trendcl); adjustedcl=nhclsmooth*multcl; offsetcl=amean0-mean(adjustedcl(t1-999:t2-999)); newcl=adjustedcl+offsetcl; newcl=newcl'; % nhlongcompose=0.3333*(newmbh+newjones'+newcl')'; % % 3. READ IN AND PROCESS PROXY TEMPERATURE RECORDS % M=8; load 'china-series1.dat' load 'itrdb-long-fixed.dat' load 'westgreen-o18.dat' load 'torny.dat' load 'chesapeake.dat' load 'mongolia-darrigo.dat' load 'dahl-jensen-gripbh1yrinterp.txt' load 'dahl-jensen-dye3bh1yrinterp.txt' % read in years x1=china_series1(:,1); x2=itrdb_long_fixed(:,1); x3=westgreen_o18(:,1); x4=torny(:,1); x5=chesapeake(:,1); x6=mongolia_darrigo(:,1); x7=dahl_jensen_gripbh1yrinterp(:,1); x8=dahl_jensen_dye3bh1yrinterp(:,1); % read in proxy values y1=china_series1(:,2); y2=itrdb_long_fixed(:,2); y3=westgreen_o18(:,2); y4=torny(:,2); y5=chesapeake(:,2); y6=mongolia_darrigo(:,2); y7=dahl_jensen_gripbh1yrinterp(:,2); y8=dahl_jensen_dye3bh1yrinterp(:,2); % Store decadal correlation of each proxy record with local available % overlapping CRU gridpoint surface temperature record (see Mann and Jones, 2003) corr(1)=0.22; corr(2)=0.52; corr(3)=0.75; corr(4)=0.32; corr(5)=0.31; corr(6)=0.40; corr(7)=0.53; corr(8)=0.52; % Estimate Area represented by each proxy record based on latitude of % record and estimated number of temperature gridpoints represented by record pi=3.14159; factor=pi/180.0; lat(1)=32.5; dof(1)=4; lat(2)=37.5; dof(2)=2; lat(3)=77; dof(3)=0.667; lat(4)=68; dof(4)=3.5; lat(5)=37.0; dof(5)=1.0; lat(6)=47; dof(6)=1; lat(7)=73; dof(7)=0.667; lat(8)=65; dof(8)=0.667; for j=1:M area(j)=dof(j)*cos(lat(j)*factor); end % determine min and max available years over all proxy records % minarray=[min(x1) min(x2) min(x3) min(x4) min(x5) min(x6) min(x7) min(x8)]; maxarray=[max(x1) max(x2) max(x3) max(x4) max(x5) max(x6) max(x7) max(x8)]; tbegin=max(minarray); tend1=min(maxarray); tend=max(maxarray); % initialize proxy data matrix notnumber = -9999; for j=1:M for i=1:minarray(j)-1 time(i)=i; mat(i,j)=notnumber; end for i=minarray(j):tend time(i)=i; end for i=minarray(j):maxarray(j) if (j==1) mat(i,j)=y1(i-minarray(j)+1); end if (j==2) mat(i,j)=y2(i-minarray(j)+1); end if (j==3) mat(i,j)=y3(i-minarray(j)+1); end if (j==4) mat(i,j)=y4(i-minarray(j)+1); end if (j==5) mat(i,j)=y5(i-minarray(j)+1); end if (j==6) mat(i,j)=y6(i-minarray(j)+1); end if (j==7) mat(i,j)=y7(i-minarray(j)+1); end if (j==8) mat(i,j)=y8(i-minarray(j)+1); end end % added in Jones and Mann (2004), extend series ending between % 1980 calibration period end and 2001 boundary by persistence of % last available value through 2001 for i=maxarray(j)+1:tend if (j==1) mat(i,j)=y1(maxarray(j)-minarray(j)+1); end if (j==2) mat(i,j)=y2(maxarray(j)-minarray(j)+1); end if (j==3) mat(i,j)=y3(maxarray(j)-minarray(j)+1); end if (j==4) mat(i,j)=y4(maxarray(j)-minarray(j)+1); end if (j==5) mat(i,j)=y5(maxarray(j)-minarray(j)+1); end if (j==6) mat(i,j)=y6(maxarray(j)-minarray(j)+1); end if (j==7) mat(i,j)=y7(maxarray(j)-minarray(j)+1); end if (j==8) mat(i,j)=y8(maxarray(j)-minarray(j)+1); end end end time=time'; data=[time mat]; % decadally lowpass of proxy series at f=0.1 cycle/year as described earlier for j=1:M unfiltered=mat(minarray(j):tend,j); filt=lowpass(unfiltered,0.1,0,0); for i=1:minarray(j)-1 filtered(i,j)=mat(i,j); end for i=minarray(j):tend filtered(i,j)=filt(i-minarray(j)+1); end end % standardize data % first remove mean from each series for j=1:M icount=0; amean(j)=0; for i=1:tend if (filtered(i,j)>notnumber) icount=icount+1; amean(j)=amean(j)+filtered(i,j); end end amean(j)=amean(j)/icount; end % now divide through by standard deviation for j=1:M icount=0; asum=0; for i=1:tend if (filtered(i,j)>notnumber) asum=asum+(filtered(i,j)-amean(j))^2; icount=icount+1; end end sd(j)=sqrt(asum/icount); for i=1:tend standardized(i,j)=filtered(i,j); if (mat(i,j)>notnumber) standardized(i,j)=(filtered(i,j)-amean(j))/sd(j); end end end % % 4. Calculate NH mean temperature reconstruction through weighted (and % unweighted) composites of the decadally-smoothed proxy indicators % % impose weighting scheme for NH mean composite for j=1:M % weighting method 1: weight each proxy series by approximate area % weighting method 2: weight each proxy series by correlation between % predictor and local gridpoint series over available overlap period % during calibration interval % weighting method 3: weight each proxy series by correlation between % predictor and NH mean series over calibration interval: % weightlong(j)=lincor(nhlong,standardized(1856:1980,j)); % weighting method 4: combine 1 and 3 % weighting method 5: combine 1 amd 2 (this is the 'standard' weighting % scheme chosen by Mann and Jones (2003) % use standard weighting scheme weight(j)=corr(j)*area(j); end % perform reconstructions based on: % (1) the 6 proxy temperature records available over interval AD 200-1980 % (2) all 8 proxy temperature records available over interval AD 553-1980 istart0=200; istart1=200; istart2=553; nseries1=0; nseries2=0; weightsum1=0; weightsum2=0; for j=1:M if (istart1>=minarray(j)) nseries1=nseries1+1; weightsum1=weightsum1+weight(j); end if (istart2>=minarray(j)) nseries2=nseries2+1; weightsum2=weightsum2+weight(j); end end % calculate composites through 1995 (too few series available after that date) % As discussed above, persistence is used to extend any series ending % between 1980 and 1995 as described by Jones and Mann (2004). tend=1995; for i=istart1:tend unweighted1(i)=0; unweighted2(i)=0; weighted1(i)=0; weighted2(i)=0; for j=1:M if (istart1>=minarray(j)) unweighted1(i)=unweighted1(i)+standardized(i,j); weighted1(i)=weighted1(i)+weight(j)*standardized(i,j); end if (istart2>=minarray(j)) unweighted2(i)=unweighted2(i)+standardized(i,j); weighted2(i)=weighted2(i)+weight(j)*standardized(i,j); end end end unweighted1=unweighted1/nseries1; unweighted2=unweighted2/nseries2; weighted1=weighted1/weightsum1; weighted2=weighted2/weightsum2; unweighted1(1:istart1-1)=0; unweighted2(1:istart2-1)=0; weighted1(1:istart1-1)=0; weighted2(1:istart2-1)=0; % scale composite to have same variance as decadally-smoothed instrumental % NH series % Mann and Jones (2003) and Jones and Mann (2004) used for this purpose % the extended (1753-1980) NH series used in: % Mann, M.E., Rutherford, S., Bradley, R.S., Hughes, M.K., Keimig, F.T., % Optimal Surface Temperature Reconstructions using Terrestrial Borehole Data, % Journal of Geophysical Research, 108 (D7), 4203, doi: 10.1029/2002JD002532, 2003. % That series has a decadal standard deviation sd=0.1123 % If instead, the 1856-2003 CRU instrumental NH mean record is used, with % a decadal standard deviation of sd=0.1446, the amplitude of the reconstruction % increases by a factor 1.29 (this scaling yields slightly lower verification % scores) load nhem-long.dat nhemlong=nhem_long(:,2); longsmooth=lowpass(nhemlong,0.10,0,0); sd0=std(longsmooth); % use weighted (rather than unweighted) composite in this case series1=weighted1; % center composites on 1856-1980 calibration period y=series1(t1:t2)'; amean1=mean(series1(t1:t2)); compseries1=series1(t1:t2)-amean1; mult1=sd0/std(compseries1); % scale composite to standard deviation of instrumental series and re-center % to have same (1961-1990) zero reference period as CRU NH instrumental % temperature record adjusted1=series1*mult1; offset1=amean0-mean(adjusted1(t1:t2)); compose1=adjusted1+offset1; compose1=compose1'; series2=weighted2; y=series2(t1:t2)'; amean2=mean(series2(t1:t2)); compseries2=series2(t1:t2)-amean2; mult2=sd0/std(compseries2); adjusted2=series2*mult2; offset2=amean0-mean(adjusted2(t1:t2)); compose2=adjusted2+offset2; compose2=compose2'; % % 5. UNCERTAINTY ESTIMATION, AND STATISTICAL VERIFICATION % % estimate uncertainty in reconstruction % nominal (white noise) unresolved calibration period variance calibvar=lincor(smoothlong,compose1(t1:t2))^2; uncalib=1-calibvar; sdunc=sd0*sqrt(uncalib); % note: this is the *nominal* white noise uncertainty in the reconstruction % a spectral analysis of the calibration residuals [as discussed briefly in % Mann and Jones, 2003] indicates that a peak at the multidecadal timescale % that exceeds the white noise average residual variance by a factor of % approximately 6. A conservative estimate of the standard error in the % reconstruction thus inflates the nominal white noise estimate "sdunc" by a % factor of sqrt(6) sdlow = sdunc*sqrt(6) % calculate long-term verification statistics for reconstruction % use composite of Mann et al (1998)/Crowley and Lowery (2000)/Jones et al (1998) % and AD 1600-1855 interval overlapcomp=nhlongcompose(1:981); % work with longer reconstruction (back to AD 200) overlaprecon=compose1(1000:1980)'; %overlaprecon=compose2(1000:1980)'; %calculate verification R^2 series11=overlaprecon(601:856); series22=overlapcomp(601:856); verifrsq=lincor(series11,series22)^2 % calculate verification RE var1=0.0; var2=0.0; var3=0.0; var4=0.0; var5=0.0; am0=0.0; % insure convention of zero mean over calibration interval for i=857:981 am0=am0+overlapcomp(i); end am0=am0/125; for i=601:856 var1=var1+(overlapcomp(i)-am0)^2; var2=var2+(overlapcomp(i)-overlaprecon(i))^2; end verifRE=1-var2/var1 -- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Gabriele Hegerl Division of Earth and Ocean Sciences, Nicholas School for the Environment and Earth Sciences, Box 90227 Duke University, Durham NC 27708 Ph: 919 684 6167, fax 684 5833 email: [4]hegerl@duke.edu, [5]http://www.env.duke.edu/faculty/bios/hegerl.html Prof. Phil Jones Climatic Research Unit Telephone +44 (0) 1603 592090 School of Environmental Sciences Fax +44 (0) 1603 507784 University of East Anglia Norwich Email p.jones@uea.ac.uk NR4 7TJ UK ---------------------------------------------------------------------------- References 1. ftp://holocene.evsc.virginia.edu/pub/mann/borehole-jgr03.pdf 2. mailto:mann@virginia.edu 3. http://www.evsc.virginia.edu/faculty/people/mann.shtml 4. mailto:hegerl@duke.edu 5. http://www.env.duke.edu/faculty/bios/hegerl.html