# Introduction to climate stress testing using foreSIGHT

#### Anjana Devanand, David McInerney, Sam Culley, Bree Bennett, Seth Westra

#### 14 Sep 2025

# 1. Introduction

A variable and changing climate presents significant challenges to the functioning and/or performance of both natural and engineered systems. Managed systems—both engineered and managed natural systems—traditionally have been designed under the assumption that future climate conditions will mirror those experienced in the past. Yet with the continuing advance of climate change, there is a need to understand how systems might perform under a range of plausible future climate conditions or, conversely, what system interventions might be required so that systems continue to achieve desired levels of performance. Given the complexity of most climate-sensitive systems, formalised approaches are required to understand likely climate impacts and evaluate the viability of adaptive measures to minimise climate vulnerability.

To this end, scenario-neutral (or ‘bottom-up’) approaches (Prudhomme et al. 2010, Brown 2011, Culley et al. 2016) are advocated as a means of rigorously stress testing a system under a range of plausible future climate conditions. These approaches treat the system’s behaviour and performance as the central concerns of the analysis, and enable better understanding of the complex climate-system relationships to support adaptation decision making. These approaches can be combined with ‘top-down’ climate impact assessment methods through the integration of projections from climate models and/or other lines of evidence. The foreSIGHT package contains functions that support both ‘bottom-up’ system stress testing, and the analysis of the implication of ‘top-down’ climate projections on system performance.

This vignette demonstrates the high-level work flow of climate ‘stress-testing’ using the foreSIGHT package. The document is intended as a quick start guide and uses a simplified case study of a rainwater tank system model to demonstrate how the functions in the package can be used together. The focus is on accessible methods for generating climate scenarios, though these simplified approaches have limitations in capturing certain features of a changing climate—such as extremes or wet–dry day patterns. For more flexible (but more complex) simulation options, consult the foreSIGHT Guide to Stochastic Simulation vignette.

## 1.1. Objectives and application areas of foreSIGHT

The objectives of foreSIGHT are to support climate impact and vulnerability assessments and the assessment of adaptation options by:

- stress testing climate-sensitive systems, including both ‘current’ system configurations, as well as potential alternative system configurations as part of the development of adaptation strategies;
- comparing the climate sensitivity of multiple system configurations to inform adaptation decision making; and
- comparing stress-testing outcomes with the results from ‘top-down’ climate impact assessments to better understand future risk for each system configuration.

The foreSIGHT modelling software adopts a rigorous quantitative approach to stress testing, that has been designed with several core assumptions in mind:

- that the system dynamics (either ‘current’ or alternative system configurations) can be represented and adequately described by a numerical system model that provides a mapping between weather/climate variables and relevant system performance measures; and
- that the system model is forced by hydroclimatic time series data.

Indeed, it is this latter feature that gives the software its name (the SIGHT in foreSIGHT stands for System Insights from the Generation of Hydroclimatic Timeseries). In particular, foreSIGHT has been designed specifically for the quantitative analysis of systems that exhibit dynamics in time, with examples of such systems including:

- environmental systems (either natural or managed) that may be resilient to individual natural hazards but become vulnerable to multiple sequential hazards or long-term structural shifts in the climate;
- water resource systems with natural (e.g. soil moisture, groundwater) and/or human-constructed (e.g. reservoirs, managed aquifer recharge) storages, for which past weather can affect current system performance;
- agricultural systems where crop outcomes (e.g. yield and various quality measures) are influenced by the weather throughout a growing season and even between seasons;
- renewable energy systems such as solar, wind and hydroelectricity and/or coupled storage solutions (e.g. pumped hydroelectricity or lithium battery systems); and
- systems that depend on one or several of the above systems, such as mining (often dependent on groundwater and/or surface water reserves), transportation (often sensitive to flooding and various other natural hazards), tourism (often highly dependent on ecosystem health) and so forth.

The focus on detailed numerical modelling and system ‘stress testing’ highlights that foreSIGHT is particularly suited to situations where the consequences of system performance degradation and/or failure as a result of climate change are likely to be significant, as well as for quantitative decision making and/or engineering design. It is assumed that a high-level (qualitative) risk assessment would have already been conducted and the outcome of that assessment is that a detailed quantitative analysis is required.

## 1.2. foreSIGHT workflow for climate stress testing

The foreSIGHT workflow is shown the diagram below, and comprises five distinct steps that are further described in the paragraphs below. A core aspect of the foreSIGHT functionality is to evaluate how the system performs under a range of plausible climate scenarios created by perturbing statistical properties of observed climate time series. The workflow involves the steps shown in the following diagram, each of which are discussed in the case study presented in Section 2. As highlighted in the previous section, at this point it is assumed that a detailed quantitative analysis of a system is required (based, for example, on the outcomes of a qualitative risk assessment) and that a numerical system model is available.

Workflow of climate ‘stress-testing’ using foreSIGHT

Step A. The process of system stress testing involves assessing how a system’s behaviour (including its ‘function’ or ‘performance’) varies as a result of plausible climatic changes. These changes are described by means of climate “attributes”, which we define as statistical measures of weather variables. Examples of attributes are annual total rainfall, annual number of wet days, and annual average temperature. In this step, the attributes that are deemed to be most relevant for a particular system are identified. These attributes are generally selected based on a priori understanding of system dynamics and likely system vulnerability. The maximum-minimum bounds of the perturbations in the selected attribute, and the type of sampling within this range are also decided. The attributes and perturbations are used to create an exposure space. The outcome of this step is a set of sampled points within an exposure space, that provide the ‘targets’ for time series generation algorithms in Step B.

Step B. This step involves generation of perturbed time series corresponding to the sampled points of target perturbations created in Step A. A reference (typically observed) time series of the relevant hydro-climate variables is required to create the perturbed time series using a selected method of perturbation.

The supported perturbation methods in foreSIGHT include

- application of scaling factors to the supplied time series, which is described in this vignette.

- use of the ‘inverse method’ of Guo et al (2018) to optimise the parameters of stochastic weather generator type models to generate time series with desired perturbed attributes. If stochastic models are used for time series generation, multiple replicates of time series that correspond to the same target can be generated to better represent stochastic (‘natural’) variability. The application of these methods is described in the foreSIGHT Guide to Stochastic Simulation vignette.

The outcome of this step is a set of perturbed time series that correspond as closely as possible to each point in the exposure space.

Step C. The perturbed time series generated in Step B are used to drive the system model and simulate system ‘performance’. The performance metrics should represent measures that are most relevant to the system under consideration, and can include a variety of economic, social and/or environmental measures. It is assumed that the performance metrics are calculated within the system model and thus represent the outputs from that model (i.e. the foreSIGHT package does not calculate the performance metrics itself). The outcome of this step is a quantitative representation of how system performance varies across the exposure space.

Step D. This step visualises the system performance metrics calculated in step C to understand the system responses to the perturbations in the selected climate attributes. If minimum or maximum threshold criteria of the performance metrics are defined, these thresholds can be used to identify instances of unsatisfactory system performance/system failure. In this step, the performance metrics are visualised in the perturbation space of the climate attributes, i,e, the axes used for visualisation are the perturbed climate attributes. Such figures are henceforth named ‘performance spaces’—these visualisations enable identification of combinations of perturbations that result in changes to system performance. In cases where the ‘stress-test’ includes multiple perturbed attributes and performance metrics, multiple visualisations of performance spaces are used to represent all combinations of attributes/metrics. If alternate climate information is available from other sources of evidence (eg: ‘top-down’ approaches), they can be superimposed on the visualisations generated in this step. Inclusion of this additional climate data may provide information about the plausibility of the perturbations in the attributes. The outcome of this step are plots of the system performance spaces/thresholds and understanding of the system responses to the climate perturbations.

Step E. This step involves analysis of alternate system configurations/policies in order to support decision making. Visualisations are created for the alternate system choices to compare their performance. The outcome of this step are plots of the performance spaces/thresholds for all system choices and understanding of the preferred choices under climate perturbations.

These five steps complete the framework of climate impact assessment using foreSIGHT.

Several parts of the foreSIGHT workflow can involve significant computational effort and thus runtimes. This is particularly relevant for Step B (generation time series that represent the target attributes) and Step C (simulating the time series through the system model), and depends on the approach to perturbation and selected system model, respectively. In relation to perturbation methods, there are multiple ways to perturb historical weather time series, from application of additive or multiplicative factors to historical weather data through to the use of sophisticated stochastic generators. To support fast run times, in this example we adopt scaling approaches of adding/subtracting constants to historical temperature and multiplying precipitation by percentages relative to historical. However, these methods has the following limitations:

- Some statistical properties such as the rainfall wet-dry patterns or extremes cannot be perturbed
- Multiple attributes cannot be perturbed in combination
- It is not possible to hold some desired attributes at historical levels while perturbing others
- The length of the generated time series cannot be longer than the supplied reference time series

The implication of the above is that key modes of system vulnerability may be missed by adopting ‘scaling’ approaches. The user should consider these points while using foreSIGHT for other case studies, with further description of the use of more sophisticated approaches for time series perturbation based on stochastic weather generators provided in the foreSIGHT Guide to Stochastic Simulation vignette.

Given the diversity of possible use-cases and heterogeneity of system models, foreSIGHT itself does not have system modelling capabilities, but instead is designed to flexibly integrate with a range of system models that are either written in R, or written in a different language but can be coupled to R via wrapper functions. Details on system model coupling to foreSIGHT are provided below.

To explore and illustrate key elements of foreSIGHT functionality, this vignette uses as an inbuilt rainwater tank system model is provided as part of the software package as described next.

# 2. Case Study - Climate ‘Stress-test’ of a Rainwater Tank System

The section demonstrates how foreSIGHT can be used to climate ‘stress-test’ the rain water tank system model under a range of plausible climate scenarios created by perturbing the attributes of the observed climate time series. This is achieved through the application of a framework involving the steps shown in the workflow in section 1.2.

The system model is a representation of a domestic rainwater tank system, which has been designed to meet both indoor (grey water) and outdoor (garden irrigation) water demands. Although this system model example is simpler than anticipated real-world usages of the foreSIGHT model, it nevertheless provides important insights associated with system sensitivities, the role of temporal dynamics and the behaviour of storages, the interaction between supply and demand, and the identification and comparison of multiple system configurations. The core functionality of this model is now described.

Schematic of the system model

A schematic representation of the rainwater tank system model is shown in the figure above. Rain falling on the roof of a house is captured and directed towards the rainwater tank. Before the rainwater is able to enter the tank, a depth of water (called the first flush) is removed from the start of each storm for water quality reasons. The water remaining after the first flush extraction flows into the rainwater tank. The amount of water supplied by the tank is calculated based on the water level in the tank. The indoor and outdoor water demands deplete the water stored in the tank. The indoor water demand is assumed to be constant throughout the year, and the outdoor water demand varies seasonally and temperature. The operation of the rain water tank system model is thus dependent upon the climate variables rainfall and temperature. The tank model simulates rainwater capture and water use processes at a daily time step using rainfall and temperature time series as input. The parameters of the model that the user should specify are:

- Area of the roof used for rain water harvesting

- Volume of the tank

- Number of people using the water

- Depth of water removed as the first flush.

These parameters can be varied for alternate system designs.

The system model estimates the performance of the rainwater tank using a range of metrics, including:

- Average Daily Deficit - the volume of average deficit in water supplied by the tank in litres
- Reliability - the fraction of days on which the full demand could be supplied

This example system model provides sufficient scope for climate stress testing using foreSIGHT. This is because the tank responds to multiple climate drivers (i.e. rainfall and temperature), and the removal of the first flush at the start of the storm means that the wet-dry pattern of the rainfall and the seasonality of the demand pattern may become important in the functioning of the tank. The system model is available as the `tankWrapper()` function in foreSIGHT. The performance metrics available in the tank model can be viewed using the `viewTankMetrics()` function.

## 2.1. Step A: Identify attributes for perturbation and create an exposure space

The function `createExpSpace()` can be used to create an exposure space. The term ‘exposure space’ refers to the perturbed values of the attributes of weather time-series that form the basis of stress testing, representing the ‘exposure’ element of a traditional climate impact assessment. For example, an exposure space may consist of multiplicative perturbations of 0.8 to 1 in annual total precipitation, and shifts of 0-1oC in average temperature.

To use the function `createExpSpace()`, the user should provide information on the following:

- which attributes should be perturbed
- the sampling ranges (minima and maxima) number and either the number of samples or sampling increment of the perturbed attributes
- type of perturbation

These are discussed in turn.

Choice of which attributes to perturb. The selection of attributes to perturb is given by the argument `attPerturb`. The choice of which attributes to select for perturbation depends upon the anticipated sensitivity of the system. Furthermore, is it possible to perturb the attributes using the available perturbation method? This is relevant as part of the choice between ‘scaling’ methods (the focus of this tutorial) and stochastic methods (covered in the Guide to Stochastic Simulation). Note that for stochastic simulation, we also need to hold/tie attributes to constrain the numerical optimization and ensure realistic time series, but this is not necessary for scaling approaches.

Number of ranges and number of samples/increments. For each perturbed attribute, the minimum-maximum values of perturbation (`attPerturbMin`, `attPerturbMax`) should be specified, as well as either the number of samples to be generated (`attPerturbSamp`) or the sampling increment (`attPerturbBy`). By default, perturbations in the attributes of variables like rainfall are specified in multiplicative units, while the perturbations in the attributes of variables like temperature are specified in additive units. The function `createExpSpace()` uses this information to create equidistant perturbation points between the minimum-maximum values. For example, if `attPerturbMin` of an attribute of rainfall is specified as 0.9, `attPerturbMax` as 1.1, and `attPerturbSamp` as 3, the perturbation points of the attribute will be 0.9, 1, and 1.1.

Sampling strategy. The function argument `attPerturbType` can be used to specify the type sampling of the exposure space. Please refer to function documentation using `?createExpSpace` to check the options available in the package. Two common sampling approaches are ‘one-at-a-time’ (OAT), and ‘regGrid’ sampling. In the case of a ‘OAT’ sampling, each attribute is perturbed one-at-a-time while holding all other attributes constant. In contrast, in a ‘regGrid’ sampling the attributes are perturbed simultaneously to create an exposure space encompassing all combination of perturbations in the selected attributes.

Attributes in foreSIGHT are specified as combinations of variables names (e.g. ‘P’ for precipitation), aggregation time step (e.g. ‘day’ for daily), stratification type (e.g. ‘all’ for all months/seasons or DJF for December-February), function name (e.g. ‘tot’ for total), and operation (e.g. ‘m’ for mean over years).

In this tutorial we focus on attributes of the following form:

- `P_day_all_tot_m` is the mean annual total rainfall

- `P_day_all_seasRatioMayAug` is the ratio of May-August rainfall to annual rainfall.

- `Temp_day_all_avg` is the average daily temperature

The definitions of climate attributes supported by foreSIGHT can be viewed using the helper function `viewAttributeDef()` available in the package.

    viewAttributeDef("P_day_all_tot_m")
    #> [1] "Mean annual rainfall"
    viewAttributeDef("P_day_all_seasRatioMayAug")
    #> [1] "Ratio of May-Aug to total rainfall"
    viewAttributeDef("Temp_day_all_avg")
    #> [1] "Average daily temperature"

Please refer to the Guide to Stochastic Simulation vignette for further details on the naming conventions for attributes.

The combinations of perturbed attributes specified by the user provide the intended sampling of the exposure space. The function `createExpSpace()` can be used to create an exposure space by specifying `attPerturb` and relevant information regarding the sampling strategy for each attribute. A typical usage of this function is shown below.

In the below function call, the attributes are perturbed one-at-a-time (OAT) as follows:

- `P_day_all_tot_m` is perturbed from 70% to 110% of the historical average (i.e. 30% reductions to 10% increase) using 9 samples between the sampling bounds

- `P_day_all_seasRatioMarAug` is perturbed from 90% to 120% of the historical average using 7 samples

- `Temp_day_all_avg` is perturbed from 0oC to 1.5oC relative to the historical average using 7 samples

An example whereby the `P_day_all_tot_m` and `P_day_all_seasRatioMarAug` attributes are perturbed jointly using ‘regGrid’ sampling is:

    # specify perturbed attributes
    attPerturb <- c("P_day_all_tot_m", "P_day_all_seasRatioMarAug") # specify perturbation type and minimum-maximum ranges of the perturbed attributes
    attPerturbType <- "regGrid"
    attPerturbSamp <- c(5, 5)
    attPerturbMin <- c(0.7, 0.9)
    attPerturbMax <- c(1., 1.2) # create the exposure space
    expSpace.regGrid <- createExpSpace(
      attPerturb = attPerturb,
      attPerturbSamp = attPerturbSamp,
      attPerturbMin = attPerturbMin,
      attPerturbMax = attPerturbMax,
      attPerturbType = attPerturbType
    )
    #> Note: There are no attributes held at historical levels
    plotExpSpace(expSpace.regGrid)

Target locations in the exposure space

The function `plotExpSpace()` above is used to plot and check the range and sampling resolution of a created 2D exposure space.

Each point in this exposure space is what is referred to as a ‘target’ location, in the sense that it describes the desired attribute values of the time series to be generated in Step B to inform the system stress test.

## 2.2. Step B: Generate perturbed time series

Having created the exposure space and the target locations, the next step is to create hydro-climate time series with these desired attribute values. The function `generateScenarios()` can be used to generate perturbed time series corresponding to all the target locations in the exposure space.

Reference climate data format

The first step is to supply the reference climate data in the appropriate format, which is a list with data for each climate variable, and the associated ‘times’.

Daily observed precipitation and temperature over the period from 2007 to 2016 obtained by combining data from multiple station locations to represent the general climate of Adelaide, South Australia is included in the demonstration, and may be be loaded using the data command.

    data("tankDat")
    lapply(tank_obs, head)
    #> $P
    #> [1] 0.0 0.0 0.0 0.0 0.0 4.5
    #>
    #> $Temp
    #> [1] 25.50 24.50 29.75 32.25 32.50 26.50
    #>
    #> $times
    #> [1] "2007-01-01 UTC" "2007-01-02 UTC" "2007-01-03 UTC" "2007-01-04 UTC"
    #> [5] "2007-01-05 UTC" "2007-01-06 UTC"

We note that foreSIGHT can use both single site and multisite data. For multisite data, the climate variable is a matrix, with columns representing the different sites, and rows the times. An example of multisite data is shown below:

    data("barossaDat")
    lapply(barossa_obs, head)
    #> $P
    #>      X23300 X23302 X23305 X23309 X23312 X23313 X23317 X23318 X23321 X23363
    #> [1,]   24.6   24.4   16.8   25.4   17.3   33.8   17.2   20.5   17.3   28.8
    #> [2,]    2.5    3.8    1.4    3.0    1.8    5.6    1.7    1.6    1.8    4.1
    #> [3,]    0.0    0.3    0.8    0.0    0.0    0.5    0.0    0.0    0.0    0.2
    #> [4,]    1.3    2.0    0.0    0.0    0.0    2.0    0.0    0.3    0.0    1.3
    #> [5,]    0.0    0.0    0.0    0.0    0.0    0.8    0.0    0.0    0.0    0.4
    #> [6,]    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
    #>      X23373 X23752 X23756
    #> [1,]   17.3   26.7   31.2
    #> [2,]    1.8    4.4    5.1
    #> [3,]    0.0    0.0    0.3
    #> [4,]    0.0    1.0    1.8
    #> [5,]    0.0    0.0    0.3
    #> [6,]    0.0    0.0    0.0
    #>
    #> $times
    #> [1] "1970-01-01 UTC" "1970-01-02 UTC" "1970-01-03 UTC" "1970-01-04 UTC"
    #> [5] "1970-01-05 UTC" "1970-01-06 UTC"

The variable `times` can have a range of time steps from daily to annual, and uses the POSIXct date-time class. `times` can be created as follows:

    # format required for "times" in reference climate data
    # specify timezone as 'UTC' (Coordinated Universal Time) to avoid issues with daylight savings
    timeStart <- as.POSIXct("2007-01-01", tz = "UTC")
    timeEnd <- as.POSIXct("2016-12-31", tz = "UTC")
    times <- seq(timeStart, timeEnd, by = "days") # use a daily time step

Calculating climate attributes for reference data

foreSIGHT contains a function that can be used to calculate the attributes of interest for climate data supplied by the user: `calculateAttributes()`. The usage of this function is shown below. The function is intended for use with the reference data:

    attSel <- c("P_day_all_tot_m", "P_day_all_seasRatioMarAug", "Temp_day_all_avg")
    tank_obs_atts <- calculateAttributes(tank_obs, attSel)
    tank_obs_atts
    #>           P_day_all_tot_m P_day_all_seasRatioMarAug          Temp_day_all_avg
    #>               449.9300000                 0.6566799                17.4369012

Note that `calculateAttributes()` can be used to calculate attributes from climate model projections, which can be used to inform sampling bounds (Step A) and projecting alternative lines of evidence on performance visualization (Step D and E).

Generating scenarios

The next step is to specify which technique to use for perturbation, which include: (1) simple/seasonal scaling, which applies an additive or multiplicative change to the user-supplied observed time series to create perturbed time series. Simple scaling applies the same factors throughout the year, while seasonal scaling uses different factors in different seasons. (2) stochastic simulation uses one of several available stochastic weather generators in combination with the inverse approach of Guo et al. (2018). This vignette focuses on the simple/seasonal scaling approach. Scaling methods are commonly used because they are fast and easy to apply and can effectively represent changes in annual or seasonal averages and totals—factors often considered key drivers of system response. However, this approach has limitations, particularly in its inability to capture changes in extremes, wet–dry sequencing, etc. For details on how to implement the more flexible stochastic simulation approach refer to the Guide to Stochastic Simulation vignette.

Typical function calls to `generateScenarios()`, using tank reference data, the OAT and regGrid exposure spaces defined above, and with the scaling approach, are shown below.

The outputs of the simulation are saved to a list object (called `sim.OAT` and `sim.regGrid`), which is then typically used as an input to the system model as part of Step C. These objects contains perturbed time series for all the target locations in the exposure space and the resulting simulation file are typically large in size. The summary (metadata) of the simulation may be obtained using the function `getSimSummary()`, which can be used as input to the plotting functions in foreSIGHT.

    # generate perturbed time series using simple scaling
    sim.OAT <- generateScenarios(
      reference = tank_obs, # input observed data
      expSpace = expSpace.OAT, # exposure space created by the user
      controlFile = "scaling"
    ) # using simple/seasonal scaling
    #> Generating replicate number 1 out of 1 replicates...
    #> Generating target number 1 out of 23 targets...
    #> Generating target number 2 out of 23 targets...
    #> Generating target number 3 out of 23 targets...
    #> Generating target number 4 out of 23 targets...
    #> Generating target number 5 out of 23 targets...
    #> Generating target number 6 out of 23 targets...
    #> Generating target number 7 out of 23 targets...
    #> Generating target number 8 out of 23 targets...
    #> Generating target number 9 out of 23 targets...
    #> Generating target number 10 out of 23 targets...
    #> Generating target number 11 out of 23 targets...
    #> Generating target number 12 out of 23 targets...
    #> Generating target number 13 out of 23 targets...
    #> Generating target number 14 out of 23 targets...
    #> Generating target number 15 out of 23 targets...
    #> Generating target number 16 out of 23 targets...
    #> Generating target number 17 out of 23 targets...
    #> Generating target number 18 out of 23 targets...
    #> Generating target number 19 out of 23 targets...
    #> Generating target number 20 out of 23 targets...
    #> Generating target number 21 out of 23 targets...
    #> Generating target number 22 out of 23 targets...
    #> Generating target number 23 out of 23 targets...
    #> Simulation completed
    # get metadata for later use
    simSummary.OAT <- getSimSummary(sim.OAT)

    # generate perturbed time series using simple scaling
    sim.regGrid <- generateScenarios(
      reference = tank_obs, # input observed data
      expSpace = expSpace.regGrid, # exposure space created by the user
      controlFile = "scaling"
    ) # using simple/seasonal scaling
    #> Generating replicate number 1 out of 1 replicates...
    #> Generating target number 1 out of 25 targets...
    #> Generating target number 2 out of 25 targets...
    #> Generating target number 3 out of 25 targets...
    #> Generating target number 4 out of 25 targets...
    #> Generating target number 5 out of 25 targets...
    #> Generating target number 6 out of 25 targets...
    #> Generating target number 7 out of 25 targets...
    #> Generating target number 8 out of 25 targets...
    #> Generating target number 9 out of 25 targets...
    #> Generating target number 10 out of 25 targets...
    #> Generating target number 11 out of 25 targets...
    #> Generating target number 12 out of 25 targets...
    #> Generating target number 13 out of 25 targets...
    #> Generating target number 14 out of 25 targets...
    #> Generating target number 15 out of 25 targets...
    #> Generating target number 16 out of 25 targets...
    #> Generating target number 17 out of 25 targets...
    #> Generating target number 18 out of 25 targets...
    #> Generating target number 19 out of 25 targets...
    #> Generating target number 20 out of 25 targets...
    #> Generating target number 21 out of 25 targets...
    #> Generating target number 22 out of 25 targets...
    #> Generating target number 23 out of 25 targets...
    #> Generating target number 24 out of 25 targets...
    #> Generating target number 25 out of 25 targets...
    #> Simulation completed

    # get metadata for later use
    simSummary.regGrid <- getSimSummary(sim.regGrid)

