The WaterNetwork component uses the America's Water Network to describe how downstream flows depend on upstream flows and withdrawals.

The America's Water Network is described in detail in its working paper. The data contained in the network and its loading are described in the sister waternet notebook.

The component definition is very simple:


In [1]:
using Mimi

@defcomp WaterNetwork begin
    gauges = Index()

    # External
    added = Parameter(index=[gauges, time], unit="1000 m^3") # Water added at node from runoff
    removed = Parameter(index=[gauges, time], unit="1000 m^3") # Water removed from node
    returned = Parameter(index=[gauges, time], unit="1000 m^3") # Water returns to a node from canals

    inflows = Variable(index=[gauges, time], unit="1000 m^3") # Sum of upstream outflows
    outflows = Variable(index=[gauges, time], unit="1000 m^3") # inflow + added - removed + returned
end


Out[1]:
WaterNetwork

All of these values are defined for each gauge, or node on the river network. They present a water balance: $$outflows_i(t) = \sum_{j \in N(i)} inflows_i(t) + added_i(t) - removed_i(t) + returned_i(t)$$

Where $N(i)$ are the neighbors of gauge $i$, where neighbors are defined as the upstream sources of flow.

Initialization

The added parameter is set to the runoff calculated by VIC and loaded in weather.jl.

In optimization mode, both the returned and removed parameters are set by the optimization to satisfy allocation demands.

In simulation mode,

  • The removed parameter is set to the contents of data/cache/counties/extraction/withdrawals.jld, if available.
  • The returned parameter is set to the contents of data/cache/counties/extraction/returns.jld, if available.

The removed and returned variables are generated by optimize-surface.jl, which finds plausible extractions to satisfy USGS demands.

Optimization

The logic for optimization is complicated by the process of matching canal-based withdrawals to gauge-based flows. However, aside from this, the entire network follows two simple principes.

  • Withdrawal applied to any gauge withdraws water from all downstream gauges.
  • Runoff added to any guage increases the allottment at all downstream gauges.

The following shows these gauge-specific allotments, and some interesting caveats.

First, the allotments take a while to calculate (which they are done by constraintoffset_waternetwork_outflows(m); in optimize.jl if need be), so they are cached in the file data/partialhouse2.jld.


In [3]:
using OptiMimi

cwro = deserialize(open("../data/partialhouse2.jld", "r"))


Out[3]:
OptiMimi.LinearProgrammingHall(:WaterNetwork,:outflows,[NaN,NaN,NaN,NaN,1363.79,1975.22,NaN,NaN,14.7898,84.0744  …  NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN])

This is a Linear Programming Hall file, used for constraint vectors.

Additional informatino on the gauges is stored in the waternet.RData file, generated outside of Julia.


In [4]:
using DataFrames

stations = read_rda("../data/waternet.RData", convertdataframes=true)["network"]


Out[4]:
collectioncolidlatlonelevnextptdist
1rivdis34440.5-124.111.0175558646.0
2rivdis34542.58-124.0635.0NANA
3rivdis34641.51-123.962.0182098794.0
4rivdis34743.58-123.5628.019439671.0
5rivdis34846.93-123.317.08278927.0
6rivdis34944.95-123.0332.0202601583.0
7rivdis35046.26-122.916.02075423709.0
8rivdis35440.7-118.21259.074621163.0
9rivdis35646.43-117.16200.0929338087.0
10rivdis35748.18-117.03610.021929927.0
11rivdis35832.73-114.6330.01354935958.0
12rivdis35939.36-112.031506.0106252570.0
13rivdis36036.86-111.58947.0154363115.0
14rivdis36138.98-110.151231.01644428841.0
15rivdis36248.11-104.46574.02179250852.0
16rivdis36347.61-104.16573.0216724412.0
17rivdis36429.83-101.38353.0128191608.0
18rivdis36527.5-99.5106.01245545890.0
19rivdis36628.03-97.859998.06502927.0
20rivdis36742.86-97.39999347.01901813636.0
21rivdis36847.93-97.08002237.0218308906.0
22rivdis36941.03-96.29999307.01791642335.0
23rivdis37029.3-96.1000120.01272983028.0
24rivdis37136.15-96.0188.015112748.0
25rivdis37240.66-95.83002276.0176524659.0
26rivdis37329.58-95.7512.01171858740.0
27rivdis37435.26-95.23001146.01471329808.0
28rivdis37538.9-94.88231.0523710977.0
29rivdis37630.43-94.8500111.01284386805.46218764059
30rivdis37730.35-94.100013.01290814759.0
&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip

Unfortunately, the allotment entries in the optimization problem are neither in the same order (they're in vertex order), nor have all of the entries that are in the water network.


In [5]:
println(nrow(stations))
println(length(cwro.f))


22619
45118

We need to use the wateridverts dictionary, which maps from gauge IDs to vertices in the network.


In [7]:
using Graphs

wateridverts = deserialize(open("../data/cache/wateridverts.jld", "r"))
wateridverts
#waternet = deserialize(open("../data/cache/waternet.jld", "r"))


Out[7]:
Dict{UTF8String,Graphs.ExVertex} with 22559 entries:
  "usgs.02456000"     => vertex [3630] "usgs.02456000"
  "junction.13805-dn" => vertex [11658] "junction.13805-dn"
  "junction.11142-up" => vertex [12918] "junction.11142-up"
  "usgs.03488000"     => vertex [5098] "usgs.03488000"
  "usgs.12113347"     => vertex [13374] "usgs.12113347"
  "junction.13981-up" => vertex [20305] "junction.13981-up"
  "junction.12762-up" => vertex [4676] "junction.12762-up"
  "usgs.0214627970"   => vertex [2418] "usgs.0214627970"
  "usgs.01396800"     => vertex [975] "usgs.01396800"
  "usgs.05532000"     => vertex [6898] "usgs.05532000"
  "junction.19865-dn" => vertex [7569] "junction.19865-dn"
  "junction.9399-dn"  => vertex [4855] "junction.9399-dn"
  "usgs.03025000"     => vertex [3794] "usgs.03025000"
  "usgs.09067200"     => vertex [10866] "usgs.09067200"
  "usgs.08261000"     => vertex [10597] "usgs.08261000"
  "usgs.04286000"     => vertex [6117] "usgs.04286000"
  "junction.22562-up" => vertex [22524] "junction.22562-up"
  "junction.7795-dn"  => vertex [16679] "junction.7795-dn"
  "usgs.11313472"     => vertex [12777] "usgs.11313472"
  "usgs.14337800"     => vertex [14571] "usgs.14337800"
  "reservoir.2566"    => vertex [17636] "reservoir.2566"
  "junction.12093-dn" => vertex [20077] "junction.12093-dn"
  "junction.14022-up" => vertex [20314] "junction.14022-up"
  "canal.09528800"    => vertex [18159] "canal.09528800"
  "usgs.11433200"     => vertex [13026] "usgs.11433200"
  ⋮                   => ⋮

And now we can collect all of the information to produce a map.


In [8]:
lats = []
lons = []
runoff = []
for ii in 1:nrow(stations)
    gaugeid = "$(stations[ii, :collection]).$(stations[ii, :colid])"
    if gaugeid in keys(wateridverts)
        runoff = [runoff; cwro.f[vertex_index(wateridverts[gaugeid])]]
        lats = [lats; stations[ii, :lat]]
        lons = [lons; stations[ii, :lon]]
    end
end

runoff[isna(runoff)] = 0


Out[8]:
0

In [44]:
using RCall

cex = convert(Vector{Float64}, runoff) / 1e9

R"library(maps)"
R"map('usa', xlim=c(-140, -60), ylim=c(20, 50))"
R"points($lons, $lats, cex=$cex)"


Out[44]:
RCall.RObject{RCall.NilSxp}
NULL