diff --git a/DESCRIPTION b/DESCRIPTION
index ebb64510c14507414eca7f7d56277461ebb689cb..d6cc52d69cffce3ba11bbfa18a864c1d43160a77 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: SWATdoctR
 Type: Package
 Title: Finding the right diagnoses and treatments for SWAT+ models
-Version: 0.1.9
+Version: 0.1.11
 Author: c(person("Svajunas", "Plunge",
              email = "svajunas_plunge@sggw.edu.pl",
              role = c("aut")),
diff --git a/NAMESPACE b/NAMESPACE
index b46f35a6376df0ff96f8c73d656a5f8c5bdf729b..07f1319b0f1581ae4fe866fddefb676f0629584c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,6 +14,7 @@ export(plot_ps)
 export(plot_qtile)
 export(plot_variable_at_harvkill)
 export(plot_water_partition)
+export(plot_waterbalance)
 export(print_avannual_qtile)
 export(print_triggered_mgt)
 export(report_mgt)
@@ -79,6 +80,7 @@ importFrom(plotly,layout)
 importFrom(plotly,plot_ly)
 importFrom(plotly,plotly_build)
 importFrom(plotly,subplot)
+importFrom(png,readPNG)
 importFrom(processx,run)
 importFrom(purrr,list_c)
 importFrom(purrr,list_rbind)
diff --git a/R/plot_waterbalance.R b/R/plot_waterbalance.R
new file mode 100644
index 0000000000000000000000000000000000000000..62b7f67a491b49bb6cb892df940068a3bf7b7fdf
--- /dev/null
+++ b/R/plot_waterbalance.R
@@ -0,0 +1,105 @@
+#' Plot the average annual water balance
+#'
+#' This function uses average annual outputs from basin_wb_aa and
+#' basin_aqu_aa and plot the simulated values of water balance components
+#' in a flow chart.
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}
+#'   as a single object.
+#'   To plot the water balance at least the output option \code{outputs = 'wb'}
+#'   must be set in \code{run_swat_verification()}.
+#'
+#' @return Returns a ggplot flow chart with the simulated average annual water
+#'   balance.
+#'
+#' @importFrom dplyr %>% distinct
+#' @importFrom png readPNG
+#' @import ggplot2
+#' @export
+#'
+plot_waterbalance <- function(sim_verify) {
+  img <- readPNG(paste0(system.file(package = "SWATdoctR"), '/extdata/swatplus_wb.png'))
+
+  wb_aa  <- round(unlist(sim_verify$basin_wb_aa[,c(8:12, 14:18, 20:23, 27:32, 34:39, 43)]), 2)
+  aqu_aa <- round(unlist(sim_verify$basin_aqu_aa[,c(8:10, 12:13, 22:24)]), 2)
+
+  ratios <- round(c(wb_aa['et'] / wb_aa['precip'],
+                    (wb_aa['surq_gen'] + wb_aa['latq'] + wb_aa['qtile'] + aqu_aa['flo']) / wb_aa['precip'],
+                    (wb_aa['surq_cha'] + wb_aa['latq_cha'] + wb_aa['qtile'] + aqu_aa['flo_cha']) / wb_aa['precip'],
+                    wb_aa['surq_gen'] + wb_aa['latq'] + wb_aa['qtile'] + aqu_aa['flo'],
+                    wb_aa['surq_cha'] + wb_aa['latq_cha'] + wb_aa['qtile'] + aqu_aa['flo_cha'],
+                    aqu_aa['flo'] / (wb_aa['surq_gen'] + wb_aa['latq'] + wb_aa['qtile'] + aqu_aa['flo'])
+  ),
+  2)
+
+  val <- c('', '', wb_aa, aqu_aa, '', '', '', '', '', '', '', ratios)
+
+  x <- c(0, 11.81,
+         2.25, 2.78, 4.69, 7.16, 7.16,
+         6.73, 0.90, 1.30, 0.90, 0.52,
+         3.80, 3.30, 3.30, 3.30, 2.60,
+         3.80, 7.16, 4.69, 0.12, 0.12,
+         9.90, 9.90, 9.90, 9.90, 9.90,
+         9.90, 3.30,
+         7.16, 0.35, 0.35, 4.44, 2.15,
+         9.90, 9.90, 9.90,
+         7.16, 7.16, 7.16, 7.16, 7.16,
+         7.16, 7.16,
+         9.90, 9.90, 9.90, 10.75, 10.75,
+         9.90)
+  y <- c(   0, 10.00,
+            8.34, 5.80, 5.12, 5.12, 3.48,
+            1.87, 7.80, 6.52, 5.68, 4.85,
+            8.37, 4.30, 3.80, 3.30, 5.12,
+            8.87, 2.85, 6.35, 9.70, 2.83,
+            5.68, 5.15, 4.61, 3.98, 3.45,
+            2.95, 2.80,
+            1.58, 2.08, 1.58, 0.43, 2.22,
+            2.12, 1.60, 1.07,
+            9.70, 9.20, 8.70, 8.20, 7.70,
+            6.90, 6.20,
+            9.20, 8.70, 8.20, 6.90, 6.20,
+            7.70)
+  a <- c( 0, 0,
+          0,   90,    0,    0,    0,
+          90,   90,   90,   90,   90,
+          0,    0,    0,    0,    0,
+          0,    0,    0,    0,    0,
+          0,    0,    0,    0,    0,
+          0,    0,
+          0,    0,    0,   90,   90,
+          0,    0,    0,
+          0,    0,    0,    0,    0,
+          0,    0,
+          0,    0,    0,    0,    0,
+          0)
+  l <- c('','',
+         'precip:\n', 'sno_fall: ', 'sno_melt: ', 'surq_gen: ', 'latq: ',
+         'perc: ', 'et: ', 'ecanopy: ', 'eplant: ', 'esoil: ',
+         'cn: ', 'sw_init: ', 'sw_final: ', 'sw_ave: ', 'sno_pack: ',
+         'pet: ', 'qtile: ', 'irr: ', 'surq_ro: ', 'latq_ro: ',
+         'surq_cha: ', 'surq_res: ', 'surq_ls: ', 'latq_cha: ', 'latq_res: ',
+         'latq_ls: ', 'sw_change: ',
+         'flo: ', 'dep_wt: ', 'stor: ', 'seep: ', 'revap: ',
+         'flo_cha: ', 'flo_res: ', 'flo_ls: ',
+         'Water balance ratios: ', 'et / precip: ', 'wyld / precip :', 'wyld_cha / precip :', 'flo / wyld: ',
+         'wyld = surq_gen + latq + qtile + flo :',
+         'wyld_cha = surq_cha + latq_cha +\n                      qtile + flo_cha :',
+         '', '', '', '', '',
+         ''
+  )
+
+  l <- paste0(l, val)
+
+  ggplot() +
+    annotation_raster(img, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, interpolate = FALSE) +
+    geom_text(aes(x = x,
+                  y = y,
+                  label = l,
+                  angle = a), vjust = 0, hjust = 0) +
+    scale_x_continuous(expand = c(0,0)) +
+    scale_y_continuous(expand = c(0,0)) +
+    coord_equal() +
+    geom_point() +
+    theme_void()
+}
diff --git a/R/run_swat_verify.R b/R/run_swat_verify.R
index ec265e271a8c2adfba4f068c4f0dfd58c7be344a..78634bc4b84d4664cdf3d5d5057b9a7ba8c4f960 100644
--- a/R/run_swat_verify.R
+++ b/R/run_swat_verify.R
@@ -73,8 +73,11 @@ run_swat_verification <- function(project_path, outputs = c('wb', 'mgt', 'plt'),
       model_output$hru_pw_day <- read_tbl('hru_pw_day.txt', run_path, 3) %>% lwr
     }
     if ('wb' %in% outputs) {
-      model_output$basin_wb_day <- read_tbl('basin_wb_day.txt', run_path, 3) %>% lwr
-      model_output$basin_pw_day <- read_tbl('basin_pw_day.txt', run_path, 3) %>% lwr
+      model_output$basin_wb_day  <- read_tbl('basin_wb_day.txt', run_path, 3) %>% lwr
+      model_output$basin_pw_day  <- read_tbl('basin_pw_day.txt', run_path, 3) %>% lwr
+      model_output$basin_wb_aa   <- read_wb_aa(run_path) %>% lwr
+      model_output$basin_aqu_aa  <- read_tbl('basin_aqu_aa.txt',  run_path, 3) %>% lwr
+      # model_output$basin_cha_aa  <- read_tbl('basin_sd_cha_aa.txt',  run_path, 3) %>% lwr
       model_output$hru_wb_aa <- read_tbl('hru_wb_aa.txt', run_path, 3) %>% lwr
       tryCatch({
         model_output$recall_yr <- read_tbl('recall_yr.txt', run_path, 3) %>% lwr
diff --git a/R/update_input_files.R b/R/update_input_files.R
index ba2733a8a6d7aad95f81a9deae9b11872cb1c67f..e038cd1fd3283a7f6f9cc74c3e3c4106a9964101 100644
--- a/R/update_input_files.R
+++ b/R/update_input_files.R
@@ -47,8 +47,10 @@ set_print_prt <- function(project_path, run_path, outputs, years_skip) {
   }
 
   if ('wb' %in% outputs) {
-    print_prt[11] <- "basin_wb                     y             n             n             n  "
+    print_prt[11] <- "basin_wb                     y             n             n             y  "
     print_prt[14] <- "basin_pw                     y             n             n             n  "
+    print_prt[15] <- "basin_aqu                    n             n             n             y  "
+    print_prt[18] <- "basin_sd_cha                 n             n             n             y  "
     print_prt[33] <- "hru_wb                       n             n             n             y  "
     print_prt[45] <- "recall                       n             n             y             y  "
   }
diff --git a/inst/extdata/swatplus_wb.png b/inst/extdata/swatplus_wb.png
new file mode 100644
index 0000000000000000000000000000000000000000..5cd774ac8b18a00f9e34697272276cde1058232a
Binary files /dev/null and b/inst/extdata/swatplus_wb.png differ
diff --git a/man/plot_waterbalance.Rd b/man/plot_waterbalance.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..a4b0b8a9f9f46e8bd980166d0bb529ad44dfa548
--- /dev/null
+++ b/man/plot_waterbalance.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_waterbalance.R
+\name{plot_waterbalance}
+\alias{plot_waterbalance}
+\title{Plot the average annual water balance}
+\usage{
+plot_waterbalance(sim_verify)
+}
+\arguments{
+\item{sim_verify}{Simulation output of the function \code{run_swat_verification()}
+as a single object.
+To plot the water balance at least the output option \code{outputs = 'wb'}
+must be set in \code{run_swat_verification()}.}
+}
+\value{
+Returns a ggplot flow chart with the simulated average annual water
+  balance.
+}
+\description{
+This function uses average annual outputs from basin_wb_aa and
+basin_aqu_aa and plot the simulated values of water balance components
+in a flow chart.
+}