CS of 70

[kdorheim]

Problem ?? The current version of Hector solves with a CS of 70 with. However Hector doeclim does not. It is unclear to me if this is an actual bug or not but with a Tgav value of only 13 in 2100 seems suspicious.

So running...

library(hector)
core2 <- newcore(system.file('input/hector_rcp85.ini', package = 'hector'))
setvar(core = core2, dates = NA, var = ECS(), values = 70, unit = 'degC' )
reset(core2)
run(core2, runtodate = 2100)
fetchvars(core2, dates = 2100, vars = c(GLOBAL_TEMP(), RF_TOTAL() ))
shutdown(core2)

Returns ....

> library(hector)
> core2 <- newcore(system.file('input/hector_rcp85.ini', package = 'hector'))
> setvar(core = core2, dates = NA, var = ECS(), values = 70, unit = 'degC' )
> reset(core2)
Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini> run(core2, runtodate = 2100)
Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   2100
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini
> fetchvars(core2, dates = 2100, vars = c(GLOBAL_TEMP(), RF_TOTAL() ))
             scenario year variable     value units
1 unnamed hector core 2100     Tgav 12.993035  degC
2 unnamed hector core 2100     Ftot  9.267068  W/m2
> shutdown(core2)
Hector core (INACTIVE)

With that kind of CS and RF in 2100 then we would expect the 2100 temp before ocean heat uptake to be like 174 degC right?

Given that temp change before ocean heat uptake would be CS/3.7 * RF =degC 

So if Hector really can solve with a CS of 70 are were confident about the 12.9 deg C? Or am I doing something wrong with the R hector code?

But when I run an old version of Hector, Hector doeclim on pic /pic/projects/GCAM/Dorheim/CMS/hector

it throws the expected error

[dorh012@constance01 hector]$ ./source/hector ./input/hector_rcp85_CStest.ini
Thu Feb  7 13:36:19 2019:NOTICE:printLogHeader: hector version 2.0
Thu Feb  7 13:36:19 2019:NOTICE:main: hector wrapper start
Thu Feb  7 13:36:19 2019:DEBUG:h_reader: h_reader created for ./input/hector_rcp85_CStest.ini
Thu Feb  7 13:36:19 2019:NOTICE:parse: Parsing ./input/hector_rcp85_CStest.ini
Thu Feb  7 13:36:19 2019:NOTICE:main: Creating and initializing the core.
Thu Feb  7 13:36:19 2019:WARNING:printLogHeader: hector version 2.0
Thu Feb  7 13:36:19 2019:NOTICE:main: Setting data in the core.
Thu Feb  7 13:36:19 2019:WARNING:setData: Disabling onelineocean
Thu Feb  7 13:36:20 2019:NOTICE:main: Adding visitors to the core.
Thu Feb  7 13:36:20 2019:DEBUG:addVisitor: Core adding a visitor
Thu Feb  7 13:36:20 2019:DEBUG:addVisitor: Core adding a visitor
Thu Feb  7 13:36:20 2019:NOTICE:main: Calling prepareToRun()
Thu Feb  7 13:36:20 2019:WARNING:prepareToRun: Disabling onelineocean
Thu Feb  7 13:36:20 2019:NOTICE:prepareToRun: Computing dependencies and re-ordering components...
Thu Feb  7 13:36:20 2019:NOTICE:prepareToRun: Preparing to run...
Thu Feb  7 13:36:20 2019:NOTICE:prepareToRun: Spinning up model...
Thu Feb  7 13:36:21 2019:NOTICE:run_spinup: Carbon model is spun up after 352 steps
Thu Feb  7 13:36:21 2019:NOTICE:prepareToRun: Model spun up after 352 steps
Thu Feb  7 13:36:21 2019:NOTICE:main: Running the core.
Thu Feb  7 13:36:21 2019:NOTICE:run: Running...
* Program exception:
msg:    Assertion failed: bad Tk value
func:   ocean_csys_run
file:   ocean_csys.cpp
ffile:  ocean_csys.cpp

line:   130

[rplzzz]

@kdorheim The difference that you are seeing between the R version of hector and the standalone version is that in the R version you are only running up to 2100, and in the standalone version you are running to 2300. The exception in ocean_csys occurs sometime after 2100.

As to the question of whether we should believe the results we are seeing with large climate sensitivity, I am skeptical. Last night I was able to get hector v 1.1 up and running again (harder than it sounds, due to copious bit rot in our old versions), and with a large climate sensitivity like 70, the code fails with that "bad Tk value" error very quickly, around 1877. Here is what I was able to capture from the logs:

Hector v 1.1, S = 3

Thu Feb  7 15:47:06 2019:DEBUG:run:  internal_Ftot=0.0330914 tgav=0.0209389 degC in 1877 
Thu Feb  7 15:47:06 2019:DEBUG:log_state: ----- State of LL box -----
Thu Feb  7 15:47:06 2019:DEBUG:log_state:    carbon = 783.16 Pg C -> 783.183 Pg C
Thu Feb  7 15:47:06 2019:DEBUG:log_state:    T=22.0196 degC, surfacebox=1, active_chemistry=1
Thu Feb  7 15:47:06 2019:DEBUG:log_state: ----- State of HL box -----
Thu Feb  7 15:47:06 2019:DEBUG:log_state:    carbon = 140.431 Pg C -> 140.433 Pg C
Thu Feb  7 15:47:06 2019:DEBUG:log_state:    T=2.01964 degC, surfacebox=1, active_chemistry=1

Hector v 1.1, S = 70

Thu Feb  7 15:37:55 2019:DEBUG:run:  internal_Ftot=0.76339 tgav=-10.2338 degC in 1877
Thu Feb  7 15:37:55 2019:DEBUG:log_state: ----- State of LL box -----
Thu Feb  7 15:37:55 2019:DEBUG:log_state:    carbon = 808.803 Pg C -> 809.533 Pg C
Thu Feb  7 15:37:55 2019:DEBUG:log_state:    T=13.0347 degC, surfacebox=1, active_chemistry=1
Thu Feb  7 15:37:55 2019:DEBUG:log_state: ----- State of HL box -----
Thu Feb  7 15:37:55 2019:DEBUG:log_state:    carbon = 144.208 Pg C -> 144.299 Pg C
Thu Feb  7 15:37:55 2019:DEBUG:log_state:    T=-6.96529 degC, surfacebox=1, active_chemistry=1

Hector current version, S = 3

Thu Feb  7 17:51:02 2019:DEBUG:run:  tgav=0.103704 degC in 1877
Thu Feb  7 17:51:02 2019:DEBUG:log_state: ----- State of LL box -----
Thu Feb  7 17:51:02 2019:DEBUG:log_state:    carbon = 782.991 Pg C -> 783.026 Pg C
Thu Feb  7 17:51:02 2019:DEBUG:log_state:    T=22.1048 degC, surfacebox=1, active_chemistry=1
Thu Feb  7 17:51:02 2019:DEBUG:log_state: ----- State of HL box -----
Thu Feb  7 17:51:02 2019:DEBUG:log_state:    carbon = 140.407 Pg C -> 140.41 Pg C
Thu Feb  7 17:51:02 2019:DEBUG:log_state:    T=2.10478 degC, surfacebox=1, active_chemistry=1

Hector current version, S = 70

Thu Feb  7 17:47:18 2019:DEBUG:run:  tgav=0.11967 degC in 1877
Thu Feb  7 17:47:18 2019:DEBUG:log_state: ----- State of LL box -----
Thu Feb  7 17:47:18 2019:DEBUG:log_state:    carbon = 782.936 Pg C -> 782.964 Pg C
Thu Feb  7 17:47:18 2019:DEBUG:log_state:    T=22.1197 degC, surfacebox=1, active_chemistry=1
Thu Feb  7 17:47:18 2019:DEBUG:log_state: ----- State of HL box -----
Thu Feb  7 17:47:18 2019:DEBUG:log_state:    carbon = 140.395 Pg C -> 140.397 Pg C
Thu Feb  7 17:47:18 2019:DEBUG:log_state:    T=2.11605 degC, surfacebox=1, active_chemistry=1

So, basically, at the time that v 1.1 is crashing due to out of bounds temperature with the high S value, v 2.1 is showing scarcely any difference between the S=3 and S=70 runs at all. I get that the version with DOECLIM should be different and in some sense better, but the effect of climate sensitivity in this case seems implausibly small to me.

That's where I left things last night. I'll continue to dig into this and update you with what I find.

[rplzzz]

@JGCRI/hector

Here is how climate sensitivity was used in v 1.1: https://github.com/JGCRI/hector/blob/12db8826e196d7883b0f4c5f6b36de290f422ce8/source/components/temperature_component.cpp#L205-L210 So, loosely speaking, all else being equal the temperature is proportional to S. That's not exactly the whole story, since heatflux is probably affected by temperature in some way, but when you increase S, the only way to have global temperature not go up proportionally is to heat up the oceans.

The use of S in the current version is rather more complicated (btw, hector currently writes "version 2.1.0" in its log files, but we do not have a v2.1.0 tag): https://github.com/JGCRI/hector/blob/fc05ec42a9d050b7d625c90ad633c7299961d843/src/temperature_component.cpp#L208-L214 This is the only place where the climate sensitivity value is used (aside: this is in prepareToRun, so it explains why you have to reset the model to get the changed value of S to take effect). It enters exclusively through the expression q2co2/S. What is q2co2? If you guessed "3.7 W/m2", give yourself a cookie.

So, S affects the constants cfs, cfl, and kls (lines 212, 213, and 214). Notably, in all three cases there is a term that is inversely proportional to S, and it is added to a term that is independent of S. So, in the limit that S goes to infinity, all three of these constants asymptote to constant values. Therefore, unless there is another sneaky use of S that I've missed, under these equations, for very large S, the effect of climate sensitivity eventually saturates.

To figure out exactly when the saturation occurs, we'll have to compute the relative values of those two terms. And, it wouldn't hurt to comb through and double check that S doesn't creep into the equations through some other avenue. I may do those things later this morning. For now, I'm going to leave it here for the rest of you to comment on. A "climate sensitivity" that eventually saturates seems a little suspect to me, but I'm just a mathematician, so this isn't really my area of expertise. I'd like to hear the rest of you weigh in on whether this is reasonable or not.

[rplzzz]

@JGCRI/hector Here are plots of the DOECLIM variables affected by S

As you can see, none of them have any behavior to speak of beyond S of about 25. In fact, as this closeup shows, most of the change happens for S < 10

In order to interpret these curves, I think we will need to do a little more digging into how these parameters are used. However, I think it is now at least clear why very large values of S don't have the effect we expect; it's because for S > 10, changing S doesn't really do anything.

Question: Do any of you have a paper or other source with the equations this code was based on? Before we spend too much more time digging into this, it would be nice to know if the equations were even programmed correctly in the first place.

[rplzzz]

All right, here is my final dispatch on this subject. The equations in question come from Kriegler's dissertation [1] , parts of which I have now read, and let me tell you, it's a barn burner. The relevant equations are in Appendix A, equations (A.19) -- (A.21). I have checked the formulae in the code against the those equations, and everything seems to be in order. Also, the figures in my last comment compare reasonably well to Figure A.6 in the dissertation, so everything appears to be working as designed.

Kriegler doesn't really comment on the saturation of the effect of climate sensitivity on the rest of the model, so it's possible that he didn't consider it remarkable. However, there doesn't seem to be a simple, intuitive explanation of why this should be so. I think one would have to recreate the derivation of equations (A.19) -- (A.21) rather carefully to pinpoint the cause.

The only remaining question is, what does this mean for our interpretation of Hector's "climate sensitivity" parameter. It clearly does not represent "the temperature change from a doubling of CO₂", but I'm wondering if "climate sensitivity" ever really means that, in any but the crudest models. The fact that Kriegler didn't consider this phenomenon particularly interesting suggests that it is well known. So, perhaps the description of climate sensitivity that is normally bandied about is only a crude approximation, valid in certain limiting cases. A final possibility is that if you analyze the equilibrium conditions closely enough, you really do see the expected temperature change, but that the very similar values of cfl/cfs (the quantities that enter into the differential equations for temperature) for large values of S mean that it will take a very long time to reach those equilibria.

In any case, here's the summary:

1. The model is working as designed.
2. Large values of S (greater than about 10) are nearly identical to one another in their effects on the short-term evolution of the system.
3. When we eventually write about anything having to do with the value of "climate sensitivity" in Hector, we need to decide whether our "climate sensitivity" (using DOECLIM) really is the same as the climate sensitivity that everyone else talks about. If it isn't, then we need to make sure we communicate that fact when we discuss it.
4. In the context of our Bayesian calibration, there probably is no problem here. The high values of S aren't causing any particular harm, and the CDF of our prior is over 0.97 by the time we get to S=25. That should be enough to regularize the Monte Carlo sampling, and that is all we really need from our prior.

Diagnosis: Hypochondriasis. There never was anything wrong with the patient. Administer a placebo and put back to work ASAP.

[1] https://publishup.uni-potsdam.de/opus4-ubp/frontdoor/deliver/index/docId/497/file/kriegler.pdf

[bpbond]

Question: Do any of you have a paper or other source with the equations this code was based on?

The idea of climate sensitivity is a lot older than that dissertation, FWIW. I think of the canonical reference as Hansen et al. 1984 (in Climate Processes and Climate Sensitivity, Geophysical Monograph 29, American Geophysical Union, Washington, DC, 1984), pp. 130–163), but the basic idea comes from much earlier feedback studies in electrical engineering.

A "climate sensitivity" that eventually saturates seems a little suspect to me 2. Large values of S (greater than about 10) are nearly identical to one another in their effects on the short-term evolution of the system.

Well, an S of 70 is ridiculously large and improbably unlikely–that is, impossible–in any physical sense. We (and folks who have studied the PDF of S like Roe & Baker 2007) only expect this approximation of what is finally an emergent property of the climate system to be valid for a range of 0-10 degC (see Figure 3 in R&B).

