Forecast-based operation of re-purposed small reservoirs for floods, farms, and (low) flows

1. Code documentation This document serves as a manual for the collection of data and programs created for this project (Figure 1). Version 2023a or later of MATLAB is required to open the data and run the functions. GIS software, such as ArcGIS or QGIS, may be needed for processing some data to use for the developed programs, but is not strictly required. Figure 1. Graphical overview of the inputs and model functions covered in this document. 1.1. Data sources Required data for the execution of the programs, as well as their sources, are listed in Table 1. These have been preprocessed and provided in the submission package. Table 1. Data included in the submission package. Data Located in file Use Source General reservoir data Reservoir operating rules (Vollstau, Dauerstau, Qcrit) allStudyReservoirs.mat Calculate reservoir operating capacity LUBW Betriebsregeln Inflow time series for selected reservoirs allStudyReservoirs.mat Determine reservoir operation mode, lowflow time series LARSIM Results from perfect knowledge models for selected reservoirs allStudyReservoirs.mat - - Forecasts for selected reservoirs [ReservoirName]_hrb_vhs.mat Calculating effects of forecasts on operation LARSIM, forecast tests provided by HYDRON Results from forecast models for selected reservoirs allVHSReservoirs.mat - - Agricultural Irrigation Demand data Geographic location of the reservoirs reservoirLoc.mat Find weather, soil, and crop data at the location of the reservoir LUBW Reservoir SHP Soil map boart1000_mask.shp Calculating agricultural demand Bundesanstalt für Geowissenschaften und Rohstoffe Crop maps crop20XX.tif Calculating agricultural demand Blickensdörfer et al. (2021) Soil parameters plant_soil_parameters.mat, soilDatabase.m Calculating agricultural demand Bodenkundliche Kartieranleitung FAO-56 crop parameters plant_soil_parameters.mat, plantDatabase.m Calculating agricultural demand Various sources (see [TABLE REF IN AG DEMAND PAPER]) Irrigation demand time series for selected reservoirs processedRawResults.mat Input for irrigation models - Climate data Daily precipitation (HYRAS-DE-PRE, v5.0) dailyP_ETRS_19912022.mat Calculating agricultural demand DWD (2022) Daily FAO-56 reference grass evapotranspiration dailyFAO_ETref_19912021.mat Calculating agricultural demand DWD (2018) Daily maximum temperature dailyTmax_19952020.mat Calculating agricultural demand DWD (2024b) Daily minimum temperature dailyTmin_19952020.mat Calculating agricultural demand DWD (2024d) Daily mean relative humidity dailyHumidity_19952020.mat Calculating agricultural demand DWD (2024c) Daily averaged wind speed @ 10 meters height dailyAvgWind_10m_19952020.mat Calculating agricultural demand DWD (2024a); upscaled from hourly 1.2. Programs An overview of the programs developed for this project and in accordance to the theoretical basis are in Table 2. The vast majority of code for this project is contained in hrb.m; however, some functions in hrb.m which have potential usefulness outside of the context of reservoirs are also included here as separate functions. Table 2. Program files included in the submission package. Name Use Input(s) Output(s) hrb.m Class definition file + reservoir operation functions See section for more details See section for more details applyAllReservoirs.m Loops through all reservoirs in a cell array See section for more details See section for more details calcExceedFlow.m Calculates the percentile exceedance flow of an inflow time series Inflow time series, desired percentile, length of moving window Percentile exceedance flow given as values for one year and for the length of the inflow time series Agricultural Demand plantDatabase.m Sets up a database of plant parameters FAO-56 dual crop plant parameters from various sources Plant parameter database for FAO56 calculations soilDatabase.m Sets up a database of soil parameters Field capacity, wilting point, and readily evaporable water from a combination of KA5 and USDA soil types Soil parameter database for FAO56 calculations compileWeatherData.mlx Sets up a weather table for use in the calcIrrigDemand funcitons Live script reads the various weather data provided in Table 1 Timetable of weather data for irrigation demand calculations calcETcact.m Calculates the FAO56 dual crop ETc,act and its corresponding irrigation demand for a given plant / soil combination and weather data. Called in calcIrrigDemand functions Plant type, soil type, weather data Table of results from the FAO-56 dual crop method generateARUs.m Finds all unique combinations of plant / soil combinations and their estimated areas within a square area around a given location. Called in calcIrrigDemand functions Desired location, desired buffer area size, crop map, soil map Table of ARUs and their corresponding areas calcIrrigDemandStatic.m Calculates the irrigation demand series for all ARUs within a buffer region around a reservoir (single crop map) Crop map, soil map, desired location, buffer area, weather data Table of irrigation demand outputs for each ARU in m3/day calcIrrigDemandDynamic.m Calculates the irrigation demand series for all ARUs within a buffer region around a reservoir (multiple crop maps) Crop map, soil map, desired location, buffer area, weather data Table of irrigation demand outputs for each ARU in m3/day processIrrigResults.m Reformats results of calcIrrigDemand for a more legible format Results of calcIrrigDemand Refined time series for each reservoir in a cell array Forecasting extractLILAVHS.m Reads in a LILA file for forecasts and processes it into a MATLAB-readable format; called in reformatLilaVHS File path of file to be processed Table containing each of the blocks of data reformatLilaVHS.m Recursively processes all LARSIM forecasts in a given catchment Folder containing all files to be processed, list of station names Table containing each of the blocks of data compileVHS.m Creates separate tables for each unique station in the table produced by reformatLilaVHS Folder containing results of reformatLilaVHS, list of station names Separate tables for each unique station Plotting resultPlotting.m Plots all of the results for the perfect-knowledge streamflow case Array of hrb objects that have been streamflowoptimized Various figures agDemandResults.m Plots all of the results for the perfect-knowledge irrigation case Array of hrb objects that have been irrigation-optimized, irrigation demand time series Various figures vhsResults.mlx Plots all of the results for the forecasted streamflow and irrigation results Array of hrb objects that have been optimized under forecasting, irrigation demand time series Various figures 1.2.1. hrb.m The file hrb.m is the class definition file and is mandatory for loading, reading, and working with the reservoir data, as it is constructed using object-oriented programming. Each reservoir is stored as a MATLAB object, which can contain a variety of different properties. The majority of developed functions (called “methods”) are stored in the hrb.m file. A further explanation of object-oriented programming in MATLAB is provided in the MATLAB documentation and will not be elaborated here. 1.2.1.1. Object Properties The object properties that can be stored are shown in Table 3, with mandatory properties listed in bold. The remaining properties are either optional for identification purposes or calculated through other functions. They are organized into several groups: - Reservoir properties store information about the reservoir. - Time series properties store time series that are unique to the reservoir, including results from the descriptive model operation. - Current operation parameters store information about the reservoir operation. - Volume parameters store information about threshold volume levels in the reservoir. - Optimization values store results from the flood-only model and streamflow-optimization model. - Agriculture demand values store results from the irrigation model. Table 3. Object properties defined in the class definition file hrb.m. Property Description Source Reservoir Properties hrbName Reservoir name User input catchment LARSIM model catchment User input hauptschluss Dam type (in-river or diversion) User input use_short Usage code (HW = flood only, HWX = flood + multipurpose) User input use_long Usage(s) as a list User input months_w Months of winter operation, if available User input days_w Day of month when winter operation begins User input measured Measured data, if exists User input fuellkurve Volume-height-area curves User input class_id_long Combination of size, dam type, usage, and inundation type User input rf Rueckhaltefaktor; also called availability factor (AF) Calculated from calcRF rf_d Originally-estimated RF (if applicable) User input; calculated based on estimates from Regionalisierung rf_d_p Rueckhaltefaktor percentile User input; calculated based on estimates from Regionalisierung; percentile based on class area_ratio Ratio of LARSIM catchment area to delineated catchment area; used to scale LARSIM discharge User input; LARSIM, catchment delineation Time Series Properties qin Area-corrected inflow timeseries to the reservoir User input; calculated by multiplying area_ratio with qin_d qin_d LARSIM inflow time series to the reservoir User input; LARSIM qd Standard minimum flow for one year; stores Q70, Q80, and Q90 Calculated from calcMinFlow based on qin qd_ts Minimum flow time series with the dimensions of qin; stores Q70, Q80, and Q90 Produced by qdTS from qin and qd v Volume time series from the descriptive model Calculated from descriptiveModel qout Discharge time series from the reservoir operating in the descriptive model Calculated from descriptiveModel ssi Standardized Streamflow Index based on qin; stores 1-month, 3-month, and 6-month accumulation periods Calculated from calcSSI Current Operation Properties Qr_d Regelabfluss (flood protection level); default value (if seasonal, assumed to be the summer value). User input Qr_w Regelabfluss (flood protection value) in the winter User input HWE Discharge-volume curve for extreme floods for the Hochwasserentlastungsanlage User input HWE_2 Discharge-volume curve for extreme floods for the secondary Hochwasserentlastungsanlage User input GA Discharge-volume curve for extreme floods through the Grundablass User input ED Fastest possible emptying duration for the operating capacity (Vv-Vd), assuming the reservoir can be emptied at Qr_d Calculated from calcED Volume Properties Vd Dauerstau volume (if seasonally operated, the summer value) User input Vd_w Dauerstau volume (winter) User input Vv Standard reservoir capacity User input Vh1 Hochwasserstauziel 1 User input Vh2 Hochwasserstauziel 2 User input Vk Crown volume User input Streamflow Optimization Values (Perfect Knowledge) qout_od Outflow time series calculated from the flood model Calculated from floodOptModel v_od Volume time series calculated from the flood model Calculated from floodOptModel pD_d Streamflow drought penalty time series based on floodOptModel; considered “default” Calculated from floodOptModel pF_d Streamflow flood penalty time series based on floodOptModel; considered “default” Calculated from floodOptModel pMonthly Total monthly default penalties across the time series Calculated from penaltyRidgePlot results Stores results from streamflow optimization under perfect knowledge Calculated from multiOptModel baseStats “Base statistics” (relating to time / volume / penalty) for flood and drought, calculated based on each model (semi-natural, floodoptimized, streamflow-optimized) Calculated from baseStats using results from multiOptModel Agriculture Demand Values (Perfect Knowledge) ag_results Stores results from irrigModel based on streamflow-optimized storage rules Calculated from irrigModel using results from multiOptModel and calcIrrigDemand ag_ben Stores the agricultural water supply benefit from irrigModel, based on ag_results Calculated from irrigModel using results from multiOptModel and calcIrrigDemand Forecasting Values vhs_summary Contains the output table of VHS optimization showing the results of all runs of the forecasting model and their performance Calculated from multiOptForecastModel vhs_irrig Contains the subset of vhs_summary that is viable for irrigation benefit Calculated from multiOptForecastModel vhs_stream Contains the subset of vhs_summary that is viable for streamflow supplementation Calculated from multiOptForecastModel 1.2.1.2. Object Methods This file also stores almost all of the functions (methods) developed for this project. The list of all methods are provided in Table 4. Not all functions will be called when running the core models, as some may either be called within other functions or were developed to provide additional visualizations or data that was ultimately not needed. The main functions that will be called when setting up and running the models are in bold. Table 4. Functions included in the class definition file hrb.m. The main functions needed are in bold, and functions with an optional plotting variable (plotFig) are indicated with *. Method Name Description Input Output hrb Initializes an hrb object hrbName, qin_d, area_ratio Updated hrb object with qin calcSSI Calculates standardized streamflow index with 1-, 3-, and 6-month accumulation periods hrb object Updated hrb object with streamflow drought indices descriptiveModel Uses Vd, Vv, Qr, HWE, HWE2, GA to model existing reservoir operation hrb object Updated hrb object with timeseries of volume and Qout plotVQout Plots the results of descriptiveModel hrb object Figure calcMinFlow* Calculates Q70, Q80, and Q90 based on qin hrb object Updated hrb object, timestep (accepts ‘daily’ or ‘hourly’) qdTS Converts the standard minimum flow for one year into a time series the hrb object, qd Updated hrb object with qd fitted to the length of qin length of qin. Called within calcMinFlow floodOptModel* Models the floodoptimized (i.e. simplified from the existing) reservoir operation. Differs from descriptiveModel in that it does not include HWE, HWE2, or GA in the calculations, and calculates “default” penalties hrb object, (one column of) qd_ts Updated hrb object with Qout and volume timeseries penalty Calculates flood and drought penalties. This function is called in both floodOptModel and compPenalty Target Qout time series (as an array), hrb object Drought and flood penalty (as an array) optModel Runs an instance of the combined operation model with a given Qr and low flow time series hrb object, (one column of), Qr value to test, qd_ts Various results relating to the operation (Qout, volume, released volume, stored volume) compPenalty Compares the “default” penalties to the penalties produced by optModel Qout (as an array), hrb object, (one column of) qd_ts Drought and flood penalties, maintained flood protection (true / false), difference in drought penalty multiOptModel* Optimization function for streamflow supplementation hrb object, number of desired runs, (one column of) qd_ts Updated hrb object with model results; contains matrix of all runs and their results, the optModel results for the most beneficial run, and the optimal Qr tested. If plotFig is active, it will plot the highest-benefit run irrigModel* Runs an instance of the combined operation model with a given Qr to hrb object, irrigation demand timeseries, Qr Updated hrb object with timeseries of model results To combine the files into a single forecast input table, we can use compileVHS.m. This function takes as inputs the path to the folder where these files are contained as well as an n by 2 matrix where the first column is the shortened station name and the second column is how it appears in ergebnis.lila: Then, after calling the function a file with the name [short name]_hrb_vhs.mat will be saved in the given folder. This file will contain a table where the forecast at each timestep is saved as a timetable: 2. Tutorials This section demonstrates tutorials on how to run and operate the programs discussed in the documentation. To run these programs, MATLAB version 2023a or newer is required, and basic familiarity with MATLAB is assumed. In this tutorial, we walk through the model set up, runs, and some visualization options with a singular reservoir (L12 Schwaigern). We also explore different options and processing tutorials for use in other regions. For some of these options, familiarity with GIS software is required. We use QGIS 3.40.2 for the examples here. 2.1. Setting up files All files and folders must be on the MATLAB path in order for them to operate. Right-click the folder and select the option Add to Path > Selected Folders and Subfolders. Load the example discharge file named schwaigern_discharge.mat. Opening the variable discharge should show the following: The discharge is the inflow time series into the reservoir and must contain a time series as the first column and the discharge as the second column. 2.2. Setting up a new hrb object The hrb object is the cornerstone of the model. We create an object using the constructor function, giving the name, inflow time series, and (optionally) an area ratio to scale the inflow. The area ratio is necessary if the inflow time series is modeled with a slightly different area than the delineated area of the reservoir and must be rescaled. The constructor function takes the name of the reservoir (taken as a string) and the discharge (as a table) as follows: A new hrb object is created. The constructor function automatically renames the columns of the inflow time series as “time” and “Q” and saves the inflow time series as qin. The original discharge time series qin_d is stored separately in case the area ratio needs to be revised. Opening the hrb object shows that most of the variables are empty: Adding additional information about the reservoir is done via dot indexing. For example, let’s add the name of the major catchment, the Neckar: The property is updated accordingly. Now, let’s add the operational properties of the reservoir. This model requires three basic operational parameters: the permanent volume (Dauerstau, here the Vd), the operational capacity (Vollstau, here Vv), and the flooding limit (Regelabfluss, or Qr_d). Do note that the capacities must be given in 1000 m3. Let’s add the following values for Schwaigern: Property Schwaigern Value Vd 0 Vv 210,000 Qr_d 3.32 We now have the basic information necessary to run the basic model. 2.3. Running the flood operation model The flood operation model first requires the calculation of the minimum discharge time series, in this case the percentile exceedance flow. This is done simply by calling the function calcMinFlow: This adds the properties qd and qd_ts. The qd is the percentile exceedance flow time series, whereas qd_ts is the qd expanded to match the length of the qin time series. This function automatically calculates the 70th, 80th, and 95th percentile exceedance flows. This can also produce a figure showing the different percentile exceedance thresholds in comparison to the discharges across all years: With these properties, the floodOptModel can be run with a selected qd_ts column: This function calculates the flood-only (default) operation and the accompanying streamflow penalty, and produces an illustrative figure: 2.3.1. Optional: Finding the percentile exceedance time series without an hrb object There may be a time in which a calculation of the percentile exceedance flow without an hrb object is desired, or when a different percentile exceedance flow other than the 70th, 80th, or 95th percentiles are desired. Should that be the case, the function calcExceedFlow can be called where the qin timeseries is in daily or hourly timesteps, the prctile is the desired percentile as a decimal, and mDays is an optional variable defining the window over which the data for the exceedance flow is collected. This is set as 30 days by default. The outputs exceedFlow and exceedFlow_ts can be reassigned to an hrb object as qd and qd_ts, if desired. 2.3.2. Optional: Running the descriptive model The flood operation model is a simplified operation model that only takes into account the operating capacity (the difference between Vv and Vd). The descriptive model is a more sophisticated status-quo model that takes into account any extra flooding capacity and operating rules in the form of a volume-discharge curve (HWE). For Schwaigern, the included HWE obtained when opening schwaigern_HWE can be added: The descriptive model can then be run: The function also returns a figure showing the inflow and outflow time series. However, the descriptive model is not intended for use with the rest of the study. 2.3.3. Optional: Getting penalty ridge plots A function for visualizing the distributions of penalty and deficit volume can be called: The resulting figure can be quite useful, as the resulting statistics are divided by month. This can allow easy visualization of when penalties are the most extreme (most negative) across the observed time period, as well as when flood penalties (if they exist) occur. Each month shows the frequency distribution of penalty or deficit volume, while the color indicates the maximum. 2.4. Running the streamflow operation model The streamflow operation model requires that the flood operation model has been run, as the streamflow model optimizes itself based on the flood operation model. The model takes as inputs the object, the number of runs desired (nQr), and the desired qd_ts: The model then runs the model nQr times, each time using a different retention flow (Qr) between the maximum of the qd_ts and the flooding threshold Qr_d. The results of the optimization are stored in the object as results, with the optimal Qr labeled as ‘Qr_o’ (in this case, 0.1963 m3/s): The results of each run are stored in the cell (1,2), where the first column shows the Qr tested, the second indicates if flood protection was maintained (as true / false), and the third indicates the total penalty benefit. The results of the actual model operation for Qr_o are stored in the cell (2,2). Alongside the outflow and volume time series (Qout and V), other operation data is included: - the flood volume retained at each time step (Vflood) - the volume released for flood protection at each time step (Vrelease_f) - the cumulative volume released for flood protection (Vrelease_f_c) - the volume stored for drought at each time step (Vstore) - the volume released for drought protection at each time step (Vrelease_d) - the cumulative volume released for drought protection (Vrelease_d_c) - the remaining volume needed to reach the drought protection level (Vneed) - the code for the operation mode at each time step (modTS). The codes represent the following modules: 1. Flood operation (flood storage) 2. Pre-flood release 3. Drought release 4. Drought storage 5. Normal operation (maintain storage) 6. Failed drought release The penalties of the Qr_o run, as well as its comparison to the default (flood-only) penalty, are stored in the cell {3,2}. The first year of comparisons are set as NaN, as we allow this first year as a warm-up: The result is a table where the rows are the crop type, the columns are the soil type, and the value indicates the number of pixels covered by this ARU: 2.5.5. Optional: Finding the FAO-56 outputs for a single ARU If the specific FAO-56 outputs (e.g. soil moisture, crop coefficients) for a single ARU are desired, this can also be done separately. After loading in plant_soil_parameters.mat and weatherData_schwaigern.mat, let’s run this for the combination of winter_wheat and Tonschluffe (tu). The soil and plant names should be given as they are listed in plant_soil_parameters.mat: The last variable—irrigPerc, which here is set to 1—is an adjustable value indicating the efficiency of the irrigation. This is the percentage of root zone depletion that is refreshed by an irrigation event. If no irrigation (i.e. a rainfed model) is desired, set irrigPerc to 0. For a less efficient irrigation setup, set irrigPerc between 0 and 1. The results will show the following: Note that this calculates all the variables for the entire length of the weather time series and for a single crop / soil combination. 2.6. Optional: Adapting irrigation demand inputs 2.6.1. Compiling the weather input table using compileWeatherData.mlx The file compileWeatherData.mlx is a live script that processes the provided weather data for a given location into a table by extracting the data from the closest raster cell in each dataset. In order to use it, first give the desired startand end-dates of the time series, as well as the location in XY coordinates: The fileName and fileLocation should be adjusted as needed. The entire script can now be run: This will produce a weather data input table identical to that in weatherData_Schwaigern.mat. 2.6.2. Using different crop maps 2.6.2.1. Processing a crop map shapefile Here, we demonstrate the process of converting a shapefile-based crop map into a raster format readable by this toolbox using the California Central Valley Classified Croplands dataset (Smith et al., 2024) to extract a cropland raster for 2005: The planting codes represent the following crops: R: Rice P: Pasture G: Grain F: Field C: Citrus / mediterranean D: Deciduous fruit and nut V: Vineyard X: Fallow and young perennial T: Truck, nursery, and berry We can add a new field to the attribute that converts the codes into numbers. We start by opening the attribute table and toggling the editing mode: Then we can open the field calculator: We can fill in a series of expressions evaluating the value in the field “2005” and assigning a number to each: Upon clicking “OK”, there should be a new column: Save the changes to the attribute table and un-toggle the editing mode. 2.6.2.2. Revising generateARUs.m When changing crop maps, these should be raster-based, where the value of each cell is a different crop type. There is typically a legend that allows this decoding in the metadata or accompanying documents. In the data from Schwieder et al. (2024), which is what was used for this study, this looks like the following: The name / value pairs from the new dataset must be updated in the subfunction convertCodes in the function generateARUs.m: The value of the cell should be input in [VALUE], and the plant name should replace [PLANT NAME]. If the demand of a plant type should not be calculated, the plant name should be given as ‘other’. It is absolutely critical that the name of the plant provided here matches the name of the crop in plants, or else the crop parameters cannot be called! 2.6.3. Revising crop parameters Crop parameters are provided based on tabulated values from Allen et al. (1998); Pereira et al. (2021a); Pereira et al. (2021b); Rallo et al. (2021) and appear in plant_ soil_parameters.mat as the structure plants. Data for each plant is stored as a structure: Each structure contains the FAO-56 dual crop parameters: 2.6.3.1. Manual revision of crop parameters Individual parameters can be changed via dot notation. For example, if we wanted to change the Lini of dry_broadbean, we could do so with the following command: The parameter would be updated accordingly upon execution: 2.6.3.2. Batch revision of crop parameters The crop parameter structure is constructed using the script plantDatabase.m. The same data for dry_broadbean (as well as all other given crops) is stored here, alongside the source of each parameter: By replacing the values as desired and running the script, the parameters will be updated. 2.6.3.3. Adding new crops New crops can be added by following the format provided in plantDatabase.m: The lines must first be un-commented and the brackets replaced with the proper information in order to run the script. All values except rd_min, the minimum rooting depth, are required. If Lini is not given, assume zero. Planting months must be given either as a single month, written in its full English name (e.g. ‘April’ or ‘September’) or two consecutive months as a cell array (e.g. {‘April’,’May’}). Once complete, the script may be run as usual. 2.6.4. Using different soil maps The provided soil parameters are based on the original study location in Germany. As such, the soil parameters provided are based on the German soil map. For use in other locations, different soil parameters will be needed. This section of the tutorial guides the processing of the global USDANRCS soil texture class map (Knoben, 2021) from raster to polygon, and the revision of crop parameters to match. 2.6.4.1. Processing a soil map raster The file is contained in the folder Global_USDA-NRCS_texture. The raster should look like this upon loading in: Each raster cell has a value from 0 to 12, where each value corresponds to a different texture class: 0: no class assigned (source data sand, silt, clay percentages all contain "no data" values) 1: Clay 2: Clay loam 3: Loam 4: Loamy sand 5: Sand 6: Sandy clay 7: Sandy clay loam 8: Sandy loam 9: Silt 10: Silty clay 11: Silty clay loam 12: Silt loam To ease the processing demands, we will first crop the map to a smaller mask area (the original study area, Baden-Württemberg): The length parameter indicates the maximum length of the value name. We set it to 15, as the names of the soil types can get long. Click on the “VALUE” header name once to select the column, then again to sort the values in ascending order. We can now enter the names of the soil textures: The soil names should not have spaces or special characters and do not have to match these exactly. They should, however, be unique and exactly match the names of the soils in the revised soil parameters. Click “Save” to save your changes and click to untoggle the editing mode. Congratulations! You have now processed another soil map for use with this toolbox. 2.6.4.2. Processing a soil map shapefile When using a different, pre-prepared soil map shapefile, two things must be checked: the name of the soil name field, and the soil names themselves. If the soil name field is incorrect, it can be changed to a readable name (“SOILNAME” or “BODART_GR”) by accessing Layer Properties > Fields > Toggle editing mode: The field name can then be changed: Remember to save changes and exit editing mode. If the soil names themselves have special characters (hyphens, underscores, and so on), these should be removed from the attributes. This can be done manually in the Attribute Table with editing mode active: Remember to save changes and exit editing mode. 2.6.4.3. Editing / adding new soils The soil parameters are structured similarly to the plant parameters. The process of revising soil parameters manually is therefore the same as with the crop parameters. For example, let’s look at the soil ss: A revision of fieldCapacity can be done using dot notation and the value is updated: Batch updates are possible by modifying and running soilDatabase.m. This script also contains a template for adding new soil types: Run the soilDatabase.m script after applying the changes to see them reflected in soils. 2.6.5. Saving changes to plants and soils In order to save the changes into the launchable file, both plants and soils should be re-saved as a new plant_soil_parameters.mat file. This can be done by selecting both variables in the Workspace section: Then by saving the file as plant_soil_parameters.mat. Congratulations! Your updated plant and soil parameters have been changed. These can now be used in the irrigation demand calculations. 2.7. Running the irrigation operation model With the irrigation time series (here, we use the irrig_demand from the static cropping example) and the results of the optimized streamflow operation model (schwaigern from the streamflow operation model), we can now run the irrigation operation model irrigModel using the optimal Qr listed in results{4,2}: This function takes the optimal Qr found in multiOptModel. The results of the irrigation operation are stored in the hrb property ag_results: These take on the same form as the results in the streamflow operation model (section 2.4). Additionally, a plot of the operation time series is produced: Note that irrigModel does not re-check the flood protection status. Thus, it is recommended to use the function multiOptIrrigModel: This function re-tests all values of Qr that were tested in the streamflow multiOptModel. If tf is 'true', then the optimal Qr is the same in both the streamflow and the irrigation models, and there is no increase in downstream flooding. Else, an error message should appear in error. Additionally, this adds the agricultural benefit (in % demand fulfillment) as ag_ben in the hrb object. However, this does not have a plotting option. It is therefore recommended to first run multiOptModel to determine the ag_ben and the flood maintenance conditions, then to use irrigModel to plot as needed. Congratulations! You have now simulated the operation for agricultural use. 2.7.1. Optional: Using a particular Qr in irrigModel If the streamflow-optimized Qr is not desired or should not be used, a Qr value can be manually given as a third input to the irrigModel function: The Qr is subject to the same rules as in the streamflow operation model (section 2.4.2). 2.7.2. Optional: Recreating figures using agDemandResults.m The figures from Ho et al. (2025) can be recreated by running the script agDemandResults.m. 2.8. Running the forecast operation model The forecast operation model is an expansion of both the flood and streamflow operation models that makes a decision at the current timestep based on available information about the future. There are two versions of the model: one for streamflow, and one for irrigation. Each usage should be optimized independently of the other. Because the formats of forecasts can be quite variable, it is difficult to create pre-processing scripts for every situation. We will instead introduce the format in which the forecasts should be given and allow the reader to prepare the data instead. However, if the forecasts were produced in LARSIM, the information in section 1.3 may be useful. 2.8.1. Setting up the forecast input table The forecast input table has two critical components: the forecast itself (stored as a timetable) and the time at which the forecast starts (stored as a datetime variable, VHSZeitpunkt). This should be sorted into a table as follows: This example table is provided in Mittelurbach.mat as the variable vhs_reservoir. The VHSZeitpunkt is used to evaluate whether a forecast is available for decisionmaking, whereas the Data stores the forecast itself. All available forecasts should be included in this table. Moreover, both the VHSZeitpunkt and Data variables should be named exactly as shown. Note that the first three columns (Station, Datenart, and Datentyp) are not mandatory. Omission of these columns will not impact the model operation. The forecast Data takes the form of a timetable, where the first column is the simulated data up to the VHSZeitpunkt and the second column is the forecast. Again, these table variables should be named exactly as shown. 2.8.2. Optional: Running a single forecast operation model using forecastModel The function forecastModel is the base forecast operation model that has options for streamflow or for irrigation usage. The modelType can take either ‘streamflow’ or ‘irrigation’. The modelObjective is the qd_ts in the case of ‘streamflow’ and the irrig in the case of ‘irrigation’. Please refer to the relevant sections for more details. The vhs_reservoir is the forecast input table discussed above. The Qr, percQ_thresh and percQ_release are inputs for calibration. We can test this operation using the provided object Mittelurbach and forecast input table vhs_reservoir (which are contained in the file Mittelurbach.mat) and sample values for percQ_thresh and percQ_release: The output modelOut is a timetable that shows the resulting time series of outflow, volume, release / storage amounts, and module codes throughout the model’s operation: However, this function does not specifically need to be called, as it is called in various other functions. 2.8.3. Optional: Running a single forecast streamflow model The first of these functions is optForecastStreamModel: The input variable plotFig is optional; however, if plotFig is given as 1 (as is above the case), the model will produce a figure: In addition to modelOut, there are three other model outputs. dPen_f returns the total change in flood penalty (compared to the flood operation model), whereas benD gives the streamflow penalty benefit. The output penalties is the time series of flood and drought penalties. This function is called in the process of optimizing the forecasting model. 2.8.4. Optional: Running a single forecast irrigation model The function optForecastIrrigModel is analogous to optForecastStreamModel: Setting the optional input variable plotFig to 1 again plots a figure. However, instead of returning the streamflow penalty benefit and penalties, this model returns the irrigation benefit. 2.8.5. Optimizing the forecast operation model Both versions of the forecast operation model can be optimized using the function multiOptForecastModel. This function takes an (perfect-knowledge streamflow and irrigation optimized) hrb object, the number of Qr, percQ_thresh, and percQ_release values to test, the irrigation demand time series, and the forecast input table. Here, we will test 3 values for each parameter: This function runs both optForecastIrrigModel and optForecastStreamModel for combinations of Qr, percQ_thresh, and percQ_release. It returns an hrb object with three new properties: • vhs_summary, which contains the model results (i.e. streamflow and irrigation benefit and total change in flood penalty) for each combination; • vhs_stream, which contains the subset of vhs_summary for which optimization for streamflow does not increase flood penalty; and • vhs_irrig, which contains the subset of vhs_summary for which optimization for irrgation does not increase flood penalty. Congratulations! You have successfully run the forecast optimization model. 2.8.6. Optional: Using operateForecasts.mlx The file operateForecasts.mlx runs through the entire forecast operation optimization sequence that was used for the original study. This can be useful for double-checking work. 2.8.1. Optional: Recreating figures using vhsResults.mlx The figures from the preprint on forecasting can be recreated by running the script vhsResults.mlx.