def simulate_water_quality(tree, cell_res, fn,
                           pct=1.0, current_cell=None, precolumbian=False):
    """
    Perform a water quality simulation by doing simulations on each of
    the cell types (leaves), then adding them together by summing the
    values of a node's subtrees and storing them at that node.

    `tree` is the (sub)tree of cell distributions that is currently
    under consideration.

    `pct` is the percentage of calculated water volume to retain.

    `cell_res` is the size of each cell/pixel in meters squared
     (used for turning inches of water into volumes of water).

    `fn` is a function that takes a cell type and a number of cells
    and returns a dictionary containing runoff, et, and inf as
    volumes.

    `current_cell` is the cell type for the present node.
    """
    # Internal node.
    if 'cell_count' in tree and 'distribution' in tree:
        n = tree['cell_count']

        # simulate subtrees
        if n != 0:
            tally = {}
            for cell, subtree in tree['distribution'].items():
                simulate_water_quality(subtree, cell_res, fn,
                                       pct, cell, precolumbian)
                subtree_ex_dist = subtree.copy()
                subtree_ex_dist.pop('distribution', None)
                tally = dict_plus(tally, subtree_ex_dist)
            tree.update(tally)  # update this node

        # effectively a leaf
        elif n == 0:
            for pol in get_pollutants():
                tree[pol] = 0.0

    # Leaf node.
    elif 'cell_count' in tree and 'distribution' not in tree:
        # the number of cells covered by this leaf
        n = tree['cell_count']

        # canonicalize the current_cell string
        split = current_cell.split(':')
        if (len(split) == 2):
            split.append('')
        if precolumbian:
            split[1] = make_precolumbian(split[1])
        current_cell = '%s:%s:%s' % tuple(split)

        # run the runoff model on this leaf
        result = fn(current_cell, n)  # runoff, et, inf
        runoff_adjustment = result['runoff-vol'] - (result['runoff-vol'] * pct)
        result['runoff-vol'] -= runoff_adjustment
        result['inf-vol'] += runoff_adjustment
        tree.update(result)

        # perform water quality calculation
        if n != 0:
            soil_type, land_use, bmp = split
            runoff_per_cell = result['runoff-vol'] / n
            liters = get_volume_of_runoff(runoff_per_cell, n, cell_res)
            for pol in get_pollutants():
                tree[pol] = get_pollutant_load(land_use, pol, liters)