FUND Model for FoS September 2017 Newsletter¶

Revised 2017-11-11 to correct for tonnes of CO2 vs tonnes of carbon.¶

In [1]:
using Mimi, Plots, PlotRecipes, StatPlots, Distributions
pyplot() # PyPlot backend is best supported
# Set include to the full path from home directory to the FUND code
include("/Users/iancameron/Documents/FUND Model/fund.jl-master/src/fund.jl")
WARNING: Method definition describe(AbstractArray) in module StatsBase at /Users/iancameron/.julia/v0.5/StatsBase/src/scalarstats.jl:560 overwritten in module DataFrames at /Users/iancameron/.julia/v0.5/DataFrames/src/abstractdataframe/abstractdataframe.jl:407.
WARNING: Method definition count() in module IterTools at deprecated.jl:49 overwritten in module Iterators at deprecated.jl:49.
WARNING: Method definition count(Number) in module IterTools at deprecated.jl:49 overwritten in module Iterators at deprecated.jl:49.
WARNING: Method definition count(Number, Number) in module IterTools at deprecated.jl:49 overwritten in module Iterators at deprecated.jl:49.
Out[1]:
getfund (generic function with 1 method)

Set various model parameters for this series of calculations

In [2]:
modelStartYear = 1950 # Built into the FUND model
modelEndYear = 3000 # Model runs for 1051 years, including modelStartYear
finalYear = 2300 # Final year for calculating SCCO2
yearsToRun = 1051 # Model runs from 1950 to 3000 inclusive
discountRate = 0.03
marginalCO2Year = 2020
numMdYears = finalYear - marginalCO2Year # No. of years to calculate marginal damages
convert1995To2016Dollars = 1.575 # Converts 1995 $US to 2016 $US
# Years of marginal CO2 release for SCCO2 calculations
marginalCO2Years = [2010, 2020, 2030, 2040, 2050]
# Added 2017-11-11: Ratio of molecular mass of CO2 to C
ratioCo2ToC = 3.664
# senarios array is the list of directory names for the 5 FUND scenarios (1x5 array) stored in scenariosPath
# NB: Each directory contains all 101 of FUND's data files, but only ch4em.csv, n2oem.csv, scenacei.csv,
# scenpgrowth and scenypcgrowth differ in each of the 5 directories. This simplifies the code for accessing
# any scenario.
scenarios = ["Base", "Image", "Merge", "Message", "MiniCAM"]
scenariosPath = "/Users/iancameron/Documents/FUND Model/fund.jl-master/scenarios"
# scenarioLabels is 1x5 matrix of labels used in various plots
scenarioLabels = ["Base" "Image" "Merge" "Message" "MiniCAM"]
# regionLegends represents FUND's index of 16 geographical regions
regionLabels = ["USA" "CAN" "WEU" "JPK" "ANZ" "EEU" "FSU" "MDE" "CAM" "LAM" "SAS" "SEA" "CHI" "MAF" "SSA" "SIS"]
Out[2]:
1×16 Array{String,2}:
 "USA"  "CAN"  "WEU"  "JPK"  "ANZ"  …  "SEA"  "CHI"  "MAF"  "SSA"  "SIS"

USA is United States; CAN is Canada; WEU is Western Europe; JPK is Japan-Korea, ANZ is Australia-New Zealand, EEU is Eastern Europe; FSU is Former Soviet Union; MDE is Middle East, CAM is Central America, LAM is Latin America; SAS is South Asia; SEA is South-east Asia, CHI is China; MAF is North Africa; SSA is Sub-Saharan Africa; SIS is Small Island States.

ECS Determination¶

FUND uses a gamma probability distribution for ECS with parameters α=6.47815626 and θ=0.547629469, with a minimum of 1.0. This gives a mode of 3.0, but because the distribution is skewed to the right, the average value is closer to 3.5 (see below).

A 2015 paper https://judithcurry.com/2014/09/24/lewis-and-curry-climate-sensitivity-uncertainty/ by Lewis and Curry, based on long-term observations, derived an ECS of 1.6. This was later reduced to 1.02.

In [3]:
using Distributions, StatPlots
α = 6.47815626
θ = 0.547629469
p = StatPlots.Gamma(α, θ)
# Plot the gamma distribution
StatPlots.plot(p, lw=2, legend=nothing, xlims=(1,9), ylims=(0,0.35),
    title="Fig. 1 - Probability Distribution for FUND's ECS",
    ylabel="Probability", xlabel="ECS, °C")
StatPlots.plot!(xticks=1.0:1.0:9.0) # Specify tick marks for x-axis in place of default values
Out[3]:
In [4]:
png("Fig_1")
In [5]:
# Calculate mean ECS based on above distribution
total = 0.0
for i = 1:100000
    total += max(1.0, rand(Gamma(α, θ)))
end
println("Average ECS = ", total / 100000.0)
Average ECS = 3.5442461262492166
In [6]:
# Use these ECS values 
baseECS = 3.5
newECS = 1.0
Out[6]:
1.0

Marginal Damage Function¶

Function calcMarginalDamages calculates marginal damages for 1 t of CO2 for all 16 regions in each year, returning a matrix extending from emissionYear+1 to modelEndYear, as well as an array of global marginal damages (sum of the 16 regions). It accepts index (in above scenarios array of directory names), emissions year and climate sensitivity as input parameters. This code is an adaptation of the marginaldamages.jl file included with FUND. It enables easy calculation of marginal damages resulting from different climate sensitivities, which is the only input parameter. The function returns arrays of regional and global marginal damages, together with the two models themselves.

In [7]:
function calcMarginalDamages(index::Int, emissionYear::Int, climateSensitivity::Float64)
    # Create models m1 and m2, using the parameters from the files in the scenarios[index] directory
    m1 = getfund(datadir=joinpath(dirname(@__FILE__), scenariosPath, scenarios[index]))
    m2 = getfund(datadir=joinpath(dirname(@__FILE__), scenariosPath, scenarios[index]))
    
    # Change the climate sensitivity from the default one supplied in the data files
    setparameter(m1, :climatedynamics, :ClimateSensitivity, climateSensitivity)
    setparameter(m2, :climatedynamics, :ClimateSensitivity, climateSensitivity)
    
    # Add a marginal CO2 component to m2
    addcomponent(m2, adder, :marginalemission, before=:climateco2cycle)
    addem = zeros(yearsToRun+1)
    addem[getindexfromyear(marginalCO2Year):getindexfromyear(marginalCO2Year)+9] = 1.0
    setparameter(m2,:marginalemission,:add,addem)
    connectparameter(m2,:marginalemission,:input,:emissions,:mco2)
    connectparameter(m2, :climateco2cycle,:mco2,:marginalemission,:output)
    
    # Run the models and get two yearsToRun x 16 matrices of damages
    run(m1)
    run(m2)
    damage1 = m1[:impactaggregation, :loss]
    damage2 = m2[:impactaggregation, :loss]
    
    # Matrix of marginal damages between m1 and m2 for each year and region
    # Division by 10^7 is as in marginaldamages.jl  code
    damageDiff = (damage2.-damage1)/10000000.0
    
    # Fill in arrays of regional and global marginal damages
    startIndex = emissionYear - modelStartYear # Starting index in damageDiff for calculating SCC
    numYears = finalYear - emissionYear
    regionalMDs = zeros(Float64, numYears, 16) # Create 16-column matrix of Float64 initialized to 0.0
    globalMDs = zeros(numYears) # Create an array of Float64 initialized to 0.0
    for i = 1:numYears
        for j = 1:16
            md = damageDiff[startIndex + i - 1, j]
            regionalMDs[i, j] = md
            globalMDs[i] += md
        end
     end
    # Return regional & global marginal damages
    return regionalMDs, globalMDs
end
Out[7]:
calcMarginalDamages (generic function with 1 method)

Function discountedValue calculates present value of an array of future values.

In [8]:
function discountedValue(numYears::Int, discountRate::Float64, values = Float64[])
    discountedSum = 0.0
    discountFactor = 1.0
    for i = 1:numYears
        discountFactor *= (1.0 + discountRate) # discountRate is a fraction, not percentage
        discountedSum += values[i] / discountFactor
    end
    return discountedSum
end
Out[8]:
discountedValue (generic function with 2 methods)
In [9]:
models = Model[] # Empty array of Models
resize!(models, 5) # Resize the array to 5
# Run the models after setting the base ECS (3.5)
for i = 1:5
    models[i] = getfund(datadir=joinpath(dirname(@__FILE__), scenariosPath, scenarios[i]))
    setparameter(models[i], :climatedynamics, :ClimateSensitivity, baseECS)
    run(models[i])
end

FUND's Component Connections¶

In [10]:
# Print a list for models[1] - the others are identical
println(models[1])
showing model component connections:
1. scenariouncertainty component
    incoming parameters:
      none
    outgoing variables:
      - pgrowth in population component
      - pgrowth in socioeconomic component
      - ypcgrowth in socioeconomic component
      - forestemm in emissions component
      - aeei in emissions component
      - acei in emissions component
      - ypcgrowth in emissions component
