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")
Out[1]:
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]:
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)
In [6]:
# Use these ECS values
baseECS = 3.5
newECS = 1.0
Out[6]:
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]:
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]:
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])
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]:
In [25]:
plot(plotTempRise, plotSeaRise,
title = ["Fig. 7 - Temperature Rise, ECS=1.0°C" "Fig. 8 - Sea-level Rise, ECS=1.0°C"],
titlefont=font(10))
Out[25]:
In [26]:
png("Figs_7_8")
In [27]:
# Calculate regional and global marginal damages for the Base scenario (1), CO2 release in 2017,
# with ECS of 3.5 (baseECS) and 1.0 (newECS)
regMDsECS1, globalMDsECS1 = calcMarginalDamages(1, marginalCO2Year, baseECS)
regMDsECS2, globalMDsECS2 = calcMarginalDamages(1, marginalCO2Year, newECS)
Out[27]:
In [28]:
# Create matrix of marginal damages and plot it
bothGlobalMDs = zeros(numMdYears, 2) # 2-D matrix of zeros
for i =1:numMdYears
# Put both sets in MDs into numMdYears x 2 matrix
# 2017-11-11: Added ratioCo2ToC in next 2 lines.
bothGlobalMDs[i, 1] = globalMDsECS1[i] * convert1995To2016Dollars / ratioCo2ToC
bothGlobalMDs[i, 2] = globalMDsECS2[i] * convert1995To2016Dollars / ratioCo2ToC
end
plotMDs = Plots.plot((marginalCO2Year+1):finalYear, bothGlobalMDs, grid=true, legend=true,
lw=2, label=["ECS = 3.5°C" "ECS = 1.0°C"],
xlims=(2000,2300), ylims=(-0.5,3.5), xticks=(2000:50:2300), yticks=(-0.5:0.5:3.5),
title="Fig. 9 - Base Scenario Global Marginal Damages for 1 t of CO2 Emitted in $marginalCO2Year",
xlabel="year", ylabel="2016 USD/year", titlefont=font(10))
Out[28]:
In [29]:
png("Fig_9")
In [30]:
# Plot regional marginal damages for 2018-2030 for regiona 1-16 to find what regions are responsible for
# sharp changes in global plot for ECS = 3.5
# 2017-11-11: Added ratioCo2ToC to following line.
plotRegMDs = Plots.plot((marginalCO2Year+1):finalYear, regMDsECS1 * convert1995To2016Dollars / ratioCo2ToC,
grid=true, legend=true, lw=2, label=regionLabels,
xlims=(2000,2300), ylims=(-0.4,1.0), xticks=(2000:50:2300), yticks=(-0.4:0.1:1.0),
title="Regional Marginal Damages for 1 t of CO2 Emitted in $marginalCO2Year, ECS = 3.5",
xlabel="year", ylabel="2016 USD/year", titlefont=font(10))
Out[30]:
In [31]:
# Above plot is not clear, so plot regional marginal damages for regions 1-8 and then 9-16
regMDs = zeros(numMdYears, 8)
regLabels = ["USA" "CAN" "WEU" "JPK" "ANZ" "EEU" "FSU" "MDE"]
for i = 1:numMdYears
for j = 1:8
# 2017-11-11: Added ratioCo2ToC to following line
regMDs[i, j] = regMDsECS1[i, j] * convert1995To2016Dollars / ratioCo2ToC
end
end
Plots.plot((marginalCO2Year+1):finalYear, regMDs, grid=true, legend=true,
lw=2, label=regLabels,
xlims=(2000,2300), ylims=(-0.3,0.4), xticks=(2000:50:2300), yticks=(-0.3:0.1:0.4),
title="Regional Marginal Damages for 1 t of CO2 Emitted in $marginalCO2Year, ECS = 3.5",
xlabel="year", ylabel="2016 USD/year", titlefont=font(10))
Out[31]:
In [32]:
# Repeat above for regions 9-16
regMDs = zeros(numMdYears, 8)
regLabels = ["CAM" "LAM" "SAS" "SEA" "CHI" "MAF" "SSA" "SIS"]
for i = 1:numMdYears
for j = 1:8
# 2017-11-11: Added ratioCo2ToC to following line
regMDs[i, j] = regMDsECS1[i, j+8] * convert1995To2016Dollars / ratioCo2ToC
end
end
Plots.plot((marginalCO2Year+1):finalYear, regMDs, grid=true, legend=true,
lw=2, label=regLabels,
xlims=(2000,2300), ylims=(-0.4,0.9), xticks=(2000:50:2300), yticks=(-0.4:0.1:0.9),
title="Regional Marginal Damages for 1 t of CO2 Emitted in $marginalCO2Year, ECS = 3.5",
xlabel="year", ylabel="2016 USD/year",
titlefont=font(10))
Out[32]:
In [33]:
# Repeat above regional plots for ECS = 1.0
regMDs = zeros(numMdYears, 8)
regLabels = ["USA" "CAN" "WEU" "JPK" "ANZ" "EEU" "FSU" "MDE"]
for i = 1:numMdYears
for j = 1:8
# 2017-11-11: Added ratioCo2ToC to following line
regMDs[i, j] = regMDsECS1[i, j] * convert1995To2016Dollars / ratioCo2ToC
end
end
Plots.plot((marginalCO2Year+1):finalYear, regMDs, grid=true, legend=true,
lw=2, label=regLabels,
xlims=(2000,2300), ylims=(-0.3,0.35), xticks=(2000:50:2300), yticks=(-0.3:0.1:0.3),
title="Regional Marginal Damages for 1 t of CO2 Emitted in $marginalCO2Year, ECS = $newECS",
xlabel="year", ylabel="2016 USD/year", titlefont=font(10))
Out[33]:
In [34]:
# Repeat above for regiona 9-16
regMDs = zeros(numMdYears, 8)
regLabels = ["CAM" "LAM" "SAS" "SEA" "CHI" "MAF" "SSA" "SIS"]
for i = 1:numMdYears
for j = 1:8
# 2017-11-11: Added ratioCo2ToC to following line
regMDs[i, j] = regMDsECS1[i, j+8] * convert1995To2016Dollars / ratioCo2ToC
end
end
Plots.plot((marginalCO2Year+1):finalYear, regMDs, grid=true, legend=true,
lw=2, label=regLabels,
xlims=(2000,2300), ylims=(-0.4,0.9), xticks=(2000:50:2300), yticks=(-0.4:0.1:0.9),
title="Regional Marginal Damages for 1 t of CO2 Emitted in $marginalCO2Year, ECS = $newECS",
xlabel="year", ylabel="2016 USD/year", titlefont=font(10))
Out[34]:
SCCO2 Calculation¶
In [35]:
# Calculate table of SCCO2 for ECS of 3.5 (baseECS) and 1.0 (new ECS) for range of CO2 emission years
# Show SCCO2 in 1995 and 2016 $US
println("SCCO2 Calculations by Year of Emission, for ECS of $baseECS °C and $newECS °C, using 3% Discount Rate")
for i = 1:5
regMDsECS1, globalMDsECS1 = calcMarginalDamages(1, marginalCO2Years[i], baseECS)
regMDsECS2, globalMDsECS2 = calcMarginalDamages(1, marginalCO2Years[i], newECS)
# 2017-11-11: Added ratioCo2ToC to following lines
scc1_1995 = discountedValue(finalYear - marginalCO2Years[i], discountRate, globalMDsECS1) / ratioCo2ToC
scc1_2016 = scc1_1995 * convert1995To2016Dollars
scc2_1995 = discountedValue(finalYear - marginalCO2Years[i], discountRate, globalMDsECS2) / ratioCo2ToC
scc2_2016 = scc2_1995 * convert1995To2016Dollars
year = marginalCO2Years[i]
println("$year $scc1_1995 $scc1_2016 $scc2_1995 $scc2_2016")
end
In [36]:
ecs = 1.0
deltaECS = 0.2
ecsValues = zeros(41)
scco2Values = zeros(41)
for i = 1:40
regMDs, globalMDs = calcMarginalDamages(1, marginalCO2Year, ecs)
scco2 = discountedValue(finalYear - marginalCO2Year, discountRate, globalMDs)
ecsValues[i] = ecs
# 2017-11-11: Added ratioCo2ToC to following line
scco2Values[i] = scco2 * convert1995To2016Dollars / ratioCo2ToC
ecs += deltaECS
end
In [37]:
Plots.plot(ecsValues, scco2Values, seriestype=:line, grid=true, legend=false, lw=2,
xlims=(1,9), ylims=(-10,20), xticks=(1:1:9), yticks=(-10:5:20),
title="Fig. 10 - SCCO2 vs ECS for Base Scenario, Year of Emission: 2020",
xlabel="ECS in °C", ylabel="2016 USD/t(CO2)", titlefont=font(10))
Out[37]:
In [38]:
png("Fig_10")
In [ ]: