class: center, middle, inverse, title-slide # Climate Modelling ### Pietro Monticone & Davide Orsenigo ### 2019-06-14 | Turin University --- # Contents - Introduction - Packages - Input Parameters - Fit from Textbook Spreadsheet - Functions - Run 0 - Run 1 - Run 2 - Run 3 & 4 - Hysteresis Cycles - Embedding Daisies in EBM # Packages ```r library(tidyverse) library(plotly) library(scales) ``` --- # Input parameters ```r S <- 1370 A <- 204 B <- 2.17 K <- 3.86 ai <- 0.62 ab <- 0.25 Temperature <- rep(0,300) Temp <- rep(0,100) Temp2 <- rep(0,180) Temp5 <- rep(0,300) J <- rep(0,100) Albedo <- rep(0,300) t <- c(1:300) Zones <- seq(-89, 89, by = 2) ``` --- ## Fit from Textbook Spreadsheet .small[ ```r zones <- c(85,75,65,55,45,35,25,15,6) coszones <- cospi(zones/180) sunWt <- c(0.5,0.531,0.624,0.77,0.892,1.021,1.12,1.189,1.219) df <- data.frame(zones,sunWt) f <- function(zones,d,n,k){ d*cos(n*zones)^2+k } Fit <- nls(sunWt ~ f(zones,d,n,k),data= df, start=list(d=1,n=0.015,k=0.3)) ggplot(df, aes(zones,sunWt))+geom_point()+geom_smooth(aes(y=f(zones,0.7768699,0.0164348,0.4617747)))+xlab("Latitude")+ylab("Solar Weight Factor") ``` <img src="ClimateModelling_slides_files/figure-html/unnamed-chunk-3-1.png" width="50%" /> ] --- # Functions .pull-left[ .small[ ```r Func <- function(x){0.7768699*cos(0.0164348*x)^2+0.4617747} gauss <- function(x,m,sd,b){ ((24+b)/(0.00798*sqrt(2*pi*sd^2)))*exp(-(x-m)^2/(2*sd^2))-b} Sun1 <- function(x,a){a*x/100} Sun2 <- function(x){1370*(sinpi((x+90)/180))^2} Sun3 <- function(x){ 1370-(((1370)/(sqrt(2*pi*1^2)))*exp(-(x-50)^2/(2*1^2)))} Sun4 <- function(x){ifelse(x==50,1370/3,1370)} Sun5 <- function(x){(1/100) * (abs(x-150)+25)} Sun6 <- function(x){1370*(1+0.1*cospi(x/180))} Incident <- function(x,y){x*y/4} Step <- function(x,c){ifelse(x<c, 0.6, 0.3)} alb <- function(x,a,b){(-exp(2.2*x+10)/(exp(2.2*x+10)+1))*(a-b)+a} ``` ] ] $$ S_1(t)= \frac{1370}{100}x $$ $$ S_2(t)= 1370(\sin^2(x+90)) $$ $$ S_3(t)= 1370\Big(1-\frac{1}{\sqrt{2\pi}}e^{-(x-50)^2/2}\Big) $$ $$ S_4(t)=1370\Big(1-\frac{1}{3}\delta(t-50)\Big) $$ $$ S_5(t)=\frac{1}{100}(|x-150|+25) $$ $$S_6(t)=1370\Big(1+\frac{1}{10} cos(x)\Big) $$ --- $$ a(T)=\frac{e^{2.2T+10}}{e^{2.2T+10}+1}(b-a)+a $$ ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-5-1.png)<!-- --> --- # Introduction ```r Ti <- gauss(Zones,0,50,31.6) cosZones <- abs(cospi(Zones/180)) SunWt <- Func(Zones) Rin <- Incident(S,SunWt) T <- gauss(Zones,0,50,31.6) a <- alb(T,ai,ab) ``` -- # Run 0 ```r for(i in c(1:300)) { Tcos <- cosZones*T Tm <- sum(Tcos)/sum(cosZones) T <- (Rin*(1-a)+K*Tm-A) / (B+K) a <- alb(T,0.62,0.25) Temperature[i] <- Tm } ``` --- .small[ ```r data1 <- data.frame(Zones,T,a,Ti) plot1 <- ggplot(data1,aes(Zones)) + geom_line(aes(y=T, colour = "Temperature",linetype="Temperature")) + geom_line(aes(y=a*25, colour = "Albedo",linetype="Albedo"))+ geom_line(aes(y=Ti, colour = "Initial Temperature",linetype="Initial Temperature"))+xlab("Latitude")+ylab("Temperature")+scale_colour_manual(name="Legend",values=c("blue","red","red"))+scale_linetype_manual(name="Legend",values=c("Temperature"=1, "Albedo"=1, "Initial Temperature"=2)) ``` ] <img src="ClimateModelling_slides_files/figure-html/unnamed-chunk-9-1.png" width="65%" /> --- .small[ ```r data2 <- data.frame(t,Temperature) plot2 <- ggplot(data2) + geom_line(aes(t, Temperature),colour = 'green')+xlab("Time")+ylab("Mean Temperature") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-11-1.png)<!-- --> --- # Run 1 ```r TEMP1 <- matrix(NA, nrow=90, ncol=100) for(j in c(1:100)){ T <- gauss(Zones,0,50,31.6) a <- alb(T,ai,ab) S <- Sun1(1370,j) Rin <- Incident(S,SunWt) for(i in c(1:300)) {Tcos <- cosZones*T Tm <- sum(Tcos)/sum(cosZones) T <- (Rin*(1-a)+K*Tm-A) / (B+K) a <- alb(T,ai,ab)} TEMP1[,j] <- T J[j] <- j Temp[j] <- Tm } ``` --- .small[ ```r data3 <- data.frame(J,Temp) plot3 <- ggplot(data3,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Non Dynamic Mean Temperature (T(t) independent from T(t-1))") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-14-1.png)<!-- --> --- .small[ ```r plot_ly(z=~TEMP1)%>% add_surface() %>% layout( title = "With Initialization", scene = list( xaxis = list(title = "S/100"), yaxis = list(title = "Latitude"), zaxis = list(title = "T") )) ```
] --- # Run 2 ```r T2 <- gauss(Zones,0,50,31.6) TEMP2 <- matrix(NA, nrow=90, ncol=180) a <- alb(T2,ai,ab) Sarr <- rep(0,180) for(j in c(1:180)){ S <- Sun2(j) Rin <- Incident(S,SunWt) for(i in c(1:300)) {Tcos <- cosZones*T2 Tm <- sum(Tcos)/sum(cosZones) T2 <- (Rin*(1-a)+K*Tm-A) / (B+K) a <- alb(T2,ai,ab) } Sarr[j] <- Sun2(j) TEMP2[,j] <- T2 J[j] <- j Temp2[j] <- Tm } ``` --- .small[ ```r data3 <- data.frame(J,Temp2) plot3 <- ggplot(data3,aes(J,Temp2))+geom_line(aes(J, Temp2),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Sinusoidal Perturbation") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-18-1.png)<!-- --> --- .small[ ```r plot_ly(z=~TEMP2) %>% add_surface() %>% layout( title = "Without Initialization", scene = list( xaxis = list(title = "S/100"), yaxis = list(title = "Latitude"), zaxis = list(title = "T") )) ```
] --- # Run 3 & 4 ```r TEMP3 <- matrix(NA, nrow=90, ncol=200) T3 <- gauss(Zones,0,50,31.6) a <- alb(T3,ai,ab) for(j in c(1:200)){ S <- Sun3(j) Rin <- Incident(S,SunWt) for(i in c(1:300)) {Tcos <- cosZones*T3 Tm <- sum(Tcos)/sum(cosZones) T3 <- (Rin*(1-a)+K*Tm-A) / (B+K) a <- alb(T3,ai,ab) } TEMP3[,j] <- T3 J[j] <- j Temp[j] <- Tm } ``` --- .small[ ```r data3 <- data.frame(J,Temp) plot3 <- ggplot(data3,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Gaussian Perturbation") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-22-1.png)<!-- --> --- .small[ ```r plot_ly(z=~TEMP3) %>% add_surface()%>% layout( title = "Without Initialization", scene = list( xaxis = list(title = "S/100"), yaxis = list(title = "Latitude"), zaxis = list(title = "T") )) ```
] --- ```r TEMP4 <- matrix(NA, nrow=90, ncol=200) T4 <- gauss(Zones,0,50,31.6) a <- alb(T4,ai,ab) for(j in c(1:200)){ S <- Sun4(j) Rin <- Incident(S,SunWt) for(i in c(1:300)) {Tcos <- cosZones*T4 Tm <- sum(Tcos)/sum(cosZones) T4 <- (Rin*(1-a)+K*Tm-A) / (B+K) a <- alb(T4,ai,ab) } TEMP4[,j] <- T4 J[j] <- j Temp[j] <- Tm } ``` --- .small[ ```r data3 <- data.frame(J,Temp) plot3 <- ggplot(data3,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Delta Perturbation") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-26-1.png)<!-- --> --- .small[ ```r plot_ly(z=~TEMP4) %>% add_surface()%>% layout( title = "Without Initialization", scene = list( xaxis = list(title = "S/100"), yaxis = list(title = "Latitude"), zaxis = list(title = "T") )) ```
] --- # Hysteresis Cycles --- .small[ ```r data <- data.frame(Sarr,Temp2) plot <- ggplot(data,aes(Sarr,Temp2))+geom_point(aes(Sarr, Temp2),colour = 'yellow')+ylab("Mean Temperature")+xlab("S_2(t)")+ggtitle("Mean temperature vs. Sinusoidally Fluctuating Incoming Radiation") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-29-1.png)<!-- --> --- ```r Temp5 <- rep(0,250) T5 <- gauss(Zones,0,50,31.6) a <- alb(T5,ai,ab) Sarr <- rep(0,250) for(j in c(0:250)){ S <- Sun5(j) Rin <- 1370*Incident(S,SunWt) for(i in c(1:60)) {Tcos <- cosZones*T5 Tm <- sum(Tcos)/sum(cosZones) T5 <- (Rin*(1-a)+K*Tm-A) / (B+K) a <- alb(T5,ai,ab) } Sarr[j] <- Sun5(j) J[j] <- j Temp5[j] <- Tm } ``` --- .small[ ```r dataw <- data.frame(Sarr,Temp5) plot5 <- ggplot(dataw,aes(Sarr,Temp5,group=1))+geom_point(aes(Sarr, Temp5),colour = 'yellow')+geom_segment(aes(x = Sarr[89], y = Temp5[89], xend = Sarr[90], yend = Temp5[90]),colour = 'yellow')+geom_segment(aes(x = Sarr[258], y = Temp5[258], xend = Sarr[259], yend = Temp5[259]),colour = 'yellow')+ylab("Mean Temperature")+xlab("S_5(t)")+ggtitle("Mean temperature vs. Linearly Fluctuating Incoming Radiation") ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-32-1.png)<!-- --> --- # Embedding Daisies in EBM --- ## Input Parameters ```r S <- 1370 A <- 204 B <- 2.17 K <- 3.86 ai <- 0.62 ab <- 0.25 aW <- 0.75 aB <- 0.25 c <- 7 k <- 0.003265*0.75 T0 <- 20 D <- 0.3 # Death Rate Zones <- seq(-89, 89, by = 2) ``` --- ## Functions ```r Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 } gauss <- function(x,m,sd,b){ ((24+b)/(0.00798*sqrt(2*pi*sd^2)))*exp(-(x-m)^2/(2*sd^2))-b } Incident <- function(x,y){x*y/4} alb <- function(x,a,b){ (-exp(1.6*x+10)/(exp(1.6*x+10)+1))*(a-b)+a} #naked ``` --- ## Initialization ```r cosZones <- abs(cospi(Zones/180)) SunWt <- Func(Zones) Rin <- Incident(S,SunWt) T <- gauss(Zones,0,50,31.6)-6 w <- rep(0.5,length(Zones)) b <- rep(0.2,length(Zones)) u <- rep(0.3,length(Zones)) a <- w*aW+b*aB+u*alb(T,0.62,0.25) Barr <- rep(0,500) Warr <- rep(0,500) Uarr <- rep(0,500) Tarr <- rep(0,500) I <- rep(0,500) TEMP <- matrix(NA, nrow=90, ncol=500) ``` --- ## EBM without Daisies ```r for(i in c(1:500)) { S <- Sun6(i) # oppure costante S <- 1370 Rin <- Incident(S,SunWt) Tcos <- cosZones*T Tm <- sum(Tcos)/sum(cosZones) T <- (Rin*(1-a)+K*Tm-A) / (B+K) TEMP[,i] <- T a <- alb(T,0.62,0.25) I[i] <- i Tarr[i] <- T[45] } ``` --- .small[ ```r plot <- ggplot(data.frame(Zones,w,b,u,T),aes(Zones))+geom_line(aes(y=w, colour = "% White"))+geom_line(aes(y=b, colour = "% Black"))+geom_line(aes(y=u, colour="% Bare Ground"))+geom_line(aes(y=T/5, colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Latitude")+ggtitle("Without Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red")) ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-38-1.png)<!-- --> --- .small[ ```r plot2 <- ggplot(data.frame(I,Barr,Warr,Uarr,Tarr),aes(I))+geom_line(aes(y=Barr,colour="% Black"))+ geom_line(aes(y=Warr,colour="% White"))+ geom_line(aes(y=Uarr, colour="% Bare Ground"))+geom_line(aes(y=Tarr/25,colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Time")+ggtitle("Without Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red")) ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-40-1.png)<!-- --> --- .small[ ```r plot_ly(z=~TEMP)%>% add_surface() %>% layout(title="Without Daisies", scene = list( xaxis = list(title = "Time"), yaxis = list(title = "Latitude"), zaxis = list(title = "T") )) ```
] --- ## EBM with Daisies ```r a <- w*aW+b*aB+u*alb(T,0.62,0.25) T <- gauss(Zones,0,50,31.6)-6 TEMP <- matrix(NA, nrow=90, ncol=500) for(i in c(1:500)) { S <- Sun6(i) Rin <- Incident(S,SunWt) Tcos <- cosZones*T Tm <- sum(Tcos)/sum(cosZones) T <- (Rin*(1-a)+K*Tm-A) / (B+K) TEMP[,i] <- T Tw <- T+c*(a-aW) Tb <- T+c*(a-aB) Fw <- 1-k*(T0-Tw)^2 Fb <- 1-k*(T0-Tb)^2 for(j in c(1:length(Zones))){ if(Fw[j]<0){Fw[j]=0} if(Fb[j]<0){Fb[j]=0} } w <- w+w*(u*Fw-D) b <- b+b*(u*Fb-D) for(j in c(1:length(Zones))){if(w[j]<0.001){w[j]=0.001} if(b[j]<0.001){b[j]=0.001} } u <- 1-w-b a <- w*aW+b*aB+u*alb(T,0.62,0.25) Barr[i] <- b[45] Warr[i] <- w[45] Uarr[i] <- u[45] I[i] <- i Tarr[i] <- T[45] } ``` --- .small[ ```r plot <- ggplot(data.frame(Zones,w,b,u,T),aes(Zones))+geom_line(aes(y=w, colour = "% White"))+geom_line(aes(y=b, colour = "% Black"))+geom_line(aes(y=u, colour="% Bare Ground"))+geom_line(aes(y=T/5, colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Latitude")+ggtitle("With Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red")) ``` ] ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-44-1.png)<!-- --> --- ```r plot2 <- ggplot(data.frame(I,Barr,Warr,Uarr,Tarr),aes(I))+geom_line(aes(y=Barr,colour="% Black"))+ geom_line(aes(y=Warr,colour="% White"))+ geom_line(aes(y=Uarr, colour="% Bare Ground"))+geom_line(aes(y=Tarr/25,colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Time")+ggtitle("With Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red")) ``` ![](ClimateModelling_slides_files/figure-html/unnamed-chunk-46-1.png)<!-- --> --- .small[ ```r plot_ly(z=~TEMP)%>% add_surface() %>% layout(title="With Daisies", scene = list( xaxis = list(title = "Time"), yaxis = list(title = "Latitude"), zaxis = list(title = "T") )) ```
]