# 1. Start a function that takes the model parameters that might change as arguments i.e albedo and emissivity
ClimateChangeModel<- function(Propn.ocean, alpha, epsilon, No.Years, StartingT.s, plot=T){
# 2. SET MODEL PARAMETERS
A = 4*pi*(6378000^2) #earths surface area
C = (A*Propn.ocean*100*1025*3.993)+(A*(1-Propn.ocean)*5*1600*0.8) # total surface heat capacity
S = 1368 #solar constant, in Watts
sigma = 5.6704e-8 # Stefan-Boltzmann constant using scientific notation
Time.int = (365/2)*24*60*60 #Time interval : the number of seconds in half a year
# 3. INITIATE VECTORS
Year = seq(from=0, to=20, by=0.5) # 6-month intervals
Rate.in = rep(S/4, length(Year)) # Rate in is a constant
Rate.out = numeric(length(Year)) # Empty Vectors
DeltaT = numeric(length(Year))
T.s = numeric(length(Year))
T.s[1] = StartingT.s # Set initial T.s
T.a = numeric(length(Year))
# 4. LOOP for every timestep
for (i in 1:length(Year)) {
# calculate T.a from T.s for the ith timestep
T.a[i] = T.s[i]/(2^0.25)
# calucute rate out for the ith timestep
Rate.out[i] = alpha*S/4 + (1-epsilon)*sigma*T.s[i]^4 + epsilon*sigma*T.a[i]^4
# Calculate DeltaT for the ith timestep
DeltaT[i] = (1/C)*(Rate.in[i]-Rate.out[i])*(A/1000)*Time.int
# If we have not reached the last timestep, Add the change in temperature to the surface
# tempeture to be used in the next timestep. If it is the last timestep, do nothing
if (i < length(Year)) {
T.s[i+1] = T.s[i] + DeltaT[i]
} else {}
}
# 5. COMBINE INTO DATA FRAME FOR PLOTTING
ClimateChange = data.frame(Year, Rate.in, Rate.out, DeltaT, T.s, T.a)
#print(round(head(ClimateChange),2))
# 6. PLOT RESULTS
if (plot==T){
with(ClimateChange, plot(T.s ~ Year, type="o", col="blue",
ylim=c(min(T.a)-10,max(T.s)+10), ylab="Temperature (K)", xlab = "Time (years)",
main = "Climate Change Model"))
lines(T.a ~ Year, type="o", col="red")
axis(1, at = seq(1,20, by=1), labels = FALSE, tcl = -0.25)
axis(2, at = seq(0,400, by=10), labels = FALSE, tcl = -0.25)
legend("topright", c(expression(T[s]), expression(T[a])), lty = c(1, 1), col = c("blue", "red"), horiz = TRUE)
}
return(ClimateChange)
}
a0.3 <- ClimateChangeModel(Propn.ocean=0.7, alpha=0.3, epsilon = 0.78, No.Years=20, StartingT.s = 288)
a0.1 <- ClimateChangeModel(Propn.ocean=0.7, alpha=0.1, epsilon = 0.78, No.Years=20, StartingT.s = 288)
The temperature increases, plateauing after ~5 years, as equilibrium has been reached and incoming and outgoing radiation are once again balanced.
A lower albedo means the reflectivity of the Earth has decreased, so there is less incoming solar radiation to be deflected and more entering the planet. This increase in thermal energy raises both atmospheric and surface temperatures.