Repository: WikiWatershed/tr-55

Issue:
Add simple water quality routine.
Barry Evans has provided an algorithm to extend TR-55 results with simple water quality predictions based on the same workflow of the existing model.

Discussion:

Patch:
diff --git a/README.md b/README.md
index 72023c1..0bda9ed 100644
--- a/README.md
+++ b/README.md
@@ -4,73 +4,174 @@ A Python implementation of TR-55.
 
 ## Functions
 
-`simulate_all_tiles` and `simulate_year` are the two functions intended for consumers of this module.
+`simulate_day`, `simulate_year`, `simulate_modifications` are the three functions most likely to be of direct interest for users of this module.
 
-### `simulate_all_tiles`
+### `simulate_day`
 
-The `simulate_all_tiles` function simulates the events of a single day.  It takes three arguments:
+The `simulate_day` function simulates the events of a single day.  It takes four arguments:
 
-   1. `parameters` which is either a `date` object (which is mapped to *precipitation* and *evapotranspiration* using sample year data) or that argument is a single day's *precipitation* and *evapotranspiration* as a tuple of two floats (in units of inches).
+   1. `day` has one of two types.  If it is an integer, it representing the day of the year (where `day == 0` represents October 15th) and is used to retrieve the precipitation and evapotranspiration for that day from the sample year table.  Alternatively, `day` can be a tuple of precipitation and evapotranspiration.
 
       If the sample year data are to be considered typical, *precipitation* should generally be between 0 and 3 inches daily, and *evapotranspiration* should be between 0 and 0.2 inches daily.
 
-   2. `tile_census` is a dictionary containing the number of appearances of each type of tile in the query area, e.g.
+   2. `tile_census` is a dictionary containing the number of appearances of each type of tile in the query area, as well as modifications, e.g.
    ```Python
-       {
-           "result": {
-               "cell_count": 147,
-               "distribution": {
-                   "d:hi_residential": 33,
-                   "c:commercial": 42,
-                   "a:deciduous_forest": 72
-               }
-           }
-       }
+    {
+        "cell_count": 147,
+        "distribution": {
+            "d:hi_residential": {"cell_count": 33},
+            "c:commercial": {"cell_count": 42},
+            "a:deciduous_forest": {"cell_count": 72}
+        },
+        "modifications": {
+            ":no_till": {
+                "cell_count": 30,
+                "distribution": {
+                    "d:hi_residential": {"cell_count": 10},
+                    "c:commercial": {"cell_count": 20}
+                }
+            },
+            "d:rock": {
+                "cell_count": 5,
+                "distribution": {
+                    "a:deciduous_forest": {"cell_count": 5}
+                }
+            }
+        }
+    }
    ```
-   Here, `"d:hi_residential"` indicates *High-Intensity Residential* land use on top of *Type D* soil.
+   Here, `"d:hi_residential"` indicates *High-Intensity Residential* land use on top of *Type D* soil, `"c:commerical"` indicates *Commercial* land use on top of *Type C* soil, and so on.
 
-   3. `pre_columbian` is a boolean argument which controls whether or not to simulate the given area under pre-Columbian conditions.  When this is set to true, all land-uses other than *water* and *wetland* are treated as *mixed forest*.
+   The `":no_till"` key and its value say that 10 units of `"d:hi_residential"` should be turned into `"d:no_till"` and 20 units of `"c:commercial"` should be turned into `"c:no_till"`.  The `"d:rock"` key and its value say that 5 units of `"a:deciduous_forest"` should be reclassified as `"d:rock"`.
 
