27 Calculating ARID
Let us build the next module, ARID, on top of what we just implemented in the load-weather-data module.
27.1 Module variables and parameters
In addition to the weather variables in our dataset, ARID requires several other variables related to soil properties and cover. We will make an exception to our code legibility rule and use the notation from the original implementation (Wallach et al. 2019).
We also introduce a way to handle multiple approaches to parameter setting. We create a chooser named experiment-type in the interface and use it to activate one of two approaches: “user-defined”, which takes the values of sliders in the interface, and “random”, which samples these values randomly within a predetermined range. We also introduce two useful procedures, parameters-check and parameters-to-default, to help us define and return to a specific “default” parameter configuration used as a point of reference.
globals
[
...
;;; default constants
MUF ; Water Uptake coefficient (mm^3.mm^-3)
WP ; Water content at wilting Point (cm^3.cm^-3)
;;;; ETr
albedo_min
albedo_max
;;;; Soil Water Balance model global parameters
WHC_min
WHC_max
DC_min
DC_max
z_min
z_max
CN_min
CN_max
...
]
...
patches-own
[
...
;;;; soil
DC ; Drainage coefficient (mm^3 mm^-3).
z ; root zone depth (mm).
CN ; Runoff curve number.
FC ; Water content at field capacity (cm^3.cm^-3)
WHC ; Water Holding Capacity of the soil (cm^3.cm^-3). Typical range from 0.05 to 0.25
ARID ; ARID index after Woli et al. 2012, ranging form 0 (no water shortage) to 1 (extreme water shortage)
WAT ; Water content in the soil profile for the rooting depth (mm)
WATp ; Volumetric Soil Water content (fraction : mm.mm-1). calculated as WAT/z
;;;; cover
albedo ; canopy reflection or albedo
netSolarRadiation ; net solar radiation discount canopy reflection or albedo
ETr ; reference evapotranspiration
]27.2 Setting module parameters
We use a randomised configuration of most of these patch variables, using hyperparameters to define the envelope of variation. Indus Village eventually escapes this initial solution by calculating these variables from various datasets and submodels. However, for simplicity, we will use it as-is. As Wallach et al. (2019), we will use two global variables as constants (MUF and WP) and create some extra procedures that will help us register and apply a default configuration of the hyperparameters (parameters-check and parameters-to-default).
to set-constants
; "constants" are variables that will not be explored as parameters
; and may be used during a simulation.
; MUF : Water Uptake coefficient (mm^3 mm^-3)
set MUF 0.096
; WP : Water content at wilting Point (cm^3.cm^-3)
set WP 0.06
end
to set-parameters
; set random seed
random-seed seed
;;; load weather input data from file
load-weather-input-data-table
parameters-check
;;; weather parameters are left with default values, but effectively ignored given that input weather is used.
set albedo_min 1E-6 + random-float 0.3
set albedo_max albedo_min + random-float 0.3
if (experiment-type = "user-defined")
[
set albedo_min par_albedo_min
set albedo_max par_albedo_max
;;; Soil Water Balance model
set WHC_min water-holding-capacity_min
set WHC_max water-holding-capacity_max
set DC_min drainage-coefficient_min
set DC_max drainage-coefficient_max
set z_min root-zone-depth_min
set z_max root-zone-depth_max
set CN_min runoff-curve_min
set CN_max runoff-curve_max
]
if (experiment-type = "random")
[
set albedo_min 1E-6 + random-float 0.3
set albedo_max albedo_min + random-float 0.3
;;; Soil Water Balance model
set WHC_min random-float 0.1
set WHC_max WHC_min + random-float 0.1
set DC_min 1E-6 + random-float 0.45
set DC_max DC_min + random-float 0.45
set z_min random-float 1000
set z_max z_min + random-float 1000
set CN_min random-float 40
set CN_max CN_min + random-float 50
]
end
to parameters-check
;;; check if values were reset to 0 (NetLogo does that from time to time...!)
;;; and set default values (assuming they are not 0)
if (par_albedo_min = 0) [ set par_albedo_min 0.1 ]
if (par_albedo_max = 0) [ set par_albedo_max 0.5 ]
if (water-holding-capacity_min = 0) [ set water-holding-capacity_min 0.05 ]
if (water-holding-capacity_max = 0) [ set water-holding-capacity_max 0.25 ]
if (drainage-coefficient_min = 0) [ set drainage-coefficient_min 0.3 ]
if (drainage-coefficient_max = 0) [ set drainage-coefficient_max 0.7 ]
if (root-zone-depth_min = 0) [ set root-zone-depth_min 200 ]
if (root-zone-depth_max = 0) [ set root-zone-depth_max 2000 ]
if (runoff-curve_min = 0) [ set runoff-curve_min 30 ]
if (runoff-curve_max = 0) [ set runoff-curve_max 80 ]
end
to parameters-to-default
;;; set parameters to a default value
set par_albedo_min 0.1
set par_albedo_max 0.5
set water-holding-capacity_min 0.05
set water-holding-capacity_max 0.25
set drainage-coefficient_min 0.3
set drainage-coefficient_max 0.7
set root-zone-depth_min 200
set root-zone-depth_max 2000
set runoff-curve_min 30
set runoff-curve_max 80
end
...
to setup-soil-water-properties
ask patchesWithElevationData
[
set albedo albedo_min + random-float (albedo_max - albedo_min)
; Water Holding Capacity of the soil (cm^3 cm^-3).
set WHC WHC_min + random-float (WHC_max - WHC_min)
; DC : Drainage coefficient (mm^3 mm^-3)
set DC DC_min + random-float (DC_max - DC_min)
; z : root zone depth (mm)
set z z_min + random (z_max - z_min)
; CN : Runoff curve number
set CN CN_min + random (CN_max - CN_max)
; FC : Water content at field capacity (cm^3.cm^-3)
set FC WP + WHC
; WAT0 : Initial Water content (mm)
set WAT z * FC
]
end27.3 Calculating Reference Evapotranspiration (ETr)
In our previous implementation of set-day-weather-from-input-data, we must now add a new final step, where netSolarRadiation and ETr (reference evapotranspiration) are set for each patch. We also need to implement a procedure to estimate ETr based on an FAO standard (Allen et al. 1998).
Allen, Richard G., Luis S. Pereira, Dirk Raes, and Martin Smith. 1998. Crop Evapotranspiration - FAO Irrigation and Drainage Paper No. 56. Rome: FAO. http://www.fao.org/3/X0490E/x0490e00.htm.
to set-day-weather-from-input-data [ dayOfYear year ]
...
ask patchesWithElevationData
[
set netSolarRadiation (1 - albedo) * solarRadiation
set ETr get-ETr
]
end
...
to-report get-ETr
;;; useful references:
;;; Suleiman A A and Hoogenboom G 2007
;;; Comparison of Priestley-Taylor and FAO-56 Penman-Monteith for Daily Reference Evapotranspiration Estimation in Georgia
;;; J. Irrig. Drain. Eng. 133 175–82 Online: http://ascelibrary.org/doi/10.1061/%28ASCE%290733-9437%282007%29133%3A2%28175%29
;;; also: Jia et al. 2013 - doi:10.4172/2168-9768.1000112
;;; Allen, R. G., Pereira, L. A., Raes, D., and Smith, M. 1998.
;;; “Crop evapotranspiration.”FAO irrigation and drainage paper 56, FAO, Rome.
;;; also: http://www.fao.org/3/X0490E/x0490e07.htm
;;; constants found in: http://www.fao.org/3/X0490E/x0490e07.htm
;;; see also r package: Evapotranspiration (consult source code)
let windSpeed 2 ; as recommended by: http://www.fao.org/3/X0490E/x0490e07.htm#estimating%20missing%20climatic%20data
;;; estimation of saturated vapour pressure (e_s) and actual vapour pressure (e_a)
let e_s (get-vapour-pressure T_max + get-vapour-pressure T_min) / 2
let e_a get-vapour-pressure T_min
; ... in absence of dew point temperature, as recommended by
; http://www.fao.org/3/X0490E/x0490e07.htm#estimating%20missing%20climatic%20data
; however, possibly min temp > dew temp under arid conditions
;;; slope of the vapor pressure-temperature curve (kPa ºC−1)
let DELTA 4098 * (get-vapour-pressure T) / (T + 237.3) ^ 2
;;; latent heat of vaporisation = 2.45 MJ.kg^-1
let lambda 2.45
;;; specific heat at constant pressure, 1.013 10-3 [MJ kg-1 °C-1]
let c_p 1.013 * 10 ^ -3
;;; ratio molecular weight of water vapour/dry air
let epsilon 0.622
;;; atmospheric pressure (kPa)
let P 101.3 * ((293 - 0.0065 * elevation) / 293) ^ 5.26
;;; psychometric constant (kPa ºC−1)
let gamma c_p * P / (epsilon * lambda)
;;; Penman-Monteith equation from: fao.org/3/X0490E/x0490e0 ; and from: weap21.org/WebHelp/Mabia_Alg ETRef.htm
; 900 and 0.34 for the grass reference; 1600 and 0.38 for the alfalfa reference
let C_n 900
let C_d 0.34
let ETr_temp (0.408 * DELTA * netSolarRadiation + gamma * (C_n / (T + 273)) * windSpeed * (e_s - e_a)) / (DELTA + gamma * (1 + C_d * windSpeed))
report ETr_temp
end
to-report get-vapour-pressure [ temp ]
report (0.6108 * exp(17.27 * temp / (temp + 237.3)))
end27.4 Implementing the main algorithm
Next, we add update-WAT, which contains the calculations that finally outputs ARID.
to update-WAT
; Soil Water Balance model
; Using the approach of:
; 'Working with dynamic crop models: Methods, tools, and examples for agriculture and enviromnent'
; Daniel Wallach, David Makowski, James W. Jones, François Brun (2006, 2014, 2019)
; Model description in p. 24-28, R code example in p. 138-144.
; see also https://github.com/cran/ZeBook/blob/master/R/watbal.model.r
; Some additional info about run off at: https://engineering.purdue.edu/mapserve/LTHIA7/documentation/scs.htm
; and at: https://en.wikipedia.org/wiki/Runoff_curve_number
; Maximum abstraction (mm; for run off)
let S 25400 / CN - 254
; Initial Abstraction (mm; for run off)
let IA 0.2 * S
; WATfc : Maximum Water content at field capacity (mm)
let WATfc FC * z
; WATwp : Water content at wilting Point (mm)
let WATwp WP * z
; Change in Water Before Drainage (Precipitation - Runoff)
let RO 0
if (RAIN > IA)
[ set RO ((RAIN - 0.2 * S) ^ 2) / (RAIN + 0.8 * S) ]
; Calculating the amount of deep drainage
let DR 0
if (WAT + RAIN - RO > WATfc)
[ set DR DC * (WAT + RAIN - RO - WATfc) ]
; Calculate rate of change of state variable WAT
; Compute maximum water uptake by plant roots on a day, RWUM
let RWUM MUF * (WAT + RAIN - RO - DR - WATwp)
; Calculate the amount of water lost through transpiration (TR)
let TR min (list RWUM ETr)
let dWAT RAIN - RO - DR - TR
set WAT WAT + dWAT
set WATp WAT / z
set ARID 0
if (TR < ETr)
[ set ARID 1 - TR / ETr ]
end27.5 Merging with spatial data module
So far, we still lack spatial data. Let us implement the necessary code to import the processed data we output from the flows module.
extensions [ csv gis ]
...
breed [ sites site ]
breed [ flowHolders flowHolder ]
...
globals
[
patchesWithElevationData
noElevationDataTag
maxElevation
width
height
;;; GIS data holders
sitesData_EMIII-MMIA
sitesData_MMIB
elevationData
riversData
...
]
sites-own
[
name
siteType
period
]
patches-own
[
elevation ; elevation above sea level [m]
flow_direction ; the numeric code for the (main) direction of flow or
; drainage within the land unit.
; Following Jenson & Domingue (1988) convention:
; NW = 64, N = 128, NE = 1,
; W = 32, <CENTRE>, E = 2,
; SW = 16, S = 8, SE = 4
flow_receive ; Boolean variable stating whether or not the land unit receives
; the flow of a neighbour.
flow_accumulation ; the amount of flow units accumulated in the land unit.
; A Flow unit is the volume of runoff water flowing from one land unit
; to another (assumed constant and without losses).
flow_accumulationState ; the state of the land unit regarding the calculation of flow
; accumulation (auxiliary variable).
isRiver
...
]
...
to import-map-with-flows
import-world "data/terrainWithFlows/BlockC_module2_flows world.csv"
;;; reduce patch size in pixels
set-patch-size 3
end
...
27.6 Visualisation
Before advancing, we can implement a display procedure, expanding it also to be able to paint patches according to the new patch variables, using the “chooser” parameter display-mode:
to refresh-view
if (display-mode = "elevation")
[
ask patchesWithElevationData [ display-elevation ]
]
if (display-mode = "albedo")
[
ask patchesWithElevationData [ display-albedo ]
]
if (display-mode = "ETr")
[
let maxETr max [ETr] of patchesWithElevationData
ask patchesWithElevationData [ display-ETr maxETr ]
]
if (display-mode = "drainage coefficient (DC)")
[
ask patchesWithElevationData [ display-DC ]
]
if (display-mode = "root zone depth (z)")
[
let maxZ max [z] of patchesWithElevationData
ask patchesWithElevationData [ display-z maxZ ]
]
if (display-mode = "runoff curve number (CN)")
[
let maxCN max [CN] of patchesWithElevationData
ask patchesWithElevationData [ display-CN maxCN ]
]
if (display-mode = "water content at field capacity (FC)")
[
let maxFC max [FC] of patchesWithElevationData
ask patchesWithElevationData [ display-FC maxFC ]
]
if (display-mode = "water holding Capacity (WHC)")
[
let maxWHC max [WHC] of patchesWithElevationData
ask patchesWithElevationData [ display-WHC maxWHC ]
]
if (display-mode = "soil water content (WATp)")
[
let maxWATp max [WATp] of patchesWithElevationData
ask patchesWithElevationData [ display-WATp maxWATp ]
]
if (display-mode = "ARID coefficient")
[
ask patchesWithElevationData [ display-arid ]
]
end
to display-elevation
let elevationGradient 100 + (155 * (elevation / maxElevation))
set pcolor rgb (elevationGradient - 100) elevationGradient 0
end
to display-albedo
set pcolor 1 + 9 * albedo
end
to display-ETr [ maxETr ]
ifelse (maxETr = 0)
[ set pcolor 25 ]
[ set pcolor 22 + 6 * (1 - ETr / maxETr) ]
end
to display-DC
set pcolor 112 + 6 * (1 - DC)
end
to display-z [ maxZ ]
set pcolor 42 + 8 * (1 - z / maxZ)
end
to display-CN [ maxCN ]
set pcolor 72 + 6 * (1 - CN / maxCN)
end
to display-FC [ maxFC ]
set pcolor 82 + 6 * (1 - FC / maxFC)
end
to display-WHC [ maxWHC ]
set pcolor 92 + 6 * (1 - WHC / maxWHC)
end
to display-WATp [ maxWATp ]
set pcolor 102 + 6 * (1 - WATp / maxWATp)
end
to display-ARID
set pcolor 12 + 6 * ARID
end27.7 Combining with flow accumulation algorithm
Finally, we finalise this module by implementing a solution that uses flow_accumulation to modulate variation in ARID as a proxy for the effect of regional hydrology on the local soil water balance. The solution uses the patch variable ARID_modifier, set during setup according to the parameter ARID-decrease-per-flow-accumulation and the local relative flow accumulation (flow_accumulation / maxFlowAccumulation). It then modifies ARID each day as a simple scalar.
patches-own
[
...
ARID_modifier ; modifier coefficient based on the relative value of flow_accumulation
]
...
to setup
clear-all
; --- loading/testing parameters -----------
import-map-with-flows ; import-world must be the first step
set-constants
set-parameters
setup-patches
; --- core procedures ----------------------
set currentYear weatherInputData_firstYear
set currentDayOfYear 1
;;; values are taken from input data
set-day-weather-from-input-data currentDayOfYear currentYear
ask patchesWithElevationData [ update-WAT ]
; --- display & output handling ------------------------
refresh-view
; -- time -------------------------------------
reset-ticks
end
...
to setup-patches
setup-soil-water-properties
setup-ARID-modifier
end
...
to setup-ARID-modifier
ask patchesWithElevationData
[
set ARID_modifier (1 - ARID-decrease-per-flow-accumulation * (flow_accumulation / maxFlowAccumulation))
]
end
...
to go
; --- core procedures -------------------------
;;; values are taken from input data
set-day-weather-from-input-data currentDayOfYear currentYear
ask patchesWithElevationData [ update-WAT modify-ARID ]
; --- output handling ------------------------
refresh-view
; -- time -------------------------------------
advance-time
tick
; --- stop conditions -------------------------
if (currentYear = weatherInputData_lastYear and currentDayOfYear = last weatherInputData_yearLengthInDays) [stop]
end
...
to modify-ARID
set ARID ARID * ARID_modifier
end27.8 Checking the milestone File (module 4)
Our solution to link flow_accumulation and ARID is undoubtedly arbitrary. Such solutions should always be temporary and prompt further research and coding excursions. For the tutorial, however, we can move forward.
Screenshot of the ‘ARID’ module (tick 100)
Screenshot of the ‘ARID’ module (tick 150)
Screenshot of the ‘ARID’ module (tick 200)
See the fully implemented version of this module: BlockC_module3_ARID.nlogo.