## 2.3. Step C: Simulate system performance

To understand how the system responds to perturbations in the climate attributes, it is necessary to run the system model for each of the simulated weather time series in order to obtain one or several measures of system performance that relate to each point in the exposure space. The perturbed time series generated in the previous section can be used to run the system model and calculate the system performance at all locations in the exposure space using the function `runSystemModel`.

`runSystemModel()` requires a perturbed simulation (`sim`) and system model function (`systemModel`) as input arguments, as well as a list of system model arguments (`systemArgs`) that need to be supplied to the system model, and a vector of performance metrics (`metrics`) that represent the system model output(s). This section describes the structure of the `systemModel` and demostrates how it can be called to calculate system performance throughout the exposure space.

System model

`simulateSystem` is the core system model function that simulates the system and calculates and returns multiple performance metrics. `systemModel` is a wrapper function that calls `simulateSystem`, which is intended to interface with `runSystemModel`. The `systemModel` function:
- receives `data, systemArgs` in the specific format
- translates/reformats the inputs to the format required by `simulateSystem` (if necessary)
- subsets (if necessary) and returns the relevant `metrics`

`systemModel` takes in arguments `data`, `systemArgs`, and `metrics`. `data` is a list of climate data with format described above. `systemArgs` is a list containing the system arguments that are required by `simulateSystem`. `metrics` is a vector of strings containing the names of the performance metrics that `systemModel` should return. It is recommended that the names of the performance metrics also include the units of the metrics. This will ensure that the units are available in the names of the performance metrics outputs created using `runSystemModel` and will be included in the legend labels of plots created using the downstream performance plotting functions in foreSIGHT.

This case study uses as the system model the rainwater tank function, `tankWrapper()` described at the start of Section 2.

The `systemModel` code for the rainwater tank system model (`tankWrapper`) is:

    tankWrapper
    #> function (data, systemArgs, metrics)
    #> {
    #>     performance <- tankPerformance(data = data, roofArea = systemArgs$roofArea,
    #>         nPeople = systemArgs$nPeople, tankVol = systemArgs$tankVol,
    #>         firstFlush = systemArgs$firstFlush, write.file = systemArgs$write.file,
    #>         fnam = systemArgs$fnam)
    #>     performanceSubset <- performance[metrics]
    #>     return(performanceSubset)
    #> }
    #> <bytecode: 0x0000026ada87b9d0>
    #> <environment: namespace:foreSIGHT>

A typical function call to `tankWrapper` is shown below, which returns the system performance metrics specified by the user.

    # Load example climate data
    data(tankDat)

    # View the metrics available for use
    tankMetrics <- viewTankMetrics()
    #> [1] "volumetric reliability (fraction)" "reliability (fraction)"
    #> [3] "system efficiency (%)"             "storage efficiency (%)"
    #> [5] "average tank storage (L)"          "average daily deficit (L)"

    # User input: system model parameters
    systemArgs <- list(
      roofArea = 205, # roof area in m2
      nPeople = 1, # number of people using water
      tankVol = 2400, # tank volume in L
      firstFlush = 2.0, # depth of water removed each event in mm
      write.file = FALSE
    ) # do not write output to file

    # performance metric chosen for reporting
    metrics <- c("average daily deficit (L)", "reliability (fraction)")

    performanceOut <- tankWrapper(data = tank_obs, systemArgs = systemArgs, metrics = metrics)
    performanceOut
    #> $`average daily deficit (L)`
    #> [1] 27.42634
    #>
    #> $`reliability (fraction)`
    #> [1] 0.8075554

    # Now try a different metric e.g. volumetric reliability
    performanceOut <- tankWrapper(data = tank_obs, systemArgs = systemArgs, metrics = tankMetrics[1])
    performanceOut
    #> $`volumetric reliability (fraction)`
    #> [1] 0.6235379

The below code template can be used to create scripts to simulated system performances in R. More advanced templates for creating wrappers that call external system models from R, or procedures that can be used to write climate inputs and read system metrics from external models that cannot be run from R, are provided in Section 3.2.

Creating wrappers for system models in R

In general, to use custom system models in R, the user should define a wrapper function `systemModel` adhering to the input-output requirements described below. The code below shows the generalised structure of the `systemModel` wrapper function.

    systemModel <- function(data, # data.frame with columns: year, month, day, *var1*, *var2* etc.
                            systemArgs, # list containing the arguments of simulateSystem
                            metrics) { # names of performance metrics (with units of the metrics)

      # convert data to format required for simulateSystem
      # Note that "reformat" is a dummy function shown here for
      # illustration
      dataforSimulateSystem <- reformat(data)

      # call simulateSystem and get system performance metrics
      # simulateSystem is the core system model function
      systemPerformance <- simulateSystem(
        data = dataforSimulateSystem,
        arg1 = systemArgs[[1]],
        arg2 = systemArgs[[2]],
        ...
      )
      # subset & return metrics (can name performance metrics
      #                          here if required)
      performanceSubset <- systemPerformance[metrics]
      return(performanceSubset)
    }

Running system model for perturbed climates