Diagnosis: Hypochondriasis

Haha!

[rplzzz]

The idea of climate sensitivity is a lot older than that dissertation

I'm aware of that. My question was about how climate sensitivity gets translated into coefficients in the differential equations for this particular model. The path from the simple equilibrium condition in the linearized energy balance equation (which is how climate sensitivity is derived) to the coefficients in a completely different set of equations in DOECLIM is long and convoluted, and it's not clear, to me, at least, that S still has the same meaning when you get to the end of that path.

Well, an S of 70 is ridiculously large and improbably unlikely–that is, impossible–in any physical sense.

That depends on what you think S means. If you plug an S of 70 into DOECLIM, you get a physically reasonable result; it's just that that result doesn't behave much differently than if you had plugged in S = 10. (Update: This is only true in the short term. If you run the model out far enough, the differences do become substantial. See further notes below)

We ... only expect this approximation of what is finally an emergent property of the climate system to be valid for a range of 0-10 degC

And that is assuming that the parameter that we are calling "climate sensitivity" in hector represents the same thing as what other people are calling "climate sensitivity". It did in versions prior to us integrating DOECLIM, but it's not clear that it still does.

Also, if people think climate sensitivity has to lie in the range 0-10, then why so much fuss over what the tails look like in our prior distributions? If we think the distribution should be truncated anyhow, then the tails are really unimportant.

[rplzzz]

BTW, @bpbond, thanks for the link to the Roe & Baker paper. It's very interesting.

After going over it, I don't believe your reading of figure 3 is correct. Although that figure cuts off at S of about 10, plugging the parameter values for curve (B) into their equation (3) gives an area under the S > 10 tail of the curve of about 0.8. Indeed, the tail is very long; the area under the S>25 tail is 0.02. Both of these figures are reasonably close to what our log-normal prior produces.

Here's a comparison of the two curves, Roe & Baker solid, and log-normal dashed.

As you can see, the upper tails aren't that different. The main difference is in the lower tail, where the log-normal prior allows much lower values than the R&B function.

Since this is supposed to be a weakly informative prior, it's worth comparing the entropy of the two distributions. It turns out the entropy of the R&B distribution is about 2.2 nats, while the log-normal (with our parameters) is about 2.6 nats. For my money that makes the log-normal the better choice (and I have some other arguments that I don't have time to go into right now), but either would serve admirably. Note also, however, that either distribution will likely have us evaluating some surprisingly high values of S, at least occasionally.

[ssmithClimate]

The conventional definition of climate sensitivity is temperature at equilibrium from CO2 doubling. In a dynamic model, that can formally only be determined by running the model out to Infinitity. Can we run hector 10,000 years to check? (Or at least several 1,000 to see what its asymptoting toward.)

[rplzzz]

I have done some runs along those lines, and I'll post them later this afternoon along with some additional thoughts I had on the subject last night. The short answer is, for low values of the sensitivity we can, but higher values will trigger the temperature out of range error sometime in the 2200s.

[cahartin]

I'm not sure we can run Hector out 1000 years since we aren't modeling the appropriate physics/chemistry needed. In particular at time series that long we need to account for circulation changes and carbonate chemistry. Typically SCMs are only valid in the less than 500 yr range.

[ssmithClimate]

To test climate sensitivity, a fixed CO2 concentration run should be used, so we won’t need to worry about the carbon cycle part. We’d be testing DOECLIM heat transport.

[rplzzz]

Here are a set of runs that go out to year 3500. The CO₂ concentration starts at 286 in 1850, ramps up to twice that value over the next 70 years, where it remains constant for the rest of the run. I zeroed out the emissions of all the other gases, except for the halocarbons, because there are so many of them.

It's actually pretty hard to run much past this point, since things start to go haywire in the carbon accounting. In fact, the S=50 run couldn't get much past 2700.

Here's what we can glean from these results.

1. Although the difference in the coefficients for large values of S is tiny, the differences accumulate over time. As a result, by the time you get out to 2700, the S=25 and S=50 runs are looking rather more different than their coefficient values might suggest.
2. The higher the value of S, the longer it takes to reach equilibrium, which is what you would expect.
3. The equilibrium temperature rise seems to be a little lower than the nominal value of S. For example, the S=3 run only makes it to about 2.7 C, and the S=10 run doesn't look like it's going to make it to 10 C.
4. Therefore, Hector's "climate sensitivity" isn't exactly the climate sensitivity as it's usually defined, but it seems to be reasonably close.
5. Finally, to bring us full circle to the question that started this little odyssey, very large values of S are viable and seem to produce results that make sense, at least in the short run.