OMS-NetZero / FAIR

Finite-amplitude Impulse Response simple climate model
https://docs.fairmodel.net
Apache License 2.0
123 stars 62 forks source link

Running idealised scenarios for non-CO2 emissions #28

Open michellecain opened 6 years ago

michellecain commented 6 years ago

Hello,

I've just started using FAIR in Oxford and I've run into an issue that Richard and I haven't been abkle to figure out.

I want to run some idealised emissions pulses using FAIR. I succeeded when emitting CO2 only. I'm using the multigas version. This is roughly what I did:

Import RCP emissions (because I was running that before and had it set up). Set the emissions all to zero except for one pulse of CO2. Then I run the following:

C,F,T = fair_scm(emissions=emissions,useMultigas=True,scale=scale,efficacy=scale)

where 'scale' is a list of all zeroes except for the first entry, which is 1, because I want to use CO2 only.

This seems to work, in that the F and T look as expected.

However, when I try and do the same for methane (ie I set a nonzero pulse of methane, and change scale so the second entry is 1 and the rest zero) then I get negative forcing and cooling that does not reflect the emission pulse.

I've inserted a set of plots from this methane pulse run. In the concentration and forcing plots, methane is in yellow. I have no idea why it goes down so much, as this is simply a result of one giant pulse that I add at the start. Perhaps it is due to some assumptions in the model that I don't know about. If you have any advice on how to fix this, or to do a pulse emission of non-CO2 gas in a better way, then please let me know.

Thanks Michelle

ch4idealised

chrisroadmap commented 6 years ago

Hi Michelle,

I set up the multi-gas version of FAIR to give realistic results when spun up from pre-industrial and historical emissions under the RCP scenarios, so it may give strange results for a pulse emission of methane.

The first thing to note is that natural, time varying emissions of CH4 and N2O are switched on by default in multigas mode and that present-day and pre-industrial emissions are not the same. This may be why you see a negative RF for zero emissions. Secondly there is a pre-industrial concentration of non-CO2 gases provided which is passed to the radiative forcing calculation. And thirdly under both the default Etminan forcing calculation and the alternative option (Myhre) the methane forcing depends also on concentrations of N2O. There may be other things at play which are also affected by methane emissions (tropospheric ozone and oxidation to CO2). You may have already done this but you may be able to find a documentation branch on the main GitHub which gives some examples (unfortunately I don't think it's exhaustive).

I'm currently travelling but if you can hold on until after Easter I'll provide the combination of options you need for a pulse emission of methane (or make code changes to enable it).

Thanks, Chris

On Thu, 29 Mar 2018, 13:19 michellecain, notifications@github.com wrote:

Hello,

I've just started using FAIR in Oxford and I've run into an issue that Richard and I haven't been abkle to figure out.

I want to run some idealised emissions pulses using FAIR. I succeeded when emitting CO2 only. I'm using the multigas version. This is roughly what I did:

Import RCP emissions (because I was running that before and had it set up). Set the emissions all to zero except for one pulse of CO2. Then I run the following:

C,F,T = fair_scm(emissions=emissions,useMultigas=True,scale=scale,efficacy=scale)

where 'scale' is a list of all zeroes except for the first entry, which is 1, because I want to use CO2 only.

This seems to work, in that the F and T look as expected.

However, when I try and do the same for methane (ie I set a nonzero pulse of methane, and change scale so the second entry is 1 and the rest zero) then I get negative forcing and cooling that does not reflect the emission pulse.

I've inserted a set of plots from this methane pulse run. In the concentration and forcing plots, methane is in yellow. I have no idea why it goes down so much, as this is simply a result of one giant pulse that I add at the start. Perhaps it is due to some assumptions in the model that I don't know about. If you have any advice on how to fix this, or to do a pulse emission of non-CO2 gas in a better way, then please let me know.

Thanks Michelle

[image: ch4idealised] https://user-images.githubusercontent.com/23382921/38088348-6a38c9be-3353-11e8-80cb-9a23e13cb274.png

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/OMS-NetZero/FAIR/issues/28, or mute the thread https://github.com/notifications/unsubscribe-auth/AGu-nFywimL4eAWOthSCVojRFHg0vd-_ks5tjNFqgaJpZM4TAPMS .

michellecain commented 6 years ago

That explains a lot, thank you! It would be great if you could suggest the best setup for running an idealised pulse of methane (or any gas, as I'd like to look at other gases too) when you are back. Thanks! Michelle

chrisroadmap commented 6 years ago

Hi Michelle,

This should work for methane. I had to tune the pre-industrial concentrations of CH4 and N2O ever so slightly as it appears that the pre-industrial concentrations and the first year of natural emissions aren't quite in balance (I think its because one or the other is actually defined as a mid-year value so is 6 months out of sync in the time series).

import numpy as np
import fair

emissions = np.zeros((300, 40))         # number of years times 40
emissions[:,0] = np.arange(1850,2150)   # column 0 is year. It shouldn't make a
    # difference here but it's safer to use one from 1850 onwards.
emissions[0,3] = 1000.                  # column 3 is methane emissions in MtCH4.
C_pi = np.zeros(31)
C_pi[0:3] = [278, 721.89409166, 272.96017722] # set pre-industrial concs of
    # gases other than CO2, CH4 and N2O to zero
# all forcing scaling factors zero except methane
scale = np.zeros(13)
scale[1] = 1.

C,F,T = fair.forward.fair_scm(emissions,
    natural=np.array([209.2492,  11.1555]),
    scale=scale,
    F_volcanic=0,                       # turn volcanic forcing off
    F_solar=0,                          # turn solar forcing off
    C_pi=C_pi,                          # use the pre-industrial concs
    )
michellecain commented 6 years ago

Thanks, this seems to work for methane for me. Is it OK to use it in a similar way for CO2, N2O and other gases, by changing the column in the emissions to match the intended species, and change the index of scale to match?

chrisroadmap commented 6 years ago

Hi Michelle,

N2O and minor gases should work fine in a similar way. Forcing index number 2 is N2O, and number 3 is a catch-all for minor gases, whereas individual concentrations of the GHGs are available. Hopefully this should be semi-clear from the docs...

CO2 should also work. In multi-gas mode the Etminan forcing routine is used which takes band overlaps with N2O into account. You can also run in CO2-only mode (set useMultigas=False, use just a one-column array of CO2 emissions for emissions, and do not supply an array of natural emissions if you do this). This will give a slightly different result for CO2 forcing however, as this forcing relationship uses the Myhre et al. (1998) logarithmic dependence on CO2 concentrations.

Hope this helps! Chris

michellecain commented 6 years ago

Hi Chris

Thanks so much for your earlier help. I have a new but related query which is hopefully straightforward to answer.

I wanted to test varying the lifetime of methane. As the lifetime is 9.3 in my input file, I thought that this following code should be identical to running without specifying the "lifetimes":

lt=np.zeros(31)
lt[1]=9.3
C_lt1,F_lt1,T_lt1 = fair.forward.fair_scm(emissions,
        natural=np.array([209.2492,  11.1555]),
        scale=scale,
        F_volcanic=0,                       # turn volcanic forcing off
        F_solar=0,                          # turn solar forcing off
        C_pi=C_pi,                          # use the pre-industrial concs
        lifetimes=lt,
    )

However, the curves look very slightly different:

lifetime9_3

And then when I vary this methane lifetime, I get some odd looking behaviour too:

lifetimes

Am I going about thing wrong? Is there a different way I should be using to vary the lifetimes. I thought this was easier than editing the lifetime input file and re-running.

thanks

Michelle

chrisroadmap commented 6 years ago

Hi Michelle,

The problem is that the non-CH4 lifetimes are set to zero in your example. This will affect the radiative forcing for N2O due to its codependence with methane, and thus temperature. Here's how to fix:

# Here's the original MWE
import numpy as np
import fair

emissions = np.zeros((300, 40))         # number of years times 40
emissions[:,0] = np.arange(1850,2150)   # column 0 is year. It shouldn't make a
    # difference here but it's safer to use one from 1850 onwards.
emissions[0,3] = 1000.                  # column 3 is methane emissions in MtCH4.
C_pi = np.zeros(31)
C_pi[0:3] = [278, 721.89409166, 272.96017722] # set pre-industrial concs of
    # gases other than CO2, CH4 and N2O to zero
# all forcing scaling factors zero except methane
scale = np.zeros(13)
scale[1] = 1.

C,F,T = fair.forward.fair_scm(emissions,
    natural=np.array([209.2492,  11.1555]),
    scale=scale,
    F_volcanic=0,                       # turn volcanic forcing off
    F_solar=0,                          # turn solar forcing off
    C_pi=C_pi,                          # use the pre-industrial concs
    )
# Here's your addition. Note the warnings...
lt=np.zeros(31)
lt[1]=9.3
C_lt1,F_lt1,T_lt1 = fair.forward.fair_scm(emissions,
        natural=np.array([209.2492,  11.1555]),
        scale=scale,
        F_volcanic=0,                       # turn volcanic forcing off
        F_solar=0,                          # turn solar forcing off
        C_pi=C_pi,                          # use the pre-industrial concs
        lifetimes=lt,
    )
C:\Users\mencsm\FAIR\fair\forward.py:378: RuntimeWarning: divide by zero encountered in double_scalars
  C[t,2] = C[t-1,2] - C[t-1,2]*(1.0 - np.exp(-1.0/lifetimes[2])) + (
C:\Users\mencsm\FAIR\fair\forward.py:384: RuntimeWarning: divide by zero encountered in true_divide
  lifetimes[3:]))) + (0.5 * (
# Plot methane concentrations - all seems to be in order, at this resolution...
import matplotlib.pyplot as pl
pl.plot(C[:,1])
pl.plot(C_lt1[:,1])
pl.ylabel('methane concentration ppb')
pl.xlabel('year');

image

# plot temperature outputs - here's the difference
pl.plot(T)
pl.plot(T_lt1)
pl.ylabel('temperature anomaly')
pl.xlabel('year')
Text(0.5,0,'year')

output_3_1

# but in your example, the lifetimes of the non-methane gases are an array of zeros.
# this overrides the default lifetimes from FAIR.

# returns 31 gases as a list in order
import fair.constants.lifetime as lifetimes
lt2 = lifetimes.aslist
lt2
[1000.0,
 9.3,
 121.0,
 50000.0,
 10000.0,
 3100.0,
 222.0,
 5.2,
 16.1,
 28.2,
 13.4,
 47.1,
 38.9,
 7.7,
 3200.0,
 45.0,
 100.0,
 85.0,
 190.0,
 1020.0,
 26.0,
 5.0,
 11.9,
 9.2,
 17.2,
 16.0,
 2.9,
 65.0,
 20.0,
 0.8,
 1.0]
# Now set the methane lifetime of this new array to 9.3
lt2[1] = 9.3
C_lt2,F_lt2,T_lt2 = fair.forward.fair_scm(emissions,
        natural=np.array([209.2492,  11.1555]),
        scale=scale,
        F_volcanic=0,                       # turn volcanic forcing off
        F_solar=0,                          # turn solar forcing off
        C_pi=C_pi,                          # use the pre-industrial concs
        lifetimes=lt2,                      # use the imported lifetimes + methane
    )

# plot temperature - they should be the same
pl.plot(T)
pl.plot(T_lt2)
pl.ylabel('temperature anomaly')
pl.xlabel('year')
Text(0.5,0,'year')

output_5_1

# Now experiment with different lifetimes for CH4 but keep other GHGs the same.
lt3 = lifetimes.aslist
for lt_ch4 in [9.3, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0]:
    lt3[1] = lt_ch4
    C_lt3,F_lt3,T_lt3 = fair.forward.fair_scm(emissions,
        natural=np.array([209.2492,  11.1555]),
        scale=scale,
        F_volcanic=0,                       # turn volcanic forcing off
        F_solar=0,                          # turn solar forcing off
        C_pi=C_pi,                          # use the pre-industrial concs
        lifetimes=lt3,                      # use the imported lifetimes + methane
    )

    # plot temperature - they should be the same
    pl.plot(T_lt3, label='lifetime %s' % lt_ch4)
pl.ylabel('temperature anomaly')
pl.xlabel('year')
pl.legend()
<matplotlib.legend.Legend at 0xfc18898>

image

Now, if you use a different lifetime to 9.3 years, you will see that the temperature continues to rise or fall after the initial perturbation, and in fact the signal of the perturbation is swamped. This is because the natural emissions are no longer in a steady state with the lifetime and the methane concentrations approach a new equilibrium that differs from 721 ppb. The easiest fix for this is to vary the natural emissions by trial and error, and run a long-ish (1000 years may be sufficient for methane depending on your lifetime) integration, analysing the concentrations of methane at the end of the run and seeing if it (a) equilibriates and (b) equals the pre-industrial value. Maybe I need to code this up as a module - I did it offline to get the initial natural emissions of CH4 and N2O based on the default lifetimes.