A typical function call to `runSystemMode` that uses the ‘OAT’ exposure space generated above is shown below.

    # define the arguments of the systemModel, here tankWrapper
    systemArgs <- list(
      roofArea = 205,
      nPeople = 1,
      tankVol = 2400,
      firstFlush = 2.0,
      write.file = FALSE
    )

    metrics <- c("average daily deficit (L)", "reliability (fraction)", "volumetric reliability (fraction)")

    # run the system model
    systemPerf.OAT <- runSystemModel(
      sim = sim.OAT, # simulation; the perturbed time series
      systemModel = tankWrapper, # the system model function
      systemArgs = systemArgs, # argument to the system model function
      metrics = metrics
    ) # selected performance metrics

    # the output contains three performance metrics for all target locations
    utils::str(systemPerf.OAT)
    #> List of 3
    #>  $ average daily deficit (L)        : num [1:23, 1] 33.2 32 31 30 29.1 ...
    #>  $ reliability (fraction)           : num [1:23, 1] 0.758 0.767 0.777 0.786 0.794 ...
    #>  $ volumetric reliability (fraction): num [1:23, 1] 0.544 0.56 0.575 0.589 0.601 ...

The output of `runSystemModel` contains the values of the performance metrics at all target locations in the exposure space.

Similarly, we can calculate system performance based on the ‘regGrid’ exposure space.

    # run the system model
    systemPerf.regGrid <- runSystemModel(
      sim = sim.regGrid, # simulation; the perturbed time series
      systemModel = tankWrapper, # the system model function
      systemArgs = systemArgs, # argument to the system model function
      metrics = metrics
    ) # selected performance metrics

    # the output contains three performance metrics for all target locations
    utils::str(systemPerf.regGrid)
    #> List of 3
    #>  $ average daily deficit (L)        : num [1:25, 1] 29.8 28 26.5 25.2 24.1 ...
    #>  $ reliability (fraction)           : num [1:25, 1] 0.769 0.785 0.797 0.811 0.821 ...
    #>  $ volumetric reliability (fraction): num [1:25, 1] 0.583 0.607 0.629 0.647 0.663 ...

## 2.4. Step D: Visualise system performance

The foreSIGHT package contains functions to plot the system performance calculated using `runSystemModel`. The output of `runSystemModel` is termed the “performance space” as it contains the system performance at each target location in the exposure space. The performance plotting functions use the system performance and the simulation summary as input arguments. There are three functions to plot performance available in foreSIGHT. Typical usages of these functions are demonstrated below using the case study data.

plotPerformanceOAT: The function plots lines showing the changes in a system performance metric with one-at-a-time (OAT) perturbations in attributes. This function is intended for use with an ‘OAT’ exposure space, assuming all other attributes are held constant (usually at their historical levels). However, if ‘OAT’ perturbations exist in a ‘regGrid’ exposure space, the function will subset these targets to create the plots. This subset can be thought of as a slice through the exposure space when the other attributes are kept at historical levels. If the exposure space does not contain attribute values at historical levels, the ‘OAT’ plots cannot be created. `plotPerformanceOAT()` will print an error to inform that there are no ‘OAT’ perturbations in the exposure space in such an instance.

For example, OAT plots may be created for the case study example in this vignette as shown below.

    # OAT perturbation plot of the first metric
    p1_2 <- plotPerformanceOAT(systemPerf.OAT, # performance metrics
      simSummary.OAT, # summary of the simulation
      metric = "volumetric reliability (fraction)"
    ) # name of the metric
    #> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
    #> ℹ Please use `linewidth` instead.
    #> ℹ The deprecated feature was likely used in the foreSIGHT package.
    #>   Please report the issue at
    #>   <https://github.com/ClimateAnalytics/foreSIGHT/issues>.
    #> This warning is displayed once every 8 hours.
    #> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
    #> generated.
    #> Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
    #> ℹ Please use the `linewidth` argument instead.
    #> ℹ The deprecated feature was likely used in the foreSIGHT package.
    #>   Please report the issue at
    #>   <https://github.com/ClimateAnalytics/foreSIGHT/issues>.
    #> This warning is displayed once every 8 hours.
    #> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
    #> generated.

Changes in performance metric: volumetric reliability (fraction)

Changes in performance metric: volumetric reliability (fraction)

plotPerformanceSpace: The function plots a heatmap and contours of a given system performance metric at multiple target locations in the exposure space. This function can be used to plot from two dimensions selected by the user. The function returns a `ggplot` object that can be used to place the plots in panels using higher dimensions, as desired by the user. In some cases, there may be a clear performance ‘threshold’, above or below which the system performance becomes undesirable and/or triggers a system ‘failure’ (for example, an agreed minimum specified level of system reliability). In this case, the user may specify the threshold value of the performance metric as an input argument, resulting in the addition of a thick contour line to the plot in order to mark this threshold in the performance space. If the exposure space has more than two perturbed attributes (i.e., a multi dimensional performance space), the user may specify the perturbed attributes to be used as the axes of the plot.

The case study has only two perturbed attributes. The performance space heatmap and contours may be plotted using the code shown below. In this case a threshold of 28 litres per day for the average maximum water deficit is used to illustrate a hypothetical case where the rainwater tank is no longer economically viable for deficit values above this threshold. This threshold assumes that the rainwater tank was sized on the threshold of economic viability using historical climate assumptions.The threshold of 28L/day corresponds to about 10% of the peak (summer) season residential water use of a single person household.

    # performance space of the first metric with a user specified threshold of 27
    p2 <- plotPerformanceSpace(systemPerf.regGrid, # performance metric
      simSummary.regGrid, # simulation summary
      metric = "average daily deficit (L)", # the name of the metric
      perfThresh = 28, # user-defined performance threshold
      perfThreshLabel = "Max Deficit"
    ) # custom label for the threshold

Performance space and user defined threshold

It is possible to add various lines of evidence to provide guidance on which parts of the exposure space are more or less plausible in a future climate. For example it is possible to superimpose projections from climate models to the performance space plotted using `plotPerformanceSpace()`. This climate data should contain values of projected changes in attributes that are used as the axes of the performance space, and which need to be developed separately from the foreSIGHT workflow. For example, one might extract relevant attribute values from a 30 year future timeslice from the relevant climate model output, potentially after downscaling, bias correction or other processing.