2. population component
    incoming parameters:
      - pgrowth from scenariouncertainty component
      - enter from impactsealevelrise component
      - leave from impactsealevelrise component
      - dead from impactdeathmorbidity component
    outgoing variables:
      - globalpopulation in socioeconomic component
      - populationin1 in socioeconomic component
      - population in socioeconomic component
      - population in emissions component
      - population in impactagriculture component
      - population in impactbiodiversity component
      - population in impactcardiovascularrespiratory component
      - population in impactcooling component
      - population in impactdiarrhoea component
      - population in impactextratropicalstorms component
      - population in impactforests component
      - population in impactheating component
      - population in impactvectorbornediseases component
      - population in impacttropicalstorms component
      - population in vslvmorb component
      - population in impactdeathmorbidity component
      - population in impactwaterresources component
      - population in impactsealevelrise component
3. geography component
    incoming parameters:
      - landloss from impactsealevelrise component
    outgoing variables:
      - area in socioeconomic component
      - area in impactsealevelrise component
4. socioeconomic component
    incoming parameters:
      - area from geography component
      - globalpopulation from population component
      - populationin1 from population component
      - population from population component
      - pgrowth from scenariouncertainty component
      - ypcgrowth from scenariouncertainty component
      - eloss from impactaggregation component
      - sloss from impactaggregation component
      - mitigationcost from emissions component
    outgoing variables:
      - income in emissions component
      - income in impactagriculture component
      - income in impactbiodiversity component
      - plus in impactcardiovascularrespiratory component
      - urbpop in impactcardiovascularrespiratory component
      - income in impactcooling component
      - income in impactdiarrhoea component
      - income in impactextratropicalstorms component
      - income in impactforests component
      - income in impactheating component
      - income in impactvectorbornediseases component
      - income in impacttropicalstorms component
      - income in vslvmorb component
      - income in impactwaterresources component
      - income in impactsealevelrise component
      - income in impactaggregation component
5. emissions component
    incoming parameters:
      - income from socioeconomic component
      - population from population component
      - forestemm from scenariouncertainty component
      - aeei from scenariouncertainty component
      - acei from scenariouncertainty component
      - ypcgrowth from scenariouncertainty component
    outgoing variables:
      - mitigationcost in socioeconomic component
      - mco2 in climateco2cycle component
      - globch4 in climatech4cycle component
      - globn2o in climaten2ocycle component
      - globsf6 in climatesf6cycle component
      - cumaeei in impactcooling component
      - cumaeei in impactheating component
6. climateco2cycle component
    incoming parameters:
      - mco2 from emissions component
      - temp from climatedynamics component
    outgoing variables:
      - acco2 in climateforcing component
      - acco2 in impactagriculture component
      - acco2 in impactextratropicalstorms component
      - acco2 in impactforests component
7. climatech4cycle component
    incoming parameters:
      - globch4 from emissions component
    outgoing variables:
      - acch4 in climateforcing component
8. climaten2ocycle component
    incoming parameters:
      - globn2o from emissions component
    outgoing variables:
      - acn2o in climateforcing component
9. climatesf6cycle component
    incoming parameters:
      - globsf6 from emissions component
    outgoing variables:
      - acsf6 in climateforcing component
10. climateforcing component
    incoming parameters:
      - acco2 from climateco2cycle component
      - acch4 from climatech4cycle component
      - acn2o from climaten2ocycle component
      - acsf6 from climatesf6cycle component
    outgoing variables:
      - radforc in climatedynamics component
11. climatedynamics component
    incoming parameters:
      - radforc from climateforcing component
    outgoing variables:
      - temp in climateco2cycle component
      - temp in climateregional component
      - temp in biodiversity component
      - temp in ocean component
12. biodiversity component
    incoming parameters:
      - temp from climatedynamics component
    outgoing variables:
      - nospecies in impactbiodiversity component
13. climateregional component
    incoming parameters:
      - inputtemp from climatedynamics component
    outgoing variables:
      - temp in impactagriculture component
      - temp in impactbiodiversity component
      - temp in impactcardiovascularrespiratory component
      - temp in impactcooling component
      - regtmp in impactdiarrhoea component
      - temp in impactforests component
      - temp in impactheating component
      - temp in impactvectorbornediseases component
      - regstmp in impacttropicalstorms component
      - temp in impactwaterresources component
14. ocean component
    incoming parameters:
      - temp from climatedynamics component
    outgoing variables:
      - sea in impactsealevelrise component
15. impactagriculture component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - temp from climateregional component
      - acco2 from climateco2cycle component
    outgoing variables:
      - agcost in impactaggregation component
16. impactbiodiversity component
    incoming parameters:
      - temp from climateregional component
      - nospecies from biodiversity component
      - income from socioeconomic component
      - population from population component
    outgoing variables:
      - species in impactaggregation component
17. impactcardiovascularrespiratory component
    incoming parameters:
      - population from population component
      - temp from climateregional component
      - plus from socioeconomic component
      - urbpop from socioeconomic component
    outgoing variables:
      - cardheat in impactdeathmorbidity component
      - cardcold in impactdeathmorbidity component
      - resp in impactdeathmorbidity component
18. impactcooling component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - temp from climateregional component
      - cumaeei from emissions component
    outgoing variables:
      - cooling in impactaggregation component
19. impactdiarrhoea component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - regtmp from climateregional component
    outgoing variables:
      - diadead in impactdeathmorbidity component
      - diasick in impactdeathmorbidity component
20. impactextratropicalstorms component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - acco2 from climateco2cycle component
    outgoing variables:
      - extratropicalstormsdead in impactdeathmorbidity component
      - extratropicalstormsdam in impactaggregation component
21. impactforests component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - temp from climateregional component
      - acco2 from climateco2cycle component
    outgoing variables:
      - forests in impactaggregation component
22. impactheating component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - temp from climateregional component
      - cumaeei from emissions component
    outgoing variables:
      - heating in impactaggregation component
23. impactvectorbornediseases component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - temp from climateregional component
    outgoing variables:
      - dengue in impactdeathmorbidity component
      - schisto in impactdeathmorbidity component
      - malaria in impactdeathmorbidity component
24. impacttropicalstorms component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - regstmp from climateregional component
    outgoing variables:
      - hurrdead in impactdeathmorbidity component
      - hurrdam in impactaggregation component
25. vslvmorb component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
    outgoing variables:
      - vsl in impactdeathmorbidity component
      - vmorb in impactdeathmorbidity component
26. impactdeathmorbidity component
    incoming parameters:
      - vsl from vslvmorb component
      - vmorb from vslvmorb component
      - population from population component
      - dengue from impactvectorbornediseases component
      - schisto from impactvectorbornediseases component
      - malaria from impactvectorbornediseases component
      - cardheat from impactcardiovascularrespiratory component
      - cardcold from impactcardiovascularrespiratory component
      - resp from impactcardiovascularrespiratory component
      - diadead from impactdiarrhoea component
      - hurrdead from impacttropicalstorms component
      - extratropicalstormsdead from impactextratropicalstorms component
      - diasick from impactdiarrhoea component
    outgoing variables:
      - dead in population component
      - deadcost in impactaggregation component
      - morbcost in impactaggregation component
27. impactwaterresources component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - temp from climateregional component
    outgoing variables:
      - water in impactaggregation component
28. impactsealevelrise component
    incoming parameters:
      - population from population component
      - income from socioeconomic component
      - sea from ocean component
      - area from geography component
    outgoing variables:
      - landloss in geography component
      - enter in population component
      - leave in population component
      - drycost in impactaggregation component
      - protcost in impactaggregation component
      - entercost in impactaggregation component
      - wetcost in impactaggregation component
      - leavecost in impactaggregation component
29. impactaggregation component
    incoming parameters:
      - income from socioeconomic component
      - heating from impactheating component
      - cooling from impactcooling component
      - agcost from impactagriculture component
      - species from impactbiodiversity component
      - water from impactwaterresources component
      - hurrdam from impacttropicalstorms component
      - extratropicalstormsdam from impactextratropicalstorms component
      - forests from impactforests component
      - drycost from impactsealevelrise component
      - protcost from impactsealevelrise component
      - entercost from impactsealevelrise component
      - deadcost from impactdeathmorbidity component
      - morbcost from impactdeathmorbidity component
      - wetcost from impactsealevelrise component
      - leavecost from impactsealevelrise component
    outgoing variables:
      - eloss in socioeconomic component
      - sloss in socioeconomic component

In [11]:
populations = zeros(finalYear - modelStartYear + 1, 5)
emissionsCO2 = zeros(finalYear - modelStartYear + 1, 5)
concentrationsCO2 = zeros(finalYear - modelStartYear + 1, 5)
globalTempRise = zeros(finalYear - modelStartYear + 1, 5)
globalSeaLevelRise = zeros(finalYear - modelStartYear + 1, 5)
for i = 1:(finalYear - modelStartYear + 1)
    for j = 1:5
        populations[i, j] = models[j][:socioeconomic, :globalpopulation][i] / 1.0e9
        # 2017-11-11: Added ratioCo2ToC to following line.
        emissionsCO2[i, j] = models[j][:climateco2cycle, :mco2][i] * ratioCo2ToC / 1000.0
        concentrationsCO2[i, j] = models[j][:climateco2cycle, :acco2][i]
        globalTempRise[i, j] = models[j][:climatedynamics, :temp][i]
        globalSeaLevelRise[i, j] = models[j][:impactsealevelrise, :sea][i]
    end
