From d1b026f0f5fcefe1022fb8e656e6cb9b8e138b4f Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 00:23:53 +0100 Subject: [PATCH 01/10] Update references.bib Added references for the OHC (Indian Ocean) tutorial --- vignettes/references.bib | 257 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 257 insertions(+) diff --git a/vignettes/references.bib b/vignettes/references.bib index 7f2bb2ee..fc4f18b1 100755 --- a/vignettes/references.bib +++ b/vignettes/references.bib @@ -782,3 +782,260 @@ @article{lee2018spatio pages={1--39}, year={2018} } + +@article{trenberth2014inbalance, + title={Earth’s Energy Imbalance}, + author={Trenberth, Kevin E and Fasullo, John T and Balmaseda, Magdalena A}, + journal={Journal of Climate}, + volume={27}, + number={9}, + pages={3129--3144}, + year={2014}, + publisher={American Meteorological Society} +} + +@article{duan2023storage, + title={Heat Storage in the Upper Indian Ocean: The Role of Wind-Driven Redistribution}, + author={Duan, Jing and Li, Yuanlong and Lin, Pengfei and Wang, Fan}, + journal={Journal of Climate}, + volume={36}, + number={7}, + pages={2221--2242}, + year={2023}, + publisher={American Meteorological Society} +} + +@article{desbruyeres2017trends, + title={Global and Full-Depth Ocean Temperature Trends during the Early Twenty-First Century from Argo and Repeat Hydrography}, + author={Desbruy{\`e}res, Damien and McDonagh, Elaine L and Thierry, Virginie}, + journal={Journal of Climate}, + volume={30}, + number={6}, + pages={1985--1997}, + year={2017}, + publisher={American Meteorological Society} +} + +@article{orsi1995acc, + title={On the meridional extent and fronts of the Antarctic Circumpolar Current }, + author={Orsi, Alejandro H and Whitworth, Thomas and Nowlin, Worth D}, + journal={Deep-Sea Research Part I}, + volume={42}, + number={5}, + pages={641--673}, + year={1995}, + publisher={Elsevier} +} + +@article{giglio2016subantarctic, + title={Subantarctic and Polar Fronts of the Antarctic Circumpolar Current and Southern Ocean Heat and Freshwater Content Variability: A View from Argo}, + author={Giglio, Donata and Johnson, Gregory C}, + journal={Journal of Physical Oceanography}, + volume={46}, + number={3}, + pages={749--768}, + year={2016}, + publisher={American Meteorological Society} +} + +@incollection{hood2024io, + author={Hood, Raleigh R and Ummenhofer, Caroline C and Philips, Helen E and Sprintall, Janet}, + title={Introduction to the Indian Ocean}, + editor={Ummenhofer, Caroline C and Hood, Raleigh R} + booktitle={The Indian Ocean and its Role in the Global Climate System}, + publisher={Elsevier}, + pages={1--31}, + year={2024}, + chapter={1} +} + +@incollection{phillips2024circulation, + author={Philips, Helen E and Menezes, Viviane V and Nagura, Motoki and McPhaden, Michael J and Vinayachandran, P N and Beal, Lisa M}, + title={Indian Ocean circulation}, + editor={Ummenhofer, Caroline C. and Hood, Raleigh R} + booktitle={The Indian Ocean and its Role in the Global Climate System}, + publisher={Elsevier}, + pages={169--203}, + year={2024}, + chapter={8} +} + +@article{huang2002eof, + title={Interannual variability in the tropical Indian Ocean}, + author={Huang, Bohua and Kinter, James L}, + journal={Journal of Geophysical Research: Oceans}, + volume={107}, + number={C11}, + year={2002}, + publisher={American Geophysical Union} +} + +@article{dandapat2021mld, + title={A numerical study on the role of atmospheric forcing on mixed layer depth variability in the Bay of Bengal using a regional ocean model}, + author={Dandapat, Sumit and Chakraborty, Arun and Kuttippurath, Jayanarayanan and Bhagawati, Chirantan and Sen, Radharani}, + journal={Ocean Dynamics}, + volume={71}, + pages={963--979}, + year={2021}, + publisher={Springer} +} + +@article{ali2018nio, + title={Dominant Modes of Upper Ocean Heat Content in the North Indian Ocean}, + author={Ali, Meer Mohammed and Singh, Neetu and Kumar, Manchikanti Suresh and Zheng, Yangxing and Bourassa, Mark and Kishtawal, Chandra Mohan and Rao, Chandu Venkateswara}, + journal={Climate}, + volume={6}, + number={3}, + year={2018}, + publisher={MDPI} +} + +@article{sandeep2018freshwater, + title={Impact of riverine freshwater forcing on the sea surface salinity simulations in the Indian Ocean}, + author={Sandeep, K K and Pant, Vimlesh and Girishkumar, M S and Rao, A D}, + journal={Journal of Marine Systems}, + volume={185}, + pages={40--58}, + year={2018}, + publisher={Elsevier} +} + +@article{dandapat2020runoff, + title={Impact of excess and deficit river runoff on Bay of Bengal upper ocean characteristics using an ocean general circulation model}, + author={Dandapat, Sumit and Gnanaseelan, C and Parekh, Anat}, + journal={Deep-Sea Research Part II: Topical Studies in Oceanography}, + volume={172}, + year={2020}, + publisher={Elsevier} +} + +@article{zhang2023bengal, + title={Surface cross-equatorial pathways of seawater from the Bay of Bengal}, + author={Zhang, Zhengbei and Wang, Jing and Hao, Jiajia and Yuan, Dongliang and Wang, Kungxiang}, + journal={Frontiers in Marine Science}, + volume={10}, + year={2023} +} + +@article{jaishree2023seasdyn, + title={Exploring the dynamics of seasonal surface features using coastal and regional ocean +community model}, + author={Jaishree, D and Ravichandran, P T and Thattal, D V}, + journal={Global Journal of Environmental Science and Management}, + volume={9}, + number={4}, + pages={741--752}, + year={2023} +} + +@article{trott2017eddies, + title={Variability of the Somali Current and eddies during the southwest monsoon regimes}, + author={Trott, Corinne B and Subrahmanyam, Bulusu and Murty, V S N}, + journal={Dynamics of Atmospheres and Oceans}, + volume={79}, + pages={43-55}, + year={2017}, + publisher={Elsevier} +} + +@incollection{roxy2020warming, + author={Roxy, M K and Gnanaseelan, Chellappan and Parekh, Anant and Chowdary, Jasti S and Singh, Shikha and Modi, Aditi and Kakatkar, Rashmi and Mohapatra, Sandeep and Dhara, Chirag and Shenoi, S C and Rajeevan, M}, + title={Indian Ocean Warming}, + editor={Krishnan, R. and Sanjay, J. and Gnanaseelan, Chellappan and Mujumdar, Milind and Kulkarni, Ashwini and Chakraborty, Supriyo} + booktitle={Assessment of Climate Change over the Indian Region}, + publisher={Springer}, + pages={191--206}, + year={2020}, + chapter={10} +} + +@article{chen2015equatorial, + title={Intraseasonal variability of upwelling in the equatorial Eastern Indian Ocean}, + author={Chen, Gengxin and Han, Weiqing and Li, Yuanlong and Wang, Dongxiao and Shinoda, Toshiaki}, + journal={Journal of Geophysical Research: Oceans}, + volume={120}, + number={11}, + pages={7598--7615}, + year={2015}, + publisher={American Geophysical Union} +} + +@article{mabarrok2023thermocline, + title={Assessment of thermocline depth bias in the Seychelles-Chagos Thermocline Ridge of the Southwestern Indian Ocean simulated by the CMIP6 models}, + author={Mubarrok, Saat and Azmiuddin, Fuad and Jang, Chan Joo}, + journal={Frontiers in Marine Science}, + volume={10}, + year={2023} +} + +@article{yokoi2008sctr, + title={Seasonal Variation of the Seychelles Dome}, + author={Yokoi, Takaaki and Tozuka, Tomoki and Yamagata, Toshio}, + journal={Journal of Climate}, + volume={21}, + number={15}, + pages={3740--3754}, + year={2008}, + publisher={American Meteorological Society} +} + +@article{mckenna2024sst, + title={Understanding Biases in Indian Ocean Seasonal SST in CMIP6 Models}, + author={McKenna, Sebastian and Santoso, Agus and Gupta, Alex Sen and Taschetto, Andr{\'e} S}, + journal={Journal of Geophysical Research: Oceans}, + volume={129}, + number={2}, + year={2024}, + publisher={American Geophysical Union} +} + +@article{vic2014mesoscale, + title={Mesoscale dynamics in the Arabian Sea and a focus on the Great Whirl life cycle: A numerical investigation using ROMS}, + author={Vic, C and Roullet, G and Carton, X and Capet, X}, + journal={Journal of Geophysical Research: Oceans}, + volume={119}, + number={9}, + pages={6422--6443}, + year={2014}, + publisher={American Geophysical Union} +} + +@article{chatterjee2019somali, + title={Annihilation of the Somali upwelling system during summer monsoon}, + author={Chatterjee, Abishek and Kumar, B Praveen and Prakash, Satya and Singh, Prerna}, + journal={Scientific Reports}, + volume={9}, + year={2019}, + publisher={Nature} +} + +@article{lhegaret2018gyres, + title={Shallow Cross-Equatorial Gyres of the Indian Ocean Driven by Seasonally Reversing Monsoon Winds}, + author={L'H{\'e}garet, Pierre and Beal, Lisa M and Elipot, shane and Laurindo, Lucas}, + journal={Journal of Geophysical Research: Oceans}, + volume={123}, + number={12}, + pages={8902--8920}, + year={2018}, + publisher={American Geophysical Union} +} + +@article{schott2009circulation, + title={Indian Ocean circulation and climate variability}, + author={Schott, Friedrich A and Xie, Shang-Ping and McCreary, Julian P}, + journal={Reviews of Geophysics}, + volume={47}, + number={1}, + year={2009}, + publisher={Wiley Online Library} +} + +@article{wu2020gyrestructure, + title={Structure and Seasonal Variation of the Indian Ocean Tropical Gyre Based on Surface Drifters}, + author={Wu, Wei and Yu, Yan and Qian, Yu-Kun and Cheng, Xuhua and Wang, Tianyu and Zhang, Lianyi and Peng, Shiqiu}, + journal={Journal of Geophysical Research: Oceans}, + volume={125}, + number={5}, + year={2020}, + publisher={American Geophysical Union} +} From 87aa348704976a3cfed3f8c4d3b609a066c2970b Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 15:41:12 +0100 Subject: [PATCH 02/10] Adding the reference formatting file --- vignettes/american-geophysical-union.csl | 703 +++++++++++++++++++++++ 1 file changed, 703 insertions(+) create mode 100644 vignettes/american-geophysical-union.csl diff --git a/vignettes/american-geophysical-union.csl b/vignettes/american-geophysical-union.csl new file mode 100644 index 00000000..3a6bf508 --- /dev/null +++ b/vignettes/american-geophysical-union.csl @@ -0,0 +1,703 @@ + + From f0fdf19a69718d7dbc8a8b1c8b195a149f930011 Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 15:41:50 +0100 Subject: [PATCH 03/10] Add the OHC tutorial --- vignettes/ohc_tutorial.Rmd | 845 +++++++++++++++++++++++++++++++++++++ 1 file changed, 845 insertions(+) create mode 100644 vignettes/ohc_tutorial.Rmd diff --git a/vignettes/ohc_tutorial.Rmd b/vignettes/ohc_tutorial.Rmd new file mode 100644 index 00000000..36fb0b05 --- /dev/null +++ b/vignettes/ohc_tutorial.Rmd @@ -0,0 +1,845 @@ +--- +title: "Mapping Ocean Heat Content from Synthetic Observations in Indian Ocean" +output: + bookdown::html_document2: + base_format: rmarkdown::html_vignette + fig_caption: yes +bibliography: "references.bib" +csl: "american-geophysical-union.csl" +link-citations: yes +author: "Anrijs Abele, Colin Morice, Man Ho Suen, Geoffrey Dawson" +vignette: > + %\VignetteIndexEntry{Mapping Ocean Heat Content from Synthetic Observations in Indian Ocean} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- +# Introduction + +Knowledge of the spatio-temporal evolution of ocean heat content (OHC) is essential +to our understanding of both past and future climate change. At global scales, ocean +heat uptake is a key element of the global energy and sea level budgets and provides +a powerful constraint on climate sensitivity and past climate forcings. At regional +scales, OHC changes are important to our understanding of climate variability and +for the development, initialization, and evaluation of seasonal-to-decadal prediction +systems. + +The world's oceans have a huge volume so even small changes in temperature can correspond +to large changes in energy. 90% of the excess energy trapped in the earth system by +anthropogenic greenhouse gas emissions goes into the oceans [@trenberth2014inbalance]. +The Indian Ocean serves as a significant reservoir of this extra energy, despite its small relatively small area [@duan2023storage]. @desbruyeres2017trends found that a half of the global top-layer (<700 m) heat uptake is in the Indian Ocean. + +## Measuring Ocean Heat Content + +To measure heat content, we measure temperatures through the whole volume of the +upper oceans, with measurements at different depths as well as at different locations. + +Historically, measurement of ocean temperature at depth have been obtained obtained using +a range of instrument deployed from ships. More recently, since the early 2000s, +autonomous robotic floats, known as Argo floats, have greatly increased the number +of available observations of the upper 2000m of the world’s oceans. These data sources +provide profiles of ocean temperatures from the surface downwards. + +## Tutorial Data Source + +This tutorial uses data from the MapEval4OceanHeat (ME4OH) project. The data +are derived from a 1/10th of a degree latitude-longitude ocean simulation that has been +sub-sampled to real world observed locations. While the synthetic profile data are +intended to be used to benchmark ocean heat mapping methods, the data are used in this +tutorial as a case study for application of 4DModeller to map ocean heat content. + +In this tutorial we will use point data of ocean heat for upper layers of the ocean +provided in the ME4OH data files (integrated temperature profiles from the surface down +to a depth of 306.25 m, multiplied by water density and heat capacity). + +The tutorial will use data from the Indian Ocean region. + +```{r, message=FALSE, warning=FALSE} +library(INLA) +library(inlabru) +library(sf) +library(sp) +library(ggplot2) +library(gridExtra) +library(scales) +``` + +## Loading Data +Here we load the Indian ocean shapefile and dataset. We smooth the boundary to make the mesh creation easier. +```{r loadData, error=TRUE, message=FALSE, warning=FALSE} +fdmr::retrieve_tutorial_data(dataset = "ohc", force_update = TRUE) + +# Load the dataset +df <- + fdmr::load_tutorial_data(dataset = "ohc", filename = "df_indian_ocean_2005.rds") + +# Load the polygon +ocean_sf <- + fdmr::load_tutorial_data(dataset = "ohc", filename = "indian_ocean.rds") + +# Deal with CRS and smooth the boundary with a buffer +ocean_crs <- fmesher::fm_crs(ocean_sf) +fmesher::fm_crs(ocean_sf) <- NULL +ocean_sf2 <- sf::st_union(sf::st_buffer(ocean_sf, 1.05)) +``` + +We plot the smoothed Indian ocean polygon to check if there are any issues. +```{r polygonPlot, error=TRUE, fig.id=TRUE, fig.cap="Smoothed polygon", fig.align='center'} +ggplot2::ggplot() + + ggplot2::geom_sf(data=ocean_sf2) +``` + +Then we convert it into a `sp` object. +```{r} +fmesher::fm_crs(ocean_sf) <- ocean_crs +ocean_sp <- sf::as_Spatial(ocean_sf) +``` + +We regularly sample some points ($20\times 20$) to put in the mesh builder function under the `fdmr` package to get a feeling of the parameter setting for the mesh. +```{r oceanPoints, error=FALSE} +ocean_pts <- sf::st_sample(ocean_sf, size = 20*20, + type = "random") + +ocean_pts_sp <- sf::as_Spatial(ocean_pts) + +colnames(ocean_pts_sp@coords) <- c("LONG", "LAT") +``` + +Here comes the mesh_builder function. We set the max edge and offset a bit larger to cater the computation. +```{r meshBuilder, eval=FALSE, error = FALSE} +fdmr::mesh_builder(spatial_data = as.data.frame(ocean_pts_sp), + max_edge = c(50,100), + offset = c(50,100), + cutoff = 1) +``` +We create the outer boundary by extending the outer boundary. +```{r boundary, error=TRUE} +ocean_bnd <- sf::st_cast( + sf::st_sf(geometry = fmesher::fm_nonconvex_hull(ocean_sf2, 5)), + "POLYGON" +) +``` + +We create the mesh for the Indian ocean with the outer boundary. +We remove Sri Lanka and keep Madagascar in the cutoff. +```{r mesh, error=TRUE, fig.id=TRUE, fig.cap="Mesh with original polygon and all datapoints", fig.align='center'} +boundary <- fmesher::fm_as_segm_list(base::list(ocean_sf2, ocean_bnd)) + +mesh <- fmesher::fm_mesh_2d_inla(boundary = boundary, + max.edge = c(5, 7), + offset = c(3, 3), + cutoff = 1) + +me4pts <- sf::st_as_sf(df, coords = c("lon", "lat")) + +ggplot2::ggplot() + ggplot2::geom_sf(data = ocean_sf2) + inlabru::gg(mesh) + inlabru::gg(data = me4pts, size=0.1) + ggplot2::theme(axis.title.x=element_blank(), axis.title.y=element_blank()) +``` + +# Model + +We considered one year of top layer density of the OHC (between 0 and 306.25 m depth) in 2005. Observations are indexed by the month of occurrence, resulting in 12 indices (equal to the number of time points). To pass the data to the function, we converted the data frame to an `sp` object. +```{r dataFormatting, error=TRUE} +# Transform latitudes and longtiudes to numeric +df$lat <- base::sapply(df$lat, as.numeric) +df$lon <- base::sapply(df$lon, as.numeric) + +# Add sp coordinates +df_f <- df +sp::coordinates(df_f) <- c('lon', 'lat') + +# Determine the group size +n_time <- base::as.integer(base::length(base::unique(df_f@data$time))) +``` +## Prior selection + +To implement the SPDE approach, we define the range un uncertainty priors. We assume that the probability of process spatial range (physically, a distance where correlation between two observations falls to 0.1) being under 25 degrees in latitude is 0.2. We also set the probability to 0.01 that the marginal standard deviation of the process exceeds 1. + +The process is assumed to be evolving temporally as an autoregressive process of the first order (AR1) with the probability of the temporal autocorrelation parameter to be over 0 at 0.9. + +```{r spde_setup, error=FALSE} +prior_range <- 25 +prior_Pgt <- 0.2 + +# Define the SPDE +ohc_spde <- INLA::inla.spde2.pcmatern(mesh, prior.range = c(prior_range, prior_Pgt), prior.sigma = c(1, 0.01)) + +# Define temporal parameter +rhoprior <- base::list(theta = base::list(prior = 'pccor1', param = c(0, 0.9))) +``` + +## Define the model formula + +We use a spatio-temporal model with only random effects and without an intercept. +```{r formula, eval=FALSE} +model_formula <- dohc_L1 ~ -1 + f( + main = coordinates, + model = ohc_spde, + group = time, + ngroup = n_time, + control.group = list( + model = 'ar1', + hyper = rhoprior + ) +) +``` + +## Fit the model + +```{r model, eval=FALSE} +bru_model <- + inlabru::bru(model_formula, + data = df_f, + family = "gaussian", + options = list( + verbose = TRUE, + bru_verbose = 4 + ) +) +``` + +# Output +The `model_viewer` of 4Dmodeller interactively shows the hyperparameters and predicted random field from the `INLA` model fit. In further sections we also plot the mean and standard deviation of the predicted field for March and October as well as the measures of temporal evolution. + +```{r modelviewer, eval=FALSE} +model_viewer(model_output = bru_model, mesh = mesh, measurement_data = df_f, data_distribution = "Gaussian") +``` + +To estimate posterior statistics from resampled model. We only use a small number of samples -- 100 (the default) -- to reduce the computational expanse. It is possible to use this function to also re-estimate the model fit with a new set of values. + +> ***NOTE:*** It takes approximately 20 min (depending on hardware) to run the prediction. You can instead run the consecutive cell to load the prediction already calculated by us. + +```{r eval=FALSE} +pred <- predict(bru_model, n.samples = 100) +``` + +If you decided to run the cell above, you do not need to run the cell below. +```{r error=TRUE} +pred <- fdmr::load_tutorial_data(dataset = "ohc", filename = "prediction_2005_L1.rds") +``` + +## Spatial fields + +First, we prepare the results for plotting. +```{r spatial_processing, error=TRUE} +# Add month to predictor +## Number of mesh nodes +mn <- mesh$n + +## Append the month index to the predictor +Predictor <- pred$Predictor +months <- data.frame(month = base::rep(base::seq(1,12,1), each = mn)) +Predictor$month <- months + +## Create a list for each month with mean values +mthpred <- list() +sdpred <- list() + +for (i in 1:12) { + mthpred[i] <- list(Predictor$mean[(1+mn*(i-1)):(mn*i)]) + sdpred[i] <- list(Predictor$sd[(1+mn*(i-1)):(mn*i)]) +} + +# Set the world map as a background +wmap <- ggplot2::map_data("world") + +# Assign the map coordinate extrema +mesh_lon_min <- base::min(mesh$loc[,1]) +mesh_lon_max <- base::max(mesh$loc[,1]) +mesh_lat_min <- base::min(mesh$loc[,2]) +mesh_lat_max <- base::max(mesh$loc[,2]) + +# Create a mask that only accepts the points over the ocean + +## Mesh node coords to sf object +mesh_sf <- sf::st_as_sf(base::data.frame(base::cbind(LAT = mesh$loc[,2], LONG = mesh$loc[,1])), coords = c("LONG", "LAT")) + +## The mask +mask_ocean <- base::unlist(sf::st_intersects(ocean_sf, mesh_sf, sparse = TRUE)) + +# Create the legend name +title <- bquote("DOHC" ~ (TJ ~ m^-2)) +title_diff <- bquote("dDOHC" ~ (TJ ~ m^-2)) +``` + +Then we plot the mean field for March and October. + +```{r meanPlots, error = FALSE, fig.id=TRUE, fig.cap="Mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +mean_median <- stats::median(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) +mean_min <- base::min(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) +mean_max <- base::max(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) + +plot_march <- + ggplot2::ggplot() + + inlabru::gg(data = mesh, + color = as.numeric(Predictor$mean[Predictor$month == 3]) + ) + + ggplot2::geom_polygon(data = wmap, + aes(x = long, y = lat, group = group), + fill = 'grey', + alpha = 1.0) + + ggplot2::coord_equal(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::ggtitle("Mean March") + + ggplot2::scale_fill_gradient2(name = title, + low = "blue", + high = "red", + midpoint = mean_median, + limits = c(mean_min, mean_max), + oob = squish + ) + + ggplot2::theme(legend.position="bottom", + legend.key.width = ggplot2::unit(1.8, "lines"), + plot.title = element_text(hjust = 0.5), + axis.title.x=element_blank(), + axis.title.y=element_blank(), + legend.title.align = 0.5 + ) + + ggplot2::guides( + fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + +plot_october <- + ggplot2::ggplot() + + inlabru::gg(data = mesh, + color = as.numeric(Predictor$mean[Predictor$month == 10]) + ) + + ggplot2::geom_polygon(data = wmap, + aes(x = long, y = lat, group = group), + fill = 'grey', + alpha = 1.0) + + ggplot2::coord_equal(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::ggtitle("Mean October") + + ggplot2::scale_fill_gradient2(name = title, + low = "blue", + high = "red", + midpoint = mean_median, + limits = c(mean_min, mean_max), + oob = squish + ) + + ggplot2::theme(legend.position = "bottom", + legend.key.width = ggplot2::unit(1.8, "lines"), + plot.title = element_text(hjust = 0.5), + axis.title.x=element_blank(), + axis.title.y=element_blank(), + legend.title.align = 0.5 + ) + + ggplot2::guides( + fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + +gridExtra::grid.arrange(plot_march, plot_october, ncol=2) +``` + +The map reproduces main features of ocean circulation, including the Subantarctic front (around 45--50$^\circ$S where transition between blue and violet occurs) [@orsi1995acc; @giglio2016subantarctic], Agulhas Current near southern African coast and Agulhas Return Current emanating from it eastwards at 40$^\circ$S [@hood2024io]. It also shows a region of increased OHC between 12$^\circ$ and 26$^\circ$S, coinciding with the subtropical gyre [@phillips2024circulation]. The maps show decreased upper-level OHC between 5$^\circ$ and 10$^\circ$S, in particular, at the Seychelles-Chagos thermocline ridge (SCTR) [@hood2024io]. Another centre of higher top-level OHC is in Arabian Sea while OHC is relatively lower in Bay of Bengal [@huang2002eof]. + +The differences between March and October are especially different in Arabian Sea. The western part of the basin and also at the Indian coast has a lower OHC in March compared to October. Similarly, in Bay of Bengal higher OHC appears in the central for October and western part in March. As expected, the heat content is relatively low due to weaker mixing during inter-monsoon [@dandapat2021mld]. The difference between March and October could also be explained by increased riverine freshwater influx as well as precipitation and cloudiness [@ali2018nio] in October compared to March [@sandeep2018freshwater], which causes a reduction in the OHC storage immediately at the plume area due to higher heat capcity of freshwater [@dandapat2020runoff] and also where freshwater is advected [@zhang2023bengal; @jaishree2023seasdyn] along the coastline. + +Similarly we show the standard deviation for March and October. +```{r stdevPlots, error = FALSE, fig.id=TRUE, fig.cap="Standard deviation field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +stdev_median <- stats::median(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) +stdev_max <- base::max(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) +n_march <- base::nrow(df[df$time == 3, ]) +n_october <- base::nrow(df[df$time == 10, ]) + +plot_march <- + ggplot2::ggplot()+ + inlabru::gg(mesh, + color = as.numeric(Predictor$sd[Predictor$month == 3])) + + ggplot2::geom_polygon(data = wmap, + aes(x = long, y = lat, group = group), + fill = 'grey', + alpha = 1.0) + + ggplot2::coord_equal(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::geom_point(data = df[df$time == 3, ], + aes(x = lon, + y = lat), + size = 0.01) + + ggplot2::ggtitle(paste0("Stdev March (n = ", toString(n_march), ")")) + + ggplot2::scale_fill_gradient2(name = title, + low = 'navyblue', + mid = 'darkmagenta', + high = 'darkorange1', + midpoint = stdev_median, + limits = c(0, stdev_max), + oob = squish + ) + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "bottom", + legend.key.width = ggplot2::unit(1.8, "lines"), + plot.title = ggplot2::element_text(hjust = 0.5), + axis.title.x = ggplot2::element_blank(), + axis.title.y = ggplot2::element_blank() + ) + + ggplot2::guides( + fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + + +plot_october <- + ggplot2::ggplot() + + inlabru::gg(mesh, + color = as.numeric(Predictor$sd[Predictor$month == 10]) + ) + + ggplot2::geom_polygon(data = wmap, + aes(x = long, y = lat, group = group), + fill = 'grey', + alpha = 1.0) + + ggplot2::coord_equal(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::geom_point(data = df[df$time == 10, ], + aes(x = lon, + y = lat), + size = 0.01) + + ggplot2::ggtitle(paste0("Stdev October (n = ", toString(n_october), ")")) + + ggplot2::scale_fill_gradient2(name = title, + low = 'navyblue', + mid = 'darkmagenta', + high = 'darkorange1', + midpoint = stdev_median, + limits = c(0, stdev_max), + oob = squish + ) + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "bottom", + legend.key.width = ggplot2::unit(1.8, "lines"), + plot.title = ggplot2::element_text(hjust = 0.5), + axis.title.x = ggplot2::element_blank(), + axis.title.y = ggplot2::element_blank() + ) + + ggplot2::guides(fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + + +gridExtra::grid.arrange(plot_march, plot_october, ncol=2) +``` + +The standard deviation is higher in areas with low number of profiles; especially, if there are no profiles in the vicinity like during March in Great Australian Bight. Despite the number of profiles being larger in October than March, the standard deviations over the whole field remained similar. This can be explained by the distribution of those extra profiles--mostly in a close proximity to other profiles. + +We also further explore the resampled posterior standard deviation. +```{r sterrPlots, error = FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} + +mc_serr_median <- stats::median(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) +mc_serr_max <- base::max(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) + +plot_march <- + ggplot2::ggplot()+ + inlabru::gg(mesh, + color = as.numeric(Predictor$mean.mc_std_err[Predictor$month == 3])) + + ggplot2::geom_polygon(data = wmap, + aes(x = long, y = lat, group = group), + fill = 'grey', + alpha = 1.0) + + ggplot2::coord_equal(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::ggtitle("Sterr of mean March") + + ggplot2::scale_fill_gradient2(name = title, + low = 'navyblue', + mid = 'darkmagenta', + high = 'darkorange1', + midpoint = mc_serr_median, + limits = c(0, mc_serr_max), + oob = squish + ) + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "bottom", + legend.key.width = ggplot2::unit(1.8, "lines"), + plot.title = ggplot2::element_text(hjust = 0.5), + axis.title.x = ggplot2::element_blank(), + axis.title.y = ggplot2::element_blank() + ) + + ggplot2::guides( + fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + + +plot_october <- + ggplot2::ggplot() + + inlabru::gg(mesh, + color = as.numeric(Predictor$mean.mc_std_err[Predictor$month == 10]) + ) + + ggplot2::geom_polygon(data = wmap, + aes(x = long, y = lat, group = group), + fill = 'grey', + alpha = 1.0) + + ggplot2::coord_equal(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::ggtitle("Sterr of mean October ") + + ggplot2::scale_fill_gradient2(name = title, + low = 'navyblue', + mid = 'darkmagenta', + high = 'darkorange1', + midpoint = mc_serr_median, + limits = c(0, mc_serr_max), + oob = squish + ) + + ggplot2::theme_bw() + + ggplot2::theme(legend.position = "bottom", + legend.key.width = ggplot2::unit(1.8, "lines"), + plot.title = ggplot2::element_text(hjust = 0.5), + axis.title.x = ggplot2::element_blank(), + axis.title.y = ggplot2::element_blank() + ) + + ggplot2::guides(fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + + +gridExtra::grid.arrange(plot_march, plot_october, ncol=2) +``` + +As one can see, both standard deviation and resampled standard error of mean are highest in the same areas. However, the magnitude is much lower. + +Finally, we evaluate the differences between synthetic observations and predictions at the locations of those observations. First, we calculate the predicted values at each location. +```{r error=TRUE} +# Define the predictions and observations for March +df_march <- df[df$time == 3, ] +df_sf_march <- sf::st_as_sf(df_march, coords = c("lon", "lat")) +pred_march <- Predictor$mean[Predictor$month == 3] + +# Evaluate at observation locations +pred_eval_march <- fmesher::fm_evaluate(mesh, pred_march, loc = df_sf_march$geometry) +pred_diff_march <- base::data.frame(diff = (df_march$dohc_L1 - pred_eval_march)) +pred_diff_march$Lat <- df_march$lat +pred_diff_march$Long <- df_march$lon +geom_diff_march <- sf::st_as_sf(pred_diff_march, coords = c("Long", "Lat")) + +# Do the same for October +df_october <- df[df$time == 10, ] +df_sf_october <- sf::st_as_sf(df_october, coords = c("lon", "lat")) +pred_october <- Predictor$mean[Predictor$month == 10] +pred_eval_october <- fmesher::fm_evaluate(mesh, pred_october, loc = df_sf_october$geometry) +pred_diff_october <- base::data.frame(diff = (df_october$dohc_L1 - pred_eval_october)) +pred_diff_october$Lat <- df_october$lat +pred_diff_october$Long <- df_october$lon +geom_diff_october <- sf::st_as_sf(pred_diff_october, coords = c("Long", "Lat")) +``` + +Then we plot these differences. +```{r differencePlots, error = FALSE, fig.id=TRUE, fig.cap="Differences between synthetic observations and predictions for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +diff_suptitle_3 <- bquote("Diff (synthetic - prediction)"~ .(toString(df_march$year[1]))/.(toString(df_march$month[1]))) +diff_suptitle_10 <- bquote("Diff (synthetic - prediction)"~ .(toString(df_october$year[1]))/.(toString(df_october$month[1]))) +min_diff <- base::min(c(pred_diff_march$diff, pred_diff_october$diff), na.rm = TRUE) +max_diff <- base::max(c(pred_diff_march$diff, pred_diff_october$diff), na.rm = TRUE) + +plot_march <- + ggplot2::ggplot() + + ggplot2::geom_polygon( + data = wmap, + aes(x = long, + y = lat, + group = group + ), + fill = 'grey', + alpha = 1.0 + ) + + ggplot2::geom_sf(data = geom_diff_march, aes(color = diff), + size = 1.0) + + ggplot2::coord_sf(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::ggtitle(diff_suptitle_3) + + ggplot2::scale_color_gradient2( + name = title_diff, + low = "blue", + high = "red", + midpoint = 0, + limits = c(min_diff, max_diff), + oob = squish + ) + + ggplot2::theme(legend.position="bottom", + legend.key.width = ggplot2::unit(2.5, "lines"), + plot.title = element_text(hjust = 0.5), + axis.title.x=element_blank(), + axis.title.y=element_blank(), + legend.title = element_text(hjust = 0.5) + ) + + ggplot2::guides( + fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + +plot_october <- + ggplot2::ggplot() + + ggplot2::geom_polygon( + data = wmap, + aes(x = long, + y = lat, + group = group + ), + fill = 'grey', + alpha = 1.0 + ) + + ggplot2::geom_sf(data = geom_diff_october, aes(color = diff), + size = 1.0) + + ggplot2::coord_sf(xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE) + + ggplot2::ggtitle(diff_suptitle_10) + + ggplot2::scale_color_gradient2( + name = title_diff, + low = "blue", + high = "red", + midpoint = 0, + limits = c(min_diff, max_diff), + oob = squish + ) + + ggplot2::theme(legend.position="bottom", + legend.key.width = ggplot2::unit(2.5, "lines"), + plot.title = element_text(hjust = 0.5), + axis.title.x=element_blank(), + axis.title.y=element_blank(), + legend.title = element_text(hjust = 0.5) + ) + + ggplot2::guides( + fill = guide_colorbar(title.position = "left", + title.vjust = 1, + title.hjust = 1 + ) + ) + +gridExtra::grid.arrange(plot_march, plot_october, ncol=2) +``` + +The differences show that there is no significant bias. The regions with larger differences are typically around zones of higher mesoscale activity, for instance, Agulhas and Agulhas Return current, Leeuwin Current near Western Australian coastline, as well as western Arabian Sea (Somali Current). There is a marked increase in differences in the Arabian Sea in October compared to March, owing to a higher number of eddies after southwest monsoon in this region [@trott2017eddies]. + +## Temporal evolution + +Next we would like to show how the field evolves temporally. The field-mean evolution of mean and standard deviation of the posterior solution is shown in the plots below. + +```{r temporalEval, error=TRUE, fig.id=TRUE, fig.cap="Area-averaged mean and standard deviation for the Indian Ocean in 2005", fig.align='center', fig.width=4, fig.height=3} +meanFunc <- function(x) { + base::mean(base::data.frame(x)[mask_ocean,]) +} + +mean_mth <- lapply(mthpred, meanFunc) +mean_std <- lapply(sdpred, meanFunc) + +list_summary <- list( + "Mean" = base::cbind(x = 1:12, y = unlist(mean_mth)), + "Stdev" = base::cbind(x = 1:12, y = unlist(mean_std)) +) + +summ <- base::data.frame(base::do.call(rbind, list_summary)) +summ$parameter <- + base::factor(base::rep(names(list_summary), + times = base::sapply(list_summary, nrow) + ), + levels = c("Mean", "Stdev") + ) + +ggplot2::ggplot(summ, ggplot2::aes(x = x, y = y)) + + ggplot2::geom_line() + + ggplot2::facet_wrap(~parameter, scales = "free") + + ggplot2::xlab('Month') + + ggplot2::ylab(title) + + ggplot2::scale_x_continuous(breaks = scales::pretty_breaks()) + + ggplot2::theme_bw() + +``` + +The plots show that the mean over ocean peaks between March and June and reaches minimum in September to October, in line with the southern hemisphere (austral) seasonal cycle. The uncertainties are the highest between November and May due to a lower number of measurements. + +To show the temporal evolution across the Indian Ocean basin, we also plotted March and October values of top-level DOHC along a transect. We selected the TOGA/WOCE (Tropical Ocean Global Atmosphere/World Ocean Circulation Experiment) transect IX12, which runs from Gulf of Aden to Freemantle in Western Australia. This transect is frequently (12--15 times a year) conducted by the commercial fleet that deploy expendable bathythermographs (XBTs) along the route. The data has been collected since 1983 until now. + +We first define the locations of this transect and select some points for further exploration. +```{r} +# The bounding box +min_lat_transect <- -34.5833 +max_lat_transect <- 11.9 +min_lon_transect <- 51.8167 +max_lon_transect <- 114.6167 + +# Set the matrix with selected points +transect_points <- + base::matrix( + data = c(56.42322, 73.55196, 86.11196, 104.95196, 8.4903437, -4.1880047, -13.4846647, -27.4296547), + nrow = 4, + ncol = 2, + dimnames = list(c(), c("lon", "lat"))) + +# Create a line +transect_matrix <- + base::matrix( + data = c(max_lon_transect, min_lon_transect, min_lat_transect, max_lat_transect), + nrow = 2, + ncol = 2, + dimnames = list(c(), c("lon", "lat"))) + +transect_line <- + sf::st_geometry( + sf::st_linestring(transect_matrix), + crs = 4326) + +# Sample the line at 75 points +tr_sample_a <- sf::st_sample(transect_line, size = 75, type = "regular") +tr_points_a <- sf::st_coordinates(tr_sample_a)[, 1:2] +colnames(tr_points_a) <- c("lon", "lat") +``` + +The location on top of the mesh is shown in the figure below. A small sample of points along the transect that we choose to further show the change in time is also displayed. + +```{r transectMap, error = FALSE, fig.id=TRUE, fig.cap="Transect IX12 on a map", fig.align="center", fig.width=6, fig.height=4} +ggplot2::ggplot() + + inlabru::gg(mesh) + + ggplot2::geom_polygon( + data = wmap, + aes(x = long, + y = lat, + group = group + ), + fill = 'grey', + alpha = 1.0 + ) + + ggplot2::coord_equal( + xlim = c(mesh_lon_min, mesh_lon_max), + ylim = c(mesh_lat_min, mesh_lat_max), + expand = FALSE + ) + + ggplot2::ggtitle("Mesh and transect IX12 with a selection of points") + + ggplot2::geom_line( + data = data.frame(transect_matrix), + aes(x = lon, y = lat), + color = "black" + ) + + ggplot2::geom_point( + data = data.frame(transect_points), + aes(x = lon, y = lat), + color = "black", + size = 2 + ) + + ggplot2::annotate( + "text", + x = transect_points[1:4, c("lon")], + y = transect_points[1:4, c("lat")] + 3, + label = sprintf("%d", 1:4) + ) + + ggplot2::theme( + plot.title = element_text(hjust = 0.5), + axis.title.x=element_blank(), + axis.title.y=element_blank() + ) +``` + +Similarly as with the difference calculations above, we can determine the predicted values at the sampled points along the transect and also at the set of four points we choose to inspect more thoroughly. +```{r error=TRUE} +# Create an empty list and matrices for evalutation +monthlist <- c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December") +eval_list <- base::vector("list", 12) +eval_point <- matrix(, nrow = 12, ncol = base::nrow(transect_points)) +std_min_point <- matrix(, nrow = 12, ncol = base::nrow(transect_points)) +std_max_point <- matrix(, nrow = 12, ncol = base::nrow(transect_points)) + +# Evaluate field mean and standard deviation for each month at locations along the transect and select points +for (i in 1:12) { + eval.part <- base::vector("list", 0) + eval.part$eval <- fmesher::fm_evaluate(mesh, as.numeric(Predictor$mean[Predictor$month == i]), loc = tr_points_a) + eval.part$std <- fmesher::fm_evaluate(mesh, as.numeric(Predictor$sd[Predictor$month == i]), loc = tr_points_a) + eval.part$ymin <- eval.part$eval - eval.part$std + eval.part$ymax <- eval.part$eval + eval.part$std + eval.part$Lat <- tr_points_a[, c("lat")] + eval.part$Long <- tr_points_a[, c("lon")] + eval.part$Month <- base::rep(monthlist[i], base::nrow(tr_points_a)) + eval_list[[i]] <- eval.part + + eval_mean <- fmesher::fm_evaluate(mesh, as.numeric(Predictor$mean[Predictor$month == i]), loc = transect_points) + std_mean <- fmesher::fm_evaluate(mesh, as.numeric(Predictor$sd[Predictor$month == i]), loc = transect_points) + eval_point[i,] <- eval_mean + std_min_point[i,] <- eval_mean - std_mean + std_max_point[i,] <- eval_mean + std_mean +} + +# Create and empty list to create timeseries at four points +point_eval <- base::vector("list", 4) + +# Add the string to determine whether the point is in northern or southern hemisphere +ns_list <- base::lapply(transect_points[, c("lat")], function(x) {if (sign(x) > 0) {"N"} else {"S"}}) + +# Create a dataframe for each of the points in list +for (i in 1:4) { + point_eval[[i]] <- + data.frame( + param = paste0(i, " (", base::round(base::abs(transect_points[i, c("lat")]), digits = 1), intToUtf8(176), ns_list[[i]], ", ", base::round(base::abs(transect_points[i, c("lon")]), digits = 1), intToUtf8(176), "E)"), + x = 1:12, + y = eval_point[, i], + ymin = std_min_point[, i], + ymax = std_max_point[, i] + ) +} +``` + +The plot below shows the predicted DOHC for the top layer in March and October (austral autumn and spring, respectively) along the transect shown in the previous plot. The line shows the mean field predictions while envelope -- the standard deviation. +```{r transectPlot, error = FALSE, fig.id=TRUE, fig.cap="Timeseries for 2005 at selected points", fig.align='center', fig.width=6, fig.height=4} +x_title <- paste0("Latitude (", intToUtf8(176), "N)") + +ggplot2::ggplot() + + ggplot2::geom_ribbon(data = base::as.data.frame(eval_list[3]), aes(x = Lat, ymin = ymin, ymax = ymax, color = Month, fill = Month, linetype = Month), alpha = 0.5) + + ggplot2::geom_line(data = base::as.data.frame(eval_list[3]), aes(x = Lat, y = eval, color = Month, linetype = Month)) + + ggplot2::geom_ribbon(data = base::as.data.frame(eval_list[10]), aes(x = Lat, ymin = ymin, ymax = ymax, color = Month, fill = Month, linetype = Month), alpha = 0.5) + + ggplot2::geom_line(data = base::as.data.frame(eval_list[10]), aes(x = Lat, y = eval, color = Month, linetype = Month)) + + ggplot2::scale_x_continuous(breaks = scales::pretty_breaks()) + + ggplot2::xlab(x_title) + + ggplot2::ylab(title) + + ggplot2::ggtitle("DOHC along IX12 in 2005") + + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +``` + +The figure shows an increase in DOHC northward from the southernmost point until 17$^\circ$S in the southern hemisphere. There the curve peaks for the first time and then rapidly decays until 12$^\circ$S. The curve has a local minimum between -13$^\circ$S and -5$^\circ$S, after which DOHC starts to increase again, peaking at 7$^\circ$N and leveling off slightly afterwards. The seasonal (autumn/spring) pattern shows that DOHC is higher in March than October in most places. Exception is the segment between 17$^\circ$S and 0$^\circ$, where curves match (descending limb) or where October values are higher (local minimum). +The whole segment is known to have a persistent cooling, which has been present throughout the year for multiple decades as shown in other studies [@roxy2020warming; @duan2023storage]. +The minimum can be explained by the equatorial upwelling (associated with enhanced vertical advection) in so-called equatorial cold tongue (this particular region known as SCTR), which is caused by negative wind stress curl due to surface wind divergence [@chen2015equatorial] between equatorial westerlies and southeasterly trade wind [@mabarrok2023thermocline]. Wind stress is regulated by the Indian monsoon due to the change in wind direction (from easterlies in boreal summer to westerlies in boreal winter) [@yokoi2008sctr]. +The segment where both curves match is shown by @mckenna2024sst to have a higher net heat flux in March while net advection brings more heat in October due to negative vertical advection in March (associated with the same process as for the "cold tongue") while the meridional heat advection is stronger in October. +Elsewhere, March OHC has higher than October OHC. In southern hemisphere, it is expected due to increased insolation during austral summer and the lag for the ocean to heat up. However, we also observe higher OHC in northern hemisphere, where maximum values occur at the end of boreal winter. This segment has similar net heat flux values because it peaks twice during the year and both peaks are of similar magnitude [@mckenna2024sst]. However, advection contributes to more heat loss in October than March. +The final (northernmost) segment lies in close proximity to land and is affected by coastal upwelling near Somali coast [@vic2014mesoscale; @chatterjee2019somali]. + +The panel below shows the mean and standard deviation at the four sample locations along the transect shown on the map earlier. Compared to the field mean, we can observe that temporal evolution is spatially variable +```{r timeseriesPlots, error = FALSE, fig.id=TRUE, fig.cap="Timeseries for 2005 at selected points", fig.align='center', fig.width=6, fig.height=6} +ggplot2::ggplot(do.call(rbind, point_eval)) + + ggplot2::geom_line(aes(x = x, y = y), color = "black", linetype = "dashed") + + ggplot2::geom_point(aes(x = x, y = y), color = 'black', size = 1.6) + + ggplot2::geom_ribbon(aes(x = x, ymin = ymin, ymax = ymax), fill = 'grey', color = 'black', alpha = 0.2) + + ggplot2::facet_wrap(~param, nrow = 2, ncol = 2, scales = "fixed") + + ggplot2::scale_x_continuous(breaks = scales::pretty_breaks()) + + ggplot2::xlab('Month') + + ggplot2::ylab(title) + + ggplot2::ggtitle("2005 timeseries along IX12") + + ggplot2::theme_bw() + + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +``` + +The northernmost point (1) has the highest DOHC values on average. There is a pronounced seasonal cycle, which has a maximum plateau (within standard deviation) from March to September and more pronounced trough reaching minimum in November--December. The plateau coincides with increased insolation when the sun is in zenith over the northern hemisphere. However, it does not produce the typical biannual maximum. The area is subject to monsoon winds, which cause seasonal Somali Current reversal (southward during boreal summer) and creation of a small gyre (Great Whirl) [@lhegaret2018gyres]. The winds reverse from southwesterly (JJA) to northeasterly (DJF), causing enhanced boreal winter surface cooling in the northern Arabian Sea due to the reduction in latent heat flux when wind brings dry and cool continental air from Asia [@schott2009circulation]. This suggests that an interplay between insolation and wind-driven surface cooling is likely causing the seasonal changes in this location. +The standard deviation envelope narrows slightly at the descending limb, likely due to a higher number of observations in the area during that time. + +Point 2 has a relatively even profile with lower average than at point 1. There is a peak in June that levels off over a couple of months until August. This could partly be caused by weaker equatorial upwelling as well as increased net heat flux as discussed earlier. This location is also near the South Equatorial Counter current, which is not present at this particular point during the boreal summer [@wu2020gyrestructure]. Thus, it is less affected by negative zonal heat advection when the top-layer OHC values peak. The standard deviation narrows near the peak, but is similar to point 1 elsewhere. + +Point 3 are another relatively flat timeseries, which has a small peak in January, leveling off until April. Within a standard deviation, the rest of values are similar. South Equatorial current is present in the vicinity of point 3 and provides some heat advection. @mckenna2024sst shows that zonal and meridional advection dampens the heat loss due to the net negative heat flux during austral winter. The standard deviation is similar for all timesteps. + +The southernmost point (4) mostly has the lowest DOHC values on average, which can be related to less solar radiation at higher latitudes. There is also a prominent seasonal cycle with a maximum in February and minimum in November. The standard deviation is similar for all months and also of similar magnitude as point 3. + +## References {-} + +
\ No newline at end of file From e5a7e70cfe3eed3ec229dc2f05d429a26f2368f3 Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 15:45:59 +0100 Subject: [PATCH 04/10] Update CHANGELOG.md --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8d3bcbc5..317917ad 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fixed `model_viewer` app to work with random effect models - [PR #317](https://github.com/4DModeller/fdmr/pull/317) - Updated `plot_map_leaflet` to feature the same selection of basemap tiles as `mapview` and fix features not appearing - [PR #310](https://github.com/4DModeller/fdmr/pull/310) +### Added +- Added ocean heat content mapping tutorial - [PR #321](https://github.com/4DModeller/fdmr/pull/321) + ## [0.2.0] - 2023-12-19 ### Fixed From 3e1912472bf641fc3e99c355f846ac592a9dd149 Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 15:48:28 +0100 Subject: [PATCH 05/10] Update _pkgdown.yml --- _pkgdown.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 4abdf5c0..7ff16348 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -45,6 +45,8 @@ navbar: href: articles/hydro.html - text: "Priors exploration" href: articles/priors.html + - text: "Ocean heat content mapping" + href: articles/ohc_tutorial.html - text: --- - text: "FAQ" - text: INLA Crash FAQ From cec8c9f1e54e58302862141677ba1a8d62cf3734 Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 15:51:06 +0100 Subject: [PATCH 06/10] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 317917ad..05380de5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,7 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated `plot_map_leaflet` to feature the same selection of basemap tiles as `mapview` and fix features not appearing - [PR #310](https://github.com/4DModeller/fdmr/pull/310) ### Added -- Added ocean heat content mapping tutorial - [PR #321](https://github.com/4DModeller/fdmr/pull/321) +- Added ocean heat content mapping tutorial - [PR #318](https://github.com/4DModeller/fdmr/pull/318) ## [0.2.0] - 2023-12-19 From 11e5d3447e4bf4a804ee1277ac94128e2f7e5c4b Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 16:01:09 +0100 Subject: [PATCH 07/10] Update references.bib --- vignettes/references.bib | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vignettes/references.bib b/vignettes/references.bib index fc4f18b1..caf74fe7 100755 --- a/vignettes/references.bib +++ b/vignettes/references.bib @@ -841,7 +841,7 @@ @article{giglio2016subantarctic @incollection{hood2024io, author={Hood, Raleigh R and Ummenhofer, Caroline C and Philips, Helen E and Sprintall, Janet}, title={Introduction to the Indian Ocean}, - editor={Ummenhofer, Caroline C and Hood, Raleigh R} + editor={Ummenhofer, Caroline C and Hood, Raleigh R}, booktitle={The Indian Ocean and its Role in the Global Climate System}, publisher={Elsevier}, pages={1--31}, @@ -852,7 +852,7 @@ @incollection{hood2024io @incollection{phillips2024circulation, author={Philips, Helen E and Menezes, Viviane V and Nagura, Motoki and McPhaden, Michael J and Vinayachandran, P N and Beal, Lisa M}, title={Indian Ocean circulation}, - editor={Ummenhofer, Caroline C. and Hood, Raleigh R} + editor={Ummenhofer, Caroline C. and Hood, Raleigh R}, booktitle={The Indian Ocean and its Role in the Global Climate System}, publisher={Elsevier}, pages={169--203}, @@ -941,7 +941,7 @@ @article{trott2017eddies @incollection{roxy2020warming, author={Roxy, M K and Gnanaseelan, Chellappan and Parekh, Anant and Chowdary, Jasti S and Singh, Shikha and Modi, Aditi and Kakatkar, Rashmi and Mohapatra, Sandeep and Dhara, Chirag and Shenoi, S C and Rajeevan, M}, title={Indian Ocean Warming}, - editor={Krishnan, R. and Sanjay, J. and Gnanaseelan, Chellappan and Mujumdar, Milind and Kulkarni, Ashwini and Chakraborty, Supriyo} + editor={Krishnan, R. and Sanjay, J. and Gnanaseelan, Chellappan and Mujumdar, Milind and Kulkarni, Ashwini and Chakraborty, Supriyo}, booktitle={Assessment of Climate Change over the Indian Region}, publisher={Springer}, pages={191--206}, From df7e974a3d53b6f36e72ad6e65c98d81fe0f1f97 Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 16:29:16 +0100 Subject: [PATCH 08/10] Disable map plots for ohc_tutorial.Rmd --- vignettes/ohc_tutorial.Rmd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/vignettes/ohc_tutorial.Rmd b/vignettes/ohc_tutorial.Rmd index 36fb0b05..dad4fdde 100644 --- a/vignettes/ohc_tutorial.Rmd +++ b/vignettes/ohc_tutorial.Rmd @@ -262,7 +262,7 @@ title_diff <- bquote("dDOHC" ~ (TJ ~ m^-2)) Then we plot the mean field for March and October. -```{r meanPlots, error = FALSE, fig.id=TRUE, fig.cap="Mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r meanPlots, error = FALSE, eval=FALSE, fig.id=TRUE, fig.cap="Mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} mean_median <- stats::median(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) mean_min <- base::min(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) mean_max <- base::max(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) @@ -343,7 +343,7 @@ The map reproduces main features of ocean circulation, including the Subantarcti The differences between March and October are especially different in Arabian Sea. The western part of the basin and also at the Indian coast has a lower OHC in March compared to October. Similarly, in Bay of Bengal higher OHC appears in the central for October and western part in March. As expected, the heat content is relatively low due to weaker mixing during inter-monsoon [@dandapat2021mld]. The difference between March and October could also be explained by increased riverine freshwater influx as well as precipitation and cloudiness [@ali2018nio] in October compared to March [@sandeep2018freshwater], which causes a reduction in the OHC storage immediately at the plume area due to higher heat capcity of freshwater [@dandapat2020runoff] and also where freshwater is advected [@zhang2023bengal; @jaishree2023seasdyn] along the coastline. Similarly we show the standard deviation for March and October. -```{r stdevPlots, error = FALSE, fig.id=TRUE, fig.cap="Standard deviation field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r stdevPlots, error = FALSE, eval=FALSE, fig.id=TRUE, fig.cap="Standard deviation field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} stdev_median <- stats::median(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) stdev_max <- base::max(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) n_march <- base::nrow(df[df$time == 3, ]) @@ -433,7 +433,7 @@ gridExtra::grid.arrange(plot_march, plot_october, ncol=2) The standard deviation is higher in areas with low number of profiles; especially, if there are no profiles in the vicinity like during March in Great Australian Bight. Despite the number of profiles being larger in October than March, the standard deviations over the whole field remained similar. This can be explained by the distribution of those extra profiles--mostly in a close proximity to other profiles. We also further explore the resampled posterior standard deviation. -```{r sterrPlots, error = FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r sterrPlots, error = FALSE, eval=FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} mc_serr_median <- stats::median(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) mc_serr_max <- base::max(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) @@ -842,4 +842,4 @@ The southernmost point (4) mostly has the lowest DOHC values on average, which c ## References {-} -
\ No newline at end of file +
From c6d8c0bcb825d34a55113c08c838867fac9902b3 Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 16:39:36 +0100 Subject: [PATCH 09/10] Include warnings in ohc_tutorial.Rmd --- vignettes/ohc_tutorial.Rmd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/vignettes/ohc_tutorial.Rmd b/vignettes/ohc_tutorial.Rmd index dad4fdde..0dc95426 100644 --- a/vignettes/ohc_tutorial.Rmd +++ b/vignettes/ohc_tutorial.Rmd @@ -219,7 +219,7 @@ pred <- fdmr::load_tutorial_data(dataset = "ohc", filename = "prediction_2005_L1 ## Spatial fields First, we prepare the results for plotting. -```{r spatial_processing, error=TRUE} +```{r spatial_processing, error=TRUE, message=FALSE, warning=FALSE} # Add month to predictor ## Number of mesh nodes mn <- mesh$n @@ -262,7 +262,7 @@ title_diff <- bquote("dDOHC" ~ (TJ ~ m^-2)) Then we plot the mean field for March and October. -```{r meanPlots, error = FALSE, eval=FALSE, fig.id=TRUE, fig.cap="Mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r meanPlots, error = FALSE, fig.id=TRUE, fig.cap="Mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} mean_median <- stats::median(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) mean_min <- base::min(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) mean_max <- base::max(c(Predictor$mean[Predictor$month == 3], Predictor$mean[Predictor$month == 10]), na.rm = TRUE) @@ -343,7 +343,7 @@ The map reproduces main features of ocean circulation, including the Subantarcti The differences between March and October are especially different in Arabian Sea. The western part of the basin and also at the Indian coast has a lower OHC in March compared to October. Similarly, in Bay of Bengal higher OHC appears in the central for October and western part in March. As expected, the heat content is relatively low due to weaker mixing during inter-monsoon [@dandapat2021mld]. The difference between March and October could also be explained by increased riverine freshwater influx as well as precipitation and cloudiness [@ali2018nio] in October compared to March [@sandeep2018freshwater], which causes a reduction in the OHC storage immediately at the plume area due to higher heat capcity of freshwater [@dandapat2020runoff] and also where freshwater is advected [@zhang2023bengal; @jaishree2023seasdyn] along the coastline. Similarly we show the standard deviation for March and October. -```{r stdevPlots, error = FALSE, eval=FALSE, fig.id=TRUE, fig.cap="Standard deviation field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r stdevPlots, error = FALSE, fig.id=TRUE, fig.cap="Standard deviation field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} stdev_median <- stats::median(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) stdev_max <- base::max(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) n_march <- base::nrow(df[df$time == 3, ]) @@ -433,7 +433,7 @@ gridExtra::grid.arrange(plot_march, plot_october, ncol=2) The standard deviation is higher in areas with low number of profiles; especially, if there are no profiles in the vicinity like during March in Great Australian Bight. Despite the number of profiles being larger in October than March, the standard deviations over the whole field remained similar. This can be explained by the distribution of those extra profiles--mostly in a close proximity to other profiles. We also further explore the resampled posterior standard deviation. -```{r sterrPlots, error = FALSE, eval=FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r sterrPlots, error = FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} mc_serr_median <- stats::median(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) mc_serr_max <- base::max(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) From b7db45fcfde9536c6698674863963851a2f0054d Mon Sep 17 00:00:00 2001 From: "Anrijs K. Abele" <57137242+aabelean@users.noreply.github.com> Date: Sat, 29 Jun 2024 16:59:50 +0100 Subject: [PATCH 10/10] Update ohc_tutorial.Rmd --- vignettes/ohc_tutorial.Rmd | 59 +++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/vignettes/ohc_tutorial.Rmd b/vignettes/ohc_tutorial.Rmd index 0dc95426..71442a67 100644 --- a/vignettes/ohc_tutorial.Rmd +++ b/vignettes/ohc_tutorial.Rmd @@ -26,7 +26,7 @@ systems. The world's oceans have a huge volume so even small changes in temperature can correspond to large changes in energy. 90% of the excess energy trapped in the earth system by anthropogenic greenhouse gas emissions goes into the oceans [@trenberth2014inbalance]. -The Indian Ocean serves as a significant reservoir of this extra energy, despite its small relatively small area [@duan2023storage]. @desbruyeres2017trends found that a half of the global top-layer (<700 m) heat uptake is in the Indian Ocean. +The Indian Ocean serves as a significant reservoir of this extra energy, despite its relatively small area [@duan2023storage]. @desbruyeres2017trends found that half of the global top-layer (<700 m) heat uptake is in the Indian Ocean. ## Measuring Ocean Heat Content @@ -47,7 +47,7 @@ sub-sampled to real world observed locations. While the synthetic profile data a intended to be used to benchmark ocean heat mapping methods, the data are used in this tutorial as a case study for application of 4DModeller to map ocean heat content. -In this tutorial we will use point data of ocean heat for upper layers of the ocean +In this tutorial, we will use point data of ocean heat for the upper layers of the ocean provided in the ME4OH data files (integrated temperature profiles from the surface down to a depth of 306.25 m, multiplied by water density and heat capacity). @@ -64,7 +64,7 @@ library(scales) ``` ## Loading Data -Here we load the Indian ocean shapefile and dataset. We smooth the boundary to make the mesh creation easier. +Here we load the Indian Ocean shapefile and dataset. We smooth the boundary to make the mesh creation easier. ```{r loadData, error=TRUE, message=FALSE, warning=FALSE} fdmr::retrieve_tutorial_data(dataset = "ohc", force_update = TRUE) @@ -82,7 +82,7 @@ fmesher::fm_crs(ocean_sf) <- NULL ocean_sf2 <- sf::st_union(sf::st_buffer(ocean_sf, 1.05)) ``` -We plot the smoothed Indian ocean polygon to check if there are any issues. +We plot the smoothed Indian Ocean polygon to check if there are any issues. ```{r polygonPlot, error=TRUE, fig.id=TRUE, fig.cap="Smoothed polygon", fig.align='center'} ggplot2::ggplot() + ggplot2::geom_sf(data=ocean_sf2) @@ -119,7 +119,7 @@ ocean_bnd <- sf::st_cast( ) ``` -We create the mesh for the Indian ocean with the outer boundary. +We create the mesh for the Indian Ocean with the outer boundary. We remove Sri Lanka and keep Madagascar in the cutoff. ```{r mesh, error=TRUE, fig.id=TRUE, fig.cap="Mesh with original polygon and all datapoints", fig.align='center'} boundary <- fmesher::fm_as_segm_list(base::list(ocean_sf2, ocean_bnd)) @@ -151,9 +151,9 @@ n_time <- base::as.integer(base::length(base::unique(df_f@data$time))) ``` ## Prior selection -To implement the SPDE approach, we define the range un uncertainty priors. We assume that the probability of process spatial range (physically, a distance where correlation between two observations falls to 0.1) being under 25 degrees in latitude is 0.2. We also set the probability to 0.01 that the marginal standard deviation of the process exceeds 1. +To implement the SPDE approach, we define the range un uncertainty priors. We assume that the probability of process spatial range (physically, a distance where the correlation between two observations falls to 0.1) being under 25 degrees in latitude is 0.2. We also set the probability to 0.01 that the marginal standard deviation of the process exceeds 1. -The process is assumed to be evolving temporally as an autoregressive process of the first order (AR1) with the probability of the temporal autocorrelation parameter to be over 0 at 0.9. +The process is assumed to be evolving temporally as an autoregressive process of the first order (AR1) with the probability of the temporal autocorrelation parameter being over 0 at 0.9. ```{r spde_setup, error=FALSE} prior_range <- 25 @@ -197,15 +197,15 @@ bru_model <- ``` # Output -The `model_viewer` of 4Dmodeller interactively shows the hyperparameters and predicted random field from the `INLA` model fit. In further sections we also plot the mean and standard deviation of the predicted field for March and October as well as the measures of temporal evolution. +The `model_viewer` of 4Dmodeller interactively shows the hyperparameters and predicted random field from the `INLA` model fit. In further sections, we also plot the mean and standard deviation of the predicted field for March and October as well as the measures of temporal evolution. ```{r modelviewer, eval=FALSE} model_viewer(model_output = bru_model, mesh = mesh, measurement_data = df_f, data_distribution = "Gaussian") ``` -To estimate posterior statistics from resampled model. We only use a small number of samples -- 100 (the default) -- to reduce the computational expanse. It is possible to use this function to also re-estimate the model fit with a new set of values. +To estimate posterior statistics from the resampled model. We only use a small number of samples -- 100 (the default) -- to reduce the computational expanse. It is possible to use this function to also re-estimate the model fit with a new set of values. -> ***NOTE:*** It takes approximately 20 min (depending on hardware) to run the prediction. You can instead run the consecutive cell to load the prediction already calculated by us. +> ***NOTE:*** It takes approximately 20 min (depending on the hardware) to run the prediction. You can instead run the consecutive cell to load the prediction already calculated by us. ```{r eval=FALSE} pred <- predict(bru_model, n.samples = 100) @@ -214,6 +214,7 @@ pred <- predict(bru_model, n.samples = 100) If you decided to run the cell above, you do not need to run the cell below. ```{r error=TRUE} pred <- fdmr::load_tutorial_data(dataset = "ohc", filename = "prediction_2005_L1.rds") +mesh <- fdmr::load_tutorial_data(dataset = "ohc", filename = "mesh.rds") ``` ## Spatial fields @@ -338,11 +339,11 @@ plot_october <- gridExtra::grid.arrange(plot_march, plot_october, ncol=2) ``` -The map reproduces main features of ocean circulation, including the Subantarctic front (around 45--50$^\circ$S where transition between blue and violet occurs) [@orsi1995acc; @giglio2016subantarctic], Agulhas Current near southern African coast and Agulhas Return Current emanating from it eastwards at 40$^\circ$S [@hood2024io]. It also shows a region of increased OHC between 12$^\circ$ and 26$^\circ$S, coinciding with the subtropical gyre [@phillips2024circulation]. The maps show decreased upper-level OHC between 5$^\circ$ and 10$^\circ$S, in particular, at the Seychelles-Chagos thermocline ridge (SCTR) [@hood2024io]. Another centre of higher top-level OHC is in Arabian Sea while OHC is relatively lower in Bay of Bengal [@huang2002eof]. +The map reproduces main features of ocean circulation, including the Subantarctic front (around 45--50$^\circ$S where a transition between blue and violet occurs) [@orsi1995acc; @giglio2016subantarctic], Agulhas Current near southern African coast and Agulhas Return Current emanating from it eastwards at 40$^\circ$S [@hood2024io]. It also shows a region of increased OHC between 12$^\circ$ and 26$^\circ$S, coinciding with the subtropical gyre [@phillips2024circulation]. The maps show decreased upper-level OHC between 5$^\circ$ and 10$^\circ$S, in particular, at the Seychelles-Chagos thermocline ridge (SCTR) [@hood2024io]. Another centre of higher top-level OHC is in the Arabian Sea while OHC is relatively lower in the Bay of Bengal [@huang2002eof]. -The differences between March and October are especially different in Arabian Sea. The western part of the basin and also at the Indian coast has a lower OHC in March compared to October. Similarly, in Bay of Bengal higher OHC appears in the central for October and western part in March. As expected, the heat content is relatively low due to weaker mixing during inter-monsoon [@dandapat2021mld]. The difference between March and October could also be explained by increased riverine freshwater influx as well as precipitation and cloudiness [@ali2018nio] in October compared to March [@sandeep2018freshwater], which causes a reduction in the OHC storage immediately at the plume area due to higher heat capcity of freshwater [@dandapat2020runoff] and also where freshwater is advected [@zhang2023bengal; @jaishree2023seasdyn] along the coastline. +The differences between March and October are especially different in the Arabian Sea. The western part of the basin and also at the Indian coast has a lower OHC in March compared to October. Similarly, in Bay of Bengal higher OHC appears in the central for October and the western part in March. As expected, the heat content is relatively low due to weaker mixing during inter-monsoon [@dandapat2021mld]. The difference between March and October could also be explained by increased riverine freshwater influx as well as precipitation and cloudiness [@ali2018nio] in October compared to March [@sandeep2018freshwater], which causes a reduction in the OHC storage immediately at the plume area due to higher heat capacity of freshwater [@dandapat2020runoff] and also where freshwater is advected [@zhang2023bengal; @jaishree2023seasdyn] along the coastline. -Similarly we show the standard deviation for March and October. +Similarly, we show the standard deviation for March and October. ```{r stdevPlots, error = FALSE, fig.id=TRUE, fig.cap="Standard deviation field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} stdev_median <- stats::median(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) stdev_max <- base::max(c(Predictor$sd[Predictor$month == 3], Predictor$sd[Predictor$month == 10]), na.rm = TRUE) @@ -430,10 +431,10 @@ plot_october <- gridExtra::grid.arrange(plot_march, plot_october, ncol=2) ``` -The standard deviation is higher in areas with low number of profiles; especially, if there are no profiles in the vicinity like during March in Great Australian Bight. Despite the number of profiles being larger in October than March, the standard deviations over the whole field remained similar. This can be explained by the distribution of those extra profiles--mostly in a close proximity to other profiles. +The standard deviation is higher in areas with a low number of profiles; especially, if there are no profiles in the vicinity like during March in the Great Australian Bight. Despite the number of profiles being larger in October than in March, the standard deviations over the whole field remained similar. This can be explained by the distribution of those extra profiles--mostly in a close proximity to other profiles. We also further explore the resampled posterior standard deviation. -```{r sterrPlots, error = FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} +```{r sterrPlots, error = FALSE, fig.id=TRUE, fig.cap="Resampled standard deviation of the mean field for L1 in 2005: (left) March, (right) October", fig.align='center', fig.width=8, fig.height=4} mc_serr_median <- stats::median(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) mc_serr_max <- base::max(c(Predictor$mean.mc_std_err[Predictor$month == 3], Predictor$mean.mc_std_err[Predictor$month == 10]), na.rm = TRUE) @@ -511,7 +512,7 @@ plot_october <- gridExtra::grid.arrange(plot_march, plot_october, ncol=2) ``` -As one can see, both standard deviation and resampled standard error of mean are highest in the same areas. However, the magnitude is much lower. +As one can see, both the standard deviation and resampled standard error of the mean are highest in the same areas. However, the magnitude is much lower. Finally, we evaluate the differences between synthetic observations and predictions at the locations of those observations. First, we calculate the predicted values at each location. ```{r error=TRUE} @@ -626,11 +627,11 @@ plot_october <- gridExtra::grid.arrange(plot_march, plot_october, ncol=2) ``` -The differences show that there is no significant bias. The regions with larger differences are typically around zones of higher mesoscale activity, for instance, Agulhas and Agulhas Return current, Leeuwin Current near Western Australian coastline, as well as western Arabian Sea (Somali Current). There is a marked increase in differences in the Arabian Sea in October compared to March, owing to a higher number of eddies after southwest monsoon in this region [@trott2017eddies]. +The differences show that there is no significant bias. The regions with larger differences are typically around zones of higher mesoscale activity, for instance, Agulhas and Agulhas Return Current, Leeuwin Current near the Western Australian coastline, as well as western Arabian Sea (Somali Current). There is a marked increase in differences in the Arabian Sea in October compared to March, owing to a higher number of eddies after the southwest monsoon in this region [@trott2017eddies]. ## Temporal evolution -Next we would like to show how the field evolves temporally. The field-mean evolution of mean and standard deviation of the posterior solution is shown in the plots below. +Next, we would like to show how the field evolves temporally. The field-mean evolution of mean and standard deviation of the posterior solution is shown in the plots below. ```{r temporalEval, error=TRUE, fig.id=TRUE, fig.cap="Area-averaged mean and standard deviation for the Indian Ocean in 2005", fig.align='center', fig.width=4, fig.height=3} meanFunc <- function(x) { @@ -663,9 +664,9 @@ ggplot2::ggplot(summ, ggplot2::aes(x = x, y = y)) + ``` -The plots show that the mean over ocean peaks between March and June and reaches minimum in September to October, in line with the southern hemisphere (austral) seasonal cycle. The uncertainties are the highest between November and May due to a lower number of measurements. +The plots show that the mean over ocean peaks between March and June and reaches a minimum in September to October, in line with the southern hemisphere (austral) seasonal cycle. The uncertainties are the highest between November and May due to a lower number of measurements. -To show the temporal evolution across the Indian Ocean basin, we also plotted March and October values of top-level DOHC along a transect. We selected the TOGA/WOCE (Tropical Ocean Global Atmosphere/World Ocean Circulation Experiment) transect IX12, which runs from Gulf of Aden to Freemantle in Western Australia. This transect is frequently (12--15 times a year) conducted by the commercial fleet that deploy expendable bathythermographs (XBTs) along the route. The data has been collected since 1983 until now. +To show the temporal evolution across the Indian Ocean basin, we also plotted March and October values of top-level DOHC along a transect. We selected the TOGA/WOCE (Tropical Ocean Global Atmosphere/World Ocean Circulation Experiment) transect IX12, which runs from the Gulf of Aden to Freemantle in Western Australia. This transect is frequently (12--15 times a year) conducted by the commercial fleet that deploys expendable bathythermographs (XBTs) along the route. The data has been collected since 1983 until now. We first define the locations of this transect and select some points for further exploration. ```{r} @@ -793,7 +794,7 @@ for (i in 1:4) { } ``` -The plot below shows the predicted DOHC for the top layer in March and October (austral autumn and spring, respectively) along the transect shown in the previous plot. The line shows the mean field predictions while envelope -- the standard deviation. +The plot below shows the predicted DOHC for the top layer in March and October (austral autumn and spring, respectively) along the transect shown in the previous plot. The line shows the mean field predictions while the envelope -- the standard deviation. ```{r transectPlot, error = FALSE, fig.id=TRUE, fig.cap="Timeseries for 2005 at selected points", fig.align='center', fig.width=6, fig.height=4} x_title <- paste0("Latitude (", intToUtf8(176), "N)") @@ -809,12 +810,12 @@ ggplot2::ggplot() + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) ``` -The figure shows an increase in DOHC northward from the southernmost point until 17$^\circ$S in the southern hemisphere. There the curve peaks for the first time and then rapidly decays until 12$^\circ$S. The curve has a local minimum between -13$^\circ$S and -5$^\circ$S, after which DOHC starts to increase again, peaking at 7$^\circ$N and leveling off slightly afterwards. The seasonal (autumn/spring) pattern shows that DOHC is higher in March than October in most places. Exception is the segment between 17$^\circ$S and 0$^\circ$, where curves match (descending limb) or where October values are higher (local minimum). +The figure shows an increase in DOHC northward from the southernmost point until 17$^\circ$S in the southern hemisphere. There the curve peaks for the first time and then rapidly decays until 12$^\circ$S. The curve has a local minimum between -13$^\circ$S and -5$^\circ$S, after which DOHC starts to increase again, peaking at 7$^\circ$N and levelling off slightly afterwards. The seasonal (autumn/spring) pattern shows that DOHC is higher in March than in October in most places. The exception is the segment between 17$^\circ$S and 0$^\circ$, where curves match (descending limb) or where October values are higher (local minimum). The whole segment is known to have a persistent cooling, which has been present throughout the year for multiple decades as shown in other studies [@roxy2020warming; @duan2023storage]. -The minimum can be explained by the equatorial upwelling (associated with enhanced vertical advection) in so-called equatorial cold tongue (this particular region known as SCTR), which is caused by negative wind stress curl due to surface wind divergence [@chen2015equatorial] between equatorial westerlies and southeasterly trade wind [@mabarrok2023thermocline]. Wind stress is regulated by the Indian monsoon due to the change in wind direction (from easterlies in boreal summer to westerlies in boreal winter) [@yokoi2008sctr]. +The minimum can be explained by the equatorial upwelling (associated with enhanced vertical advection) in the so-called equatorial cold tongue (this particular region is known as SCTR), which is caused by negative wind stress curl due to surface wind divergence [@chen2015equatorial] between equatorial westerlies and southeasterly trade wind [@mabarrok2023thermocline]. Wind stress is regulated by the Indian monsoon due to the change in wind direction (from easterlies in boreal summer to westerlies in boreal winter) [@yokoi2008sctr]. The segment where both curves match is shown by @mckenna2024sst to have a higher net heat flux in March while net advection brings more heat in October due to negative vertical advection in March (associated with the same process as for the "cold tongue") while the meridional heat advection is stronger in October. -Elsewhere, March OHC has higher than October OHC. In southern hemisphere, it is expected due to increased insolation during austral summer and the lag for the ocean to heat up. However, we also observe higher OHC in northern hemisphere, where maximum values occur at the end of boreal winter. This segment has similar net heat flux values because it peaks twice during the year and both peaks are of similar magnitude [@mckenna2024sst]. However, advection contributes to more heat loss in October than March. -The final (northernmost) segment lies in close proximity to land and is affected by coastal upwelling near Somali coast [@vic2014mesoscale; @chatterjee2019somali]. +Elsewhere, March OHC is higher than October OHC. In the southern hemisphere, it is expected due to increased insolation during the austral summer and the lag for the ocean to heat up. However, we also observe higher OHC in the northern hemisphere, where maximum values occur at the end of boreal winter. This segment has similar net heat flux values because it peaks twice during the year and both peaks are of similar magnitude [@mckenna2024sst]. However, advection contributes to more heat loss in October than in March. +The final (northernmost) segment lies close to land and is affected by coastal upwelling near the Somali coast [@vic2014mesoscale; @chatterjee2019somali]. The panel below shows the mean and standard deviation at the four sample locations along the transect shown on the map earlier. Compared to the field mean, we can observe that temporal evolution is spatially variable ```{r timeseriesPlots, error = FALSE, fig.id=TRUE, fig.cap="Timeseries for 2005 at selected points", fig.align='center', fig.width=6, fig.height=6} @@ -831,14 +832,14 @@ ggplot2::ggplot(do.call(rbind, point_eval)) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) ``` -The northernmost point (1) has the highest DOHC values on average. There is a pronounced seasonal cycle, which has a maximum plateau (within standard deviation) from March to September and more pronounced trough reaching minimum in November--December. The plateau coincides with increased insolation when the sun is in zenith over the northern hemisphere. However, it does not produce the typical biannual maximum. The area is subject to monsoon winds, which cause seasonal Somali Current reversal (southward during boreal summer) and creation of a small gyre (Great Whirl) [@lhegaret2018gyres]. The winds reverse from southwesterly (JJA) to northeasterly (DJF), causing enhanced boreal winter surface cooling in the northern Arabian Sea due to the reduction in latent heat flux when wind brings dry and cool continental air from Asia [@schott2009circulation]. This suggests that an interplay between insolation and wind-driven surface cooling is likely causing the seasonal changes in this location. +The northernmost point (1) has the highest DOHC values on average. There is a pronounced seasonal cycle, which has a maximum plateau (within standard deviation) from March to September and a more pronounced trough reaching a minimum in November--December. The plateau coincides with increased insolation when the sun is in zenith over the northern hemisphere. However, it does not produce the typical biannual maximum. The area is subject to monsoon winds, which cause seasonal Somali Current reversal (southward during boreal summer) and the creation of a small gyre (Great Whirl) [@lhegaret2018gyres]. The winds reverse from southwesterly (JJA) to northeasterly (DJF), causing enhanced boreal winter surface cooling in the northern Arabian Sea due to the reduction in latent heat flux when the wind brings dry and cool continental air from Asia [@schott2009circulation]. This suggests that an interplay between insolation and wind-driven surface cooling is likely causing the seasonal changes in this location. The standard deviation envelope narrows slightly at the descending limb, likely due to a higher number of observations in the area during that time. -Point 2 has a relatively even profile with lower average than at point 1. There is a peak in June that levels off over a couple of months until August. This could partly be caused by weaker equatorial upwelling as well as increased net heat flux as discussed earlier. This location is also near the South Equatorial Counter current, which is not present at this particular point during the boreal summer [@wu2020gyrestructure]. Thus, it is less affected by negative zonal heat advection when the top-layer OHC values peak. The standard deviation narrows near the peak, but is similar to point 1 elsewhere. +Point 2 has a relatively even profile with a lower average than at point 1. There is a peak in June that levels off over a couple of months until August. This could partly be caused by weaker equatorial upwelling as well as increased net heat flux as discussed earlier. This location is also near the South Equatorial Counter current, which is not present at this particular point during the boreal summer [@wu2020gyrestructure]. Thus, it is less affected by negative zonal heat advection when the top-layer OHC values peak. The standard deviation narrows near the peak, but is similar to point 1 elsewhere. -Point 3 are another relatively flat timeseries, which has a small peak in January, leveling off until April. Within a standard deviation, the rest of values are similar. South Equatorial current is present in the vicinity of point 3 and provides some heat advection. @mckenna2024sst shows that zonal and meridional advection dampens the heat loss due to the net negative heat flux during austral winter. The standard deviation is similar for all timesteps. +Point 3 are another relatively flat timeseries, which has a small peak in January, levelling off until April. Within a standard deviation, the rest of the values are similar. South Equatorial current is present in the vicinity of point 3 and provides some heat advection. @mckenna2024sst shows that zonal and meridional advection dampens the heat loss due to the net negative heat flux during austral winter. The standard deviation is similar for all timesteps. -The southernmost point (4) mostly has the lowest DOHC values on average, which can be related to less solar radiation at higher latitudes. There is also a prominent seasonal cycle with a maximum in February and minimum in November. The standard deviation is similar for all months and also of similar magnitude as point 3. +The southernmost point (4) mostly has the lowest DOHC values on average, which can be related to less solar radiation at higher latitudes. There is also a prominent seasonal cycle with a maximum in February and a minimum in November. The standard deviation is similar for all months and also of similar magnitude as point 3. ## References {-}