One may also elect to use the climate model simulations (potentially after downscaling, bias correction or other processing) as inputs to the system model to generate new performance values corresponding to each projection time series, and in this case it is possible to plot the performance values corresponding to the climate model simulations as coloured data points in plots created using `plotPerformanceSpace`, using the same colour scale. The code below provides an example, in which the system performance obtained from the time series of two simulations from each of three climate models are illustrated.

    # load an example alternate climate data available in the package
    data("egClimData")

    # expected format of egClimData - the column names should match attribute and performance metric names!
    head(egClimData)
    #>   P_day_all_tot_m P_day_all_seasRatioMarAug P_day_all_P99 Temp_day_all_avg
    #> 1       0.7472544                  1.139591     0.9098522         1.273604
    #> 2       0.7436642                  1.151786     0.9086538         1.267123
    #> 3       0.7148945                  1.109015     0.9174129         1.256753
    #> 4       0.8123943                  1.020873     0.9697000         1.222056
    #> 5       0.8173473                  1.047943     0.9218819         1.214068
    #> 6       0.8232569                  1.056295     0.9411360         1.219123
    #>       Name Avg. Deficit
    #> 1  canesm2     30.54369
    #> 2  canesm2     33.90540
    #> 3  canesm2     32.96717
    #> 4 cnrm.cm5     27.71773
    #> 5 cnrm.cm5     27.72935
    #> 6 cnrm.cm5     28.26993

    # performance space overlaid with alternate climate data
    p3 <- plotPerformanceSpace(systemPerf.regGrid[1], # performance metric
      simSummary.regGrid, # simulation summary
      perfThresh = 28, # user-defined performance threshold
      perfThreshLabel = "Max Deficit", # custom label for the threshold
      climData = egClimData
    ) # alternate climate data

Performance spaces overlaid with alternate climate data points

plotPerformanceSpaceMulti: The third function available in foreSIGHT for plotting system performance is the joint presentation of multiple system performance metrics to facilitate decision making. The function plots contours showing the number of performance metric thresholds exceeded in the performance space. The user should specify the minimum or maximum thresholds of each performance metric as input arguments for calculation.

Using the case study data, the number of thresholds exceeded may be plotted as shown below.

    # plot number of performance thresholds exceeded
    p4 <- plotPerformanceSpaceMulti(systemPerf.regGrid, # 3 performance metrics
      simSummary.regGrid, # simulation summary
      perfThreshMin = c(NA, 0.8, 0.6), # min thresholds for each metric
      #       use NA if not applicable
      perfThreshMax = c(28, NA, NA), # max thresholds for each metric
      climData = egClimData
    )

## 2.5 Step E: Evaluate system options

The process of climate ‘stress-testing’ is often undertaken to facilitate decisions involving choices between multiple system configurations or operating policies. Step E of the process involves comparison of the results from these alternate choices under climate perturbations. Steps C and D are repeated for the multiple configurations/policies and the results are used to compare the performance spaces and number of thresholds exceeded to identify choices that exhibit better performances under changes in climate attributes. Figures created using the plotting functions `plotPerformanceSpace()` and `plotPerformanceSpaceMulti()` discussed in Step D can be used for the comparison.

To illustrate such a comparison in this case study, consider the following alternate configuration of a rain water tank. The original configuration is referred to as “System A” and this alternate configuration is referred to as “System B” henceforth. System B is alternate rain water tank design that harvests water from the some roof area (`roofArea`), but has a higher tank volume (`tankVol`), compared to System A. The generated perturbed scenarios are used to simulate the performance of this alternate system (System B) in the example below.

    # plot number of performance thresholds exceeded
    systemArgsB <- list(
      roofArea = 205,
      nPeople = 1,
      tankVol = 2600,
      firstFlush = 2.0,
      write.file = FALSE
    )

    # run the system model
    systemPerfB.regGrid <- runSystemModel(
      sim = sim.regGrid, # simulation; the perturbed time series
      systemModel = tankWrapper, # the system model function
      systemArgs = systemArgsB, # argument to the system model function
      metrics = metrics
    ) # selected performance metrics

The performance spaces of the average daily deficit (L) for the two system configurations can be compared by generating the figures shown in Step D for system B, using the code below.

    data("egClimData")

    p5 <- plotPerformanceSpace(systemPerfB.regGrid[1], # performance metric
      simSummary.regGrid, # simulation summary
      perfThresh = 28, # user-defined performance threshold
      perfThreshLabel = "System B\n Max Deficit", # custom label for the threshold
      climData = egClimData, # alternate climate data
      colLim = c(18, 38)
    )

    p6 <- plotPerformanceSpaceMulti(systemPerfB.regGrid, # 3 performance metrics
      simSummary.regGrid, # simulation summary
      perfThreshMin = c(NA, 0.8, 0.6), # min thresholds for each metric
      #       use NA if not applicable
      perfThreshMax = c(28, NA, NA), # max thresholds for each metric
      climData = egClimData
    ) # alternate climate data

Performance spaces for two alternate rain water tank configurations

In addition foreSIGHT contains a function named `plotOptions()` that can be used to plot the differences in the performance metrics of two system options and the shift the performance threshold contours. The below code provides an example.The comparison shows that the threshold values of maximum deficit (28 litres) are exceeded in a fewer perturbed scenarios for systemB (lower area of the performance space), compared to systemA. Thus, systemB exhibits better performance in terms of this performance metric.

    data("egClimData")

    p7 <- plotOptions(
      performanceOpt1 = systemPerf.regGrid[1], # performance metrics of option 1
      performanceOpt2 = systemPerfB.regGrid[1], # performance metrics of option 2
      sim = simSummary.regGrid, # simulation metadata
      opt1Label = "System A", # label of option 1
      opt2Label = "System B", # label of option 2
      titleText = "Avg Deficit: System B - System A", # plot title
      perfThresh = 28, # threshold value of the metric
      perfThreshLabel = "Max. Deficit (28L)"
    ) # other climate data

Performance spaces for two alternate rain water tank configurations

Thus, in this case study the results of the ‘stress-test’ indicate that system B is preferable as it should operate satisfactorily across a wider range of conditions, including the drier climate projected by the alternate climate data.

This concludes the workflow of climate ‘stress-testing’ a system using the foreSIGHT package. These steps may be repeated to compare the performance of other system configurations or operating policies.

# 3. Advanced usage

## 3.1 Stochastic simulation

Exploring changes to some climate attributes—such as extremes, wet-dry day sequencing, and interannual variability—cannot be captured using simple scaling methods. In such cases, stochastic simulation provides a more flexible approach.

For these cases, foreSIGHT uses the ‘inverse method’ of Guo et al (2018) to optimise the parameters of stochastic weather generator type models to generate time series with desired perturbed attributes.

Using foreSIGHT’s inverse modelling framework, users can generate climate time series that match specified target attributes. This allows for the exploration of more complex and realistic climate change scenarios.

For detailed instructions, refer to the foreSIGHT Guide to Stochastic Simulation vignette.

## 3.2 Coupling with external system models

In some cases, the user may be interested in using system models that are not available in R.

Calling an external system model from foreSIGHT

If the system model can be called from the command line, then a wrapper can be used to

1. Write climate data to system model input files

2. Executue the system model via a command line call (e.g. using the base-R `system()` function)

3. Read system model output and calculate system metrics.

A template wrapper for this approach is shown below:

    systemModel <- function(data, # data.frame with columns: year, month, day, *var1*, *var2* etc.
                            systemArgs, # list containing the arguments of simulateSystem
                            metrics) { # names of performance metrics (with units of the metrics)

      # Note this template has dummy functions"writeInputData", "systemModelCommand", "readOutputData" and "calculatePerformance" which are used to illuistrate use, but are not available

      # convert data to format required for csv climate inputs used in external system model
      writeInputData(data, systemModelInputFiles)

      # use the "system()" function to run external system commands from within R
      # (i.e. communicate directly with the command line or shell of your operating system).
      # The dummy "systemModelCommand" represents the command line call to the system model, with arguments arg1, arg2, etc.
      system(paste(
        systemModelCommand,
        " arg1=", systemArgs[[1]],
        " arg2=", systemArgs[[2]],
        ...
      ))

      # Read output files from the system model.
      systemOutput <- readOutputData(outputFiles)

      # Calculate performance metrics from system model output.
      systemPerformance <- calculatePerformance(systemOutput)

      # subset & return metrics (can name performance metrics
      #                          here if required)
      performanceSubset <- systemPerformance[metrics]
      return(performanceSubset)
    }

Writing and reading foreSIGHT data files for an external system model

If the system model cannot be from the command line, the user can

1. write scenarios generated in Step B to CSV files to be used in other programming environments,

2. execute the system model externally, looping over all replicates and targets

3. read system performance metrics calculated in other languages and saved in CSV files into an R workspace.

The templates can be modified by the user for the system models of their interest.

Example code to write scenarios to CSV files for external system models:

    # Example code to write the generated perturbed time series to csv files which may be used to run
    # system models in other software environments/programming languages
    # ========================================================================================================

    repNames <- names(sim[grep("Rep", names(sim))]) # replicate names
    tarNames <- names(sim[[repNames[1]]]) # target names

    nRep <- length(repNames)
    nTar <- length(tarNames)
    varNames <- c("P", "Temp") # variable names

    # create list for climate data
    clim <- list()

    # loop over replicates and targets
    for (r in 1:nRep) {
      for (t in 1:nTar) {
        clim$times <- sim[["simDates"]] # dates of the simulation, will add variables later
        for (v in varNames) {
          # add variables
          clim[[v]] <- sim[[repNames[r]]][[tarNames[t]]][[v]][["sim"]]
        }
        # write output from given replicate/target to file (dummy function "writeInputData")
        outCSVFile <- paste0("Scenario_Rep", r, "_Tar", t, ".csv") # name the csv file as desired
        writeInputData(clim, outCSVFile)
      }
    }

    # Scenario_*.csv files can be used to run external system models

Example code to read system performance saved in CSV files to an `R workspace`

    # Example code to read the system performances calculated using the generated time series into an
    # R workspace. It is assumed that the system performances calculated in another software environment
    # are saved in separate files for each scenario (named by replicate and target numbers)
    # ========================================================================================================

    metrics <- c("performance metric 1 (%)", "performance metric 2 (fraction)") # metric names
    repNames <- names(sim[grep("Rep", names(sim))]) # replicate names
    tarNames <- names(sim[[repNames[1]]]) # target names
    nRep <- length(repNames)
    nTar <- length(tarNames)

    systemPerformance <- list() # initialised in the format that runSystemModel would return
    for (m in 1:length(metrics)) {
      systemPerformance[[metrics[m]]] <- matrix(NA, nrow = nTar, ncol = nRep)
    }

    # read from files containing metric values from an external system model
    for (r in 1:nRep) {
      for (t in 1:nTar) {
        # name of the csv file containing the system performance, header is the name of the metric
        inCSVFile <- paste0("SystemPerformances_Rep", r, "_Tar", t, ".csv")
        # check.names = FALSE is useful if the metric names contain brackets around the units
        systemPerfIn <- readMetrics(clim, outCSVFile)
        for (m in 1:length(metrics)) {
          systemPerformance[[metrics[m]]][t, r] <- systemPerfIn[[metrics[m]]]
        }
      }
    }

    # systemPerformance can be used as an input to plotPerformanceSpace, plotPerformanceOAT, and plotPerformanceMulti

# 4. Conclusions

This vignette demonstrates how foreSIGHT may be used for climate impact assessment of an environmental system. The example code provides an introduction to how the functions work together for application of a ‘stress-testing’ framework.

While the climate impact analysis has been presented as a sequential work flow here for demonstration, there are important aspects to be considered. ‘Stress testing’ typically involves iterative procedures to understand the sensitivities of the system to changes in various hydro-climate attributes, to choose the most relevant attributes for perturbation, and to compare the performance of alternate system configurations. foreSIGHT offers multiple options for an advanced user at each stage of this ‘stress-testing’ process. There are nuances surrounding the selection of these choices especially when more complex stochastic models are used to generate perturbed time series to assess real-world systems. The reader is directed to the Guide for Stochastic Simulation vignette for details on using these methods.

# 5. References

- Brown, C. (2011) Decision-scaling for robust planning and policy under climate uncertainty, World Resour. Rep., World Resour. Inst., Washington D.C.
- Culley, S., S. Noble, A. Yates, M. Timbs, S. Westra, H. R. Maier, M. Giuliani, and A. Castelletti (2016), A bottom-up approach to identifying the maximum operational adaptive capacity of water resource systems to a changing climate, Water Resour. Res., 52, 6751-6768, .
- Guo, D., S. Westra, and H. R. Maier (2018), An inverse approach to perturb historical rainfall data for scenario-neutral climate impact studies, J. Hydrol., 556, 877-890, .
- Prudhomme, C., R. L. Wilby, S. Crooks, A. L. Kay, and N. S. Reynard (2010), Scenario-neutral approach to climate change impact studies: Application to flood risk, J. Hydrol., 390, 198-209, .