This post is intended to describe my conceptual idea of one might relate the constants Nick, Arthur and Tamino obtain by regressing time series for the earth’s surface to forcings from a GISS temp file. It’s a bit fuzzy, and may contain typos/ algebra errors etc. My goal is to illustrate the basic outline of the analysis, and reduce the number of scraps of paper on my desk.
Though, the analysis may contain both typos and algebra errors, I think it will explain the idea. I’m still not 100% sure I will be able to get this to the point of solving for all the parameters…. but I think I have a sufficient number of equations to solve for all unknowns. Given the nature of the post, it is likely to be rife with notational defects, and will contain the minimum amount of explanation. I’m trying to maintain the notation I used in the previous post on the two box model. The post will be periodically updated as I find typos (which are surely exist) and algebra errors, which may also exist.
Regression
If I understand Nick & Tamino’s explanation of how they find the best fit parameters to fit global average surface temperature from the historic record, they both assume that for a two-box model with two time constants, τ+ and τ–, temperature can be obtained from the total forcing applied to the box by integrating as follows:
Ts(t)= ws+
∫ ds { exp[-(t-s)/τ+] w+ Cs GT(s) + [exp[-(t-s)/τ–] w– Cs GT(s) }
To(t)=To,0+
∫ ds { exp[-(t-s)/τ+] w+r+ Cs GT(s) + [exp[-(t-s)/τ–] w–r– Cs/ GT(s) } +
where
- In both equations t and s both represent time with the limits for integration from s=0 to s=t,
- Cs represents the heat capacity per unit area for the “surface” box,
- τ+ and τ– are time constants for the two-box model; their values are assumed to be known;
- w+ and w– are weighting factors and have dimensions [temperature * area / energy.]
- To,0 are ws the temperature of the ocean and surface boxes at equilibrium when the forcings are zero; (the odd symbol for the surface is selected because it also a fitting parameter in Tamino’s procedure as eluciated by Nick Stokes.).
- (1, r±) are the eigenvectors for the τ± solutions respectively and can be obtained algebraically when the parameters for the box model are fully specified. (That is: We have an equation relating r+ to the parameters for the box model and another relating r–).
- GT is related to the total force FT applied to the two boxes and GT is related to the forcings applied to each box, and will be defined further below.
- Additional terms representing the decay of the initial temperature to equilibrium value, but Nick’s method zero pads with zero temperatures prior to t=0. This means the solution method contains the hidden assumption that the temperatures in both boxes were near equilibrium at time=0.
The total forcing applied to the earth’s climate, FT, is the sum of the forcing added to the ocean and surface boxes.
FT = Fs+ Fo
Nothing about the curve fitting exercise proposed by Tamino constraints the forces applied to the boxes any further than indicated with (2)
For this analysis, we will also assume that Fs = x FT where x is a coefficient of proportionality. Note that x=1 indicates all of the heat is added to the top box; x=0 indicates all heat is added to the lower box. For values of x outside the range 0≤x ≤1, one box will be heated when the other is cooled. (This violates my physical intuition about how global warming would work, but we will defer that discussion until a later time.)
The quantity “GT ” represents the forcing normalized by heat capacity, and for the purpose of this presentation will be defined as
GT = FT/Cs; Gs = Fs/Cs and Go = Fo/Co
Using (3)
Go = (1-x) GT (Cs/Co)
Before proceeding, it is worth recognizing the products w+ are related to the two of the regression coefficients obtained when Nick fits the temperature data for box s to the total forcing; ws is the third regression parameter obtained from the fit; in otherwords, for the purpose of the current analyses, these quantities have known values. (I need to sort out what exactly R does to determine the relationship because Arthur’s discussion suggests the regression coefficients have dimensions (Temperature * Area/ Power); either is fine, but I’m not sufficiently familiar with the packages in R to know precisely how their constants relate to mine. With luck, someone will tell me.)
Leibniz Rule
We will now apply Leibniz rule to obtain express the time derivative of temperatures Ts and Ts in terms parameters contained in (1) and (2). With luck, I will not screw this up. Here goes:
dTs/dt =
∫ ds { exp[-(t-s)/τ+] (-w+/τ+) Cs GT(s) + [exp[-(t-s)/τ–] (-w–/τ–) CsGT(s)}
+ (w+ Cs GT + w– Cs GT)
Recall, that in Monday’s post, we saw a different expression for the time dTs/dt which was imposed by applying conservation of energy to the “surface” box. This equation required
dTs/dt = – (αs+γs)Ts + γsTo + Gs
Having two different equations for dTs/dt probably looks like a mess at this point.
However it’s actually the key to determing the magnitude of the parameters for the box model based on the regression constants. If we recognize that the post dTs/dt as expressed in (1c from ref. must be equal to (6) above, we obtain four, count them, four equations!
How? If we substitute in the equations for Ts and To from (1) and (2) and think, we will recognize that to ensure dTs/dt for both equations match for any and all values of ‘t’, then
- the coefficients on the terms multiplying exp[-(t-s)/τ–] must match in both equations. That gives us 1 equation.
- the coefficients on the terms multiplying exp[-(t-s)/τ+] must match in both equations. That gives us a second equation.
Decomposing the forcing into a constant plus a time varying term G = GT,ave+ G’T(t),the constants on in the two expression must match and- The coefficients multiplying the remaining time varying forcings GT(t) must match.
[Update: The decomposition into mean and time varying was a boo-boo.]
So, for example, setting the terms that do not vary with time equal, we obtain:
(w++ w–) Cs GT,ave =
– (αs+γs) ws+ γs To,0 + x GT,ave
The terms on the left hand side are obtained from (6); those on the right hand side are obtained from (1c from ref.)
Possibly more interesting, x, the fraction of energy total energy going into the top box is setting the terms multiplying G’T(t) equal
Although we do not yet know the magnitude of all parameters, equation 8a suggests if know the magnitude of (w++ w–) from the regression and have an estimate for the heat capacity per unit area, we can estimate the fraction of heat going into the two boxes for any given set of fitting parameters. (Or course, it’s even nice to solve for the heat capacity consistent with the solution. But, even simply having an estimate permits us to do a partial reality check.)
I’ll defer computation of the estimate so that the non-math averse can find algebraic errors (which I fear I must have made.) I’m also trying to sort out the relationship between my constants and the regression coefficients reported in comments of the previous post.
Can we do more
Obviously, we can get two more equations by balacing term multiplying the exponentials. Also, we can apply the same procedure to the equations for dTo/dt to obtain 4 more equations giving us 8 equations to restricting the value of the parameters (αs, αo, γs=γo(Cs/Co), x, r+, r–, Cs, Cs/Co, To,0). Recollecting that the solution to for the eigenvectors for two box problem gives us algebraic relations for τ+ and τ+, we find we have 10 equations and 9 unknowns, so I suspect I am overlooking an unknown (or some of the equations will turn out to be linearly dependent on one another. Update: Or, Fs= x FT is overly restrictive and I need a relationship with 2 parameters!)
In the meantime, I think this illustrates:
- Why I think I need to know all three of the coefficients from the regression; the reason is all three appear in the set of equations required to solve for the parameters. (The reason can be explained physically btw.)
- Roughly, how one could (if they make no algebra mistakes) obtain the parameters for the two box model from the regression parameters for the air temperature.
Is there more?
As some of you may have guessed, there is more to come. However, it’s going to be slow, because, as you can see, I will be doing the algebra to obtain a set of equations to solve for the unknowns. The current post being posted today to isolate the discussion of the solution method from the solution for the ocean temperature, to show why I think I have enough information to compute the ocean tempeatures, to explain why I am asking for certain parameters in comments and (even more importantly) to keep from losing this much of the analysis under the piles of stray algebraic computations on my desk! (It’s a mess.)
Also… I wouldn’t mind if people happened to note notational errors. Typing equations in html can be a Pain In The A**.
Meanwhile, more will be coming. (In fact, I may have given others enough to figure out bits of the solution.)
Lucia, let’s start at the beginning with equations 1 and 2. What Nick and Tamino and I did was fit to *one* global mean surface temperature (the GISS temp record). Not two. So – what do you believe is the relationship between the T that we fit to and the two temperatures T_s and T_o in equations 1 and 2? There is no T_s or T_o in the fitting process.
Yes, I agree that any underlying physical two-box model has two boxes that may have physically distinct temperatures. But my understanding is that the only piece of the two-box model actually used in the fitting process is the result that each of the temperatures (T_s and T_o) becomes linearly dependent on two exponentially smoothed versions of the forcing, with the two different time constants and associated coefficients. Any linear combination of T_s and T_o will have the same property. So GMST could be T_s, it could be T_o, it could be some linear combination of the two. But without a *second* observational temperature record, the fit doesn’t tell you enough to split back to separate underlying two-box temperatures.
Oh – on the question of the relationship between the “sensitivity” coefficients coming out of the R analysis and the w+, w- here, the exponential smoothing you use in eq 1 is not normalized – you have to divide by
∫_0^&infty; exp[-t/τ] dt = τ
(in Nick’s R script this is done by dividing ‘e’ by sum(e) in the definition of expsmooth) – i.e. there’s another factor of τ_+ or τ_- in your equations if you want to use the properly normalized w+, w-.
Arthur–
I agree you/ Tamino/ Nick fit one temperature not two when computing your fitting parameters. The second equation represents the temperature for the box for which you did not have data. This represents no problem to my analysis.
Of course. That’s consistent with the analysis in this post.
However, since no one has specifically specified the boxes, we are free to decree that the “surface box” is whatever whose temperature can be represented by GISS. Then the rest of the climate system is in the “other” box.
We can leave that somewhat undefined,and find the parameters matching that system. After wards, we can decide whether any partitioning of the earth’s climate system makes sense relative to the parameters.
If none exists, that’s a problem.
None sense.
The regression was done on “a box”, which we will call “box s” using GISS Temp. So, GISSTemp goes with”box s”. We may not know what’s in box s. But whatever it is, the rest of the climate system is in “box o”. So, we can compute the representative temperature of “box o”.
Arthur–
Thanks on the tau question. Also, I’ve been scrolling through comments– I find a lot of numbers and graphs, but few units. What units are you getting for w?
I guess I could back it out, but I don’t know if the time steps are in months, years, seconds etc. I don’t know if the force is in W/m^2, or KW/m^2 etc.
I recognize Tamino’s created a crude physical model reflecting basic properties of GISS E. He feeds in the very forcing data that made GISS E’s 20th century simulation look so nice and obtains a temperature time series response that makes his model look almost as nice as GISS E. Now isn’t that amazing?
Wonder who measured black carbon forcing and sulphur aerosol forcing in the stratosphere 100 years ago .. wonder why land use forcing abruptly stopped around 1990 after 100 years of steady decline. Wonder why so many forcing and trends came to a sudden stall around the millenium.
Lucia,
I think there are some things to note here. Arthur noted that we’re really only dealing with Eq 1 here. Eq 2 may be useful, but needs to be connected somehow.
In Eq (1) (and 2) I think the sign of the exponent is wrong; it should be exp[(t-s)… The limits on the integral are from -ê to t (which is why I had to specify zero padding pre 1880).
Apart from the sign (the consequences of which can be corrected just by changing the sign of the Ï„’s), I think your Liebniz rule integration is correct.
As Arthur mentioned, there’s a need to normalise the smoother – you want to be sure that the smoothed version of 1 is 1. That ensures that the fitted coefficients do add to make the sensitivity.
Your equation that follows (6) shouldn’t be different – the convolved exponential appears just because that’s the solution of the inhomogeneous de which that equation represents. Your Liebniz rule differentiation should be just verifying this.
That’s as far as I’ve got at the moment. I have to leave the keyboard for a while (it’s early morning here) but hope to get back to it. It’s a very interesting post.
The aerosol forcing smells of circularity. Unless someone can show me where they were measured independently, directly or by proxy, the possibility remains that they were tuned to make Model E fit the observed temperature anomalies.
As far as looking for algebra errors, at this point my eyes glaze over when I look at that much math.
The time steps are years. The forcing is W/m2. The anomaly table is in centidegrees, but it’s converted to degrees when read into the script so the weighting factors are in units of Km2/W and the intercept unit is K.
Lucia,
This might prove useful in the future: http://wordpress.org/extend/plugins/latex/
Nick–
In a two box model, equation 1 and 2 are connected. Both equations exist because both boxes exist. Period. If the functional form for the regression is motivated by a two-box model, you can’t get around that.
The concept of a two-box model is that a system consists of two relatively homogeneous regions with relatively homogeneous properties that are coupled thermodynamically. One can argue wither or not this will represent a realistic approximation of the atmosphere. One can predict the two-box model will fail to represent the earth’s climate. But this is the physics of the two box model and that’s the models whose solution is represented by the form Tamino wants’ to fit.
Now, as for the “linear combination” business.
In the two box model a measurement devise like a thermometer is either “in” one box or or it’s “in” the other. So, you can’t measure a “linear combination” of temperature. Your measurements either represent on box or the represent the other box. The measurement cannot be a “linear combination” of the temperature from both boxes. That is cuckoo-nonsense.
It is true that every linear combination of two temperatures can be written as integrals similar to (1) and (2), but that’s irrelevant to fact that a measurement is in one box or the other. In a two-box model, the surface temperature either represents the temperature in one box or the other the other box. It can’t be a “linear combination” of the two temperatures!
Now I admit some people are going to say: Well two boxes are too simple to use for this purpose. Maybe they are right. Maybe you need a more complex model. Maybe you need an EBM with an upwelling diffusion model for the ocean and something a two layer atmosphere.
If you want a diffusive model or upwelling diffusive model, with the air temperature at some “point”, propose that model. If you want more boxes to deal with temperature gradients do that. Each of these have different solutions for the temperature as a function of time. You want to regress on those models, use those functional forms.
On the integral: I think my sign is the one that makes
dT/dt = -1/tau T for a 1 d case. (I’ll check. I’m tempted to buy Mathematica for all the stupid signs I’m going to be dealing with!)
The lower limit of the integral should be zero if we specify the ws as the temperature at time t=0. The integrals contribute zero there. The integral prior to t=0 should just get absorbed in the constant.
That equation has nothing to do with Leibnitz rule. It’s conservation of energy for the box and precedes the general solutin (1).
Yes. (1) and (6) have to agree. That’s how you solve for the parameters.
Yes. We want the intergral of 1 to come out. Other than that, I have no idea what you are saying. I’m not normalizing anything and haven’t suggested anyone normalize anything.
I don’t know what a function call you used in R does because I don’t speak R and I’m not finding units in the comments in the previous discussion. (That doesn’t mean they aren’t there– there are just a lot of comments.)
Duplicate of Carrick’s comment on other thread:
===
I’m gonna preface this by saying that Nick did a great job putting together the framework for analyzing this problem, and this isn’t intended as a “gotcha”.
However, I am pretty sure there is an issue with how the version of code that Nick Stokes originally supplied is performing the exponential smoothing.
What he did is appropriate for periodic or probably even statistically stationary input. However, this data is far from it, because there is essentially no driving at or before t=1880 and very large driving as we approach t=2010. For this case, it is very hard for frequency-domain convolutions such as that implemented by R to not generate significant end effects.
Here is a comparison of doing it the way I consider “correct” and Nick’s original method:
2-Box Model Fit Results
The red line is from Nicks’ version and the blue is mine. The “bop-down” around 1880 is real, that is what the forcings really do. Here is a comparison of the GISS forcings (green), my filtered forcings (blue) and the R-convoluted filter (red)
GISS Forcings and 30-year filter
Note that the two filters agree very well after around 1930 (e.g. 30 years after the start of the time series). Before that there is a significant divergence between the two curves, which is reflected in the above 2-box reconstruction results.
About the only major difference other than the difference in the exponential filter is I’m fitting to the monthly data, and Nick is using the annual version. I don’t have 1-month forcing data, so I interpolated the forcing curves for that.
The exponential smoothing is really just the numerical solution to T'(t) + T(t)/tau = F(t).
After people have had a chance to comment on this if we sign off this is a better approach, I’ll post my code and my data files (I’ll provide it to anybody who wants in any case, just trying not to spam this website).
Physically the basic problem with the convolution method (and I think this is “well known”=a few people know it) is that for non-stationary problems it does not give an accurate solution of the original ODE over the entire domain of the solution space.
DeWitt Payne is correct on units. The source temperature data is GISS annual (J-D) land+ocean anomaly (divided by 100 to get degrees C), so one data point per year there. The source on forcings is the sum of all the GISS forcings from the table Tamino linked to, provided in units of W/m^2, with one data point per year for forcings as well.
The fit is to model temperature (degrees C) as a linear combination of different smoothings of the forcings plus an offset. The linear coefficients then naturally have units (as in Tamino’s original post, and repeated inconsistently by the rest of us) of degrees C per W/m^2. The offset of course is just a temperature shift. It would probably be better to just call that offset T_s0 rather than the odd “w_s” notation you have adopted here – that’s what it means, it’s units are temperature.
The offset (T_s0 or your w_s) does have the physical meaning of the model’s expected temperature anomaly value when all forcings are set to 0.
By the way, I did provide all three parameters in comment #18610 on the other thread (w1, w2, and “intercept”) – together with the standard error and other statistical properties of the fit from the R summary.
Thanks Arthur. Sorry I didn’t see the extra parameter. That was one crowded comment what will all the other surrounding stuff.
Yes. The w_s is temperature as I said. I’m leaving it as w_s to distinguish it from the unknowns. That makes it easier for me to not mess up my algebra.
Arthur–
That’s what I thought must be the case. The only difficulty is I don’t ‘speak R’ so I wanted to be sure.
Carrick – not sure what you’re arguing for here – Nick used the “open” flag in his call to ‘convolve’ in R, which zero-pads the input forcing for times earlier than 1880 (before the start of the data). The “circular” flag would have made the assumption that the input data wraps around, but that was not used.
Convolution’s very straightforward though, I suppose one could just write a function to do the sum directly. Is that what you did? I’m surprised it’s so different… perhaps an issue with the normalization for the early period?
Lucia,
OK, I’m no longer blogging in ‘jamas, and my head is clearer. You’re right about the sign in the exponent – I’d forgotten the – sign in the homogeneous equation.
The lower limit of integration is interesting, and may relate to Carrick’s comment, which I haven’t got on top of yet. Of course, any fixed lower limit gives a solution of the inhomogeneous ode; choosing the limit determines the boundary condition. If your s=0 corresponds to 1880, then it’s the same as mine, with my zero padding pre 1880. There is a conceptual difference; your w_s is thought of as the 1880 estimated temp anomaly, which is probably why you attached more importance to it than Arthur, who saw it as an offset to the anomaly base.
There would also be a practical difference if someone else used different padding.
My comment about normalising relates to your query about units. Of course your (1) is correct – normalising only scales the w± coefs. But if normalised, the smoother (my expsmooth) maps 1 to 1, but also 1 quatloo to 1 quatloo, so the smoothing with convolution is unit-neutral. The w± can be interpreted as sensitivities.
You clearly have a role in mind for eq (2), which I’ll be interested to see – my comment was just that my R model was restricted to eq (1).
Arthur – I think this is not so much a question of the convolution, but of the windowing (after all, the forcing prior to 1880 is not “really” zero). I am confident that Nick understands this process and this is just a minor oversight on what looks like nice work.
Nick–
On the sign– I have the advantage of having at least doodled that down a few times.
Yes– I have a role in mind for (2); I need both equations to solve for all the parameters. Once I have them, it’s possible to solve for T2.
I’m pondering the convolution and also Arthur’s interpretation of the constants. It’s true that if your taking the integrals back to the big bang, the constant is the value of the temperature when the forcing is zero. But, in that case, the constants should work out to zero. So, any non-zero magnitude should just be chalked up to measurement uncertainty, finite padding, statistical imprecision etc.
In my representation, those constants need not be zero because the climate might not have been in perfect equilibrium even in 1880.
I do now understand the mechanics of what you did better.
err.
R loaded.
gistemp.txt loaded ( annual to 2003)
gissforc.txt loaded ( annual to 2003)
Nick script loaded.
Executed error free.
no output
hmm should I just RTFM
Nick. might help if you post your txt files AS USED.
lucia,
Which is why I think there is an implied third box in the model as programmed in the R script and possibly by Tamino as well. The third box contains the thermometer. As you so eloquently point out , you can’t do a linear combination of forcings to obtain a temperature with a true two box model. As far as I can tell, there’s no energy exchange between the boxes in the script either, which also is not consistent with a true two box model.
I’d also like to see the response surface for the fit over a range of t1 from 0.5 to 10 and t2 from 5 to 100 using something like the F statistic for the fit at each point. From testing a few points manually while optimizing, it appears there’s a big ridge that gives highly significant fits over a wide range of time constants. For example, the optimum I found was t1=1.1 and t2=22 with an F statistic of 290.3 (2 and 120 DF) and an adjusted R2 of 0.8258. But t1=10 and t2=100 has an F statistic of 280.5 and an adjusted R2 of 0.8208. Obviously the weights change a lot. In the first case, most of the weight is on t2 while in the second case the weights are about equal.
steven mosher (Comment#18664) August 26th, 2009 at 4:41 pm ,
Edit out the gaps, gaps with intermediate headers and headers in both files. Then change the line skip number in each read statement to 1 to eliminate the bad data in gistemp in 1880. That worked for me.
Lucia,
About the only major difference other than the difference in the exponential filter is I’m fitting to the monthly data, and Nick is using the annual version. I don’t have 1-month forcing data, so I interpolated the forcing curves for that.
Dang it. you beat me to it. One problem is that you should probably not interpolate the volcanic forcings. You basically would see a spike in the forcing (eruption) following by a decay as ash ect settled out. Now the monhtly forcing data exists, I’ve been searching but came up empty.
Steve,
there are some minor text issues, which Arthur summarised in a comment:
in gistemp, you need to rub out intermediate headings, and fix some stars in 1880 (just 0’s OK).
In gissforc, there’s a line of —– at the bottom which needs rubbing out. (or better, add nlines=124 to the scan call, as done with gistemp)
R is not too intrusive with error messages, but they are there.
The data files are a bit long to post here. I think Arthur may have them at his site. If it’s a problem, I can email, or maybe post at the CA BB.
I’ll be away for few hours – when I get back I’ll post a revised code with some fixes and more informative variable names.
lucia,
Should Cs be Co in equation 2?
DeWitt Payne (Comment#18667)
I did worse than that dewitt. I preprocessed the gisstemp and forcing files so that they were just vectors of annual numbers, summing the forcings and eliminating the monthly temps all together. ( for my first pass at getting an R program to work)
Then ( stupid me) I loaded nicks script and presto! no errors, blank screen. I failed to see that in the first steps of his script he is parsing the gisstemp and summing the forcing files. Anyways,
I’m late and have to run. Be back at it tonight.
R.. crap one more language to learn. Guess I put my python self teaching on hold
Arthur:
I’m thinking the normalization is the issue here, which is something that is tied intrinsically into the convolution method.
I implemented it as a solution to the ODE in the time domain, which basically has the form (Mathematica notation):
T_n = N_n * Sum[ F_(n-k) x^k, {k,0,n}]
= N_n [F_n + F_(n-1) x + F_(n-2) x^2 + …
… + F_1 x^(n-1) + F_0 x^n]
where T_n = T(t_n) is temperature, F_n = F(t_n), t_n = n/fs, and where x = exp(-1/(fs*tau)), fs is the sampling rate and n = 0, 1, … , nmax.
This of course is the discretized version of the Green’s Function solution of the ODE T'(t) + T(t)/tau = F(t)
If you break this out by value of n, you see the pattern emerge:
T_0 = N_0 F_0
T_1 = N_1 ( F_0 * x + F_1)
T_2 = N_2 ( F_0 * x^2 + F_1 * x + F_2)
…
T_n = N_n (F_0 * x^n + F_1 * x^(n-1) + … + F_n )
The issue here is that the solution from t_0 to t_n shouldn’t depend on the total duration of the forcing available to be integrated for times t > t_n.
In other words, if the values from T_0 tof T_n you are obtaining are changing as the total number of points available nmax (nmax > n) changes, there’s a problem with your solution.
And that is easy to show in this case for the convolution method. I’ve gone from smoothing nmax=100,200,400,600 and 1550 points to illustrate the problem:
Effect of varying data size on expsmooth
I think this problem is intrinsic to convolutions in general, whether they are circular, padded or reflected, or some other end-point treatment is applied.
Anyway, that said what is the “proper” normalization?
What I am using is something similar to Nick’s normalization:
N_n = 1/Sum[x^k, {k,0,n}] = [1 + x + x^2 + … + x^n]^(-1)
This has the property of making T_n independent of the choice of nmax and E[T_n] = E[F_n] in the limit that F-> constant (E[…] is the expectation operator of course).
Since the filter is low-pass, we would like the value at f=0 (constant F) to be preserved by the normalization choice. This does that and also fixes the problem with the filtered value at T_n unphysically varying depending on the choice of nmax.
I also think this is the appropriate choice here, because we want the coefficients we are fitting to to have units of sensitivity.
Comments?
In cases the changes don’t take, some errata:
T_n = N_n (F_0 * x^n + F_1 * x^(n-1) + … + F_n )
… smoothing nmax=100,200,400,600 and 1550 points …
E[…] is the expectation operator
What I like about Lucia’s presentations are that they provide helpful background and a “lead you by the hand” approach that shows the assumptions and work leading up to each step.
That and she doesn’t assume you are an “ignorati” if you don’t understand it right away.
Carrick,
The standard exponential smooth (S(0)=X(0), S(t)=alpha*X(t)+(1-alpha)*S(t-1) is not sensitive to the number of data points smoothed. Your smoothing looks like it overweights the early points compared to later points. There’s also no need to zero pad anything. Just because you can do something more elegantly doesn’t necessarily make it better.
Carrick,
Here’s a test for your smoother. Create two 124 point data files. In the first one set the value of zero for the first 5 points and 1 thereafter. In the second put 50 zeros first and then 1 thereafter. Run your smoother on both files and tell me if they have the same values at points 6 on for the file and points 51 on for the second. I’ll be very surprised if they do.
Wolfgang Flamme (Comment#18648) August 26th, 2009 at 2:25 pm
“I recognize Tamino’s created a crude physical model reflecting basic properties of GISS E. He feeds in the very forcing data that made GISS E’s 20th century simulation look so nice and obtains a temperature time series response that makes his model look almost as nice as GISS E. Now isn’t that amazing?”
No so amazing; GIGO.
“Wonder who measured black carbon forcing and sulphur aerosol forcing in the stratosphere 100 years ago .. wonder why land use forcing abruptly stopped around 1990 after 100 years of steady decline. Wonder why so many forcing and trends came to a sudden stall around the millenium.”
Perhaps the folks at GISS in their last lives made the last century measurements of sulfate and black carbon aerosols.
I personally wonder why the GISS forcing data ends in 2003 instead of 2008 or 2009? Might it have to do with not having any volcanic forcings to offset greenhouse gas forcings?
I took the time to look over the GISS forcings curve through 2003. The volcanic forcings are consistently high throughout. In addition, it appears that between 25% and 30% of radiative forcing from greenhouse gases is assumed to have been ‘canceled’ by non-volcanic (human generated) aerosols in the GISS data. I understand now how Tamino/Arthur can come up with such high diagnosed climate sensitivities from their curve fit. But, really, it is quite a stretch.
Where is the GISS 2003 to 2009 forcing data? How does one explain no ocean heat accumulation since 2003 without volcanoes?
As everyone can probably see from the pingback above, I’ve started a climate blog, and that’s my commentary on this whole mess.
Cheers, and here’s to bringing some actually hypothesis testing into this debate, huh? 🙂
Some have asked about the non-volcano Aerosols forcing.
There may have been some measurements to back some of this up at some point, but for the most part, they are just pulled out of a hat (just quoting a recent statement by Hansen about this) (and he was not kidding).
The total direct and indirect Aerosols effect is -1.8 watts/m2 in 2003 compared to the net overall number of +1.92 watts/m2
There is obviously no annual measurements, just a smoothed function of some sort.
http://img58.imageshack.us/img58/855/modelaerosolsforcingp.png
http://data.giss.nasa.gov/modelforce/RadF.txt
Lucia (#18654)
You’re starting to sound like my good friend Gerhard Kramm there…
What is an average? A linear combination of the values being averaged (with the coefficients normalized so the average of a constant is that constant). GMST is an average of surface temperature measurements from around the world, some of which may be in one of the “boxes” in question, some in the other.
Moreover, the fit here was inspired by two-box solutions, but wasn’t based on any one two-box model (Tamino’s declaration of “atmosphere” vs “ocean” notwithstanding). I am pretty sure a given fit corresponds to an infinite collection of possible physical two-box models.
But your post here is apparently intending to prove otherwise – I just don’t see how you’re getting there yet.
Arthur Smith (Comment#18682)-That’s ridiculous. If the “long time constant” doesn’t correspond to the ocean response but to some part of the surface, then there is no justification for the long response time. I’m not sure that means that some combination of the ocean response and the atmospheric response doesn’t represent the surface temperature, I don’t know enough. But it’s absurd to say that stations measuring temperature are in different boxes! It’s a one dimensional model! And the stations are (supposed) to be measuring the same thing, not different parts of a system with two heat capacities.
Arthur
1) Of course the surface temperautre is a linear combination of temperatures measured on the surface. But do you really think this is the same as a linear combination of temperatures for the deep ocean and the atmosphere?
2) Even if you say yes, it is trivial to extend this to that bizzarre notion. In fact, it solves a difficulty. But, I intend to do the first cut assuming measurements at the surface are all in the same box rather than some surface measurements being in hte “deep ocean” box and others in the “atmosphere”.
3)
The solution method is appplies to all possiblephysical two box models. The regression identifies the most specific two-box model that best explains the data. You do get this, right?
If the two box two-box models that best explains the data does not match any partitioning of the earth’s climate system the method has a problem.
Do you need someone to draw Ven diagrams for you? We can show that it’s possible for one infinite class ofsolutions to not overlap a second infinite class of solutions, and both of those to be embedded in a third infinite class of all possible solutions.
Lucia:
But I don’t think it does identify “the most specific two-box model”, even given that we’re picking specific time-constants to start with. The fit gives you two coefficients for surface temperature in relation to the two time-constants, but that’s it (aside from the third parameter which sets the temperature scale). My understanding is that, beyond temperature scale and the two time-constants, the two-box model requires 3 parameters, of which the fit constrains only two.
I think what you’re trying to show here is that the two-box model (with the 9 parameters you’ve identified) is over-constrained once you have such a fit. So that’s the question, whether the number of corresponding models is overconstrained, exactly one, or an infinite class of models. I believe it’s an infinite class.
DeWitt Palmer:
Urg.
You’re absolutely right of course! … My normalization fails time-translational invariance. Real systems have transients when they get turned on, but as Lucia points out that’s irrelevant: We know the forcings were nonzero before 1880….
That said, the “right” normalization for both Nick Stokes and my versions is to simply take the t->-infinity limit (I think Lucia was kind of hitting on this idea above):
x = exp(-1/(tau*fs))
N_infty = (1 + x + x^2 + x^3 + x^4 + …)^-1 = 1 – x
[However, we need to be careful to evaluate 1-x using expm1(…) because x -> 1.]
Here are three different implementations now. Note my version is now functionally identical to yours, other than some shifting around of notation (I had implemented the summation recursively).
Here’s the codes (all give numerically identical results now):
ns_expsmooth<-function(v,tau,nlen){ # Nick Stokes corrected version
n=1:nlen
x = exp(-1/tau)
e=(n-1)/tau
e=exp(-rev(e))
norm = -expm1(-1/tau);
e=e*norm;
u=convolve(v,e,type="open")[n]
u
}
dwp_expsmooth<-function(v,tau,nlen) { # DeWitt Palmer's suggestion
alpha = exp(-1/tau)
alphap = -expm1(-1/tau)
u = rep(0, nlen)
y = 0;
for (i in 1:nlen) {
y = alpha * y + alphap * v[i]
u[i] = y
}
u
}
clt_expsmooth<-function(v,tau,nlen) { # CLT's solution
x = exp(-1/tau);
u = rep(0, nlen);
norm = -expm1(-1/tau);
sum = 0;
for (i in 1:nlen) {
sum = sum * x + v[i];
u[i] = sum*norm;
}
u
}
Andrew_FL and Lucia – I wasn’t aware that anybody had determined the two boxes had to be “deep ocean” and “atmosphere”. I thought the question was whether there existed a partitioning of Earth’s climate system that could correspond to the fitted two boxes, without constraining from the start that one was “deep ocean” and the other “atmosphere”. In fact, that partitioning doesn’t make sense because it leaves out ocean surface and below-surface land, which ought to be accounted for somewhere. Maybe the best representation is to have part of the ocean surface (say, the polar regions) in the slow box along with deeper ocean, and part of the ocean surface in the fast box. That would naturally lead to GMST being some weighted average of the two box temperatures, rather than corresponding directly to one or the other.
But if it’s a necessary simplification to assume GMST = Ts here, I suppose that’s ok to move forward, I don’t think it should substantively change the math.
Arthur Smith:
Well there’s Tamino:
The purpose of these articles are to critically examine Tamino’s papers. This seems to be an important element.
If we can show that the time constants are unphysical for a two-box model with the interpretations as represented by Tamino, that also really means we can’t use them to estimate sensitivity.
Nor does the T necessarily have anything to do with surface temperature, nor can we use the GISS forcings as inputs to this model without modification, since they apply to atmospheric forcings, not to some combination of ocean and atmopheric dynamics.
I agree with the substance of your comments btw, but I really do think you are just digging a deeper hole for Tamino’s two-box model to fit into.
Steven M (and others interested)
I’ve uploaded a more readable version of the code, plus the data files as used, in a zip at the CA BB site (where I’m registered as Pliny).
This is all too confusing, said Alice to the Queen; I agree with the lass, I don’t know whether I’m counting angels dancing or devils darting among the reeds; one would think that there exists calibrated thermometers located in suitable areas and read by sensible men of sensible means and with the use of an ordinary lead pencil they could make a bit of sense overall of the whole thing.
But here we stand, not sure of where we’ve been or where it is we must go and nobody has even packed the first bag.
Carrick #18686
Just a comment on normalising. I prefer the way I did it because it is more general and foolproof. You want to scale so that the result of applying the operator to 1 is 1? OK, then divide by the result of applying it to 1. Works for anything.
Arthur Smith (Comment#18687)-Well, I’m no expert but my impression is that the ocean broadly responds more slowly to forcings than the atmosphere.
Oh yeah, and the two box model doesn’t have an Arctic. It’s one dimensional.
Nick Stokes (Comment#18689)-That explains a whole lot.
Nick, I showed in a previous comment where your normalization led to problems.
Comments?
Nick, here’s another comparison:
Effect of normalization
The green line is smoothing a segment (100 points @ 12 points/year) using your old normalization, the blue is using the new normalization [replace 1/sum(e) with (-exp1m(-1/d)]. The third is using the new normalization and using the full range of data (1550 points in this case).
Carrick, I see what you mean, but I still prefer my version. I’ve convolved with a finite exponential with 124 terms. You could say it is truncated beyond that. This isn’t an issue for the convolution, because you don’t run out of numbers within the data space.
So the normalisation issue is, if you have an infinite set of 1’s in the future, what operator do you apply. I’ve assumed in effect that it is this truncated set; you’ve assumed that you go on making new terms.
Fortunately the difference will usually be small. The new terms start at exp(-4) with tau=30 and get smaller.
So what’s right? Yours, in theory. But if it matters, then the lack of pre-1880 data will be a much bigger part of the problem. Indeed, this is so – I blithely quoted to Bill Illis above the result of applying his suggested tau=100 yrs, but then the result of the slow decay and no pre-1880 data will make the result very doubtful.
I still prefer the robustness of just dividing bythe result of the application to 1 of the operator you have. True, there is ambiguity, but as I said, if that is a problem, there are worse ones.
You could have a fixed x, but shouldn’t there also be heat transfer between the two boxes proportional to the difference in temperature? We can assume equilibrium between the two boxes before 1880, which is why it’s important to calculate the change in temperature of the ocean box relative to the surface box. Or is it already there and I missed it.
I’ve modified the script to plot from 1880 instead of 1901 and was surprised to see that according to GISS, there was almost no drop in global temperature associated with the Krakatoa eruption. Considering that’s the largest of the volcanic aerosol forcings in the file, it’s somewhat surprising.
Nick, did you even look at the graphs? Clearly the normalization matters; and your method simply isn’t a robust way to normalize.
Also we aren’t “dividing by 1”, we’re multiplying by
norm = 1 — exp(—1/(fs*tau))= —expm1(—1/(fs*tau))
For tau = 30 years, norm is a small number, which is why one has to use —expm1(…) rather than 1—exp(…) to evaluate it.
No I don’t.
Using Lucia’s notation from above, the Green’s function for
T'(t) + T(t)/Ï„ = F(t)
is given by Θ(t—s) exp[—(t—s)/τ].
And the ODE response is just
T(t) = ∫ ds Θ(t—s) exp[—(t—s)/τ] F(s)
After applying Euler’s approximation, the series cuts off at s = t, and you end up with a closed form solution for the normalization (which does not depend on the window size) as I showed above.
Lucia, just a couple of sentences that could use revision for clarity:
“For values of x outside the range 0≤x ≤0, one box will be heated when the other is cooled. (This violates my physical intuition about how global warming would work, but we will defer that discussion until a later time.)”
“Obviously, we can get the two more equations by balacing term multipying the exponentials.”
Additionally, I’ll second Zeke’s comment about considering a Latex plugin of some sort.
DeWitt Payne, there should be a lot of things, but the goal of this programming exercise so far was simply to replicate Tamino’s modeling.
If you have explicit suggestions for what to do next, I’m sure there are people here who would listen.
Carrick,
I did look at your plots, but I must admit to not really understanding them. It looks to me as if you have gone to monthly data, but kept the 124 length for the exponential segment? nlen=124? It should be the full length of the temp data. In my new version, that is made explicit.
The normalisation doesn’t affect the fitted curves. It’s only a multiplier of the sensitivities.
My point about the normalisation I’ve used being the least of problems is this. Going back to annual, the fractional error in normalising is about equal to the first term ignored – exp(-124/tau). That’s fairly small even if tau=30. But if you consider a point at about the middle of the data range, the data at 1880 cuts out when the exponential is about exp(-62/tau). That implies much bigger errors at the smoothing stage.
Anyway, as I say, your version is theoretically slightly better, by that factor (1-exp(-124/tau))
Nick, I went to monthly data but that isn’t relevant to the issue here. I was noticing some issues for large tau, which I think are explained by your normalization choice.
The bottom line is for this particular application, the normalization doesn’t matter much, as long as you are using the full data set and tau isn’t too big (the ratio of tau to duration of the data set, tmax, is what matters really). As tau/tmax gets large, the truncation of your series makes a difference and affects the answers.
For the sake of comparison, what I did was truncate the series at 100, 200, 400, 600 and 1555 points, and apply smoothings for each of these. I’ve modified the calling parameters slightly, so you’d need to adjust these before trying them.
I’m including a code snippet that illustrates what I did more fully, for the case of nmax=100 (first 100 points of the time series for the total forcings), comparing the two normalizations.
You should see something like this.
First source files:
The source file for the normalization test is here:
ns_test.r
which needs this:
giss.forcing.spl.txt
My latest version of Nick’s code is here: stokestd.r
giss.temp.txt
Finally, here’s the actual code snippet (ns_test.r):
expsmooth<-function(v,tau,nlen){ # u is the smoothed v with decay exp(-1/d)
n=1:nlen
e=(n-1)/tau
e=exp(-rev(e))
e=e/sum(e)
u=convolve(v,e,type="open")[n]
u
}
ns_expsmooth<-function(v,tau,nlen){ # Nick Stokes corrected version
n=1:nlen
x = exp(-1/tau)
e=(n-1)/tau
e=exp(-rev(e))
norm = -expm1(-1/tau);
e=e*norm;
u=convolve(v,e,type="open")[n]
u
}
fin <- matrix(scan("giss.forcing.spl.txt", 0), ncol=2, byrow=TRUE)
t = fin[,1]
forcings = fin[,2]
t100 = t[1:100]
f100 = forcings[1:100]
w100 = expsmooth(f100, 12*30, length(f100)) # NS Norm
w100f = ns_expsmooth(f100, 12*30, length(f100)); # fixed norm
wns = ns_expsmooth(forcings, 12*30, length(forcings)); # 1555 pts
plot(t,forcings, type='l', col ='gray')
lines(t, w100, col='red')
lines(t, w100_ns, col='blue')
lines(t, wns, col='blue')
Carrick,
That’s not my reading of lucia’s goal. If that were the case then there’s no point in determining the ocean box temperature at all and we’re already done. Besides, a linear combination of the temperature of the ocean box and the surface box with a fixed energy split ratio, which is what Nick Stokes’ script does, is not a realistic two box model of the planet. Combining temperature anomalies from different boxes to fit the observed surface temperature anomaly is not at all the same as averaging measured surface temperatures from different locations. Those measurements are sampling from one inhomogeneous box, not two.
lucia stated above:
Tamino claims to have made a proper two box energy balance model with different time constants (amounting to different heat capacities) for the boxes. We should take him at his word rather than just cobbling up something that produces similar results and assuming that’s what he did too. We may find in the end that he did indeed combine the temperatures from the two boxes, but until we have a proper model which doesn’t combine temperatures to compare against, we don’t know.
Or to put it more simply, we haven’t proved we’ve replicated Tamino’s modeling yet.
Hi all,
Just a few random thoughts to share after looking through the system of ODEs, the pictures, and also the R code (in rough translation to Matlab).
1. If one filters the RF only the first box and then allow heat into the second (long timescale) box only by diffusion, then the fundamental solution approach (convolution by the heat kernel) has some issues. I think Lucia was pointing this out earlier. If you assume the “atmosphere” is always radiating OLR to space and has a very short adjustment timescale, then why should the heat wait around to filter into a long timescale box over a period of 30 years or whatever? Why doesn’t it just radiate away to space in ~1 year, so that hardly any of it makes it into the 30 year box before it escapes? (It’s just like an electric current flowing through the path of least resistance).
2. Pure diffusive balance won’t work anyway if the 2nd box has no other inputs or outputs. The bottom box has lower temperature than the surface box, which is an untenable equilibrium state for heat diffusion. Thus there must be some other heat sink. The physical process be vertical advection (the overturning circulation) which balances the downward heat diffusion.
3. Applying the forcing to both boxes equally and simultaneously is slightly bizarre. Since net RF is an energy flux (per area), one needs to conserve energy (the First Law!) by dividing the forcing, but how exactly this should be done is a good question. Note: a super thin mixed layer might just complicate the model, since then the solar flux is definitely not absorbed by the upper layer alone.
4. Of course two time constants generates a better fit than one if you are playing the one-side-exponential game. A longish choice of constant creates a phase lag relative to the forcing, in which case you have just created a second basis function with some orthogonal information (eyeballing it, it seems to be around 1/4 phase for the major swings).
5. Yup, we’re just going to have to truncate the exponential tail for the convolutions at some reasonable cutoff. It ought not pose a very serious problem for reasonable timescale filters, since this is really quite rough estimation.
Hopefully, we aren’t going to try inverting the output of an IIR filter that’s a large fraction of the total length of the record!
6. I guess I’m not the only one to notice that some of the forcings are not very orthogonal (pick and choose your difference of large terms!)
Comments welcome, these are just random musings (hopefully in keeping with the spirit of the webpage directory name!).
Oliver
Nick Stokes (Comment#18689)
throws errors. You left header txt in the data files.
/Users/mosher/Desktop/gissforc.txt:1:27: unexpected symbol
1: Global Mean
^
> source(“/Users/mosher/Desktop/gistemp.txt”)
Context is ‘ GLOBAL Land’Error in source(“/Users/mosher/Desktop/gistemp.txt”) :
/Users/mosher/Desktop/gistemp.txt:1:21: unexpected symbol
1: GLOBAL Land
Nick Stokes, here is a comparison of the sensitivities as you hold taus = 0.09 years and vary tau0 from 5 to 100 years. By tau0=100 years, I find a 25% error using the original normalization.
Comparison of Climate Sensitivity, original and revised normalizations
Again, I appreciate your work in setting this up originally. I hope you don’t think I’m picking on you.
If I didn’t think it mattered short/long term, I wouldn’t be pestering you about it.
DeWitt Payne:
Agreed, but that was the goal of the software, or at least my participation in it.
First step is to understand what has gone on before.
Second is to replicate it.
Third (next) is to improve upon it.
Well, up to the SOI part, I think we have replicated it. We can’t demonstrate this because Tamino is a little vague in what he has done, and all we can do is show we get quantitatively similar results.
I would prefer solving the 2-box problem (even if it isn’t really physical), but we lack values for Fo. I’m not sure how you convert SOE into climate forcing units. Perhaps somebody who hasn’t been banned could ask him that.
Again, I’m all eyes, especially if you write down an analytic model. Qualitative ideas are useful, but harder to implement and test.
Oliver,
1. Because the atmosphere radiates towards the ground as well as towards space. MODTRAN says that if you double the CO2 (280 to 560 ppmv) for the 1976 atmosphere with clear sky, the surface sees an additional 3.2 W/m2 while 2.83 W/m2 less flux is radiated to space. And that’s before the atmosphere or the surface warms up.
2. The other heat sinks are at the poles for the ocean. Convective heat transfer is almost certainly more important than diffusion.
3. I think it’s bizarre as well.
Carrick,
I’m working on it.
Carrick,
I ran your ns_test code, and I can see the problem. You’ve truncated the data period to `100 months (about 8 years) and you are smoothing with a 30 year decay period. Yes, the normalising error is then large – as a fraction about exp(-100/360). The result is about 1/(1-exp(100/360) ie about 4 times what it should be.
But the “smoothing” is also very bad, and really meaningless. Almost all the data you need for any reasonable estimate is missing.
For any reasonable application of this theory, you should have data at least twice as long as your longest decay.
Nick:
Um.
As I’ve shown the smoothing is identical to what you get with a longer window as long as you use the proper normalization. (Blue and green lines are identical for domain values of green line.)
So it’s not “very bad”, it works just fine there.
One last point then I’m done kicking this around and ready to write it up as notes and move on:
If you have (fs=sampling rate, tau = time constant):
x = exp(-1/(fs*tau))
N = [1 + x + x^2 + x^nmax]^-1
this can be explicitly written as:
N = (1-x)/(1-x^(nmax+1))
The relative error in using this normalization is then
N/N0 = 1/(1-x^(nmax+1)) ~= 1/(1 – exp(-tmax/tau)) > 1
where tmax is the duration of the data file.
Note this has the effect of inflating the magnitude of filtered forcing associated with tau, which in turn has the effect of reducing the effective climate sensitivity as tau is made larger.
Comparison of Climate Sensitivity, original and revised normalizations
Again, there is absolutely nothing wrong with using long time constants here. The only defect is the wrong use of the normalization. As I’ve shown, my recursive form is a solution to the discretized Green’s Function solution to the diffusion equation.
Increasing the domain of the solution beyond what it already is (e.g., making the time series duration longer) theoretically should have no effect on the values integrated to any time t. It does if you use the wrong normalization. (I rejected my first attempt at normalization too as people may recall when it was pointed out it wasn’t time-translation invariant.)
If you disagree with me, you can put your last word in, but I would recommend anybody who uses your code to correct this error. At the least there is no point in leaving in mistakes, even if you are trying to argue they don’t matter. Without a lot of work you can’t know “they don’t matter.”
Steven,
There is meant to be text there. The scan() commands have a skip parameter which says how many lines (of header) to skip.
I’m surprised you got errors, as I had run the code. R is tricky, though, because it keeps defaults in the workspace, which can make it hard to debug. I’ll try to run it in a new environment just from the zip file.
Carrick:
“I’m not sure how you convert SOE into climate forcing units. ”
I don’t thinl you have to. The fitting process takes care of that. You don’t add the SOI coefficient in to get the sensitivity.
On the normalisation, again I’m not opposed to using your version. I’m really just emphasising that there’s only a discrepancy when the decay is long relative to the forcing data length, when smoothing will be vary bad anyway.
Nick Stokes:
I disagree the smoothing is “very bad” : how can identical results = very bad? But if we agree to fix the normalization, let’s move on.
Carrick #18713. No, the smoothing is very bad. Of course your green is identical to the longer blue stretch; it’s algebraically the same. But exponential smoothing is always bad where the data is short relative to the decay. The true result would be heavily dependent on what happened before 1880. That is why Tamino showed results after 1900
OK I understand now.
In my parlance, what you mean isn’t smoothing per se (that works fine as I pointed out), but the accuracy of the output that comes out of the filter as a result of the missing data prior to 1880. There is no similar problem near the 2004 endpoint because of the one-sided nature of the filter of course…
I’ll see if I can quantify this. I use exponential filters enough that getting details like the normalization perfect and things like that is useful.
Arthur,
Where would you be IF Lucia banned you here?
You see what I am saying?
Arthur Smith (Comment#18410) August 24th, 2009 at 9:59 am
Lucia, I am disappointed in your attitude here. It may get lots of cheers from Watts and friends, but you know what, science is a cooperative endeavor, not mud-wrestling.
Lucia has magically changed some how from this ” science is a cooperative endeavor, not mud-wrestling”
to now?
Lucia, I do not think “THEY” see this as you and I and A.W. and Jeff ID. and many others!
Lack of growth, skin,I.Q.? I am not getting it as of now 2 years of AGW, 30+ years of watching life, they can do anything they want to, but US? only if “THEY” say it is O.K.
Thank you for this thread/blog.
Tim
Nick Stokes (Comment#18714)
Re downloaded. ran from zip.
errors.
success.
Steven, if you want to run my version (basically an extension of Nick’s) the links are in Comment#18704 above.
SOI – I’ve created a version which fits SOI as a third variable. It matches Tamino’s plot visually. The sensitivity is 0.657, a bit less that Tamino’s 0.69. I’m not sure over what range of years T does the fit – I used the full range 1880-2003.
There’s a zip file at the CA BB site. It includes a jpg file of the result.
Unfortunately, as Steve found, the previous zip file had a bug in reading the gistemp.txt file. Unfortunately, I intended to paste the nlines=124 arg from the gistemp read to the gissforc.txt read, but I cut instead of copied. I’ve updated to a new twoboxR2.zip, which is the same except for this fix.
I’ve also updated the new versions to include Carrick’s amendment for the more accurate normalisation. It does make some difference – in the SOI case, the sensitivity increased from 0.657 to 0.666 C m2/W
Arthur–
Math must not be your strong suit. Based on what I’ve done, I can solve 10 equations for 10 unknowns. Initially, I thought I might have 1 too many equations (which, as you know, can resolve itself because it can always turn out something is not linearly independent.) However, I can resolve that by further constraining the relation between the two forcings OR…. by using the rather ridiculous notion of surface temperature being a linear combination of ocean and atmospheric temperature. Or the some equations might be linear combos.
This is solveable. I can get the parameters. I will.
Arthur–
The math assumes there are two boxes. When we get the solution, we can figure out whether we can identify a partitioning that corresponds to any partitioning of the planet earth (or even any planet in the universe. If any of the parameters that must be positive are negative, that solution cannot apply to anything.)
Actually, it’s not. I just did that for this post. I an assume GMST = y Ts + (1-y) Ts and solve the set also.
lucia,
I’ve thought of a possible way to make a linear combination of the two box temperatures rational: they’re land and ocean rather than ocean and atmosphere. Temperatures over land react much faster to a change in conditions than the ocean so it makes sense that the time constant would be shorter. That would also fix the energy split to 30:70 fast:slow.
Dewitt–
Yes. If some of the thermometer are in “the ocean” and some in “the atmosphere”, we can justify the linear combination. (Of course, if we knew which were which, we could create an “ocean” series and an “atmosphere” series and do the fit using two seperate time series. The mechanics might be different from Tamino or Nick’s, but it could be done. )
Still, I can easily say create a temperature that is
T=(y)Ts + (1-y)To, say the parameters are fit to that, do the algebra, get 10 equations with 10 unknowns, and then construct the corresponding Ts and To, This is totally do-able. In fact, it resolves the problem of the extra equation.
We can then discuss whether the value of “y” makes sense. It should not, for example fall outside 0≤ y ≤ 1. Right?
lucia,
I can’t see how a linear combination of temperatures to construct an average could possibly have a negative weight, which would be the case if y was greater than 1 or less than zero, and be physical. That sort of thing was one of the reasons Jeff Id and Ryan O rejected Steig’s Antarctic temperature reconstruction.
Ok… I’m trying to confirm I’ve sorted out the units.
This comment: Arthur Smith (Comment#18610) August 25th, 2009 at 7:19 pm
Arthur provided the following:
Intercept = -0.072 C (Which happens to be statistically significant.)
W1 = 0.5549 m^2* K /Watt or units (m^2 K s/Joules ) for τ1 0.05 years
W2 = 0.046996 m^2 K/ Watt for τ1 20 years.
(I’m assuming I correctly sorted the long and short time constants. If not, this will affect the following calculations.)
my lower case w’s relate to his by
w= W/τ and so have units (m^2 K s/joules) / (s) = m^2 K /joules.
So, converting τ into seconds and dividing, I get
w1 = 3.60812E-07 m^2 K /joules
and
w2 = 7.63954E-11 m^2 K /joules
Now, while I have not solved for all 10 unknowns, I can, in principle, still suggest the Cs should be “reasonable” for the atmosphere. (Otherwise, if the solution for 10 unknowns results in an unreasonable solution, we must laugh.) If so, I can estimate the fraction of heat, x, Taminos solution shoves into the atmosphere and how much it shoves into the ocean.
If we believe the top box is the atmosphere and we think GMST corresponds to the top box, we find the heat capacity of hte top box is roughly Cs1.80E+07 Joules/(m^2K)
x = (w1+w2)Cs = 6.50E+00
So, for each unit heat added to earth’s climate roughly 6.5 units are added to the atmosphere “box”. This means 5.5 units are removed from the ocean “box” (by some mystery process. Perhaps CO2 causes the rocks in the ocean to absorb energy.)
Before we proceed to any interpretation, or any further analysis:
Have I interpreted the units correctly?
Lucia (#18730) – yes, I believe you have the units correct, however w1 is the long-time-constant series and w2 the short one (see the expsmooth calls in the R code) so you have the two reversed on that. Translating to your units, that means the coefficients for the GMST fit to the two forcing time-series become:
w1 = 8.79×10^-10 m^2 K /joules (long time-constant)
w2 = 2.98×10^-8 m^2 K/joules (short time-constant)
Generally I ignore the traditional flurry of insults here, but was it really necessary to go with “Math must not be your strong suit” in comment #18725? You really are sounding more and more like Dr. Kramm. Not a good sign!
I believe my mathematical intuition is somewhat on the high side though I know many professional mathematicians and mathematical physicists who have far more expertise than I do. The question of over-determination or under-determination (essentially the issue of dimensionality of the space in question) is one of the very central things that intuition is useful for – and my intuition feels pretty strong on this case. But forge ahead, and we’ll see.
Arthur–
As you know, I ignore most of your insults.
As for the answer to your question: Yes. It was necessary to say that. I suspect that you can understand the equations in the post. You are here in comments, as I assume you do understand, it appears you must not wish to read the equations, the argument or the discussion in the post.
In the post I have 10 equations and 9 unknowns. As discussed in comments, this can easily be resolved by either
a) using your notion that GMST is a linear combination of surface and ocean temperatures. This will introduce 1 more unknown, “y” resulting in 10 equations and 10 unknowns, or
b) insisting that the GMST is one of the two box temperatures and suggesting that the forcing applied to the two boxes is xGT+b. or
c) it could resolve itself if it turns out that the 10 equations are not linearly independent. If so we specify more physical parameters.
In the post, I have shown that the solution for the three parameters does not appear underconstrained. But if it is… so?
Arthur–
For what it’s worth, if the solution is overconstrained, that would be a big problem for the two-box model. BIG. I do not think the problem is over constrained, I have never suggested that, and it is not my intention to show it.
However, if someone did show this, it would mean the two-box model is nuts. I am confident no one will show it is nuts in this way.
Also, if you are going to get riled about someone suggesting your math ain’t all that, you should think about your own demeanor. Consider dropping your expressions of “disappointment”, complaints about which “side” my blog is on, what other analysis might be more “useful” for me to pursue. All are judgmental, insulting and the final one would, in this specific instance, mean that for some reason, you think we all understand that questioning whether or not Tamino’s box model corresponds to physical reality is somehow not “useful”. Obviously, we have a difference of opinion on that.
I’m generally not going to call you on each one of these because I know that you tend to resort to this when your technical argument is poor or non-existent. The fact that I do not call you on them does not mean I do not notice the insults. I find it better suits me to focus my response on the poverty of your technical argument and not diffuse it with compliants about your veiled insults to me, my blog visitor or your other various and sundry problems.
Wow, Arthur really opens himself to charges of hypocrisy. In addition to what Lucia quotes above, here are some examples of Arthur’s “politeness” to the host:
– “science is not mud-wrestling”
– “your assertions otherwise indicate either some deep lack of understanding of what he’s talking about, or, well, something else”
– “Maybe you’ll learn something from being banned by Tamino, but I don’t see much sign of it in this thread yet”
And on his own blog, until he corrected his error, he accused Lucia of “groping”.
And then he cries foul when Lucia responds in kind. Nice going, Arthur.
lucia (Comment#18733)
Best. Post. Ever.
Ok. So then for the assumptions above,
x=0.56 which is not implausible on its face. We would need to check whether this makes sense for an atmosphere and ocean– but at least it’s not 6.
See… checking for plausibility does not always hurt. My point is, and has always been, that this sort of checking must be done.
Lucia – claiming that the parameters correspond to no physical solution is by definition saying that the space of physical corresponding two-box models is empty, ie. the problem is over-constrained (once you include all the equations and inequalities relevant for physical reality). That seemed to be your implication and goal here. Maybe I’m wrong on that.
On the number of equations – you have only written out 2 of them (eq. 7 and 8 above). You’ve asserted that there are 10 independent ones, but they haven’t been written out yet, so it’s hard to tell. Write ’em up, and we’ll see. Unless you’re expecting the rest of us to do all the derivations from this point on?
Arthur
First: The question is whether the parameters correspond to a box model that can represent the earth’s climate system. The solution is not overconstrained, but may correspond to krypton.
Moreover, even if I had made the claim embedded in your system, your response is nonsense.
The solution could correspond to something that violates the 2nd law of thermodynamics. Assuming that you mean “overconstrained” in a mathematical sense, that result would not imply the problem is “overconstrained”. It would simply mean that the solution to the not-over constrained mathematical problem did not correspond to a physical model for a system governed by the 2nd law of thermodynamics.
If, while discussing a the solution to a math problem, you meant ‘over constrained’ in some obscure non-mathematical sense, you should define what you mean by that word. (Or find a better word that says what you mean.)
There are 10. The 2 are examples.
If you understand that this is an eignevalue problem and you understand how those two are obtained, you can easily discover the additional 10. I even explain how they are obtained:
1) 4 are applied as discussed in bullets 1-4. (I showed how you get 2)
2) Four more are applied by doing the corresponding solution for conservation of energy for the other box.
3) Two are obtained from the equation for the eigenvectors.
2+4+4 = 10.
I’m not writing out all the equations in this post. I am deferring that. But if you have a general understanding of the mathematics of the eigenvalue problem, you understand that I have 10 equations. If you don’t have a general understanding, you don’t.
Those who understand the math don’t need me to write them out, and those who don’t understand the math won’t benefit from the specific equations at this point.
Now, I plan to buy Mathematica and solve these puppies.
Arthur Smith:
When you solve the inverse problem for a physical problem, one thing you should always do is verify the solution is physical. E.g., no negative time constants, no negative densities, conservation of energy and so forth.
I would say Lucia’s goal is to test this, not necessarily to just “falsify” it.
It appears you climate types are unaware that testing for physical constraints is part of normal-operating procedure with inverse solutions.
Am I wrong?
lucia (#18736) – x=0.56 implies you are making an assumption that C_s = 1.8×10^7 J/K m^2 (assuming you calculated that based on your equation 8a above) – where’d that value for C_s come from?
Lucia:
And even more to the point, if they still don’t understand the motivation of what you are doing giving them the equations wouldn’t fix that.
Perhaps Arthur Smith could explain how having all 10 equations will fix his lack of understanding of pedagogical issues.
Arthur–
The post on the 2nd law of thermodynamics which is linked, and motivated using smaller time constants for the “short time constant”.
It’s normal in engineering. From the beginning, I’ve said it needs to be tested. I didn’t say my intuition said it would fail. I don’t know if it will fail. It might pass.
When I asked Tamino if he tested, I wanted to know if he tested. He said he has. Possibly, Arthur, who seems to be representing himself knowing the thoughts in Tamino’s mind (rather than simply reading what Tamino said) should explain how Tamino tested.
Every physical constraint resolves to a mathematical equation (first-law-type) or inequality (second-law-type). You have used some of these physical constraints already in formulating these equations in the first place – after all the two-box model is a physics-based model. You have mentioned some of the inequality-type constraints as well (x and y have to be between 0 and 1).
The way I think of problems like this is in the space of parameter values – here you have 9 or 10 parameters, so a 9-or-10-dimensional space. Every independent equation slices a dimension off that space, defining a surface of one lower dimension. Every inequality chops the space into a region that satisfies it and one that doesn’t (the equality version of the inequality defines a space of one dimension lower again).
In the end the question is what is the dimension of the space remaining after you’ve applied all the appropriate constraints from equations and inequalities. If the remaining space is completely empty, the problem is over-determined – no solutions. If it’s just one point, then there’s a unique solution. If the remaining space has a dimension of 1 or more, the problem is underdetermined and there are an infinite collection of possible solutions. Which is what I pretty strongly believe to be the case here.
That intuition is based on long experience of going through such mathematical exercises during my years as a math/physics/theoretical physics researcher: there’s not enough incoming information here to sufficiently constrain the problem as far as I can see. What you will most likely find is that, of the 10 proposed equations, some will be simple rearrangements of one or more of the others, and you don’t have enough to find a unique solution. If I get a bit of time I’ll even have a go at working through it and see what the problem is myself. I don’t think you need Mathematica for this one (or that it will help much).
The solution domain will not be empty.
If it were, that would mean the two-box model corresponds to nothing and would present a problem for Tamino.
I’ve already outlined that there is, given a rather straightforward assumption on the distribution of forces. That’s discussed between equations (3) and (4). There are other possibilities, and I would welcome discussion of them because, they are unphysical.
The purpose of Mathematica is simply to avoid algebra errors. That’s it.
I’m just representing what I understood from reading Tamino’s “brilliant, and elegant” post and some of the commentary. What I understood when Tamino said he “tested” was that he checked that the coefficients (W1 and W2 above) were both positive – since there are some choices of time constants for which you get an unphysical negative response. I can’t see what else there was to test.
Arthur–
Evidently you can’t.
It seems to me what Tamino has done (in the end) is write down the equations:
Eq. (1): T(t) = a0 + c_S F_S(t) + c_L F_L(t)
where F_S(t) and F_L(t) are exponentially filtered versions of the total forcings with constants tau_S and tau_L, and he’s shown this does a good job of fitting to the GISS time series.
Now it’s the case he’s tried to justify these particular processes based on a two-box model. But that’s just an interpretation of a mathematical model. I’m not as good at mind reading Tamino’s brain (I lack Arthur Smith’s slith abilities for that), but I would think that Tamino started with the idea of a fast and a slow process and tried drumming up a physical model to help justify it.
I happen to think there’s reason to suspicious about the physicality of the two-box model, so it is not implausible that Lucia’s exercise is going to show violations of physical constraints.
That doesn’t mean that the two-process model is flawed, just that the interpretation given was wrong.
This sort of thing happens. In physics, when people started playing with second quantization, they discovered their probability distributions could be negative, and this was a major stumbling block until they realized their interpretation was wrong: The quantities weren’t probability distributions but rather creation and annihilation operators (the value was positive if there were more particles than antiparticles and negative if more antiparticles than particles).
As a phenomenological equation, Eq. (1) encapsulates some importance aspects of variability in the atmosphere. I suspect one could eschew dispersion-type equations and go with any ODE that had a low-pass character and get similar results.
I’ve promised some code today for my group so I don’t have more time to play, but I may another obviously non-physical filter as an alternative to the exponential filter this evening and see how far I get.
Lucia (#18742)
in the original post you had C_s as 1.34×10^7 J/m^2 K which I believe gives x = 0.41, not 0.56.
However, C_s is one of the 10 parameters you’re trying to determine here, right? So we don’t know that’s the right value for x until you’ve solved the whole thing.
Arthur–
Yes, C_s is one of the parameters I would get from the equations. So, no, I can’t know it until I solve the 10 equations:
This is why I wrote
If the solution for the 10 unknowns results in a value of Cs that is not “reasonable” for the earth’s climate, that would indicate that Tamino’s solution is not supported by a two box model that corresponds to the earth.
But, for the time being, at least it looks like if the solution for the 10 parameters does give back a “reasonable” value for Cs, then the value of “x” will be reasonable.
This point is in favor of the two box model mapping into physical reality. So, even though you don’t believe I can solve the equations, you should be happy. 🙂
As Carrick suggest, even if my solution for the 10 parameters ends up not mapping into physical reality, it might be possible to find an entirely different physical model that supports Tamino’s solution. If that is found and verified to work, then you can attach physical significance to the regression parameters in your fit.
Carrick is contemplating other models that might support Tamino’s general notions about the curve fit.
But even now, it’s possible that the solution for the 10 eqn will show that Tamino’s model is just fine. I never suggested that it must be wrong. I have only said it must be checked.
Arthur–
Why do you think one negative constant would be unphysical?
Arthur Smith:
I thought I detected an odor of applied mathematics. Not meant as a slight here, but I guessed from some of your statements, you didn’t have experimental physics training.
Tell us why you think a negative climate sensitivity for a forcing that includes anthropogenic CO2 emissions would be physical?
Negative greenhouse gas effect?
Sorry, I missed that section of your comment (#18730) where you also got a result x = 6.5, when I read it the first time; I must have seen the error in the first section and responded based on that without reading the rest.
I’m still not sure why C_s is now 1.80 when it was 1.34 the other day?
Lucia:
I’m not certain it’s impossible, but I wouldn’t ordinarily expect part of a system to cool down when heat flux increases. Something like that was once proposed to happen to parts of Europe if the Atlantic thermohaline circulation shut down thanks to warming, so I suppose there are models where it could be possible, but it just seems very unlikely.
Arthur–
It should be 1.34, I probably just read quickly and mistyped in my spreadsheet. The radiative heat transfer coefficient is 1.8 e108 and I misread the digit. Cs should be the same in both the comment and the post. (I don’t proof read comments as carefully as entire post. Happens.). The result “that’s reasonable” still holds.
Carrick–
The sensitivity is the sum of the two weightings. I want to know why Arthur things both weights must be positive.
If 100% of the forcing were added to the lower box, and the system is at equilibrium before we start forcing, one weighting could be positive and one negative for the upper box. This is required because the dT/dtime in the upper box ==0 until after the lower box warms. However, d^2T/dt^2 is positive– becuase it’s temperature picks up after the lower box warms.
If both time constants are positive, only way to get these two conditions is for one weight to be positive and one negative. (I’m suspect the weight on the term containing the positive time constants should be the positive one– but I’m not even sure about that. I know the sum should be positive.)
If we don’t know what’s in the boxes, or what GMST corresponds to, how do we know both weightings are positive? The sum should be– but why each individually?
Arthur–
Fair enough.
But.. I don’t think we can say that result violates physics without having some notion what’s in the boxes or what linear combination of temperature corresponds to GMST. (That is, if we buy into the notion that GMST is a linear combination of temperatures.)
So, I don’t think that constitutes “a check”. Positive sensitivity, on the other hand, would seem to be a requirement.
Arthur:
Generally I ignore the traditional flurry of insults here, but was it really necessary to go with “Math must not be your strong suit†in comment #18725? You really are sounding more and more like Dr. Kramm. Not a good sign!
Given your earlier comments, what of this? Are you the only person qualified to throw turds, or will there be a retraction?
Ok, on the significance of the intercept value and T_s0, T_o0
First, the assumption T = y T_s + (1 – y) T_o isn’t quite general enough – you could have an offset as well – the observed (fitted) T really should be
T = c + y T_s + (1 – y) T_o
which gives you another parameter (c) to fit.
However, looking at the two-box model in the previous post on this, the equations require that, if the forcing term F is a constant value of 0, then both T_s and T_o must be constant at zero as well. I.e. The zero-anomaly for the temperatures in the two-box model as presented is exactly zero.
So Lucia’s w_s and T_o_0 must be zero as well. This removes two parameters from the list of parameters to be fit – and it also makes eq. 7 and eq. 8 identical, so we’ve lost 1 equation for each derivative term (2 equations total).
I.e. I think that leaves us with 9 variables and 8 equations.
We call it “egg throwing” Terry. Arthur Smith would have fewer enemies on this blog if he cut out the passive aggressive nonsense.
I’ve rethought the issues and agree you could have one of the parameters negative as long as the sum is positive. It has to do with how the low-pass filtering distributes the forcings differently between F_S and F_L (mainly due to the low-pass character of F_CO2 and F_sulfates, F_L is mostly just measuring the anthropogenic effects).
Anybody feel like backing out the climate sensitivity to CO2? I.e., what does this model fit predict for delta T for a doubling of CO2?
Passive aggressive? I’m just trying to hold people to a high standard of truth. As far as I can tell the discussion has moved considerably closer to real scientific questions than the “stupid Tamino, snort snort, Lucia’s wonderful” level it was at the other day.
Carrick – doubling CO2 adds (probably a little under) 4 W/m^2 forcing, so just multiply the sensitivity numbers here by 4 to roughly get the temperature change for doubling CO2. I.e. a result of 0.6 C m^2/W is 2.4 C per doubling.
Arthur–
In the appropriate baseline, this is true. However, it’s not necessarily true in the baseline for your measurements. 🙂
Arthur–
Also, for the general solution using the correct method, those don’t have to be zero even in the ‘proper’ baseline. They are initial values. The world does not have to be at equilibrium in 1880 even if Nicks solution method assumes this is so.
I’m just trying to hold people to a high standard of truth.
George Boole is spinning in his grave at that comment…
Arthur, whether you like it or not, you’ve been exhibiting passive aggressive behavior in your commenting style at times. I’m willing to leave it at that. Do we really need to explore it further?
That’s what I’m trying to do too.-
Thanks, that makes perfect sense.
Holy canoli! I’ve got some reading to do! This conversation just goes on and on.
I’m probably putting my foot in my mouth again, but my interpretation of w1 and w2, which are the box temperatures, is that w1=x*S*F(tau1, t) and w2=(1-x)*S*F(tau2, t) (if I have my nomenclature correct) where S, the climate sensitivity, is w1+w2. If x is greater than 1 then heat must flow from box 2 to box 1 to conserve energy if FT(t greater than zero) is positive which means that F(tau1,t) and F(tau2,t) are also positive so F(tau1,t) is greater than F(tau2,t) when tau2 is greater than tau1 and both are positive . But that can’t happen because that would mean that energy would somehow flow from box 2 to box 1 as the temperature of box 1 increases and the temperature of box 2 decreases ((1-x)is less than zero). That looks like a Second Law violation to me. Or am I, as usual, missing something here.
Ok, bear with me a bit here, I’m going to try reformulating Lucia’s notation and get down to the actual list of equations in question.
First of all, start by setting values for τ+ and τ- which satisfy eq. 4 from the previous post:
Eq. 1: -1/τ(+,-) = (1/2) (-b +- sqrt(b^2 – 4c))
where
b = (αo + αs) + (γo + γs)
c= (αo αs) + (αs γo +αo γs)
(substitute those b and c here, we’ll re-use c).
So Eq. 1 represents 2 equations relating the 4 α and γ parameters (and our predetermined τs).
Now define
F+(t) = ∫_-infinity^t exp[- (t-s)/τ+] F(s)/τ+ ds
and F-(t) similarly (F(s) is total forcing – F_T in Lucia’s notes here).
Let T(t) be the fit to the observational temperature with coefficients a1, a2, and a3:
Eq. 2: T(t) = a1 + a2 F+(t) + a3 F-(t)
This fit of T(t) to the forcing with two exponential smoothings indicates an underlying physical two-box model; denote the two components by (s) and (o) (tentatively “surface” vs “ocean”) that respectively have temperatures (Ts and To), heat capacities (Cs and Co) and a heat exchange term (β). Then
Eq. 3: T(t) = c + y Ts(t) + (1-y) To(t)
introducing two new parameters, c and y. Ts and To are similarly linear combinations of F+ and F-, but with a zero constant term:
Ts(t) = w+s F+(t) + w-s F-(t)
To(t) = w+o F+(t) + w-o F-(t)
which defines 4 more parameters in the respective weights.
From these definitions and Eq. 2 and 3 we get
Eq. 4: a1 + a2 F+(t) + a3 F-(t) = c + (y w+s + (1-y) w+o) F+(t) + (y w-s + (1-y) w-o) F-(t)
or
Eq. 5: a1 = c
(1 equation fixing 1 parameter)
Eq. 6: a2 = y w+s + (1-y) w+o
Eq. 7: a3 = y w-s + (1-y) w-o
(2 more equations relating the 5 parameters y, and the w’s to the fitting coefficients).
Now we have the differential equations that Ts and To satisfy, dependent on the split (x, 1-x) in forcing. First note the derivative of F(+,-) above is simple:
dF+(t)/dt = -F+(t)/τ+ + F(t)
and similarly for F-(t). Then
Eq. 8: dTs/dt = – w+s/τ+ F+(t) – w-s/τ- F-(t) + (w+s + w-s) F(t)
and similarly for To. The original two-box ODE’s were:
dTs/dt = -(αs + γs) Ts(t) + γs To(t) + x F(t)/Cs
dTo/dt = γo Ts(t) – (αo + γo) To(t) + (1-x) F(t)/Co
so substituting Eq. 8 and the original equations for Ts and To and equating the terms mutliplying F+(t), F-(t) and F(t), respectively, we find:
Eq. 9: -w+s/τ+ = – (αs + γs) w+s + γs w+o
Eq. 10: -w-s/τ- = -(αs + γs)w-s + γs w-o
Eq. 11: x = Cs (w+s + w-s)
(note Eq. 11 is the same as Lucia’s Eq. 7 and 8 above)
Eq. 12: -w+o/τ+ = γo w+s – (αo + γo) w+o
Eq. 13: -w-o/τ- = γo w-s – (αo + γo) w-o
Eq. 14: 1-x = Co (w+o + w-o)
which gives us 6 equations relating the 4 w’s, 2 α’s, 2 γ’s Cs, Co, and x (11 parameters) and our pre-determined τ’s.
The γ’s are actually related to one another through:
Eq. 15: Cs γs = Co γo
(from the definition of γ in terms of the heat exchange rate β and heat capacity) and note the constraint that each γ >= 0
Additional constraints are 0 < x < 1 and 0 < y < 1, and that the Cs and Co should be reasonable for some physical partitioning of Earth’s climate system.
I think that’s it. So by my count here we have 13 parameters (4 w’s, 2 α’s, 2 γ’s, c, x, y, Cs, Co) and 12 equations (2 from Eq. 1, Eq. 5-7, Eq. 9-14, 15) and at least 3 independent inequalities (for x, y, and γ – in principle the constraints on Cs and Co could also be expressed as inequalities). I am not entirely convinced the 12 equations are independent (in particular I think there’s some interdependency between Eq. 1 and 9,10,12,13 but haven’t had a chance to confirm that).
Other questions are:
* does the solution space for any input parameters τ(+,-), a1, a2, a3 include a nonempty subspace satisfying the first 3 inequalities?
* do the solutions from fitting the observed temperature record include a nonempty subspace where Cs and Co are reasonable physical values for the Earth?
Sound about right?
Oh yes, of course, my Eq. 9 and 12 are just a version of the eigenvalue equation for 1/τ+, and eq. 10 and 13 the same for 1/τ- – that is there are only 4 independent equations between the 2 of Eq. 1 and Eq 9,10,12,13, rather than 6. So I think we’re left with 10 equations in 13 unknowns, or a 3-dimensional solution space. It still could be empty for some input parameters once you add in the inequalities though, so I suppose worth checking on that too.
Dewitt–
There is a requirement between w1 and w2. Since the comments/post contain both large and small case ws.
Arthur–
On a notational note: I realize the blog doesn’t have latex installed and my equations are difficult to read. That said, it would at least help if you used html to create subscripts and superscripts and spaces. When you define new variables and write out things like, w+s I don’t know if that’s supposed to be w + s or w+s.
I don’t know what you are trying to do, but if I create a temperature series based on a liner combination of T=(y)Ts + (1-y) To I get 10 linearly independent equations and 10 unknowns. Generally speaking, it appears you are creating 4 extra equations by creating a 3rd ODE based on the summed variable, and creating a new solution based on the summed variable and then more or less discussing the solution method I outlined.
This will not result in equations that are linearly independent. You also seem to be defining some extra new variables for the solutions in the process. This can all be simplified– but it’s better to think of it more simply in the first place.
On specific details,
1) there is no reason to zero out the constant term in Ts or To.
2) You are failing to account partition the time varying forcing into an average plus anomaly portion. This is screwing up your solution.So, you are going astray sometime before your equation 5. Since what follows depends on that, it’s unlikely to be correct.
What a great thread. Let’s see would you see such a thread ( and such restraint from the host) on…
1. RC? Nope
2. Eli Rabbit? Nope
3. Open Mind? Nope
4. In it for the gold? Nope.
5. Atmoz, probably. ( always liked that kid)
In the end if lucia finishes Tamino’s work for him and proves that the model doesnt violate any physical laws, where will we be?
Basically, I think you have a simplified effects model for climate change. Note, the forcings all depend on prior modeling so Tamino circular reasoning ( this is not computer models) is still circular. But what you do have ( which is what I liked about LUMPY ) is a simple effects model for doing quick and dirty estimates of various future scenarios. Now if we really want to have some fun with this model ( using SOI is a cheat) we should go get our hands on the forcing files for the paleo GCM runs. I’ve seen bits and pieces of them around
“doubling CO2 adds (probably a little under) 4 W/m^2 forcing, so just multiply the sensitivity numbers here by 4 to roughly get the temperature change for doubling CO2. I.e. a result of 0.6 C m^2/W is 2.4 C per doubling.”
According to the IPCC, doubling CO2 (everything else equal) produces a forcing of 5.35 * Ln(2) = 3.71 watts/M^2, or 2.22 degrees per doubling at equilibrium for a sensitivity of 0.6. With a slow box lag of 20 years, this equals 1.41C increase 20 years after an instantaneous doubling, and 1.92C increase 40years after an instantaneous doubling.
It occurred to me that what we are really trying to do here is Kalman filtering. Or at least the system identification part of it. Looking at sec 2 of that link, we have a state variable – a vector of box temperatures, and a control vector (the forcing), and an observation model (just surface temp).
If you scroll down in that link, the Kalman-Bucy filter describes the differential equation version, which is similar to what Lucia has been describing.
Does the model include solar fluctuations?
http://www.scienceblog.com/cms/small-fluctuations-solar-activity-large-influence-climate-24510.html#content
Lucia – whenever I have used subs or sups commenting in your blog before, they always disappear in the final display of the comment. Let me try again though – in my handwritten notation, the w+s above etc. are w+s – superscripts referring to the “box”, subscript referring to the ‘+’ or ‘-‘ eigenstate.
Ok, that didn’t work. Please explain how to get html subscripts and superscripts in comments here!
On the need for the reference box temperatures To0, Ts0 to be identically zero – this is seen from the original ODE’s. The forcing F(t) can be any function of time. Let F(t) = 0 for all t. Then the solutions resolve to Ts(t) = Ts0, To(t) = To0 (i.e. constant). However, the ODE’s (eq. 1 and 2 of previous post) reduce to:
Eq. 1: 0 = -(Cs/τs) Ts0 + β (Ts0 – To0)
Eq. 2: 0 = -(Co/τo) To0 – β (Ts0 – To0)
Solving for Ts0 we get
0 = To0 * (β – Co/τo – β^2/(β – Cs/τs))
so either To0 = 0 (and Ts0 = 0) or there’s a singular relationship between the heat capacities and heat exchange term β
WordPress filters out <sub> and <sup> by default from comments (along with <img> and host of other html commands).
I believe it’s possible to change this behavior, but it’s not something I’ve ever done myself.
I found a WordPress code support page http://support.wordpress.com/code/ that suggested using [sourcecode] tags, but no luck.
[sourcecode language=’css’]
w+s
[/sourcecode]
I just hit the tip jar in appreciation of all of this great work.
Diogenes could have retired if he had found Lucia.
x
Ok, I put up a post with proper subscripts and superscripts and a tiny bit more explanation here.
Arthur Smith (Comment#18743)
re: “dimension of 1 or more, the problem is underdetermined and there are an infinite collection of possible solutions.”
As I recall, there is an infinite number of solutions for EACH undetermined dimension.
brid (Comment#18734) August 27th, 2009 at 9:50 am
nice 1
18733!
Thank you
Joel (Comment#18781)
I found her first on CA
Nick Stokes (Comment#18771)
Yup. We used Kalmen filtering to do all sorts of interesting things like tracking targets and missiles by using vehicle dynamics equations.
Arthur said: “I can tell the discussion has moved considerably closer to real scientific questions than the “stupid Tamino, snort snort, Lucia’s wonderful†level it was at the other day.”
Indeed, Lucia is wonderful. Can you imagine this thread going on at Tamino’s, where one can be ridiculed and banned for a simple and real scientific question about parameters of a model?
Lucia,
I think Kalman filtering really does do what is needed to fit a full two-box model, but it is complicated. It would go like this:
In their terms, you have an equation
x_k = F x_{k-1} + B u_{k-1} + w
x is the vector of temperatures (T_s, T_o), u is a vector of forcings (probably two terms both equal to GISS forcing, but you could split them). In your terms, F is I + A, where A is your transfer matrix. B might be 2×2 diagonal. You’ll have to start with an estimate for these, and also the variance of w.
Then there is an observation process. You can probably assume no noise, and H=(1,0), so Z = H x = T_s
So then there is a rather heavy but explicit recurrence relation, their Predict and Update, which lets you step forward from 1880, getting an optimal estimate for x, and hence T_o. That will yield the optimised R2 for T_s.
Then you have to revisit the estimates for A, B and variance of w, with a nonlinear Newton process for fitting these numbers (repeating the recurrence) to get the best value of R2.
Arthur-
Oddly… no. Remember those are anomaly equations. That equation is true when you subtract off the steady state solution for a particular applied forcing. Non zero answers only mean you didn’t magically happen to know the correct baseline for both F and T that corresponds to the “true” anomalies.
I thought what you thought two days ago. Carrick sorted me out.
On the subscripts– I’ll look for a plugin. Comments plugins only suck CPU when you submit. Display ones suck CPU every time the page displays.
Arthur–
Ahh!!! But that discussion helped!
A baseline exists in which both To and Ts are zero. So, we can set those to zero provided we remember the solutions are in a different baseline as the data! Your substitution is fine. It will make equations of the type (7) trivially true in the “right” baseline. So, we are down to 8 equations and 8 unknowns.
Thanks.
Lucia:
This is what I would say too. To=Ts by definition in equilibrium of course.
Can we calculate the change in ocean heat content from the model fit at this point? Would it be the integral over time of the difference between the smoothed forcing Fo and the total forcing times the splitting factor?
I’ve posted to show the solution for the parameter space can be expressed as 9 equations with 9 unknowns.
I’ll be closing comments here, please comment at
The third post in the series.
In answer here:
Dewitt:
I can’t because I haven’t solved the 9 equations with 9 unknowns.
Currently, it looks non-linear, but I may be able to sort that out.
Generally– There are lots of good comments above, but I’m afraid I can’t keep up with all of them. I have various projects (blog and real life) and blog-wise, I am focused on finding the range of parameters for the two box problem.
I think it’s great that other people are discussing other possible models that might be better (and I’m reading them). But even though I think some may be more promising than the two-well-mixed boxes, I want to finish thinking about that problem before jumping off into other (quite likely better) models.