Climate change comments

Information from previous work (SRK climate change and baseline)

Climate change values were estimated based on L-moment methodology, GEV distribution was selected.

Baseline records were cut from 1985 to 2014.

Code

#information from baseline hydrology
#records cut to match historical baseline under climate change (1985 to 2014)
#site.cmip6.db

site.cmip6.db<-arrow::read_parquet(
  file = here::here("outputs","site_records.par"))


#defined for site based on SRK
site.cmip6.db %>% 
  dplyr::select(Index,prcp) %>% 
  tk_zoo() %>% 
  plot(main="Local precipitation 1985 to 2014 based on SRK",
       ylab="mm/day")

Code

#Values considered for the internal vectors
FA.val<-  c(0.9999, 0.9995, 0.999, 0.998, 0.995, 0.99, 0.98, 
    0.96, 0.95, 0.9, 0.8, 0.5708, 0.5, 0.3333, 0.2, 0.1, 
    0.05, 0.04, 0.02, 0.01, 0.005, 0.002, 0.001, 5e-04, 
    1e-04)

#based on equation implemented from Hershfield (WMO) 
daily_24<-1.122


#values estimated with GEV and L-moments.
#prcp_max_rc
prcp_max_rc<-arrow::read_parquet(
  file = here::here("outputs","prcp_max_cmip6_records.par"))



#extend the table for IDF
prcp_max_rc_tbl<-apply(prcp_max_rc,1,FUN=function(x){
  
  data.frame( 
    model=x$model,
    scenario=x$scenario,
    period=x$period,
    days=x$name,
    prob=FA.val,
    proj=x$value,
    hist=x$historical,
    rc=x$rc %>% unlist()
  )
}) %>% 
  do.call(rbind,.) %>% 
  dplyr::filter(prob>=0.5) %>% 
  mutate(rp=(1/(1-prob)) %>% round(1))

As an example, 6 GCM samples were selected. In the figure below the maximum information per year is presented (historical and projected); from the GCM can be seen the ~1:2 is increased as Rx1Day (presented on SRK Climdex Appendix and Stantec bibliographical review); while more extreme events are decreasing.

Code

site_cmip6_eval<-arrow::read_parquet(
  file = here::here("outputs","site_cmip6_prcp.par"))

set.seed(12)

gcm_sample<-sample(site_cmip6_eval$model %>% 
                     unique(),size = 8)

site_cmip6_eval %>% 
  dplyr::filter(model%in%gcm_sample) %>%
  group_by(model,year,scenario,period) %>% 
  summarize(pr=max(pr)) %>% 
  data.frame() %>% 
  ggplot(data=.,aes(x=year,y=pr,fill=scenario))+
  geom_bar(stat = "Identity",position="dodge")+
  facet_grid(model~period,scales="free")+
  theme_bw()+
  labs(title="Maximum daily precipitation per year comparison between GCM historical / projections",
       y="max prcp - year [mm]") +
  theme(legend.position="bottom")

Information from Stantec and estimate SRK values

Based on AIC, AICc, BIC, ADC, the best distribution should be LN, follow by GUMBEL. To be consistent with typical climate change approaches, GEV was selected.

Also GEV is normally selected for climate change analyzes. To observe differences all of them are included.

Code

#defined for the complete period by Stantec
stantec_tbl<-data.frame(rp=c(10,25,50,100,200,500,1000,2000,10000),
                        stantec=c(65,77,87,97,108,123,135,148,181))


pandoc.table(stantec_tbl, caption="results based on stantec")

results based on stantec

rp	stantec
10	65
25	77
50	87
100	97
200	108
500	123
1000	135
2000	148
10000	181

Code

#max value by 1.13 site records changed from daily to 24 hrs.
site_max_val<-site.cmip6.db %>% 
  dplyr::select(Index,prcp) %>% 
  mutate(yr=year(Index)) %>%
  group_by(yr) %>% 
  summarize(prcp=max(prcp)*daily_24) %>% 
  pull(prcp) 

srk_all_tbl<-FA(a = site_max_val)$report%>% 
  dplyr::filter(Prob>=0.5) %>% 
  mutate(rp=(1/(1-Prob)) %>% round(1)) %>% 
  arrange(rp) %>%
  dplyr::select(rp,NORM,LN,GEV, GUMBEL,P3)

Tested distributions: [1] NORM LN GEV GUMBEL P3
———————— Chosen distributions: AIC AICc BIC ADC
LN LN LN LN
whose Maximum-Likelihood parameters are: NORM parameters of log(x): 3.909016 0.2501186

Code

srk_tbl<-srk_all_tbl %>% 
  dplyr::select(rp,GEV) %>% 
  dplyr::rename("srk"="GEV")


pander::emphasize.cols(3)

pandoc.table(srk_all_tbl,
             caption="results based on SRK, on report GEV was selected",
             round=1)

results based on SRK, on report GEV was selected

rp	NORM	LN	GEV	GUMBEL	P3
2	51.4	49.7	49.6	49.2	49.6
2.3	53.7	52	52	51.4	52
5	62.3	61.7	61.7	61.2	61.8
10	67.9	69.1	69.2	69.2	69.2
20	72.6	76	76.1	76.9	76
25	74	78.1	78.2	79.3	78
50	77.9	84.5	84.6	86.8	84.2
100	81.4	90.8	90.7	94.2	90
200	84.6	97	96.5	101.6	95.7
500	88.5	105	103.8	111.4	102.9
1000	91.2	111.1	109.1	118.8	108.2
2000	93.8	117.1	114.1	126.1	113.3
10000	99.3	131.2	125	143.3	125

Return periods and distributions

The information from Stantec and SRK is presented together. Stantec in red, while SRK in cyan; in the case of SRK, it is also presented other distributions.

In the case of SRK, the frequency analysis was obtained from the site records from 1985 to 2014 (implemented on the baseline), and matches the extend used on the climate change baseline.

Code

max_hist_version<-right_join(stantec_tbl,
                             srk_tbl,by="rp") %>% 
  arrange(rp)

max_hist_version_long1<-right_join(stantec_tbl,
                                   srk_tbl,by="rp") %>% 
  arrange(rp) %>% 
  reshape2::melt(.,id.vars=c("rp"))

max_hist_version_long2<-right_join(stantec_tbl,
                                   srk_all_tbl,by="rp") %>% 
  arrange(rp) %>% 
  reshape2::melt(.,id.vars=c("rp"))

tibble(max_val=site_max_val) %>% 
  arrange(-max_val) %>% 
  mutate(n=length(site_max_val),
         i=1:n) %>% 
  mutate(pos=(i-0.4)/(n+0.2)) %>% 
  mutate(rp=1/(pos)) %>% 
  ggplot(data=.,aes(x=rp,y=max_val))+
  geom_line(data=max_hist_version_long2,
            aes(x=rp,y=value,group=variable,col=variable))+
  geom_line(data=max_hist_version_long2 %>% 
              dplyr::filter(variable%in%c("stantec","GEV")),
            aes(x=rp,y=value,group=variable,col=variable),
            size=2,linetype="dashed")+
  geom_point(data=max_hist_version_long2,
             aes(x=rp,y=value,group=variable,col=variable))+
  geom_text_repel(data=max_hist_version_long2,
                  aes(x=rp,y=value,label=round(value,0),
                      group=variable,col=variable))+
  
  geom_point(size=2)+
  scale_x_log10()+
  labs(col="sources:",x="return period[yrs]",
       y="max precipitation at site mm/24hrs",
       title="Frequency analysis results compared with site observations (SRK, 2022)",
       subtitle="records from 1985 to 2014",
       caption="note: dashed lines are Stantec and SRK (GEV)")+
  theme_bw()+
  theme(legend.position="bottom")

Code

pandoc.table(max_hist_version,
             caption="max 24 hrs precipitation - values based on stantec/srk",
             missing="-",round=1)

max 24 hrs precipitation - values based on stantec/srk

rp	stantec	srk
2	-	49.6
2.3	-	52
5	-	61.7
10	65	69.2
20	-	76.1
25	77	78.2
50	87	84.6
100	97	90.7
200	108	96.5
500	123	103.8
1000	135	109.1
2000	148	114.1
10000	181	125

Comparison and evaluation of GCM results

These are the results from all the GCM applied at the site, the graph present the information from day 1 to day 30; while the table is just with 1 day.

Code

#define the median values
prcp_max_rc_high_tbl<-prcp_max_rc_tbl %>% 
  group_by(days,scenario,period,rp) %>% 
  summarise(rc=median(rc))

#just for the median
ggplot(data=prcp_max_rc_high_tbl,aes(x=rp,y=rc,col=days))+
  geom_point()+
  geom_line()+
  geom_hline(yintercept = 0)+
  facet_grid(scenario~period)+
  scale_x_log10()+
  scale_y_continuous(labels=scales::percent)+
  labs(y="rate of change",x="return periods [yrs]",
       title="Rate of change by GCM median values")+
  theme_bw()+
  theme(legend.position="bottom")

Code

#select the information for ssp370 and 1 day  
prcp_max_rc_median_tbl<-prcp_max_rc_high_tbl %>% 
  dplyr::filter(scenario=="ssp370",days=="day_01") %>% 
  reshape2::dcast(data = .,formula = scenario+rp~period,value.var = "rc")

pandoc.table(prcp_max_rc_median_tbl,round=2,
             caption="Median values from GCMs for ssp370 and 1 day [as percentage / 1 is 100%]")

Median values from GCMs for ssp370 and 1 day [as percentage / 1 is 100%]

scenario	rp	2040s	2060s	2080s
ssp370	2	0.16	0.24	0.29
ssp370	2.3	0.14	0.22	0.29
ssp370	5	0.08	0.14	0.19
ssp370	10	0.03	0.07	0.1
ssp370	20	-0.02	0.04	0.04
ssp370	25	-0.04	0.01	0.03
ssp370	50	-0.06	-0.06	0.01
ssp370	100	-0.13	-0.08	-0.03
ssp370	200	-0.17	-0.13	-0.07
ssp370	500	-0.22	-0.21	-0.14
ssp370	1000	-0.26	-0.25	-0.18
ssp370	2000	-0.3	-0.29	-0.22
ssp370	10000	-0.38	-0.37	-0.29

In the figure below the information is applied to Stantec values as well as SRK values; while SRK is also included but the selected distribution is GEV.

In the case of SRK estimations, because the “rate of change” of the historical values are increasing but in a small rate as the median of the rate of change for the GCM the resultant values are not monotonic. While if the values from Stantec are used, in this case, the historical are increasing in a higher rate; then climate change projections present still a increasing trend.

Code

#combine the rate of change with stantec values
prcp_max_prj<-left_join(srk_tbl,prcp_max_rc_median_tbl,by="rp") %>%
  right_join(stantec_tbl,.,by="rp") %>% 
  reshape2::melt(data=.,id.vars=c("rp","srk","stantec","scenario")) %>% 
  mutate(stantec=stantec*(value+1),
         srk=srk*(value+1)) %>% 
  # dplyr::filter(complete.cases(.)) %>%
  dplyr::rename(period="variable",rc="value") %>% 
  reshape2::melt(measure.vars=c("stantec","srk"))

ggplot(data=prcp_max_prj,
       aes(x=rp,y=value,col=period,
           group=interaction(period,variable)))+
  geom_point()+
  geom_line()+
  geom_text_repel(aes(label=round(value,0)))+
  geom_point(data=max_hist_version_long1,aes(x=rp,y=value),inherit.aes = F)+
  geom_line(data=max_hist_version_long1,aes(x=rp,y=value),inherit.aes = F,col="black")+
  scale_x_log10()+
  facet_wrap(~variable,ncol = 1)+
  theme_bw()+
  annotation_logticks(sides = "b")+
  labs(x="return periods [yrs]",y="max precipitation 24 hrs [mm]",
       title="Comparison of max precipitation historical and climate change",
       subtitle="Values evaluated based on SSP370 - median")+
  theme(legend.position="bottom")

CSV outputs

Code

data.table::fwrite(site.cmip6.db,
                   file = here::here("outputs","site_records.csv"))

data.table::fwrite(prcp_max_rc,
                   file = here::here("outputs",
                                     "prcp_max_cmip6_records.csv"))

data.table::fwrite(prcp_max_rc_tbl,
                   file = here::here("outputs",
                                     "prcp_max_cmip6_records_ext.csv"))

data.table::fwrite(site_cmip6_eval,
                   file = here::here("outputs",
                                     "site_cmip6_prcp.par.csv"))