FUND Model for FoS September 2017 Newsletter¶
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]
# 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
emissionsCO2[i, j] = models[j][:climateco2cycle, :mco2][i] / 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,32), xticks=(1950:50:2300), yticks=(0:10:30),
xlabel="year", ylabel="Gt/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]: