#R script for statistical analysis of the article "Local ship speed reduction effect on black carbon emissions measured at remote marine station"

#import packages used in the analysis
library(ggplot2)
library(ggpubr)
library(cowplot)
library(data.table)
library(dplyr)

#Define constants
#Fuel carbon content of different marine fuels
FCC_HFO = 0.8493
FCC_MGO = 0.8744

#Molar masses of C and CO2
MCO2 = 44.009
MC = 12.011

#Import dataset for sample of engine loads
shipload_df <- read.csv("shipload_df.csv", sep=";")

#Plot
pdf("Figure_3.pdf",width = 6.7,height = 3.9)
ggplot(data = shipload_df, aes(Speed,Load))+
    geom_line(aes(color=Ship),linewidth=1.1)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(x="Vessel speed (kts)",y="Main engine load")
dev.off()

#Import dataset with plume data
uto_dataset <- read.csv("uto_dataset.csv", sep=";")
uto_dataset$EGCS = ifelse(uto_dataset$scrubber == 1,"YES","NO")

#Calculate EFbc (g BC / kg fuel)
uto_dataset$EFbc = ifelse(uto_dataset$scrubber == 1,(uto_dataset$areaeBC / (uto_dataset$areaCO2/1000)) * (MCO2/MC) * FCC_HFO,(uto_dataset$areaeBC / (uto_dataset$areaCO2/1000)) * (MCO2/MC) * FCC_MGO)

#Calculate mean and SD (g BC / kg fuel) of all dataset
mean(uto_dataset$EFbc)
sd(uto_dataset$EFbc)

#Calculate loading factor to define ships in ballast and laden condition
uto_dataset$lfactor = uto_dataset$draught / uto_dataset$Ddraught
uto_dataset$Condition = ifelse(uto_dataset$lfactor < 0.9,"BALLAST","LADEN")
uto_Ballast = subset(uto_dataset,uto_dataset$Condition == "BALLAST")
uto_Laden = subset(uto_dataset,uto_dataset$Condition == "LADEN")

#Plot Figure 4
pdf("Figure_4.pdf",width = 6.7,height = 4.3)
ggplot(data = uto_dataset,aes(draught,Ddraught))+
    scale_shape_manual(values=c(9:21))+
    geom_point(aes(shape=shipType,color=Condition),size=4)+
    xlim(2,11)+ylim(2,11)+
    labs(x="Actual draught (AIS)",y="Design draught (ship database)")+
    geom_abline(slope = 1,intercept = 0,linetype="dashed")+
    theme_minimal_grid()+
    guides(color=guide_legend(nrow=2,byrow = TRUE),shape=guide_legend(nrow=2,byrow = TRUE))+
    theme(legend.position = "top",
          text = element_text(size=12),legend.title = element_blank())
dev.off()

#Define categories for EGCS and service speed
uto_EGCS = subset(uto_dataset,uto_dataset$EGCS == "YES")
uto_noEGCS = subset(uto_dataset,uto_dataset$EGCS == "NO")
uto_fast = subset(uto_dataset,uto_dataset$serviceSpeed > 15)
uto_slow = subset(uto_dataset,uto_dataset$serviceSpeed >= 15)
uto_dataset$Category = ifelse(uto_dataset$serviceSpeed > 15,"FAST","SLOW")

#Make boxplot for EGCS category
EGCS = ggplot(data=uto_dataset,aes(EGCS,log10(EFbc),fill=EGCS))+
    geom_jitter(aes(color=EGCS))+
    geom_boxplot(outlier.shape = NA)+
    stat_summary(fun="mean",geom="point",shape=17,size=3,show.legend = FALSE)+
    ylab(bquote(log[10]~(g~BC~g~fuel^-1)))+
    theme_minimal_grid()+
    theme(legend.position = "top",
          axis.title.x = element_blank(),
          axis.text.x = element_blank())

#Calculate mean and SD for EGCS category and perform Welch test
mean(uto_EGCS$EFbc)
sd(uto_EGCS$EFbc)
mean(uto_noEGCS$EFbc)
sd(uto_noEGCS$EFbc)
t.test(uto_EGCS$EFbc,uto_noEGCS$EFbc)

#Make boxplot for loading condition category
Condition = ggplot(data=uto_dataset,aes(Condition,log10(EFbc),fill=Condition))+
    geom_jitter(aes(color=Condition))+
    geom_boxplot(outlier.shape = NA)+
    stat_summary(fun="mean",geom="point",shape=17,size=3,show.legend = FALSE)+
    ylab(bquote(log[10]~(g~BC~g~fuel^-1)))+
    theme_minimal_grid()+
    theme(legend.position = "top",
          legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.text.x = element_blank())

#Calculate mean and SD for loading condition category and perform Welch test
mean(uto_Ballast$EFbc)
sd(uto_Ballast$EFbc)
mean(uto_Laden$EFbc)
sd(uto_Laden$EFbc)
t.test(uto_Ballast$EFbc,uto_Laden$EFbc)

#Make boxplot for service speed category
Speed = ggplot(data=uto_dataset,aes(Category,log10(EFbc),fill=Category))+
    geom_jitter(aes(color=Category))+
    geom_boxplot(outlier.shape = NA)+
    stat_summary(fun="mean",geom="point",shape=17,size=3,show.legend = FALSE)+
    ylab(bquote(log[10]~(g~BC~g~fuel^-1)))+
    theme_minimal_grid()+
    theme(legend.position = "top",
          legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.text.x = element_blank())

uto_fast = subset(uto_dataset,uto_dataset$Category =="FAST")
uto_slow = subset(uto_dataset,uto_dataset$Category =="SLOW")

#Calculate mean and SD for service speed category and perform Welch test
mean(uto_fast$EFbc)
sd(uto_fast$EFbc)
mean(uto_slow$EFbc)
sd(uto_slow$EFbc)
t.test(uto_fast$EFbc,uto_slow$EFbc)

Type = ggplot(data=uto_dataset)+
    geom_jitter(aes(shipType,log10(EFbc),color=shipType))+
    geom_boxplot(aes(shipType,log10(EFbc),fill=shipType),outlier.shape = NA,show.legend = FALSE)+
    stat_summary(aes(shipType,log10(EFbc)),fun="mean",geom="point",shape=17,size=3,show.legend = FALSE)+
    ylab(bquote(log[10]~(g~BC~g~fuel^-1)))+
    theme_minimal_grid()+
    theme(legend.text = element_text(size=8),
          legend.justification = "left",
          axis.title.y = element_text(size = 14),
          axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          legend.title = element_blank(),
          legend.position = "top",
          legend.spacing.x = unit(0.1,"cm"))+
    guides(color=guide_legend(nrow = 2,byrow = TRUE))+
    scale_color_discrete(labels = c("BUL","FIS","GEN","OTH","CRU","ROP","ROR","TAC","TAP","TUG"))

pdf("Figure_5.pdf",width = 6.7,height = 5.9)
plot_grid(EGCS,Condition,Type,Speed)
dev.off()

#Calculate means for each vessel type
DT = data.table(uto_dataset[c(15,34)])
means_df = data.frame(DT[,sapply(.SD,function(x) list(mean=mean(x),sd=sd(x))),by=shipType])
means_df = means_df[order(means_df$shipType),]
names(means_df)[2] = "EFbc"
names(means_df)[3] = "SD"

#Correlation tests for speed
cor.test(uto_dataset$EFbc,uto_dataset$shipSpeed)
cor.test(uto_EGCS$EFbc,uto_EGCS$shipSpeed)
cor.test(uto_noEGCS$EFbc,uto_noEGCS$shipSpeed)
cor.test(uto_fast$EFbc,uto_fast$shipSpeed)
cor.test(uto_slow$EFbc,uto_slow$shipSpeed)
cor.test(uto_Ballast$EFbc,uto_Ballast$shipSpeed)
cor.test(uto_Laden$EFbc,uto_Laden$shipSpeed)

#Correlation tests for engine load
cor.test(uto_EGCS$load,uto_EGCS$EFbc)
cor.test(uto_EGCS$Sload,uto_EGCS$EFbc)

#Calculate linear regressions for unadjusted and weather adjusted loads
EGCS_lin_load = lm(data = uto_EGCS, formula = EFbc~load)
summary(EGCS_lin_load)

EGCS_lin_Sload = lm(data = uto_EGCS, formula = EFbc~Sload)
summary(EGCS_lin_Sload)


#Analyse EGCS and no EGCS vessels separately
#Building models
uto_EGCS$logEFbc = log10(uto_EGCS$EFbc)

EGCS_model_log_poly = lm(data = uto_EGCS,formula = logEFbc~load+I(load^2))
summary(EGCS_model_log_poly)

EGCS_model_log_poly = lm(data = uto_EGCS,formula = logEFbc~Sload+I(Sload^2))
summary(EGCS_model_log_poly)

colcruise = "#1B9E77"
colropax = "#D95F02"
colroro1 = "#7570B3"
colroro2 = "#E7298A"
coltanker = "#66A61E"

ship_colors = c("#1B9E77","#D95F02","#7570B3","#E7298A","#66A61E")

#Plot obtained models
EGCS_load = ggplot(data=uto_EGCS,aes(Sload,logEFbc))+
    geom_point(aes(shape=shipType,colour=shipType),size=4)+
    stat_smooth(method = lm,formula = y~x+I(x^2),color="red",linetype="dashed")+
    stat_regline_equation(size=3.5,label.x = 0.25,formula = y~x+I(x^2))+
    stat_regline_equation(size=3.5,label.x = 0.59,label.y = -0.1,formula = y~poly(x,2),aes(label=..adj.rr.label..))+
    ylab(bquote(log[10]~(g~BC~kg~fuel^-1)))+
    xlab("Weather adjusted main engine load")+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=10))+
    scale_shape_manual(name="",
                       labels=c("CRUISE","ROPAX","RORO"),
                       values = c(15:17))+
    scale_colour_manual(name="",
                        labels=c("CRUISE","ROPAX","RORO"),
                        values = c(colcruise,colropax,colroro1))

Sload = seq(0,0.8,by=0.01)
new_df = data.frame(Sload)
prediction = predict(EGCS_model_log_poly,newdata = new_df,interval = "prediction")
new_df[c("pred","lwr","upr")] = prediction

EGCS_load_lin = ggplot(data=new_df)+
    geom_ribbon(aes(x=Sload,ymax=10^upr,ymin=10^lwr),alpha=0.2)+
    geom_point(data=uto_EGCS,aes(load,EFbc,shape=shipType,color=shipType,group=shipType),size=3)+
    geom_line(aes(Sload,10^pred),linetype="dashed",linewidth=1,color="red")+
    theme_minimal_grid()+
    stat_cor(label.x = 0.3,data=uto_EGCS,aes(load,EFbc),size=3.5,digits = 3)+
    theme(legend.position = "top",
          text = element_text(size=10))+
    ylab(bquote(g~BC~kg~fuel^-1))+
    xlab("Weather adjusted main engine load")+
    scale_shape_manual(name="",
                       labels=c("CRUISE","ROPAX","RORO"),
                       values = c(15:17))+
    scale_colour_manual(name="",
                        labels=c("CRUISE","ROPAX","RORO"),
                        values = c(colcruise,colropax,colroro1))

pdf("Figure_6.pdf",width = 6.7,height = 5.12)
plot_grid(EGCS_load,EGCS_load_lin)
dev.off()

#Model BC output on 5 different vessels
roro1 <- read.csv("roro1.csv", sep=";")
roro1 = subset(roro1,roro1$Load < 0.8)
roro1$fuel = roro1$fuelcons / 1000 * roro1$Power
roro1$CO2 = roro1$fuel * 3114
roro1$CO2_nm = roro1$CO2 / roro1$Speed
roro1$BCprod = 10^(2.10498*roro1$Load^2 -2.84832*roro1$Load -0.02357)
roro1$BC = roro1$fuel * roro1$BCprod
roro1$BC_nm = roro1$BC / roro1$Speed
roro1$BC_GWP100 = roro1$BC * 680
roro1$BC_GWP100_nm = roro1$BC_GWP100 / roro1$Speed
roro1$BC_GWP20 = roro1$BC * 2200
roro1$BC_GWP20_nm = roro1$BC_GWP20 / roro1$Speed
roro1$CO2e_100 = roro1$CO2 + roro1$BC_GWP100
roro1$CO2e_20 = roro1$CO2 + roro1$BC_GWP20
roro1$CO2e_100_nm = roro1$CO2e_100 / roro1$Speed
roro1$CO2e_20_nm = roro1$CO2e_20 / roro1$Speed

ropax <- read.csv("ropax.csv", sep=";")
ropax = subset(ropax,ropax$Load < 0.8)
ropax$fuel = ropax$fuelcons / 1000 * ropax$Power #kg/kWh * kW = kg/h
ropax$CO2 = ropax$fuel * 3114 #kg/h /1000 * gCO2 / g fuel = gCO2 / h
ropax$CO2_nm = ropax$CO2 / ropax$Speed
ropax$BCprod = 10^(2.10498*ropax$Load^2 -2.84832*ropax$Load -0.02357) #g BC / g fuel
ropax$BC = ropax$fuel * ropax$BCprod # kg fuel/h * g BC/kg fuel = g BC/h
ropax$BC_nm = ropax$BC / ropax$Speed
ropax$BC_GWP100 = ropax$BC * 680 # g CO2e / h
ropax$BC_GWP100_nm = ropax$BC_GWP100 / ropax$Speed
ropax$BC_GWP20 = ropax$BC * 2200 # g CO2e / h
ropax$BC_GWP20_nm = ropax$BC_GWP20 / ropax$Speed
ropax$CO2e_100 = ropax$CO2 + ropax$BC_GWP100
ropax$CO2e_20 = ropax$CO2 + ropax$BC_GWP20
ropax$CO2e_100_nm = ropax$CO2e_100 / ropax$Speed
ropax$CO2e_20_nm = ropax$CO2e_20 / ropax$Speed

roro2 <- read.csv("roro2.csv", sep=";")
roro2 = subset(roro2,roro2$Load < 0.8)
roro2$fuel = roro2$fuelcons / 1000 * roro2$Power #kg/kWh * kW = kg/h
roro2$CO2 = roro2$fuel * 3114 #kg/h /1000 * gCO2 / g fuel = gCO2 / h
roro2$CO2_nm = roro2$CO2 / roro2$Speed
roro2$BCprod = 10^(2.10498*roro2$Load^2 -2.84832*roro2$Load -0.02357) #g BC / g fuel
roro2$BC = roro2$fuel * roro2$BCprod # kg fuel/h * g BC/kg fuel = g BC/h
roro2$BC_nm = roro2$BC / roro2$Speed
roro2$BC_GWP100 = roro2$BC * 680 # g CO2e / h
roro2$BC_GWP100_nm = roro2$BC_GWP100 / roro2$Speed
roro2$BC_GWP20 = roro2$BC * 2200 # g CO2e / h
roro2$BC_GWP20_nm = roro2$BC_GWP20 / roro2$Speed
roro2$CO2e_100 = roro2$CO2 + roro2$BC_GWP100
roro2$CO2e_20 = roro2$CO2 + roro2$BC_GWP20
roro2$CO2e_100_nm = roro2$CO2e_100 / roro2$Speed
roro2$CO2e_20_nm = roro2$CO2e_20 / roro2$Speed

cruise <- read.csv("cruise.csv", sep=";")
cruise = subset(cruise,cruise$Load < 0.8)
cruise$fuel = cruise$fuelcons / 1000 * cruise$Power #kg/kWh * kW = kg/h
cruise$CO2 = cruise$fuel * 3114 #kg/h /1000 * gCO2 / g fuel = gCO2 / h
cruise$CO2_nm = cruise$CO2 / cruise$Speed
cruise$BCprod = 10^(2.10498*cruise$Load^2 -2.84832*cruise$Load -0.02357) #g BC / g fuel
cruise$BC = cruise$fuel * cruise$BCprod # kg fuel/h * g BC/kg fuel = g BC/h
cruise$BC_nm = cruise$BC / cruise$Speed
cruise$BC_GWP100 = cruise$BC * 680 # g CO2e / h
cruise$BC_GWP100_nm = cruise$BC_GWP100 / cruise$Speed
cruise$BC_GWP20 = cruise$BC * 2200 # g CO2e / h
cruise$BC_GWP20_nm = cruise$BC_GWP20 / cruise$Speed
cruise$CO2e_100 = cruise$CO2 + cruise$BC_GWP100
cruise$CO2e_20 = cruise$CO2 + cruise$BC_GWP20
cruise$CO2e_100_nm = cruise$CO2e_100 / cruise$Speed
cruise$CO2e_20_nm = cruise$CO2e_20 / cruise$Speed

tanker <- read.csv("tanker.csv", sep=";")
tanker = subset(tanker,tanker$Load < 0.8)
tanker = subset(tanker,tanker$Speed<21.1)
tanker$fuel = tanker$fuelcons / 1000 * tanker$Power #kg/kWh * kW = kg/h
tanker$CO2 = tanker$fuel * 3114 #kg/h /1000 * gCO2 / g fuel = gCO2 / h
tanker$CO2_nm = tanker$CO2 / tanker$Speed
tanker$BCprod = 10^(2.10498*tanker$Load^2 -2.84832*tanker$Load -0.02357) #g BC / g fuel
tanker$BC = tanker$fuel * tanker$BCprod # kg fuel/h * g BC/kg fuel = g BC/h
tanker$BC_nm = tanker$BC / tanker$Speed
tanker$BC_GWP100 = tanker$BC * 680 # g CO2e / h
tanker$BC_GWP100_nm = tanker$BC_GWP100 / tanker$Speed
tanker$BC_GWP20 = tanker$BC * 2200 # g CO2e / h
tanker$BC_GWP20_nm = tanker$BC_GWP20 / tanker$Speed
tanker$CO2e_100 = tanker$CO2 + tanker$BC_GWP100
tanker$CO2e_20 = tanker$CO2 + tanker$BC_GWP20
tanker$CO2e_100_nm = tanker$CO2e_100 / tanker$Speed
tanker$CO2e_20_nm = tanker$CO2e_20 / tanker$Speed

model_df = data.frame(Speed = c(roro1$Speed,roro2$Speed,ropax$Speed,cruise$Speed,tanker$Speed),
                      CO2 = c(roro1$CO2,roro2$CO2,ropax$CO2,cruise$CO2,tanker$CO2),
                      CO2_nm = c(roro1$CO2_nm,roro2$CO2_nm,ropax$CO2_nm,cruise$CO2_nm,tanker$CO2_nm),
                      CO2e_100 = c(roro1$CO2e_100,roro2$CO2e_100,ropax$CO2e_100,cruise$CO2e_100,tanker$CO2e_100),
                      CO2e_20 = c(roro1$CO2e_20,roro2$CO2e_20,ropax$CO2e_20,cruise$CO2e_20,tanker$CO2e_20),
                      CO2e_100_nm = c(roro1$CO2e_100_nm,roro2$CO2e_100_nm,ropax$CO2e_100_nm,cruise$CO2e_100_nm,tanker$CO2e_100_nm),
                      CO2e_20_nm = c(roro1$CO2e_20_nm,roro2$CO2e_20_nm,ropax$CO2e_20_nm,cruise$CO2e_20_nm,tanker$CO2e_20_nm),
                      BC = c(roro1$BC,roro2$BC,ropax$BC,cruise$BC,tanker$BC),
                      BC_nm = c(roro1$BC_nm,roro2$BC_nm,ropax$BC_nm,cruise$BC_nm,tanker$BC_nm),
                      BC_GWP100 = c(roro1$BC_GWP100,roro2$BC_GWP100,ropax$BC_GWP100,cruise$BC_GWP100,tanker$BC_GWP100),
                      BC_GWP20 = c(roro1$BC_GWP20,roro2$BC_GWP20,ropax$BC_GWP20,cruise$BC_GWP20,tanker$BC_GWP20),
                      BC_GWP100_nm = c(roro1$BC_GWP100_nm,roro2$BC_GWP100_nm,ropax$BC_GWP100_nm,cruise$BC_GWP100_nm,tanker$BC_GWP100_nm),
                      BC_GWP20_nm = c(roro1$BC_GWP20_nm,roro2$BC_GWP20_nm,ropax$BC_GWP20_nm,cruise$BC_GWP20_nm,tanker$BC_GWP20_nm),
                      Fuel = c(roro1$fuelcons,roro2$fuelcons,ropax$fuelcons,cruise$fuelcons,tanker$fuelcons),
                      Sspeed = c(roro1$dSpeed,roro2$dSpeed,ropax$dSpeed,cruise$dSpeed,tanker$dSpeed),
                      Ship = c(rep("Roro1",length(roro1$Speed)),
                               rep("Roro2",length(roro2$Speed)),
                               rep("Ropax",length(ropax$Speed)),
                               rep("Cruise",length(cruise$Speed)),
                               rep("Tanker",length(tanker$Speed))))

Sspeed_roro1 = 20
Sspeed_ropax = 18.5
Sspeed_roro2 = 20
Sspeed_cruise = 20
Sspeed_tanker = 15.3

Fuel = ggplot(data = model_df,aes(Speed,Fuel,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(Fuel~consumption~g~kWh^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

BC = ggplot(data = model_df,aes(Speed,BC,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(BC~g~h^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

BC_nm = ggplot(data = model_df,aes(Speed,BC_nm,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(BC~g~NM^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

CO2e_100 = ggplot(data = model_df,aes(Speed,CO2e_100,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(GWP[100]~CO[2]~e~g~h^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

CO2e_100_nm = ggplot(data = model_df,aes(Speed,CO2e_100_nm,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(GWP[100]~CO[2]~e~g~NM^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

CO2e_20 = ggplot(data = model_df,aes(Speed,CO2e_20,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(GWP[20]~CO[2]~e~g~h^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

CO2e_20_nm = ggplot(data = model_df,aes(Speed,CO2e_20_nm,color=Ship))+
    scale_color_manual(values = ship_colors)+
    geom_line(linewidth=0.8)+
    theme_minimal_grid()+
    theme(legend.position = "top",
          text = element_text(size=12))+
    labs(y=bquote(GWP[20]~CO[2]~e~g~NM^-1))+
    geom_vline(xintercept = Sspeed_roro1,color=colroro1,linetype="dashed")+
    geom_vline(xintercept = Sspeed_ropax,color=colropax,linetype="dashed")+
    geom_vline(xintercept = Sspeed_tanker,color=coltanker,linetype="dashed")+
    guides(color=guide_legend(title = "",nrow = 2,byrow = TRUE))

pdf("Figure_7.pdf",width = 6.7,height = 5.5)
plot_grid(BC,BC_nm,CO2e_100,CO2e_20)
dev.off()