end
In [12]:
plotGlobalPop = Plots.plot(modelStartYear:finalYear, populations, grid=true, legend=true, 
    lw=2, label=scenarioLabels, 
    xlims=(1950,2300), ylims=(2,12), xticks=(1950:50:2300), yticks=(2.0:1.0:11.0),
    title="Fig. 2 - Global Population for 5 Scenarios", xlabel="year", ylabel="billions",
    titlefont=font(11))
Out[12]:
In [13]:
# Save Fig. 2 as a .png file in the same directory as this notebook.
# If the file Fig_2 already exists, it is overwritten.
png("Fig_2")
In [14]:
plotCO2Emissions = Plots.plot(modelStartYear:finalYear, emissionsCO2, grid=true, legend=true,
    lw=2, label=scenarioLabels,
    xlims=(1950,2300), ylims=(0,120), xticks=(1950:50:2300), yticks=(0:10:120),
    xlabel="year", ylabel="Gt(CO2)/year", guidefont=font(9))
Out[14]:
In [15]:
plotCO2Concentrations = Plots.plot(modelStartYear:finalYear, concentrationsCO2, grid=true, legend=true, 
    lw=2, label=scenarioLabels,
    xlims=(1950,2300), ylims=(300,1400), xticks=(1950:50:2300), yticks=(0:200:1300),
    xlabel="year", ylabel="parts per million", guidefont=font(9))
Out[15]:
In [16]:
plot(plotCO2Emissions, plotCO2Concentrations,
    title = ["Fig. 3 - CO2 Emissions" "Fig. 4 - CO2 Concentration"],
    titlefont=font(10))
Out[16]:
In [17]:
png("Figs_3_4")
In [18]:
plotTempRise = Plots.plot(modelStartYear:finalYear, globalTempRise, grid=true, legend=true,
    lw=2, label=scenarioLabels,
    xlims=(1950,2300), ylims=(0,10), xticks=(1950:50:2300), yticks=(0:2:10),
    xlabel="year", ylabel="°C", guidefont=font(9))
Out[18]:
In [19]:
plotSeaRise = Plots.plot(modelStartYear:finalYear, globalSeaLevelRise, grid=true, legend=true, 
    lw=2, label=scenarioLabels,
    xlims=(1950,2300), ylims=(0,5.5), xticks=(1950:50:2300), yticks=(0:1:5),
    xlabel="year", ylabel="metres", guidefont=font(9))
Out[19]:
In [20]:
plot(plotTempRise, plotSeaRise,
    title = ["Fig. 5 - Temperature Rise, ECS=3.5°C" "Fig. 6 - Sea-level Rise, ECS=3.5°C"],
    titlefont=font(10))
Out[20]:
In [21]:
png("Figs_5_6")
In [22]:
# Change the models' climate sensitivies to 1.0, rerun them and extract the temperature and sea-level rises.
for i = 1:5
    setparameter(models[i], :climatedynamics, :ClimateSensitivity, newECS)
    run(models[i])
end
globalTempRise = zeros(finalYear - modelStartYear + 1, 5)
globalSeaLevelRise = zeros(finalYear - modelStartYear + 1, 5)
for i = 1:(finalYear - modelStartYear + 1)
    for j = 1:5
        globalTempRise[i, j] = models[j][:climatedynamics, :temp][i]
        globalSeaLevelRise[i, j] = models[j][:impactsealevelrise, :sea][i]
    end
end
In [23]:
plotTempRise = Plots.plot(modelStartYear:finalYear, globalTempRise, grid=true, legend=true,
    lw=2, label=scenarioLabels,
    xlims=(1950,2300), ylims=(0,3.0), xticks=(1950:50:2300), yticks=(0:0.5:3.0),
    xlabel="year", ylabel="°C", guidefont=font(9))
Out[23]:
In [24]:
plotSeaRise = Plots.plot(modelStartYear:finalYear, globalSeaLevelRise, grid=true, legend=true, 
    lw=2, label=scenarioLabels,
    xlims=(1950,2300), ylims=(0,2.0), xticks=(1950:50:2300), yticks=(0:0.5:2),
    xlabel="year", ylabel="metres", guidefont=font(9))
Out[24]: