Modification of Earth system model to include cloud and moisture feedbacks

In Section 6.8, a model of the atmosphere was developed and the Exercises investigated the effects of making changes to the parameterised albedo and emissivity. Here, you will build upon that model by considering how changes in atmospheric water vapour may affect the modelled temperatures.

If the temperature of the atmosphere warms, the carrying capacity of water in the air increases and the evaporation of water increases. This could have two effects.

Firstly, as water vapour is a very powerful greenhouse gas, the emissivity of the atmosphere increases; i.e., the greenhouse effect increases. This is a positive feedback because a warming climate could produce a stronger greenhouse effect that will warm the planet further allowing more water vapour that will in turn produce more warming, etc. If the feedback is strong enough, this could cause a runaway greenhouse effect where the ocean could boil. Atmospheric emissivity at atmospheric temperature Ta can be modelled by the linear relation

\[ \begin{align*} \epsilon_a &= \epsilon_{ac}(T_a-T_{a0})\ + \epsilon_{ao} \tag{1} \end{align*} \] where εac is the emissivity variation coefficient, and εao is the emissivity at the reference atmospheric temperature Ta0. Note: emissivity values are only valid between 0 and 1.

The second feedback due to water vapour could be that the larger amount of water vapour could produce more clouds that could increase the planetary albedo. This would reflect more sunlight back to space thus cooling the planet. This feedback is a negative feedback. If warming is caused by some other factor, more clouds will limit the warming. Also, if there was some other effect to cool the planet, the lower albedo, due to fewer clouds, would tend to heat the planet. It is a negative feedback because the original perturbation is made smaller. Albedo at atmospheric temperature Ta can be modelled by the linear relation

\[ \begin{align*} \alpha &= \alpha_{c}(T_a-T_{a0})\ + \alpha_{o} \tag{2} \end{align*} \] where αc is the albedo variation coefficient, Ta0 is a reference atmospheric temperature and αo is the albedo at the reference temperature. Note: albedo values are only valid between 0 and 1.

Initialisation

  1. Use initial rough values of Ta0 = 240 K, αc = 0.01 K-1, αo = 0.3, εac = 0.01 K-1 and εao = 0.78 to plot (using R) how atmospheric and surface temperature vary through a 20-year period. Do the values converge to an asymptote? If so, comment on how the values compare with the model from Section 6.8 (with constant emissivity and albedo).
# 1. Start a function that takes the model parameters that might change as arguments i.e albedo and emissivity

ClimateChangeModel<- function(Propn.ocean, No.Years, StartingT.s, T.a0, alpha.c, alpha.o, epsilon.ac, epsilon.ao, 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
      # T.a0 = 240
      # epsilon.ac = 0.01 # emissivity variation coefficient, in K^-1
      # epsilon.ao = 0.78
      # alpha.c = 0.01 # albedo variation coefficient, in K^-1
      # alpha.o = 0.3
      
      # 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))
      alpha = numeric(length(Year))
      epsilon = 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)
        # calculate alpha for the ith timestep
        alpha[i] = alpha.c * (T.a[i] - T.a0) + alpha.o
        # print(alpha[i])
        # calculate epsilon for the ith timestep
        epsilon[i] = epsilon.ac * (T.a[i] - T.a0) + epsilon.ao
        # print(epsilon[i])
        # calculate rate out for the ith timestep
        Rate.out[i] = alpha[i]*S/4 + (1-epsilon[i])*sigma*T.s[i]^4 + epsilon[i]*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, alpha, epsilon)
      # print(round((ClimateChange),5))
      # 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)
}

test <- ClimateChangeModel(Propn.ocean=0.7, No.Years=20, StartingT.s = 288, T.a0 = 240, alpha.c = 0.01, alpha.o = 0.3, epsilon.ac = 0.01, epsilon.ao = 0.78)

no_variation.test <- ClimateChangeModel(Propn.ocean=0.7, No.Years=20, StartingT.s = 288, T.a0 = 240, alpha.c = 0, alpha.o = 0.3, epsilon.ac = 0, epsilon.ao = 0.78)

Part A

Explore the model using various values of αc and εac. First, increment (up and down) the albedo variation coefficient, including use of negative values, whilst holding the emissivity variation coefficient constant. Test if the albedo and/or atmospheric emissivity can be made to go outside the valid range; record temperature values where convergence occurs.

  1. With emissivity variation constant (= 0.01 K-1), determine the maximum and minimum values of the albedo variation coefficient for which valid values of albedo and emissivity occur.
results <- data.frame(alpha.c = numeric(0), alpha_out_of_range = logical(0), epsilon_out_of_range = logical(0))

for (alpha.c in seq(-1, 1, 0.01)) {
  result <- ClimateChangeModel(
    Propn.ocean = 0.7, 
    No.Years = 20, 
    StartingT.s = 288, 
    T.a0 = 240, 
    alpha.c = alpha.c, 
    alpha.o = 0.3, 
    epsilon.ac = 0.01, 
    epsilon.ao = 0.78,
    plot = FALSE
  )
  
  alpha_out_of_range <- any(result$alpha < 0 | result$alpha > 1)
  epsilon_out_of_range <- any(result$epsilon < 0 | result$epsilon > 1)
  
  results <- rbind(results, data.frame(alpha.c, alpha_out_of_range, epsilon_out_of_range))
}

# print(results)

valid_range <- round(results$alpha.c[!results$epsilon_out_of_range & !results$alpha_out_of_range],3)
paste("The maximum and minimum values of the albedo variation coefficient for which valid values of albedo and emissivity occur is", max(valid_range), "and", min(valid_range), "respectively.")
## [1] "The maximum and minimum values of the albedo variation coefficient for which valid values of albedo and emissivity occur is 0.12 and 0 respectively."

Second, vary the emissivity variation coefficient, including use of negative values, whilst holding the albedo variation coefficient constant (repeat the tests; record values).

  1. With albedo variation constant (= 0.01 K-1), determine the maximum and minimum values of the emissivity variation coefficient for which valid values of albedo and emissivity occur.
results <- data.frame(epsilon.ac = numeric(0), alpha_out_of_range = logical(0), epsilon_out_of_range = logical(0))

for (epsilon.ac in seq(-1, 1, 0.01)) {
  result <- ClimateChangeModel(
    Propn.ocean = 0.7, 
    No.Years = 20, 
    StartingT.s = 288, 
    T.a0 = 240, 
    alpha.c = 0.01, 
    alpha.o = 0.3, 
    epsilon.ac = epsilon.ac, 
    epsilon.ao = 0.78,
    plot = FALSE
  )
  
  alpha_out_of_range <- any(result$alpha < 0 | result$alpha > 1)
  epsilon_out_of_range <- any(result$epsilon < 0 | result$epsilon > 1)
  
  results <- rbind(results, data.frame(epsilon.ac, alpha_out_of_range, epsilon_out_of_range))
}

# print(results)

valid_range <- round(results$epsilon.ac[!results$epsilon_out_of_range & !results$alpha_out_of_range],3)
paste("The maximum and minimum values of the emissivity variation coefficient for which valid values of albedo and emissivity occur is", max(valid_range), "and", min(valid_range), "respectively.")
## [1] "The maximum and minimum values of the emissivity variation coefficient for which valid values of albedo and emissivity occur is 0.02 and -0.16 respectively."

Third, vary both parameters (guided by the values found when varying only one parameter at a time) and repeat the tests/record values.

  1. Replicate the plot in #1 for key parameter values that have interesting outcomes and that produce values at/near limits of the valid range of albedo and emissivity (see previous notes on valid ranges). You should aim to present 4-8 scenarios. Plot atmospheric temperature from all scenarios on a single plot (using different colours for each scenario) and surface temperature from all scenarios on a second plot (using matching colours).
Year = seq(0, 20, 0.5)
alpha.c_params <- c(0.12,0,0.01,0.01,0.03,0.04,0.01)
epsilon.ac_params <- c(0.01,0.01,-0.16,0.02,-0.01,-0.03,-0.05)

plot(Year, numeric(length(Year)), ylim=c(283.5,292), ylab = "Surface temperature (K)", xlab = "Time (years)", 
     main = "Surface Temperatures", type="l", col = "white")
axis(1, at = seq(1,20, by=1), labels = FALSE, tcl = -0.25)

colours = rainbow(length(alpha.c_params))
t.s_legend <- data.frame(epsilon.ac = epsilon.ac_params, alpha.c = alpha.c_params)

for (i in 1:8) {
  surface_temp <- ClimateChangeModel(Propn.ocean = 0.7, No.Years = 20, StartingT.s = 288, T.a0 = 240, alpha.o = 0.3, epsilon.ao = 0.78, plot = FALSE, epsilon.ac = epsilon.ac_params[i], alpha.c = alpha.c_params[i])
  
  lines(surface_temp$Year, surface_temp$T.s, type="o", col = colours[i])
}

legend("topright", legend = paste(t.s_legend$epsilon.ac, t.s_legend$alpha.c), col = colours, lty = 1, 
       title = "         ε.ac   α.c", cex = 0.8)

plot(Year, numeric(length(Year)), ylim=c(238,246), ylab = "Atmospheric temperature (K)", xlab = "Time (years)", 
     main = "Atmospheric Temperatures", type="l", col = "white")
axis(1, at = seq(1,20, by=1), labels = FALSE, tcl = -0.25)

for (i in 1:8) {
  surface_temp <- ClimateChangeModel(Propn.ocean = 0.7, No.Years = 20, StartingT.s = 288, T.a0 = 240, alpha.o = 0.3, epsilon.ao = 0.78, plot = FALSE, epsilon.ac = epsilon.ac_params[i], alpha.c = alpha.c_params[i])
  
  lines(surface_temp$Year, surface_temp$T.a, type="o", col = colours[i])
}

legend("topright", legend = paste(t.s_legend$epsilon.ac, t.s_legend$alpha.c), col = colours, lty = 1, 
       title = "         ε.ac   α.c", cex = 0.8)

  1. By considering outcomes for the first and second experiments, describe the relative sensitivity of the temperatures to changes in each of the parameters. Sensitivity is a measure of how much change in the outcome occurs with change to the input value (e.g., zero change in outcome would indicate zero sensitivity).

As αc is negative feedback and εac is positive feedback, we would expect increasing αc decreasing εac to have a roughly similar effect on temperature. When αc and εac are at their extremes (0.12 and -0.16 respectively), we see a greater decrease in temperature for αc, indicating a greater sensitivity towards it. So, the albedo feedback mechanism could influence temperature more, hence, changes in αc have a greater impact on outcome and the model is more sensitive to changes in αc compared to εac.

  1. Create a table that summarises outcomes for several variation coefficient values (albedo across the table, emissivity down), reporting the temperature asymptote value (if an asymptote occurs), otherwise describing the behaviour of the temperature (including noting whether one or both parameters become invalid). Include the maximum/minimum values you found in #2 and #3, as well as 0.01 K-1 on each table axis, and insert an appropriate number of other values on each axis to demonstrate any patterns (aim for around seven values for each parameter).
Table 1: asymptotic and oscillatory properties of surface temperature when both αc and εac are varied. Numbers represent final surface temperature in Kelvin, x means temperature diverges, i.e. does not asymptote, * means αc and/or εac is invalid over the 20 years, and † means the surface temperature oscillates between two values over time.
εac αc = 0.00 αc = 0.01 αc = 0.04 αc = 0.06 αc = 0.08 αc = 0.10 αc = 0.12
0.02 x* 288.899 286.267 285.981 285.838† 285.753† 285.673†
0.01 291.480 287.606 286.162 285.933 285.811† 285.735† 284.788†
0.00 288.433 287.024 286.080 285.892 285.786† 285.719† x*†
-0.01 287.435 286.689 286.015 285.857 285.765† 285.702† x*†
-0.04 286.433 286.200† 285.878† 285.778† 285.713† x*† x*†
-0.10 285.925† 285.858† 285.732† 284.023† x*† x*† x*†
-0.16 285.755† 285.723† x*† x*† x*† x*† x*†
  1. Comment on the plots (#4) and the range of variation coefficient values that produce valid outcomes (#6).

Generally, as the variation coefficients get more extreme, the function tends to oscillate. This can be interpreted as the feedback mechanisms of albedo and emissivity having a larger effect compared to other determinants of temperature. At a certain point, the variation coefficients overtake other properties entirely, and temperature diverges. It’s worth noting that some cases on the plots in #4 appear to have reached an asymptote (e.g. εac = -0.16, αc = 0.01), but are just oscillating by miniscule amounts.

When changing the variation coefficents in the direction opposite to their feedback loop, εac negative and αc positive, in moderate amounts, we see lower temperatures from less energy emitted and more energy reflected in a down right through the table. The reverse occurs when changing them in the direction of their feedback loops, (i.e. εac = 0.02, αc = 0.01 and εac = 0.01, αc = 0).

Part B

In Equation 1, replace the atmospheric temperature, Ta, with surface temperature, Ts (including resetting the reference temperature to Ts0 = 288 K).

  1. Redo steps #4-7 and describe the effect for various values of the variation coefficients.
## [1] "The maximum and minimum values of the albedo variation coefficient for which valid values of albedo and emissivity occur is 0.12 and 0 respectively."
## [1] "The maximum and minimum values of the emissivity variation coefficient for which valid values of albedo and emissivity occur is 0.03 and -0.16 respectively."

Step 4B Step 5B

The range of αc is identical, but εac maximum and minimums are slightly different. The key difference compared to part A is εa appears to follow a positive feedback loop as well. This can be seen when εac = 0.02, αc = 0.01, where the curve is decreasing instead of increasing to an asymptote. Additionally, some curves are diverging (e.g. εac = -0.16, αc = 0.01), despite having all valid values, meaning the values likely become invalid after the 20 years.

Step 6B

Table 2: asymptotic and oscillatory properties of surface temperature when both αc and εac are varied. Numbers represent final surface temperature in Kelvin, x means temperature diverges, i.e. does not asymptote, * means αc and/or εac is invalid over the 20 years, and † means the surface temperature oscillates between two values over time.
εac αc = 0.00 αc = 0.01 αc = 0.04 αc = 0.06 αc = 0.08 αc = 0.10 αc = 0.12
0.02 x* 288.899 286.267 285.981† 285.838† 285.231† 285.263†
0.01 289.012 287.606 286.162 285.933† 285.811† 285.576† 284.979†
0.00 288.433 287.024 286.080 285.892† 285.786† 285.719† x*†
-0.01 288.272 287.257 286.300 286.071† 285.936† 285.842† x*†
-0.04 288.129 287.568 286.737† 286.464† 286.282† x*† x*†
-0.10 288.063† 285.858† 285.732† x† x*† x*† x*†
-0.16 288.042† x† x*† x*† x*† x*† x*†

Step 7B

Similar trends as observed in #4, with the notable case of (e.g. εac = -0.16, αc = 0.01), where the curve diverges, but has valid values, although they likely will become invalid after the set time. Oscillating curves generally occur with lower values of αc as well, at 0.06 here compared to 0.04 in part A.

Part C

Replace the linear function of Equation 1 with a non-linear form for emissivity.

\[ \begin{align*} \epsilon_{a} &= ke_s =k\ 6.1094e^{\frac{17.625\ \cdot\ (T_a - 273.15)}{243.04\ +\ (T_a - 273.15)}} \end{align*} \] (the August-Roche-Magnus formula; see https://en.wikipedia.org/wiki/Clausius-Clapeyron_relation).

  1. Determine the value of k where εa = 0.78 at Ta = 240 K. Apply this constant value of k to determine the temperature evolution and describe the effect for various values of the albedo variation coefficient (i.e., repeat steps #4-7).

\[ \begin{align*} \epsilon_{a} &=k\ 6.1094e^{\frac{17.625\ \cdot\ (T_a - 273.15)}{243.04\ +\ (T_a - 273.15)}} \\ 0.78 &=k\ 6.1094e^{\frac{17.625\ \cdot\ (240 - 273.15)}{243.04\ +\ (240 - 273.15)}} \\ 0.78 &=k\cdot 0.37766... \\ k &= 2.06555... \approx2.066 \end{align*} \]

## [1] "The maximum and minimum values of the albedo variation coefficient for which valid values of albedo and emissivity occur is 0.15 and 0.05 respectively."

Step 4C Step 5C

The range of αc is notably different here, between 0.05 to 0.15, potentially due to the more rapid rate of change of εac, so higher values are required to compensate and stablise the curve. But overall, the curves are more predictable as we are only increasing one variable (εac is no longer a parameter).

Step 6C

Table 3: asymptotic and oscillatory properties of surface temperature when αc is varied and epsilon is determined using the August-Roche-Magnus formula. Numbers represent final surface temperature in Kelvin and † means the surface temperature oscillates between two values over time.
αc = 0.05 αc = 0.07 αc = 0.09 αc = 0.11 αc = 0.13 αc = 0.15
287.855 286.325 286.000 285.847† 285.758† 285.693†

Step 7C

Again, a more predictable and simple table as only εac is being altered. We can see the point where oscillation occurs at αc = 0.11. All albedo and epsilon values are valid here, as they are within the limits, and we don’t see the interaction between the changing of both loops at once.

Appendices

Code for Step 6A

alpha.c_params <- c(0,0.01,0.04,0.06,0.08,0.10,0.12)
epsilon.ac_params <- c(0.02,0.01,0,-0.01,-0.04,-0.10,-0.16)

results.vals <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.vals) <- c("alpha.c", as.character(epsilon.ac_params))
results.vals$alpha.c <- alpha.c_params

results.oscillates <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.oscillates) <- c("alpha.c", as.character(epsilon.ac_params))
results.oscillates$alpha.c <- alpha.c_params

results.valid <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.valid) <- c("alpha.c", as.character(epsilon.ac_params))
results.valid$alpha.c <- alpha.c_params

for (i in 1:length(epsilon.ac_params)) {
  for (j in 1:length(alpha.c_params)) {
    result <- ClimateChangeModel(
      Propn.ocean = 0.7, 
      No.Years = 20, 
      StartingT.s = 288, 
      T.a0 = 240, 
      alpha.c = alpha.c_params[j], 
      alpha.o = 0.3, 
      epsilon.ac = epsilon.ac_params[i], 
      epsilon.ao = 0.78,
      plot = FALSE
    )
    
    temp.diff <- result$DeltaT
    is_asymptotic <- TRUE
    is_oscillating <- FALSE
    
        if (abs(temp.diff[3]) >= abs(temp.diff[2])) {
          is_asymptotic <- FALSE
        }
    
    
    for (year in 2:length(temp.diff) - 1) { 
    if (!is.na(temp.diff[year]) && !is.na(temp.diff[year + 1]) && sign(temp.diff[year]) != sign(temp.diff[year + 1])) {
       is_oscillating <- TRUE #checks for sign differences between temp diffs over years
       break
      }
    }
    
    asymptote.val <- 0
    
    if (is_asymptotic & !is_oscillating) {
      asymptote.val <- round(result$T.s[length(result$T.s)-1], 3)
    }

    if (is_asymptotic & is_oscillating) {
      asymptote.val <- round(mean(result$T.s[length(result$T.s)-1], result$T.s[length(result$T.s)-2]), 3)
    }

    if (!is_asymptotic) {
      asymptote.val <- NA
    }
    
    results.vals[j, i + 1] <- asymptote.val
    results.oscillates[j, i + 1] <- is_oscillating
    
    is.valid <- TRUE
    alpha_out_of_range <- any(result$alpha < 0 | result$alpha > 1)
    epsilon_out_of_range <- any(result$epsilon < 0 | result$epsilon > 1)
  
    if (alpha_out_of_range | epsilon_out_of_range) {
      is.valid <- FALSE
    }
  
    results.valid[j, i + 1] <- is.valid
  }
}
library(knitr)
kable(results.vals)
kable(results.oscillates)
kable(results.valid)

Model used for Question 8

# 1. Start a function that takes the model parameters that might change as arguments i.e albedo and emissivity

ClimateChangeModelB<- function(Propn.ocean, No.Years, StartingT.s, T.a0, alpha.c, alpha.o, epsilon.ac, epsilon.ao, 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
      # T.a0 = 240
      # epsilon.ac = 0.01 # emissivity variation coefficient, in K^-1
      # epsilon.ao = 0.78
      # alpha.c = 0.01 # albedo variation coefficient, in K^-1
      # alpha.o = 0.3
      
      # 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))
      alpha = numeric(length(Year))
      epsilon = 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)
        # calculate alpha for the ith timestep
        alpha[i] = alpha.c * (T.a[i] - T.a0) + alpha.o
        # print(alpha[i])
        # calculate epsilon for the ith timestep
        epsilon[i] = epsilon.ac * (T.s[i] - StartingT.s) + epsilon.ao
        # print(epsilon[i])
        # calculate rate out for the ith timestep
        Rate.out[i] = alpha[i]*S/4 + (1-epsilon[i])*sigma*T.s[i]^4 + epsilon[i]*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, alpha, epsilon)
      # print(round((ClimateChange),5))
      # 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)
}

Code for Step 6B

alpha.c_params <- c(0,0.01,0.04,0.06,0.08,0.10,0.12)
epsilon.ac_params <- c(0.03,0.01,0,-0.01,-0.04,-0.10,-0.16)

results.vals <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.vals) <- c("alpha.c", as.character(epsilon.ac_params))
results.vals$alpha.c <- alpha.c_params

results.oscillates <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.oscillates) <- c("alpha.c", as.character(epsilon.ac_params))
results.oscillates$alpha.c <- alpha.c_params

results.valid <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.valid) <- c("alpha.c", as.character(epsilon.ac_params))
results.valid$alpha.c <- alpha.c_params

for (i in 1:length(epsilon.ac_params)) {
  for (j in 1:length(alpha.c_params)) {
    result <- ClimateChangeModelB(
      Propn.ocean = 0.7, 
      No.Years = 20, 
      StartingT.s = 288, 
      T.a0 = 240, 
      alpha.c = alpha.c_params[j], 
      alpha.o = 0.3, 
      epsilon.ac = epsilon.ac_params[i], 
      epsilon.ao = 0.78,
      plot = FALSE
    )
    
    temp.diff <- result$DeltaT
    is_asymptotic <- TRUE
    is_oscillating <- FALSE
    
        if (abs(temp.diff[3]) >= abs(temp.diff[2])) {
          is_asymptotic <- FALSE
        }
    
    
    for (year in 2:length(temp.diff) - 1) { 
    if (!is.na(temp.diff[year]) && !is.na(temp.diff[year + 1]) && sign(temp.diff[year]) != sign(temp.diff[year + 1])) {
       is_oscillating <- TRUE #checks for sign differences between temp diffs over years
       break
      }
    }
    
    asymptote.val <- 0
    
    if (is_asymptotic & !is_oscillating) {
      asymptote.val <- round(result$T.s[length(result$T.s)-1], 3)
    }

    if (is_asymptotic & is_oscillating) {
      asymptote.val <- round(mean(result$T.s[length(result$T.s)-1], result$T.s[length(result$T.s)-2]), 3)
    }

    if (!is_asymptotic) {
      asymptote.val <- NA
    }
    
    results.vals[j, i + 1] <- asymptote.val
    results.oscillates[j, i + 1] <- is_oscillating
    
    is.valid <- TRUE
    alpha_out_of_range <- any(result$alpha < 0 | result$alpha > 1)
    epsilon_out_of_range <- any(result$epsilon < 0 | result$epsilon > 1)
  
    if (alpha_out_of_range | epsilon_out_of_range) {
      is.valid <- FALSE
    }
  
    results.valid[j, i + 1] <- is.valid
  }
}
library(knitr)
kable(results.vals)
kable(results.oscillates)
kable(results.valid)

Model used for Question 9

# 1. Start a function that takes the model parameters that might change as arguments i.e albedo and emissivity

ClimateChangeModelC<- function(Propn.ocean, No.Years, StartingT.s, T.a0, alpha.c, alpha.o, epsilon.ao, 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
      # T.a0 = 240
      # epsilon.ac = 0.01 # emissivity variation coefficient, in K^-1
      # epsilon.ao = 0.78
      # alpha.c = 0.01 # albedo variation coefficient, in K^-1
      # alpha.o = 0.3
      
      # 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))
      alpha = numeric(length(Year))
      epsilon = numeric(length(Year))
      T.s[1] = StartingT.s                     # Set initial T.s
      T.a = numeric(length(Year))
      k = epsilon.ao / (6.1094 * exp((17.625 * (T.a0 - 273.15)) / (243.04 + (T.a0 - 273.15))))
      # 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)
        # calculate alpha for the ith timestep
        alpha[i] = alpha.c * (T.a[i] - T.a0) + alpha.o
        # print(alpha[i])
        # calculate epsilon for the ith timestep
        epsilon[i] = k * 6.1094 * exp((17.625 * (T.a[i] - 273.15)) / (243.04 + (T.a[i] - 273.15)))
        # print(epsilon[i])
        # calculate rate out for the ith timestep
        Rate.out[i] = alpha[i]*S/4 + (1-epsilon[i])*sigma*T.s[i]^4 + epsilon[i]*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, alpha, epsilon)
      # print(round((ClimateChange),5))
      # 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)
}

Code for Step 6C

alpha.c_params <- c(0.05,0.07,0.09,0.11,0.13,0.15)
epsilon.ac_params <- c(0)

results.vals <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.vals) <- c("alpha.c", as.character(epsilon.ac_params))
results.vals$alpha.c <- alpha.c_params

results.oscillates <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.oscillates) <- c("alpha.c", as.character(epsilon.ac_params))
results.oscillates$alpha.c <- alpha.c_params

results.valid <- data.frame(matrix(ncol = length(epsilon.ac_params) + 1, nrow = length(alpha.c_params)))
names(results.valid) <- c("alpha.c", as.character(epsilon.ac_params))
results.valid$alpha.c <- alpha.c_params

for (i in 1:length(epsilon.ac_params)) {
  for (j in 1:length(alpha.c_params)) {
    result <- ClimateChangeModelC(
      Propn.ocean = 0.7, 
      No.Years = 20, 
      StartingT.s = 288, 
      T.a0 = 240, 
      alpha.c = alpha.c_params[j], 
      alpha.o = 0.3, 
      epsilon.ao = 0.78,
      plot = FALSE
    )
    
    temp.diff <- result$DeltaT
    is_asymptotic <- TRUE
    is_oscillating <- FALSE
    
        if (abs(temp.diff[3]) >= abs(temp.diff[2])) {
          is_asymptotic <- FALSE
        }
    
    
    for (year in 2:length(temp.diff) - 1) { 
    if (!is.na(temp.diff[year]) && !is.na(temp.diff[year + 1]) && sign(temp.diff[year]) != sign(temp.diff[year + 1])) {
       is_oscillating <- TRUE #checks for sign differences between temp diffs over years
       break
      }
    }
    
    asymptote.val <- 0
    
    if (is_asymptotic & !is_oscillating) {
      asymptote.val <- round(result$T.s[length(result$T.s)-1], 3)
    }

    if (is_asymptotic & is_oscillating) {
      asymptote.val <- round(mean(result$T.s[length(result$T.s)-1], result$T.s[length(result$T.s)-2]), 3)
    }

    if (!is_asymptotic) {
      asymptote.val <- NA
    }
    
    results.vals[j, i + 1] <- asymptote.val
    results.oscillates[j, i + 1] <- is_oscillating
    
    is.valid <- TRUE
    alpha_out_of_range <- any(result$alpha < 0 | result$alpha > 1)
    epsilon_out_of_range <- any(result$epsilon < 0 | result$epsilon > 1)
  
    if (alpha_out_of_range | epsilon_out_of_range) {
      is.valid <- FALSE
    }
  
    results.valid[j, i + 1] <- is.valid
  }
}
library(knitr)
kable(results.vals)
kable(results.oscillates)
kable(results.valid)