-The algorithm used is close to TR-55 (the algorithm found in [the USDA's Technical Release 55, revised 1986](http://www.cpesc.org/reference/tr55.pdf)), but with a few differences.  The main differences are:
+   3. `subst` is a substitution that should be performed on the census.  It is used to implement the BMPs/reclassifications mentioned above.  When `subst` is a colon followed by a BMP, then the BMP is superimposed on all of the tiles in the census.  When `subst` is a soil type and a land use, then all tiles in the census are treated as if they have that soil type and land use.
+
+   4. `pre_columbian` is a boolean argument which controls whether or not to simulate the given area under pre-Columbian conditions.  When it set to true, all land-uses other than *water* and *wetland* are treated as *mixed forest*.
+
+The algorithm implemented by this code is close to TR-55, the algorithm found in [the USDA's Technical Release 55, revised 1986](http://www.cpesc.org/reference/tr55.pdf), but with a few differences.  The main differences are:
 
    * it applies the TR-55 algorithm on a tile-by-tile basis and aggregates the results rather than running the TR-55 algorithm once over the entire area
    * it uses the *Pitt Small Storm Hydrology Model* on tiles which are of a "built-type" when precipitation is two inches or less
 
-and there are numerous other small differences, as well.
+and there are numerous other small differences.
 
-`simulate_all_tiles` returns a triple of *runoff*, *evapotranspiration*, and *infiltration* (all in units of inches).
+`simulate_day` returns a dictionary with `runoff`, `et`, and `inf` keys.  All three of these numbers are in units of inches.
 
 ### `simulate_year`
 
-The `simulate_year` function takes two parameters, `tile_census` and `pre_columbian` as described above, and simulates an entire year using a year's worth of sample precipitation and evapotranspiration data.  It returns a triple of *runoff*, *evapotranspiration*, and *infiltration* for the given area over the entire sample year.
+The `simulate_year` and simulates an entire year using one year of sample precipitation and evapotranspiration data.  It returns a dictionary of with `runoff`, `et`, and `inf` keys, and also keys which map associated with various pollutant loads.
+
+The function takes four parameters.  `tile_census`, `subst`, and `pre_columbian` as described above.  The `cell_res` argument gives the resolution in meters of the data.
+
+### `simulate_modifications`
+
+This function takes three parameters, `tile_census`, `cell_res`, and `pre_columbian`, all as described above.  The output of this function is similar to that of `simulate_year`, except it shows the results both with and without modifications.
 
 ## Usage Example
 
 The following program:
 ```Python
-from tr55.model import simulate_year
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
+import pprint
+
+from tr55.model import simulate_modifications
 
 tiles = {
-    "result": {
-        "cell_count": 147,
-        "distribution": {
-            "d:hi_residential": 33,
-            "c:commercial": 42,
-            "a:deciduous_forest": 72
+    "cell_count": 147,
+    "distribution": {
+        "d:hi_residential": {"cell_count": 33},
+        "c:commercial": {"cell_count": 42},
+        "a:deciduous_forest": {"cell_count": 72}
+    },
+    "modifications": {
+        ":no_till": {
+            "cell_count": 30,
+            "distribution": {
+                "d:hi_residential": {"cell_count": 10},
+                "c:commercial": {"cell_count": 20}
+            }
+        },
+        "d:rock": {
+            "cell_count": 5,
+            "distribution": {
+                "a:deciduous_forest": {"cell_count": 5}
+            }
         }
     }
 }
 
-(q, et, inf) = simulate_year(tiles)
-print("%f inches of runoff" % q)
-print("%f inches of et" % et)
-print("%f inches of infiltration" % inf)
+pprint.pprint(simulate_modifications(tiles))
 ```
 should produce the following output:
-```
-17.614212 inches of runoff
-15.167862 inches of et
-20.245580 inches of infiltration
+```Python
+{u'modified': {u'bod': 1155.5832690023512,
+               u'cell_count': 147,
+               u'distribution': {u'a:deciduous_forest': {u'bod': 1.9676947491691286,
+                                                         u'cell_count': 72,
+                                                         u'et': 24.675262500000017,
+                                                         u'inf': 34.53218880673186,
+                                                         u'runoff': 0.7999320266014792,
+                                                         u'tn': 0.006128885284297285,
+                                                         u'tp': 0.00019354374581991428,
+                                                         u'tss': 1.258034347829443},
+                                 u'c:commercial': {u'bod': 619.3321622220828,
+                                                   u'cell_count': 42,
+                                                   u'et': 17.42526,
+                                                   u'inf': 16.198443763932524,
+                                                   u'runoff': 21.23295052178176,
+                                                   u'tn': 4.545098932436254,
+                                                   u'tp': 0.7192244464514511,
+                                                   u'tss': 126.3137934080361},
+                                 u'd:hi_residential': {u'bod': 534.2834120310993,
+                                                       u'cell_count': 33,
+                                                       u'et': 15.08352545454543,
+                                                       u'inf': 13.967458770273502,
+                                                       u'runoff': 24.498158107332266,
+                                                       u'tn': 3.101560485095788,
+                                                       u'tp': 0.5206999354540374,
+                                                       u'tss': 63.842339912190674}},
+               u'et': 20.450586122448982,
+               u'inf': 24.67740388835977,
+               u'runoff': 11.957947247429287,
+               u'tn': 7.652788302816338,
+               u'tp': 1.2401179256513084,
+               u'tss': 191.41416766805622},
+ u'unmodified': {u'bod': 1762.555509187117,
+                 u'cell_count': 147,
+                 u'distribution': {u'a:deciduous_forest': {u'bod': 0.0,
+                                                           u'cell_count': 72,
+                                                           u'et': 26.51670000000007,
+                                                           u'inf': 34.91810000000001,
+                                                           u'runoff': 0.0,
+                                                           u'tn': 0.0,
+                                                           u'tp': 0.0,
+                                                           u'tss': 0.0},
+                                   u'c:commercial': {u'bod': 1061.0138827829708,
+                                                     u'cell_count': 42,
+                                                     u'et': 2.272860000000005,
+                                                     u'inf': 4.430079770137537,
+                                                     u'runoff': 36.375400229862464,
+                                                     u'tn': 7.786472849455671,
+                                                     u'tp': 1.2321451541995787,
+                                                     u'tss': 216.39549270630107},
+                                   u'd:hi_residential': {u'bod': 701.5416264041463,
+                                                         u'cell_count': 33,
+                                                         u'et': 6.818579999999993,
+                                                         u'inf': 8.361629213022574,
+                                                         u'runoff': 32.167342828759615,
+                                                         u'tn': 4.072508593956273,
+                                                         u'tp': 0.6837058223430238,
+                                                         u'tss': 83.8282790872751}},
+                 u'et': 15.167861632653095,
+                 u'inf': 20.24558036990151,
+                 u'runoff': 17.614211721110824,
+                 u'tn': 11.858981443411944,
+                 u'tp': 1.9158509765426026,
+                 u'tss': 300.2237717935762}}
 ```
 
 ## Testing
diff --git a/tr55/model.py b/tr55/model.py
index 4e4cfc9..3ec4b0e 100644
--- a/tr55/model.py
+++ b/tr55/model.py
@@ -1,3 +1,8 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
 """
 TR-55 Model Implementation
 
@@ -10,32 +15,34 @@ and variables used in this program are as follows:
  * `init_abs` is Ia, the initial abstraction, another form of infiltration
 """
 
-from datetime import date, timedelta
-from tr55.tablelookup import days_in_sample_year
-from tr55.tablelookup import lookup_et, lookup_p, lookup_cn
-from tr55.tablelookup import lookup_bmp_infiltration, is_bmp, is_built_type
+import sys
+
+from tr55.tablelookup import lookup_pet, lookup_cn, lookup_bmp_infiltration, \
+    is_bmp, is_built_type, precolumbian, get_pollutants
+from tr55.water_quality import get_volume_of_runoff, get_pollutant_load
+from tr55.operations import dict_minus, dict_plus
 
 
 def runoff_pitt(precip, land_use):
     """
-    The Pitt Small Storm Hydrology method.  This comes from Table D in
-    the 2010/12/27 document.  The output is a runoff value in inches.
+    The Pitt Small Storm Hydrology method.  The output is a runoff
+    value in inches.
     """
-    co1 = +3.638858398e-2
-    co2 = -1.243464039e-1
-    co3 = +1.295682223e-1
-    co4 = +9.375868043e-1
-    co5 = -2.235170859e-2
-    co6 = +0.170228067e0
-    co7 = -3.971810782e-1
-    co8 = +3.887275538e-1
-    co9 = -2.289321859e-2
-    pr4 = pow(precip, 4)
-    pr3 = pow(precip, 3)
-    pr2 = pow(precip, 2)
-
-    impervious = (co1 * pr3) + (co2 * pr2) + (co3 * precip) + co4
-    urb_grass = (co5 * pr4) + (co6 * pr3) + (co7 * pr2) + (co8 * precip) + co9
+    c1 = +3.638858398e-2
+    c2 = -1.243464039e-1
+    c3 = +1.295682223e-1
+    c4 = +9.375868043e-1
+    c5 = -2.235170859e-2
+    c6 = +0.170228067e0
+    c7 = -3.971810782e-1
+    c8 = +3.887275538e-1
+    c9 = -2.289321859e-2
+    p4 = pow(precip, 4)
+    p3 = pow(precip, 3)
+    p2 = pow(precip, 2)
+
+    impervious = (c1 * p3) + (c2 * p2) + (c3 * precip) + c4
+    urb_grass = (c5 * p4) + (c6 * p3) + (c7 * p2) + (c8 * precip) + c9
 
     runoff_vals = {
         'water':          impervious,
@@ -47,8 +54,9 @@ def runoff_pitt(precip, land_use):
         'transportation': impervious,
         'urban_grass':     urb_grass
     }
+
     if land_use not in runoff_vals:
-        raise Exception('Land use %s not a built-type' % land_use)
+        raise Exception('Land use %s not a built-type.' % land_use)
     else:
         return min(runoff_vals[land_use], precip)
 
@@ -71,7 +79,7 @@ def runoff_nrcs(precip, evaptrans, soil_type, land_use):
     """
     curve_number = lookup_cn(soil_type, land_use)
     if nrcs_cutoff(precip, curve_number):
-        return 0
+        return 0.0
     potential_retention = (1000.0 / curve_number) - 10
     initial_abs = 0.2 * potential_retention
     precip_minus_initial_abs = precip - initial_abs
@@ -81,55 +89,43 @@ def runoff_nrcs(precip, evaptrans, soil_type, land_use):
     return min(runoff, precip - evaptrans)
 
 
-def simulate_tile(parameters, tile_string, pre_columbian=False):
+def simulate_tile(parameters, tile_string):
     """
-    Simulate a tile on a given day using the method given in the
-    flowchart 2011_06_16_Stroud_model_diagram_revised.PNG and as
-    revised by various emails.
+    Simulate a single tile with some particular precipitation and
+    evapotranspiration.
 
-    The first argument can be one of two types.  It can either be a
-    date object, in which case the precipitation and
-    evapotranspiration are looked up from the sample year table.
-    Alternatively, those two values can be supplied directly via this
-    argument as a tuple.
+    The first argument contains the precipitation and
+    evapotranspiration as a tuple.
 
     The second argument is a string which contains a soil type and
-    land use separted by a colon.
-
-    The third argument is a boolean which is true if pre-Columbian
-    circumstances are to be simulated and false otherwise.
+    land use separated by a colon.
 
-    The return value is a triple of runoff, evapotranspiration, and
+    The return value is a dictionary of runoff, evapotranspiration, and
     infiltration.
     """
-    tile_string = tile_string.lower();
-    soil_type, land_use = tile_string.split(':')
-
-    pre_columbian_land_uses = set([
-        'water',
-        'woody_wetland',
-        'herbaceous_wetland'
-    ])
-
-    if pre_columbian:
-        if land_use not in pre_columbian_land_uses:
-            land_use = 'mixed_forest'
-
-    if type(parameters) is date:
-        precip = lookup_p(parameters)  # precipitation
-        evaptrans = lookup_et(parameters, land_use)  # evapotranspiration
-    elif type(parameters) is tuple:
+    if type(parameters) is tuple:
         precip, evaptrans = parameters
     else:
-        raise Exception('First argument must be a date or a (P,ET) pair')
+        raise Exception('First argument must be a (P,ET) pair')
+
+    tile_string = tile_string.lower()
+    soil_type, land_use = tile_string.split(':')
 
     if precip == 0.0:
-        return (0.0, evaptrans, 0.0)
+        return {
+            'runoff': 0.0,
+            'et': evaptrans,
+            'inf': 0.0,
+        }
 
     if is_bmp(land_use) and land_use != 'rain_garden':
         inf = lookup_bmp_infiltration(soil_type, land_use)  # infiltration
         runoff = precip - (evaptrans + inf)  # runoff
-        return (runoff, evaptrans, inf)  # Q, ET, Inf.
+        return {
+            'runoff': runoff,
+            'et': evaptrans,
+            'inf': inf,
+        }
     elif land_use == 'rain_garden':
         # Here, return a mixture of 20% ideal rain garden and 80% high
         # intensity residential.
@@ -137,9 +133,11 @@ def simulate_tile(parameters, tile_string, pre_columbian=False):
         runoff = precip - (evaptrans + inf)
         hi_res_tile = soil_type + ':hi_residential'
         hi_res = simulate_tile((precip, evaptrans), hi_res_tile)
-        return (0.2 * runoff + 0.8 * hi_res[0],
-                0.2 * evaptrans + 0.8 * hi_res[1],
-                0.2 * inf + 0.8 * hi_res[2])
+        return {
+            'runoff': 0.2 * runoff + 0.8 * hi_res[0],
+            'et': 0.2 * evaptrans + 0.8 * hi_res[1],
+            'inf': 0.2 * inf + 0.8 * hi_res[2],
+        }
 
     if is_built_type(land_use) and precip <= 2.0:
         runoff = runoff_pitt(precip, land_use)
@@ -150,75 +148,176 @@ def simulate_tile(parameters, tile_string, pre_columbian=False):
     else:
         runoff = runoff_nrcs(precip, evaptrans, soil_type, land_use)
     inf = precip - (evaptrans + runoff)
-    return (runoff, evaptrans, max(inf, 0.0))
 
+    return {
+        'runoff': runoff,
+        'et': evaptrans,
+        'inf': max(inf, 0.0),
+    }
 
-def simulate_all_tiles(parameters, tile_census, pre_columbian=False):
+
+def simulate_day(day, tile_census, subst=None, pre_columbian=False):
     """
     Simulate each tile for one day and return the overall results.
 
-    The first argument is either a day or a P,ET double (as in simulate_tile).
+    The first argument is an integer representing the day of the year
+    (day == 0 represents October 15th) or a precip, et pair.
 
-    The second argument is a dictionary (presumably converted from
-    JSON) that gives the number of each type of tile in the query
-    polygon.
+    The second argument is a dictionary that gives a census of all of the
+    tiles in the query area.
 
-    The output is a runoff, evapotranspiration, infiltration triple
-    which is an average of those produced by all of the tiles.
-    """
-    if 'result' not in tile_census:
-        raise Exception('No "result" key')
-    elif 'cell_count' not in tile_census['result']:
-        raise Exception('No "result.cell_count" key')
-    elif 'distribution' not in tile_census['result']:
-        raise Exception('No "result.distribution" key')
+    The third argument is the tile string substitution to be applied
+    to each tile.  This is used for simulating BMPs and
+    reclassifications.
 
-    global_count = tile_census['result']['cell_count']
+    The fourth argument is a boolean which is true if pre-Columbian
+    circumstances are to be simulated and false otherwise.  When this
+    argument is true, all land uses except for water and wetland are
+    transformed to mixed forest.
 
-    def simulate(tile_string, local_count):
-        """
-        A local helper function which captures various values.
-        """
-        return [(x * local_count) / global_count
-                for x in simulate_tile(parameters, tile_string, pre_columbian)]
+    The output is a runoff, evapotranspiration, infiltration dictionary.
+    """
+    if 'cell_count' not in tile_census:
+        raise Exception('No "cell_count" key.')
+    elif 'distribution' not in tile_census:
+        raise Exception('No "distribution" key.')
 
-    def simulate_plus(left, right):
+    def simulate(tile, n):
         """
-        Add two simulations.
+        Return the three values in units of inch-tiles instead of inches.
+        Inches are an inconvenient unit to work with.
         """
-        (a, b, c) = left
-        (x, y, z) = right
-        return (a+x, b+y, c+z)
+        (soil_type, land_use) = tile.split(':')
+
+        # If a substitution has been supplied, apply it.
+        if sys.version_info.major == 3:
+            types = (str)
+        else:
+            types = (str, unicode)
+        if isinstance(subst, types) and subst.find(':') >= 0:
+            (subst_soil_type, subst_land_use) = subst.split(':')
+            if len(subst_soil_type) > 0:
+                soil_type = subst_soil_type
+            land_use = subst_land_use
+
+        # If a Pre-Columbian simulation has been requested, change the
+        # land use type.
+        if pre_columbian:
+            land_use = precolumbian(land_use)
+
+        # Retrieve the precipitation and evapotranspiration, then run
+        # the simulation.
+        if isinstance(day, int):
+            parameters = lookup_pet(day, land_use)
+        else:
+            parameters = day
+
+        retval = simulate_tile(parameters, soil_type + ':' + land_use)
+        for key in retval.keys():
+            retval[key] *= n
+        return retval
+
+    results = {tile: simulate(tile, d['cell_count'])
+               for tile, d in tile_census['distribution'].items()}
+
+    return {
+        'distribution': results
+    }
 
-    results = [simulate(tile, n)
-               for tile, n in tile_census['result']['distribution'].items()]
 
-    area_sum = (0.0, 0.0, 0.0)
-    for result in results:
-        area_sum = simulate_plus(area_sum, result)
+def simulate_year(tile_census, cell_res=10, subst=None, pre_columbian=False):
+    """
+    Simulate an entire year, including water quality.
 
-    return area_sum
+    The first argument is a tile census as described in `simulate_day`.
 
+    The second argument is the tile resolution in meters.
 
-def simulate_year(tile_census, pre_columbian=False):
-    """
-    Simulate an entire year.
-    """
-    year = [date(1, 1, 1) + timedelta(days=i)
-            for i in range(days_in_sample_year())]
-    simulated_year = [simulate_all_tiles(day, tile_census, pre_columbian)
-                      for day in year]
+    The third parameter is the tile string substitution to be performed.
 
-    def day_plus(one, two):
-        """
-        Add two simulated days.
-        """
-        (runoff1, et1, inf1) = one
-        (runoff2, et2, inf2) = two
-        return (runoff1 + runoff2, et1 + et2, inf1 + inf2)
+    The fourth parameter is as in `simulate_day`.
+    """
+    # perform a daily simulation for each day of the year
+    simulated_year = [simulate_day(day, tile_census, subst, pre_columbian)
+                      for day in range(365)]
 
-    year_sum = (0.0, 0.0, 0.0)
+    # add those simulations together
+    year_sum = {}
     for day in simulated_year:
-        year_sum = day_plus(year_sum, day)
+        year_sum = dict_plus(year_sum, day)
 
-    return year_sum
+    # add the results into the census
+    retval = dict_plus(year_sum, tile_census)
+
+    # get the land use after the modification.
+    if isinstance(subst, str):
+        use_after = subst.split(':')[1]
+    else:
+        user_after = None  # noqa
+
+    # perform water quality computations
+    for (pair, result) in retval['distribution'].items():
+        use_before = pair.split(':')[1]
+        runoff = result['runoff']
+        n = result['cell_count']
+        liters = get_volume_of_runoff(runoff / n, n, cell_res)
+        for pol in get_pollutants():
+            # Try to compute the pollutant load with the land use
+            # after the modifications.  If that is impossible (there
+            # is no modification, the use is modified to a BMP, et
+            # cetera) then use the previous land use.
+            try:
+                load = get_pollutant_load(use_after, pol, liters)
+            except:
+                load = get_pollutant_load(use_before, pol, liters)
+            result[pol] = load
+
+    return retval
+
+
+def simulate_modifications(tile_census, cell_res=10, pre_columbian=False):
+    """
+    Simulate an entire year, including effects of modifications.
+    """
+    def rescale(result):
+        if isinstance(result, dict):
+            if 'cell_count' in result:
+                n = float(result['cell_count'])
+                result['runoff'] /= n
+                result['et'] /= n
+                result['inf'] /= n
+            for (key, val) in result.items():
+                rescale(val)
+
+    def tally_subresults(result):
+        retval = result.copy()
+        total = {}
+        for (pair, subresult) in result['distribution'].items():
+            total = dict_plus(total, subresult)
+        retval.update(total)
+        retval.pop('modifications', None)
+        rescale(retval)
+        return retval
+
+    # compute unmodified results
+    orig_result = simulate_year(tile_census, cell_res, None, pre_columbian)
+
+    # Compute the census ex-modifications and the subresults for each
+    # of the modifications.
+    mod_census = tile_census.copy()
+    subresults = []
+    if 'modifications' in tile_census:
+        for (subst, census) in tile_census['modifications'].items():
+            mod_census = dict_minus(mod_census, census)
+            subresult = simulate_year(census, cell_res, subst, pre_columbian)
+            subresults.append(subresult)
+
+    # Compute the result with modifications taken into account.
+    mod_result = simulate_year(mod_census, cell_res, None, pre_columbian)
+    for subresult in subresults:
+        mod_result = dict_plus(mod_result, subresult)
+
+    return {
+        'unmodified': tally_subresults(orig_result),
+        'modified': tally_subresults(mod_result)
+    }
diff --git a/tr55/operations.py b/tr55/operations.py
new file mode 100644
index 0000000..ae66328
--- /dev/null
+++ b/tr55/operations.py
@@ -0,0 +1,66 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
+import sys
+
+
+def tandem_walk(op, pred, left, right):
+    """
+    Walk two similarly-structured dictionaries in tandem, performing
+    the given operation when both halves satisfy the given predicate.
+    Otherwise, if both haves are dictionaries, recurse.
+    """
+    if pred(left) and pred(right):
+        return op(left, right)
+    elif isinstance(left, dict) and isinstance(right, dict):
+        retval = left.copy()
+        left_set = set(left.keys())
+        right_set = set(right.keys())
+        intersection = left_set & right_set
+        difference = right_set - left_set
+        for key in intersection:
+            retval[key] = tandem_walk(op, pred, left[key], right[key])
+        for key in difference:
+            retval[key] = right[key]
+        return retval
+
+
+def isnumber(obj):
+    """
+    Is obj a number?
+    """
+    if sys.version_info.major == 3:
+        types = (int, float)
+    else:
+        types = (int, long, float)
+    return isinstance(obj, types)
+
+
+def plus(x, y):
+    """
+    Sum of two numbers.
+    """
+    return x + y
+
+
+def minus(x, y):
+    """
+    Difference of two numbers.
+    """
+    return x - y
+
+
+def dict_plus(left, right):
+    """
+    Sum of two similarly-structured dictionaries.
+    """
+    return tandem_walk(plus, isnumber, left, right)
+
+
+def dict_minus(left, right):
+    """
+    Difference of two similarly-structured dictionaries.
+    """
+    return tandem_walk(minus, isnumber, left, right)
diff --git a/tr55/tablelookup.py b/tr55/tablelookup.py
index 562b55c..eb5115e 100644
--- a/tr55/tablelookup.py
+++ b/tr55/tablelookup.py
@@ -1,80 +1,42 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
 """
 Various routines to do table lookups.
 """
 
-from datetime import date
-from tr55.tables import SAMPLE_YEAR, BMPS, BUILT_TYPES, LAND_USE_VALUES
+from tr55.tables import SAMPLE_YEAR, BMPS, BUILT_TYPES, LAND_USE_VALUES, \
+    PRE_COLUMBIAN_LAND_USES, POLLUTANTS, POLLUTION_LOADS
 
 
 def lookup_ki(land_use):
     """
     Lookup the landuse coefficient.
     """
-    if land_use not in LAND_USE_VALUES or 'ki' not in LAND_USE_VALUES[land_use]:
+    if land_use not in LAND_USE_VALUES \
+       or 'ki' not in LAND_USE_VALUES[land_use]:
         raise Exception('Unknown land use: %s' % land_use)
     else:
         return LAND_USE_VALUES[land_use]['ki']
 
 
-def lookup_et(simulation_day, land_use):
+def lookup_pet(day, land_use):
     """
     Lookup/compute evapotranspiration from the tables.
     """
-    fixed_year = SAMPLE_YEAR['year_start'].year
-    simulation_day = date(fixed_year, simulation_day.month, simulation_day.day)
-
-    # Compute $ET_{\max}$ based on the time of year
-    grow_start = SAMPLE_YEAR['growing_start']
-    grow_end = SAMPLE_YEAR['growing_end']
-    grow_et = SAMPLE_YEAR['growing_etmax']
-    nongrow_et = SAMPLE_YEAR['nongrowing_etmax']
-    if grow_start <= simulation_day and simulation_day <= grow_end:
-        et_max = grow_et
-    else:
-        et_max = nongrow_et
-
-    # Compute the landuse coefficient
-    landuse_coefficient = LAND_USE_VALUES[land_use]['ki']
-
-    # Report $ET$, the evapotranspiration
-    return et_max * landuse_coefficient
-
-
-def days_in_sample_year():
-    """
-    Return the number of days in the sample year.
-    """
-    return SAMPLE_YEAR['days_per_year']
-
-
-def lookup_p(simulation_day):
-    """
-    Lookup percipitation from the SAMPLE_YEAR table.
-    """
-    fixed_year = SAMPLE_YEAR['year_start'].year
-    simulation_day = date(fixed_year, simulation_day.month, simulation_day.day)
-    seconds_per_day = 60 * 60 * 24
-    days_per_year = SAMPLE_YEAR['days_per_year']
-    day_delta = simulation_day - SAMPLE_YEAR['year_start']
-    seconds_from_start = day_delta.total_seconds()
-    days_from_start = int(seconds_from_start / seconds_per_day)
-
-    days = (days_from_start + days_per_year) % days_per_year
-    # This is a bit unattractive.  If it turns out to be a performance
-    # problem (which is unlikely), an interval tree can be used or the
-    # SAMPLE_YEAR table itself can be reorganized.
-    for consecutive_days, precipitation in SAMPLE_YEAR['precipitation']:
-        if 0 <= days and days < consecutive_days:
-            return precipitation
-        days -= consecutive_days
-    raise Exception('No percipitation data for %s' % simulation_day)
+    (precip, et_max) = SAMPLE_YEAR[day]
+    ki = LAND_USE_VALUES[land_use]['ki']
+    return (precip, et_max * ki)
 
 
 def lookup_bmp_infiltration(soil_type, bmp):
     """
     Lookup the amount of infiltration causes by a particular BMP.
     """
-    if bmp not in LAND_USE_VALUES or 'infiltration' not in LAND_USE_VALUES[bmp]:
+    if bmp not in LAND_USE_VALUES \
+       or 'infiltration' not in LAND_USE_VALUES[bmp]:
         raise Exception('%s not a BMP' % bmp)
     elif soil_type not in LAND_USE_VALUES[bmp]['infiltration']:
         raise Exception('BMP %s incompatible with soil %s' % (bmp, soil_type))
@@ -88,7 +50,8 @@ def lookup_cn(soil_type, land_use):
     """
     if land_use not in LAND_USE_VALUES:
         raise Exception('Unknown land use %s' % land_use)
-    elif 'cn' not in LAND_USE_VALUES[land_use] or soil_type not in LAND_USE_VALUES[land_use]['cn']:
+    elif 'cn' not in LAND_USE_VALUES[land_use] \
+         or soil_type not in LAND_USE_VALUES[land_use]['cn']:
         raise Exception('Unknown soil type %s' % soil_type)
     else:
         return LAND_USE_VALUES[land_use]['cn'][soil_type]
@@ -106,3 +69,48 @@ def is_built_type(land_use):
     Test to see if the land use is a "built type".
     """
     return land_use in BUILT_TYPES
+
+
+def precolumbian(land_use):
+    """
+    Project the given land use to a Pre-Columbian one.
+    """
+    if not land_use in PRE_COLUMBIAN_LAND_USES:
+        return 'mixed_forest'
+    else:
+        return land_use
+
+
+def lookup_load(nlcd_class, pollutant):
+    """
+    Get the Event Mean Concentration of `pollutant` for land use
+    class `nlcd_class`
+    """
+    if pollutant not in ['tn', 'tp', 'bod', 'tss']:
+        raise Exception('Unknown pollutant type: %s' % pollutant)
+
+    if nlcd_class not in POLLUTION_LOADS:
+        raise Exception('Unknown NLCD class: %s' % nlcd_class)
+
+    return POLLUTION_LOADS[nlcd_class][pollutant]
+
+
+def lookup_nlcd(land_use):
+    """
+    Get the NLCD number for a particular human-readable land use.
+    """
+    if land_use not in LAND_USE_VALUES:
+        raise Exception('Unknown land use type: %s' % land_use)
+
+    if 'nlcd' not in LAND_USE_VALUES[land_use]:
+        raise Exception('Land use type %s does not have an NLCD class defined',
+                        land_use)
+
+    return LAND_USE_VALUES[land_use]['nlcd']
+
+
+def get_pollutants():
+    """
+    Return the list of pollutants.
+    """
+    return POLLUTANTS
diff --git a/tr55/tables.py b/tr55/tables.py
index f94909a..a3ef230 100644
--- a/tr55/tables.py
+++ b/tr55/tables.py
@@ -1,146 +1,416 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
 """
 TR-55 tables
 """
 
-from datetime import date
-
-# The 'year_start' key points to the date on which this year-long
-# dataset starts.  The 'growing_start' and 'growing_end' keys point to
-# dates giving the respective start and end of growing season.  The
-# values in the array associated with the 'precipitation' key are
-# tuples of a number of consecutive days, and the amount of
-# precipitation during that run of days.
-SAMPLE_YEAR = {
-    'year_start': date(1, 10, 15),    # Sample year starts on 10/15
-    'days_per_year': 365,             # The number of days in the sample year
-
-    'growing_start': date(1, 4, 15),  # Growing season starts on 4/15
-    'growing_end': date(1, 10, 14),   # The last day of growing season is 10/14
-    'growing_etmax': 0.207,           # Max. e/t during growing season
-    'nongrowing_etmax': 0.0,          # ditto non-growing season
+# The values in the array are tuples of inches of precipitation and
+# evapotranspiration.  October 15th is at index 0.
+SAMPLE_YEAR = [
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.0, 0.00),
+    (0.01, 0.0),
+    (0.01, 0.0),
+    (0.01, 0.0),
+    (0.01, 0.0),
+    (0.01, 0.0),
+    (0.01, 0.0),
+    (0.01, 0.0),
+    (0.02, 0.0),
+    (0.02, 0.0),
+    (0.02, 0.0),
+    (0.02, 0.0),
+    (0.03, 0.0),
+    (0.03, 0.0),
+    (0.03, 0.0),
+    (0.04, 0.0),
+    (0.04, 0.0),
+    (0.05, 0.0),
+    (0.05, 0.0),
+    (0.06, 0.0),
+    (0.06, 0.0),
+    (0.07, 0.0),
+    (0.07, 0.0),
+    (0.08, 0.0),
+    (0.09, 0.0),
+    (0.09, 0.0),
+    (0.10, 0.0),
+    (0.11, 0.0),
+    (0.11, 0.0),
+    (0.12, 0.0),
+    (0.14, 0.0),
+    (0.15, 0.0),
+    (0.16, 0.0),
+    (0.17, 0.0),
+    (0.18, 0.0),
+    (0.20, 0.0),
+    (0.21, 0.0),
+    (0.23, 0.0),
+    (0.25, 0.0),
+    (0.26, 0.0),
+    (0.28, 0.0),
+    (0.30, 0.0),
+    (0.33, 0.0),
+    (0.34, 0.0),
+    (0.37, 0.0),
+    (0.38, 0.0),
+    (0.41, 0.0),
+    (0.44, 0.0),
+    (0.48, 0.0),
+    (0.52, 0.0),
+    (0.55, 0.0),
+    (0.59, 0.0),
+    (0.63, 0.0),
+    (0.67, 0.0),
+    (0.72, 0.0),
+    (0.77, 0.0),
+    (0.82, 0.0),
+    (0.89, 0.0),
+    (0.98, 0.0),
+    (1.09, 0.0),
+    (1.22, 0.0),
+    (1.42, 0.0),
+    (1.92, 0.0),
+    (2.71, 0.207),
+    (1.88, 0.207),
+    (1.57, 0.207),
+    (1.33, 0.207),
+    (1.18, 0.207),
+    (1.06, 0.207),
+    (0.95, 0.207),
+    (0.86, 0.207),
+    (0.76, 0.207),
+    (0.69, 0.207),
+    (0.63, 0.207),
+    (0.57, 0.207),
+    (0.54, 0.207),
+    (0.50, 0.207),
+    (0.46, 0.207),
+    (0.43, 0.207),
+    (0.40, 0.207),
+    (0.37, 0.207),
+    (0.34, 0.207),
+    (0.33, 0.207),
+    (0.30, 0.207),
+    (0.28, 0.207),
+    (0.26, 0.207),
+    (0.25, 0.207),
+    (0.23, 0.207),
+    (0.22, 0.207),
+    (0.20, 0.207),
+    (0.18, 0.207),
+    (0.17, 0.207),
+    (0.16, 0.207),
+    (0.15, 0.207),
+    (0.13, 0.207),
+    (0.12, 0.207),
+    (0.11, 0.207),
+    (0.10, 0.207),
+    (0.10, 0.207),
+    (0.09, 0.207),
+    (0.08, 0.207),
+    (0.08, 0.207),
+    (0.07, 0.207),
+    (0.07, 0.207),
+    (0.06, 0.207),
+    (0.05, 0.207),
+    (0.05, 0.207),
+    (0.04, 0.207),
+    (0.04, 0.207),
+    (0.04, 0.207),
+    (0.03, 0.207),
+    (0.03, 0.207),
+    (0.03, 0.207),
+    (0.02, 0.207),
+    (0.02, 0.207),
+    (0.02, 0.207),
+    (0.02, 0.207),
+    (0.01, 0.207),
+    (0.01, 0.207),
+    (0.01, 0.207),
+    (0.01, 0.207),
+    (0.01, 0.207),
+    (0.01, 0.207),
+    (0.01, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207),
+    (0.00, 0.207)
+]
 
-    'precipitation': [
-        (120, 0.00),
-        (7, 0.01),
-        (4, 0.02),
-        (3, 0.03),
-        (2, 0.04),
-        (2, 0.05),
-        (2, 0.06),
-        (2, 0.07),
-        (1, 0.08),
-        (2, 0.09),
-        (1, 0.10),
-        (2, 0.11),
-        (1, 0.12),
-        (1, 0.14),
-        (1, 0.15),
-        (1, 0.16),
-        (1, 0.17),
-        (1, 0.18),
-        (1, 0.20),
-        (1, 0.21),
-        (1, 0.23),
-        (1, 0.25),
-        (1, 0.26),
-        (1, 0.28),
-        (1, 0.30),
-        (1, 0.33),
-        (1, 0.34),
-        (1, 0.37),
-        (1, 0.38),
-        (1, 0.41),
-        (1, 0.44),
-        (1, 0.48),
-        (1, 0.52),
-        (1, 0.55),
-        (1, 0.59),
-        (1, 0.63),
-        (1, 0.67),
-        (1, 0.72),
-        (1, 0.77),
-        (1, 0.82),
-        (1, 0.89),
-        (1, 0.98),
-        (1, 1.09),
-        (1, 1.22),
-        (1, 1.42),
-        (1, 1.92),
-
-        (1, 2.71),
-        (1, 1.88),
-        (1, 1.57),
-        (1, 1.33),
-        (1, 1.18),
-        (1, 1.06),
-        (1, 0.95),
-        (1, 0.86),
-        (1, 0.76),
-        (1, 0.69),
-        (1, 0.63),
-        (1, 0.57),
-        (1, 0.54),
-        (1, 0.50),
-        (1, 0.46),
-        (1, 0.43),
-        (1, 0.40),
-        (1, 0.37),
-        (1, 0.34),
-        (1, 0.33),
-        (1, 0.30),
-        (1, 0.28),
-        (1, 0.26),
-        (1, 0.25),
-        (1, 0.23),
-        (1, 0.22),
-        (1, 0.20),
-        (1, 0.18),
-        (1, 0.17),
-        (1, 0.16),
-        (1, 0.15),
-        (1, 0.13),
-        (1, 0.12),
-        (1, 0.11),
-        (2, 0.10),
-        (1, 0.09),
-        (2, 0.08),
-        (2, 0.07),
-        (1, 0.06),
-        (2, 0.05),
-        (3, 0.04),
-        (3, 0.03),
-        (4, 0.02),
-        (7, 0.01),
-        (122, 0.00)
-    ]
-}
 
+# Land use descriptions may map to the same nlcd code and have the same values
+# but could be handled differently in the model execution (ex, commercial,
+# industrial, transportation). Describes the NLCD class value, the landscape
+# factor (ki) and the Curve Numbers for each hydrologic soil group for that use
+# type.
+# NOTE: Missing NLCD types 12 & 52 (plus all Alaska only types (51, 72-74)
 LAND_USE_VALUES = {
-    'water': {'ki': 0.0, 'cn': {'a': 100, 'b': 100, 'c': 100, 'd': 100}}, 
-    'li_residential': {'ki': 0.42, 'cn': {'a': 51, 'b': 68, 'c': 79, 'd': 84}},
-    'hi_residential': {'ki': 0.18, 'cn': {'a': 77, 'b': 85, 'c': 90, 'd': 92}},
-    'commercial': {'ki': 0.06, 'cn': {'a': 89, 'b': 92, 'c': 94, 'd': 95}},
-    'industrial': {'ki': 0.06, 'cn': {'a': 89, 'b': 92, 'c': 94, 'd': 95}},
-    'transportation': {'ki': 0.06, 'cn': {'a': 89, 'b': 92, 'c': 94, 'd': 95}},
-    'rock': {'ki': 0.0, 'cn': {'a': 77, 'b': 86, 'c': 86, 'd': 91}},
-    'sand': {'ki': 0.0, 'cn': {'a': 77, 'b': 86, 'c': 86, 'd': 91}},
-    'clay': {'ki': 0.0, 'cn': {'a': 77, 'b': 86, 'c': 86, 'd': 91}},
-    'deciduous_forest': {'ki': 0.7, 'cn': {'a': 30, 'b': 55, 'c': 70, 'd': 77}},
-    'evergreen_forest': {'ki': 0.7, 'cn': {'a': 30, 'b': 55, 'c': 70, 'd': 77}},
-    'mixed_forest': {'ki': 0.7, 'cn': {'a': 30, 'b': 55, 'c': 70, 'd': 77}},
-    'grassland': {'ki': 0.6, 'cn': {'a': 30, 'b': 58, 'c': 71, 'd': 78}},
-    'pasture': {'ki': 0.6, 'cn': {'a': 39, 'b': 61, 'c': 74, 'd': 80}},
-    'hay': {'ki': 0.6, 'cn': {'a': 39, 'b': 61, 'c': 74, 'd': 80}},
-    'row_crop': {'ki': 0.9, 'cn': {'a': 67, 'b': 78, 'c': 85, 'd': 89}},
-    'urban_grass': {'ki': 0.7, 'cn': {'a': 68, 'b': 79, 'c': 86, 'd': 89}},
-    'woody_wetland': {'ki': 1, 'cn': {'a': 98, 'b': 98, 'c': 98, 'd': 98}},
-    'herbaceous_wetland': {'ki': 1, 'cn': {'a': 98, 'b': 98, 'c': 98, 'd': 98}},
-    'green_roof': {'ki': 0.4, 'infiltration': {'a': 1.6, 'b': 1.6, 'c': 1.6, 'd': 1.6}},
-    'porous_paving': {'ki': 0.0, 'infiltration': {'a': 7.73, 'b': 4.13, 'c': 1.73}},
-    'rain_garden': {'ki': 0.08, 'infiltration': {'a': 1.2, 'b': 0.6, 'c': 0.2}},
-    'infiltration_trench': {'ki': 0.0, 'infiltration': {'a': 2.4, 'b': 1.8, 'c': 1.4}},
-    'cluster_housing': {'ki': 0.42},
-    'no_till': {'ki': 0.9, 'cn': {'a': 57, 'b': 73, 'c': 82, 'd': 86}}
+    'water':                {'nlcd': 11, 'ki': 0.0, 'cn': {'a': 100, 'b': 100, 'c': 100, 'd': 100}},  # noqa
+    'li_residential':       {'nlcd': 22, 'ki': 0.42, 'cn': {'a': 51, 'b': 68, 'c': 79, 'd': 84}},  # noqa
+    'hi_residential':       {'nlcd': 23, 'ki': 0.18, 'cn': {'a': 77, 'b': 85, 'c': 90, 'd': 92}},  # noqa
+    'commercial':           {'nlcd': 24, 'ki': 0.06, 'cn': {'a': 89, 'b': 92, 'c': 94, 'd': 95}},  # noqa
+    'industrial':           {'nlcd': 24, 'ki': 0.06, 'cn': {'a': 89, 'b': 92, 'c': 94, 'd': 95}},  # noqa
+    'transportation':       {'nlcd': 24, 'ki': 0.06, 'cn': {'a': 89, 'b': 92, 'c': 94, 'd': 95}},  # noqa
+    'rock':                 {'nlcd': 31, 'ki': 0.0, 'cn': {'a': 77, 'b': 86, 'c': 86, 'd': 91}},  # noqa
+    'sand':                 {'nlcd': 31, 'ki': 0.0, 'cn': {'a': 77, 'b': 86, 'c': 86, 'd': 91}},  # noqa
+    'clay':                 {'nlcd': 31, 'ki': 0.0, 'cn': {'a': 77, 'b': 86, 'c': 86, 'd': 91}},  # noqa
+    'deciduous_forest':     {'nlcd': 41, 'ki': 0.7, 'cn': {'a': 30, 'b': 55, 'c': 70, 'd': 77}},  # noqa
+    'evergreen_forest':     {'nlcd': 42, 'ki': 0.7, 'cn': {'a': 30, 'b': 55, 'c': 70, 'd': 77}},  # noqa
+    'mixed_forest':         {'nlcd': 43, 'ki': 0.7, 'cn': {'a': 30, 'b': 55, 'c': 70, 'd': 77}},  # noqa
+    'grassland':            {'nlcd': 71, 'ki': 0.6, 'cn': {'a': 30, 'b': 58, 'c': 71, 'd': 78}},  # noqa
+    'pasture':              {'nlcd': 81, 'ki': 0.6, 'cn': {'a': 39, 'b': 61, 'c': 74, 'd': 80}},  # noqa
+    'hay':                  {'nlcd': 81, 'ki': 0.6, 'cn': {'a': 39, 'b': 61, 'c': 74, 'd': 80}},  # noqa
+    'row_crop':             {'nlcd': 82, 'ki': 0.9, 'cn': {'a': 67, 'b': 78, 'c': 85, 'd': 89}},  # noqa
+    'urban_grass':          {'nlcd': 21, 'ki': 0.7, 'cn': {'a': 68, 'b': 79, 'c': 86, 'd': 89}},  # noqa
+    'woody_wetland':        {'nlcd': 90, 'ki': 1, 'cn': {'a': 98, 'b': 98, 'c': 98, 'd': 98}},  # noqa
+    'herbaceous_wetland':   {'nlcd': 95, 'ki': 1, 'cn': {'a': 98, 'b': 98, 'c': 98, 'd': 98}},  # noqa
+
+    'green_roof':           {'ki': 0.4, 'infiltration': {'a': 1.6, 'b': 1.6, 'c': 1.6, 'd': 1.6}},  # noqa
+    'porous_paving':        {'ki': 0.0, 'infiltration': {'a': 7.73, 'b': 4.13, 'c': 1.73}},  # noqa
+    'rain_garden':          {'ki': 0.08, 'infiltration': {'a': 1.2, 'b': 0.6, 'c': 0.2}},  # noqa
+    'infiltration_trench':  {'ki': 0.0, 'infiltration': {'a': 2.4, 'b': 1.8, 'c': 1.4}},  # noqa
+    'cluster_housing':      {'ki': 0.42},
+    'no_till':              {'ki': 0.9, 'cn': {'a': 57, 'b': 73, 'c': 82, 'd': 86}}  # noqa
 }
 
 # The set of best management practices that we know about.  The
@@ -151,4 +421,49 @@ BMPS = set(['green_roof', 'porous_paving',
 
 # The set of "built" land uses
 BUILT_TYPES = set(['li_residential', 'hi_residential', 'cluster_housing',
-                   'commercial', 'industrial', 'transportation', 'urban_grass'])
+                   'commercial', 'industrial', 'transportation',
+                   'urban_grass'])
+
+PRE_COLUMBIAN_LAND_USES = set([
+    'water',
+    'woody_wetland',
+    'herbaceous_wetland'
+])
+
+# The set of pollutants that we are concerned with.
+POLLUTANTS = set(['tn', 'tp', 'bod', 'tss'])
+
+# Event mean concentrations (mg/l) by pollutant and NLCD type
+# tn: Total Nitrogen, tp: Total Phosphorus,
+# bod: Biochemical Oxygen Demand, tss: Total Suspended Solids
+POLLUTION_LOADS = {
+    11: {'tn': 0,    'tp': 0,     'bod': 0,     'tss': 0},
+    12: {'tn': 0,    'tp': 0,     'bod': 0,     'tss': 0},
+    21: {'tn': 2.8,  'tp': 0.62,  'bod': 61,    'tss': 155.6},
+    22: {'tn': 4.15, 'tp': 0.8,   'bod': 309,   'tss': 147.1},
+    23: {'tn': 6.85, 'tp': 1.15,  'bod': 1180,  'tss': 141.0},
+    24: {'tn': 9.1,  'tp': 1.44,  'bod': 1240,  'tss': 252.9},
+    31: {'tn': 0.1,  'tp': 0.01,  'bod': 1320,  'tss': 10},
+    32: {'tn': 0.1,  'tp': 0.01,  'bod': 30,    'tss': 10},
+    41: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 39},
+    42: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 39},
+    43: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 39},
+    51: {'tn': 0,    'tp': 0,     'bod': 0,     'tss': 0},
+    52: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 39},
+    71: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 47},
+    72: {'tn': 0,    'tp': 0,     'bod': 0,     'tss': 0},
+    73: {'tn': 0,    'tp': 0,     'bod': 0,     'tss': 0},
+    74: {'tn': 0,    'tp': 0,     'bod': 0,     'tss': 0},
+    81: {'tn': 23.0, 'tp': 3.0,   'bod': 1000,  'tss': 500},
+    82: {'tn': 23.0, 'tp': 3.0,   'bod': 1000,  'tss': 1000},
+    90: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    91: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    92: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    93: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    94: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    95: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    96: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    97: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    98: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0},
+    99: {'tn': 0.19, 'tp': 0.006, 'bod': 61,    'tss': 0}
+}
diff --git a/tr55/water_quality.py b/tr55/water_quality.py
new file mode 100644
index 0000000..1abae85
--- /dev/null
+++ b/tr55/water_quality.py
@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
+from tr55.tablelookup import lookup_load, lookup_nlcd
+
+
+def get_volume_of_runoff(runoff, cell_count, cell_resolution):
+    """
+    Calculate the volume of runoff over the entire modeled area
+
+    Args:
+        runoff (number): Q from TR55, averaged amount of runoff over a number
+        of cells.
+
+        cell_count (integer): The number of cells included in the area
+
+        tile_resolution (integer): The size in meters that a cell represents
+
+    Returns:
+        The volume of runoff liters in of the total area of interest
+    """
+
+    # Runoff is in inches, so convert to meters which is the units for the cell
+    # area and compute the meter-cells in the group.  Multiply the resolution
+    # of the cell to get the runoff volume in cubic meters.
+    inch_to_meter = 0.0254
+
+    runoff_m = runoff * inch_to_meter
+    meter_cells = runoff_m * cell_count
+    volume_cubic_meters = meter_cells * cell_resolution  # XXX
+
+    liters = volume_cubic_meters * 1000
+
+    return liters
+
+
+def get_pollutant_load(use_type, pollutant, runoff_liters):
+    """
+    Calculate the pollutant load over a particular land use type given an
+    amount of runoff generated on that area and an event mean concentration
+    of the pollutant.  Returns the pollutant load in lbs.
+    """
+    mg_per_kg = 1000000
+    lbs_per_kg = 2.205
+
+    nlcd = lookup_nlcd(use_type)
+    emc = lookup_load(nlcd, pollutant)
+
+    load_mg_l = emc * runoff_liters
+
+    return (load_mg_l / mg_per_kg) * lbs_per_kg


Test patch:
diff --git a/test/test_model.py b/test/test_model.py
index ac66ae5..201b339 100644
--- a/test/test_model.py
+++ b/test/test_model.py
@@ -1,25 +1,28 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
 """
 Model test set
 """
 
 import unittest
-from datetime import date
 from math import sqrt
-from tr55.tablelookup import lookup_ki, is_built_type
-from tr55.model import runoff_nrcs, simulate_tile, simulate_all_tiles
-from tr55.model import simulate_year
+
+from tr55.tablelookup import lookup_ki
+from tr55.model import runoff_nrcs, simulate_tile, simulate_day
 
 # These data are taken directly from Table 2-1 of the revised (1986)
 # TR-55 report.  The data in the PS array are various precipitation
 # levels, and each respective CNx array is the calculated runoff for
 # that particular curve number with the given level of precipitation
 # corresponding to that in PS.
-PS = [1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0]
-CN55 = [0.000, 0.000, 0.000, 0.000, 0.000, 0.020, 0.080, 0.190, 0.350, 0.530, 0.740, 0.980, 1.520, 2.120, 2.780, 3.490, 4.230, 5.000, 5.790, 6.610, 7.440, 8.290]
-CN70 = [0.000, 0.030, 0.060, 0.110, 0.170, 0.240, 0.460, 0.710, 1.010, 1.330, 1.670, 2.040, 2.810, 3.620, 4.460, 5.330, 6.220, 7.130, 8.050, 8.980, 9.910, 10.85]
-CN80 = [0.080, 0.150, 0.240, 0.340, 0.440, 0.560, 0.890, 1.250, 1.640, 2.040, 2.460, 2.890, 3.780, 4.690, 5.630, 6.570, 7.520, 8.480, 9.450, 10.42, 11.39, 12.37]
-CN90 = [0.320, 0.460, 0.610, 0.760, 0.930, 1.090, 1.530, 1.980, 2.450, 2.920, 3.400, 3.880, 4.850, 5.820, 6.810, 7.790, 8.780, 9.770, 10.76, 11.76, 12.75, 13.74]
-
+PS = [1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0]  # noqa
+CN55 = [0.000, 0.000, 0.000, 0.000, 0.000, 0.020, 0.080, 0.190, 0.350, 0.530, 0.740, 0.980, 1.520, 2.120, 2.780, 3.490, 4.230, 5.000, 5.790, 6.610, 7.440, 8.290]  # noqa
+CN70 = [0.000, 0.030, 0.060, 0.110, 0.170, 0.240, 0.460, 0.710, 1.010, 1.330, 1.670, 2.040, 2.810, 3.620, 4.460, 5.330, 6.220, 7.130, 8.050, 8.980, 9.910, 10.85]  # noqa
+CN80 = [0.080, 0.150, 0.240, 0.340, 0.440, 0.560, 0.890, 1.250, 1.640, 2.040, 2.460, 2.890, 3.780, 4.690, 5.630, 6.570, 7.520, 8.480, 9.450, 10.42, 11.39, 12.37]  # noqa
+CN90 = [0.320, 0.460, 0.610, 0.760, 0.930, 1.090, 1.530, 1.980, 2.450, 2.920, 3.400, 3.880, 4.850, 5.820, 6.810, 7.790, 8.780, 9.770, 10.76, 11.76, 12.75, 13.74]  # noqa
 # INPUT and OUTPUT are data that were emailed to Azavea in a
 # spreadsheet for testing the TR-55 model implementation.
 INPUT = [
@@ -592,23 +595,26 @@ OUTPUT = [
 def simulate(precip, tile_string):
     soil_type, land_use = tile_string.split(':')
     ki = lookup_ki(land_use)
-    return simulate_tile((precip, 0.209 * ki), tile_string)
+    return simulate_tile((precip, 0.207 * ki), tile_string)
+
 
 def average(l):
     return reduce(lambda x, y: x + y, l) / len(l)
 
+
 class TestModel(unittest.TestCase):
     """
-    Model test set
+    Model test set.
     """
     def test_nrcs(self):
         """
         Test the implementation of the runoff equation.
         """
-        # This pair has CN=55 in Table C of the 2010/12/27 memo
+        # This pair has CN=55
         runoffs = [round(runoff_nrcs(precip, 0.0, 'b', 'deciduous_forest'), 2)
                    for precip in PS]
-        self.assertEqual(runoffs[4:], CN55[4:])  # Low curve number and low P cause too-high runoff
+        # Low curve number and low P cause too-high runoff
+        self.assertEqual(runoffs[4:], CN55[4:])
 
         # This pair has CN=70
         runoffs = [round(runoff_nrcs(precip, 0.0, 'c', 'deciduous_forest'), 2)
@@ -625,104 +631,174 @@ class TestModel(unittest.TestCase):
                    for precip in PS]
         self.assertEqual(runoffs, CN90)
 
-    def test_simulate_tile_horizontal(self):
+    def test_simulate_tile_1(self):
         """
-        Test the one-day simulation using sample input/output.  The number
-        0.04 is not very meaningful, this test just attempts to give
-        you some idea about the mean error of the three quantities
-        relative to precipitation.
+        Test the tile simulation using sample input/output.
+
+        The number 0.04 is not very meaningful, this test just
+        attempts to give some idea about the mean error of the three
+        quantities relative to precipitation.
         """
         def similar(incoming, expected):
             precip, tile_string = incoming
             results = simulate(precip, tile_string)
-            me = average(map(lambda x, y: abs(x - y) / precip, results, expected))
+            results = (results['runoff'], results['et'], results['inf'])
+            lam = lambda x, y: abs(x - y) / precip
+            me = average(map(lam, results, expected))
             # Precipitation levels <= 2 inches are known to be
             # problematic.  It is unclear why the 'rock' type is
             # giving trouble on soil types C and D.
-            if precip > 2 and tile_string != 'c:rock' and tile_string != 'd:rock':
+            if precip > 2 and tile_string != 'c:rock' \
+               and tile_string != 'd:rock':
                 self.assertTrue(me < 0.04, tile_string + ' ' + str(me))
         map(similar, INPUT, OUTPUT)
 
-    def test_simulate_tiles_vertical(self):
+    def test_simulate_tile_2(self):
         """
-        Test the RMSE of the runoff levels produced by the one-day
-        simulation against values sample input/output.  The number
-        0.13 is not very meaningful, this test just attempts to show
-        to deviation.
+        Test the RMSE of the runoff levels produced by the tile simulation
+        against values sample input/output.  The number 0.13 is not
+        very meaningful, this test just attempts to put a bound on the
+        deviation.
         """
-        results = [simulate(precip, tile_string)[0] / precip
+        results = [simulate(precip, tile_string)['runoff'] / precip
                    for precip, tile_string in INPUT
-                   if precip > 2 and tile_string != 'c:rock' and tile_string != 'd:rock']
+                   if precip > 2 and tile_string != 'c:rock' and
+                   tile_string != 'd:rock']
         expected = [OUTPUT[i][0] / INPUT[i][0]
                     for i in range(len(INPUT))
-                    if INPUT[i][0] > 2 and INPUT[i][1] != 'c:rock' and INPUT[i][1] != 'd:rock']
-        rmse = sqrt(average(map(lambda x, y: pow((x - y), 2), results, expected)))
+                    if INPUT[i][0] > 2 and INPUT[i][1] != 'c:rock' and
+                    INPUT[i][1] != 'd:rock']
+        lam = lambda x, y: pow((x - y), 2)
+        rmse = sqrt(average(map(lam, results, expected)))
         self.assertTrue(rmse < 0.13)
 
-    def test_simulate_all_tiles(self):
+    def test_simulate_day_invalid(self):
         """
-        Test the tile-by-tile simulation.
+        Test the daily simulation with bad responses.
         """
-        # Test invalid responses
+        non_response1 = {
+            "cell_count": 1
+        }
+
         non_response2 = {
-            "result": {  # No "distribution" key
-                "cell_count": 1
+            "distribution": {}
+        }
+
+        self.assertRaises(Exception, simulate_day, (0, non_response1))
+        self.assertRaises(Exception, simulate_day, (0, non_response2))
+
+    def test_simulate_day_valid(self):
+        """
+        Test the daily simulation with valid responses.
+        """
+        response1 = {
+            "cell_count": 2,
+            "distribution": {
+                "a:pasture": 1,
+                "c:rock": 1
             }
         }
-        non_response3 = {
-            "result": {  # No "cell_count" key
-                "distribution": {}
+
+        response2 = {
+            "cell_count": 20,
+            "distribution": {
+                "a:pasture": 10,
+                "c:rock": 10
             }
         }
-        self.assertRaises(Exception, simulate_all_tiles, (date.today(), non_response2))
-        self.assertRaises(Exception, simulate_all_tiles, (date.today(), non_response3))
 
-        # Test valid responses
+        self.assertEqual(
+            simulate_day(0, response1),
+            simulate_day(0, response2))
+
+    def test_simulate_day_precolumbian(self):
+        """
+        Test the daily simulation in Pre-Columbian times.
+        """
         response1 = {
-            "result": {
-                "cell_count": 2,
-                "distribution": {
-                    "a:pasture": 1,
-                    "c:rock": 1
-                }
+            "cell_count": 1,
+            "distribution": {
+                "d:hi_residential": 1
             }
         }
+
         response2 = {
-            "result": {
-                "cell_count": 20,
-                "distribution": {
-                    "a:pasture": 10,
-                    "c:rock": 10
-                }
+            "cell_count": 10,
+            "distribution": {
+                "d:pasture": 10
             }
         }
-        map(self.assertAlmostEqual,
-            simulate_all_tiles(date.today(), response1),
-            simulate_all_tiles(date.today(), response2))
-
-        # Test Pre-Columbian calculation
-        response3 = {
-            "result": {
-                "cell_count": 1,
-                "distribution": {
-                    "d:hi_residential": 1
-                }
+
+        result1 = simulate_day(182, response1)
+        result2 = simulate_day(182, response2)
+        self.assertNotEqual(
+            result1['distribution']['d:hi_residential'],
+            result2['distribution']['d:pasture'])
+
+        result3 = simulate_day(182, response1, pre_columbian=True)
+        result4 = simulate_day(182, response2, pre_columbian=True)
+        self.assertEqual(
+            result3['distribution']['d:hi_residential'],
+            result4['distribution']['d:pasture'])
+
+    def test_simulate_day_subst_1(self):
+        """
+        Test the daily simulation with BMP substitution.
+        """
+        response1 = {
+            "cell_count": 1,
+            "distribution": {
+                "d:hi_residential": 1
             }
         }
-        response4 = {
-            "result": {
-                "cell_count": 10,
-                "distribution": {
-                    "d:pasture": 10
-                }
+
+        response2 = {
+            "cell_count": 1,
+            "distribution": {
+                "d:no_till": 1
             }
         }
-        map(self.assertNotEqual,
-            simulate_all_tiles(date(1, 4, 15), response3),
-            simulate_all_tiles(date(1, 4, 15), response4))
-        map(self.assertEqual,
-            simulate_all_tiles(date(1, 4, 15), response3, True),
-            simulate_all_tiles(date(1, 4, 15), response4, True))
+
+        result1 = simulate_day(182, response1)
+        result2 = simulate_day(182, response2)
+        self.assertNotEqual(
+            result1['distribution']['d:hi_residential'],
+            result2['distribution']['d:no_till'])
+
+        result1 = simulate_day(182, response1, subst=':no_till')
+        self.assertEqual(
+            result1['distribution']['d:hi_residential'],
+            result2['distribution']['d:no_till'])
+
+    def test_simulate_day_subst_2(self):
+        """
+        Test the daily simulation with reclassification.
+        """
+        response1 = {
+            "cell_count": 1,
+            "distribution": {
+                "d:mixed_forest": 1
+            }
+        }
+
+        response2 = {
+            "cell_count": 1,
+            "distribution": {
+                "a:rock": 1
+            }
+        }
+
+        result1 = simulate_day(182, response1)
+        result2 = simulate_day(182, response2)
+        self.assertNotEqual(
+            result1['distribution']['d:mixed_forest'],
+            result2['distribution']['a:rock'])
+
+        result1 = simulate_day(182, response1, subst='a:rock')
+        self.assertEqual(
+            result1['distribution']['d:mixed_forest'],
+            result2['distribution']['a:rock'])
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/test_operations.py b/test/test_operations.py
new file mode 100644
index 0000000..5c5a77e
--- /dev/null
+++ b/test/test_operations.py
@@ -0,0 +1,38 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
+"""
+Operation tests.
+"""
+
+import unittest
+
+from tr55.operations import dict_plus, dict_minus
+
+
+class TestOperations(unittest.TestCase):
+    """
+    Dictionary operation test set.
+    """
+    def test_plus(self):
+        """
+        Test dictionary addition.
+        """
+        a = {'x': {'y': {'z': {'a': 1, 'b': 3, 'c': 13}}}}
+        b = {'x': {'y': {'z': {'a': 1, 'b': 5, 'c': 21}}}}
+        c = {'x': {'y': {'z': {'a': 2, 'b': 8, 'c': 34}}}}
+        self.assertEqual(dict_plus(a, b), c)
+
+    def test_minus(self):
+        """
+        Test dictionary subtraction.
+        """
+        a = {'x': {'y': {'z': {'a': 34, 'b': 144, 'c': 610}}}}
+        b = {'x': {'y': {'z': {'a': 21, 'b': 89, 'c': 377}}}}
+        c = {'x': {'y': {'z': {'a': 13, 'b': 55, 'c': 233}}}}
+        self.assertEqual(dict_minus(a, b), c)
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_tablelookup.py b/test/test_tablelookup.py
index f1b1e83..b9398b6 100644
--- a/test/test_tablelookup.py
+++ b/test/test_tablelookup.py
@@ -1,36 +1,32 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
 """
-Table Lookup test set
+Table Lookup test set.
 """
 
 import unittest
-from datetime import date
-from tr55.tables import LAND_USE_VALUES
-from tr55.tablelookup import lookup_et, lookup_p, lookup_bmp_infiltration, lookup_cn
+
+from tr55.tablelookup import lookup_pet, lookup_bmp_infiltration, lookup_cn
 
 
 class TestTablelookups(unittest.TestCase):
     """
     Table Lookup test set
     """
-    def test_lookup_p(self):
+    def test_lookup_pet(self):
         """
         Do some spot-checks on the SampleYear data.
         """
-        self.assertEqual(lookup_p(date(1, 10, 15)), 0.0)
-        self.assertEqual(lookup_p(date(1, 2, 12)), 0.01)
-        self.assertEqual(lookup_p(date(1, 2, 15)), 0.01)
-        self.assertEqual(lookup_p(date(1, 2, 19)), 0.02)
-        self.assertEqual(lookup_p(date(1, 10, 14)), 0.0)
-
-    def test_lookup_et(self):
-        """
-        Do some spot-checks on the data from Table A.
-        """
-        self.assertEqual(lookup_et(date(1, 4, 15), 'woody_wetland'), 0.207)
-        self.assertEqual(lookup_et(date(1, 10, 14), 'woody_wetland'), 0.207)
-        self.assertTrue(lookup_et(date(1, 6, 15), 'commercial') > 0.0)
-        self.assertEqual(lookup_et(date(1, 4, 14), 'woody_wetland'), 0.0)
-        self.assertEqual(lookup_et(date(1, 10, 15), 'woody_wetland'), 0.0)
+        self.assertEqual(lookup_pet(0, 'water')[0], 0.0)
+        self.assertEqual(lookup_pet(52, 'hi_residential')[1], 0.0)
+        self.assertEqual(lookup_pet(104, 'green_roof')[1], 0.0)
+        self.assertEqual(lookup_pet(156, 'cluster_housing')[0], 0.23)
+        self.assertEqual(lookup_pet(208, 'grassland')[0], 0.2)
+        self.assertEqual(lookup_pet(260, 'woody_wetland')[1], 0.207)
+        self.assertEqual(lookup_pet(352, 'hay')[1], 0.207 * 0.6)
 
     def test_lookup_bmp_infiltration(self):
         """
@@ -39,7 +35,7 @@ class TestTablelookups(unittest.TestCase):
         self.assertEqual(lookup_bmp_infiltration('d', 'green_roof'), 1.6)
         self.assertEqual(lookup_bmp_infiltration('c', 'porous_paving'), 1.73)
         self.assertEqual(lookup_bmp_infiltration('b', 'rain_garden'), 0.6)
-        self.assertEqual(lookup_bmp_infiltration('a', 'infiltration_trench'), 2.4)
+        self.assertEqual(lookup_bmp_infiltration('a', 'infiltration_trench'), 2.4)  # noqa
 
     def test_lookup_cn(self):
         """
diff --git a/test/test_water_quality.py b/test/test_water_quality.py
new file mode 100644
index 0000000..f3cf6b7
--- /dev/null
+++ b/test/test_water_quality.py
@@ -0,0 +1,44 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+
+import unittest
+from tr55.water_quality import get_volume_of_runoff, get_pollutant_load
+from tr55.tables import POLLUTION_LOADS
+
+
+class TestWaterQuality(unittest.TestCase):
+
+    def test_volume(self):
+        """
+        Test the water volume computation.
+        """
+        cell_res = 30  # meters
+        cell_count = 100
+        runoff = 0.4  # inches
+        liters = get_volume_of_runoff(runoff, cell_count, cell_res)
+
+        self.assertEquals(30480, liters,
+                          "Did not calculate the correct runoff volume")
+
+    def test_load(self):
+        """
+        Test the pollutant load computation.
+        """
+        nlcd = 24
+        pollutant = 'tn'
+        emc = POLLUTION_LOADS[nlcd][pollutant]
+        runoff_liters = 1000
+
+        expected = ((runoff_liters * emc) / 1000000) * 2.205
+        load = get_pollutant_load('commercial', pollutant, runoff_liters)
+
+        self.assertEquals(expected, load)
+
+    def test_load_bad_nlcd(self):
+        self.assertRaises(Exception, get_pollutant_load, 'asdf', 'tn', 1000)
+
+    def test_load_bad_pollutant(self):
+        self.assertRaises(Exception, get_pollutant_load, 'commercial',
+                          'asdf', 1000)