SURFACE WATER FLOW IN LANDSCAPE MODELS:
1. EVERGLADES CASE STUDY

Abstract

Many landscape models require extensive computational effort using a large array of grid cells that represent the landscape. The number of spatial cells may be in the thousands and millions, while the ecological component run in each of the cells to account for landscape dynamics is often process based and fairly complex. To compensate for the increased computational complexity of the model there is a tendency to simplify the hydrologic component that fluxes material horizontally across the landscape. Instead of full scale hydrologic models based on stable implicit schemes, computationally simpler explicit algorithms are incorporated and run with quite large time steps. As a result some fairly inadequate behavior may be observed, especially when the temporal and spatial steps are modified without due care. We illustrate these problems with a series of runs performed using the Everglades Landscape Model (Southern Florida, USA), that covers an area of more than 10,000 km2. Several algorithms for hydrologic fluxing are compared in terms of their computational complexity and stability. We argue that a compromise can be drawn by supplementing the explicit modeling scheme with a series of additional checks and conditions that provide for model stability, and with some empirical assumptions that allow the model to operate over a sufficiently large range of temporal and spatial scales.

Introduction

In landscape models hydrology plays an important driving role. Nevertheless there are no off-the-shelf universal models that can be easily adapted for a wide range of applications. Significant effort is needed to tune existing models to the specifics of the landscape and the goals of the study. This is because hydrologic modeling as a part of larger landscape modeling imparts certain conditions on the methods used. Being part of a more complicated modeling structure, the hydrologic module is required to be simple enough to run within the framework of the integrated physico-ecological model. As a result some hydrologic details need be sacrificed to make the whole task more feasible, and these details may differ from one application to another, depending upon the sizes of the study area, the physical characteristics of the slope and surface, and the goals and priorities of the modeling effort.

Another trade-off specific to landscape scale models is the coarser spatial and temporal resolution that they usually employ in contrast to the classic hydrologic methods that were primarily developed for small scale well sampled, replicated and controlled systems. While most of the methods and equations in classic hydrology have been developed for sizes on the order of meters and times on the order of hours and less, this is hardly the resolution that landscape models can afford. As a result discrete approximations of the essentially continuous hydrologic processes become a source of potential problems. In place of a continuous movement of water and constituents over the area we need to deal with essentially discrete motions, when large volumes of material are moved over large distances on relatively rare occasions.

Among the simplified approaches to surface water fluxing, the prevailing ones are based on the kinematic wave approximation of the St. Venant's equations [Beven and Wood, 1993]. These are solved numerically, with the implicit scheme bearing all the advantages of being stable for any chosen time and space intervals [Greco and Panattoni, 1977]. This approach is most common in modeling linear flood routing in rivers and canals. However in models of 2-dimensional overland flow the implicit scheme is usually substituted by a computationally simpler explicit one. For example, this approach was adopted in the SHE model [Abbott, Bathurst et al., 1986b], where channel flow was modeled by an implicit scheme, but the overland flow scheme used an explicit procedure. Because the solution of the implicit method is essentially based on the boundary conditions, the advantages of the explicit approach become especially obvious when we need to consider a region that has complex boundaries or if it is further subdivided into subregions by some impermeable borders (levees) or if it has other kinds of disturbances like a network of canals, that is not attached to the underlying elevation data.

On the other hand, with the obvious advantage of computational simplicity, when running the explicit scheme we should realize that there is always a risk of destabilizing the model, especially if we modify its temporal or spatial resolution. This was the situation we encountered in modeling hydrologic fluxes within the framework of the Everglades Landscape Model (ELM). The problem was further complicated by the relatively large size of the area to be considered. While ANSWERS [Beasley and Huggins, 1980] has been recommended for watershed sizes less than 100 km², and TOPMODEL was tested for catchments of up to 10 km² [Beven, 1984], in our case we were looking at an area orders of magnitude larger than that, roughly 10,000 km². Moreover the hydrologic component was to be embedded in a fairly detailed process-based simulation model of terrestrial and aquatic biota. This necessarily implied certain restrictions and conditions on the methods to be used.

Hydrologic Module in the Everglades Landscape Model

The Everglades Landscape Model (ELM) is designed to simulate the landscape scale vegetation response to simulated hydrology, water quality and fire. The ELM has a fine scale grid that divides the landscape into 10,000 1 km ² cells. A fairly complex General Ecological Model [Fitz et al., 1995] is run in each of the cells. This imposed certain restrictions on the time step to be used in the model. Our intent was to keep the time step at the order of magnitude of 1 day.

Overland sheet flow and groundwater movements are some of the principal fluxes of water across the landscape and they regulate the plant and animal community structure. A large network of canals, levees, and associated flow control structures are used for water management in the region. The canals are a significant water transport mechanism in the Everglades, moving large quantities of water further and more rapidly than the comparatively slow overland flow. The interaction of canals with overland sheet flow is described elsewhere.

Within each raster cell the unit model simulates a variety of hydrologic processes and parameters [Fitz et al., 1995], including the following:

1. Transpiration associated with plant growth, physiology and relative humidity.
2. Evaporation using pan evaporation estimates and pan coefficients.
3. Rainfall based on precipitation data interpolated over 9 stations.
4. Seepage of water from that stored above the sediment/soil surface into that stored in sediment pore space (either in unsaturated or saturated storage).
5. Roughness coefficient that depends on dynamic simulation of plant biomass, numeric density, and plant morphology.

After the water head in each raster cell is modified due to vertical fluxes, the surface water movement between the raster cells (and associated transport of constituents) is calculated. The simplest way to do this is to look at the flow between two adjacent cells iand i+1 as flow in an open channel and use the so called slope-area method [Boyer, 1964], which is a kinematic wave approximation of St. Venant's momentum equation. The flux (m³/d) in this case is described by the empirical Manning's equation for overland flow:

F = A R^2/3 G^1/2/ M,

where A is the cross-sectional area of water flux, R is the hydraulic radius and G is the slope of the energy gradient.

Geometry of a 1-dimensional flow. Horizontal fluxing is driven by the head difference.
The water head is determined as a sum of elevation E_i and the surface water D_i.

In discrete terms used in the model, assuming that the cells are squares of the same area S, and that the energy gradient is determined by the head difference over the grid size, the equation can be rewritten as:

F = sgn(H_i-H_i+1) |H_i-H_i+1|^1/2 D^5/3 S^1/4/ M, (1)

where F is the flux of water (m³/d) between cells, H_i and H_i+1 are the hydraulic heads (m) of the cells, which result from adding the surface water depth (D_i) to the cell elevation (E_i); Sis the area (m²) of the surface water in the cell (assuming that both cells have the same area); D = (D_i + D_i+1) / 2 is the average depth of the surface water (m) in the two interacting cells, and M = (M_i + M_i+1) / 2 is the average Manning's coefficient of surface roughness. The algorithm based on this flow calculation we will call "proper Manning's algorithm" and refer to it as Method N1.

The Manning's roughness coefficient in a cell is a function of the sediment type and the interaction of the vegetation height/density and water depth:

M = M_max - |(M_max - M_min) (2^1-D/Hm- 1)| (2)

The M = M_min roughness is the minimum Manning's coefficient for a vegetation-free cell, the M = M_max is the maximum roughness associated with the (dynamic) vegetation density in the cell, and Hm is the (dynamic) height of the macrophytes in the cell. This function returns a roughness coefficient whose value ranges from a vegetation-free minimum to a maximum at the point of full plant immersion [Petryk, Asce et al., 1975]. As water depth increases over that of the macrophyte height, the roughness decreases to an asymptote at the baseline sediment roughness [Nalluri and Judy, 1989].

St. Venant's continuity equation:
A/t = -F/x

in its discrete form provides for a means to link a cell to the neighboring ones:

D_i(t+dt) = D_i(t) + (F_i-1(t) - F_i(t)) / dt / S (3)

This equation works well as long as the outflux from the cell does not exceed the sum of the influxes and the available water volume in a cell. Moreover it would be quite unrealistic if because of the fluxing, the water head in the recipient cell exceeds the water head in the donor cell. These events obviously will become more likely as the time step dt grows. We may illustrate this on a series of test runs on a cross-section of cells. We assume an initially homogeneous distribution of water over the cells except for a 1 m higher head in the first cell. With dt = 0.005 the resulting distribution is fairly smooth, however already with this small time step we may observe slight oscillations, where the water levels out in the cells.

A.

B.

C.
Distribution of water over the cells in time with different time steps in the model. A.- t = 0.005; B. t = 0.01; C. t = 0.05 day.
As the time step is increased to dt = 0.01 and further to dt = 0.05 the oscillations increase and their amplitude tends to increase with time, eventually destabilizing the model.

"Flip-flop" checkerboard patterns observed in the model when time step is large. The checkerboard pattern appears first in the areas with the lowest water heads, the oscillations do not tend to dampen out, they tend to grow eventually destabilizing the model.

In the 2-dimensional case this instability results in, what we call, a "checker-board" pattern, the situation described by Chow et al. [Chow, 1988] as an accumulation or piling up of water, when the time step is larger than the time needed for a wave to travel the distance of the cell width. This is the Courant condition which is necessary (but not sufficient) for stability of explicit schemes:

dx / dt >= F/S.

This condition implies that the most natural way to make the model stable would be to reduce the time step, when we need to run the model at a finer spatial resolution. However this is in conflict with our initial desire to run the model at larger time steps to decrease the amount of CPU time needed, even at the expense of some precision.

To mitigate this problem another condition is added to the model, which effectively slows down the flux rate, preventing further water fluxing once the heads in two adjacent cells are balanced. Equation (1) is replaced by:

F = sgn(H_i-H_i+1) min[ S |H_i-H_i+1|/2, |H_i-H_i+1|^1/2 D^5/3 S^1/4 dt / M]/dt,(4)

In this way we insure that, whenever a flux between cells is calculated, it is not going to move more water than there is to equilibrate the water levels in the two cells. We call this the "cell-equilibrating algorithm" and refer to it as Method N2 in what follows.

With this restriction we can increase the time step quite significantly while still keeping the model stable.

Test results with the cell equilibrating algorithm, dt = 0.5 day. Model stability has very much increased, but the propagation rate went down.
Here the propagation is calculated with a time step of dt = 0.5, which is 10 times greater than the maximal time step that was used above, yet there are no observed instabilities. However the water propagation slows down dramatically . Actually the leveling out of water across the cells becomes asymptotic, very slowly tending to the equilibrium, that was so quickly achieved above. When the water head differences are not so large (less than 1m), the two methods prove to be fairly close. The cell equilibrating method seems to be especially efficient when the head differences are small. In our case it rapidly accommodates the pulses, while Method N1 results in significant instabilities unless the time step is made extremely small (~0.001 day).

Once the interaction between two cells is determined, we need to figure out how all the cells can be linked spatially over the landscape. This immediately causes some controversy, because from the 1-dimensional channel flow assumption we need to expand to the 2-dimensional water flow and there may be a variety of options to generalize the pairwise interaction described above over the multitude of cells in the landscape. In ELM we first applied the simplest pattern, when each cell was interacting with its four immediate neighbors, assuming that the flows in the North-South and East-West directions are independent and simultaneous. This approach was applied to the whole model area, which was scanned twice: first there was a loop to calculate all the potential fluxes VF_i,j in the NS and HF_i,j in the EW directions, based on the existing water heads, surface water depths and vegetation in the cells; next during the second loop the water depths D_ij were recalculated applying equation (3) modified for four fluxes. This approach, though requiring extra memory to store the two intermediate flow arrays, has the advantage of being independent of the order in which the cells are scanned.

Fluxing water simultaneously with several cells creates additional potential for destabilizing the scheme. In fact, if each flow is calculated independently, with no account of the head change due to the flux that has been already computed, one can easily drain the cell below zero. This is especially likely to occur if we keep the time steps sufficiently large. To avoid this a series of additional checks are implemented when calculating the fluxes for each cell. First the four potential fluxes HF_i-1,j, HF_ij, VF_i,j-1, and VF_ij are calculated. The variable OUT is calculated as a sum of all potential flows out of the cell and the variable IN is the sum of all potential flows into the cell. Next the potential new head D = D_ij + (IN - OUT) dt/S is compared to the minimal head in the interacting adjacent four cells. If it is lower:
D - min (D_i-1,j; D_i+1,j; D_i,j-1; D_i,j+1) < 0 (5)
it means that at some point, if the potential fluxes were to be realized, water would have to move against the gradient, further draining the donor cell. In order to avoid this unrealistic and destabilizing situation, all the fluxes out of the cell, that make the variable OUT above, are proportionally decreased so that condition (5) would no longer hold.

After some adjustments in the rate parameters the method produced fairly reasonable results even at sufficiently large time steps (dt = 0.5 day).

The CALM case-study

It is still quite a large job to calibrate the whole landscape model running it over the almost 10,000 cells of the full ELM model. A one year run takes almost 2.5 hours on a Sparc 10 Sun workstation. Therefore it was decided to cut out a relatively small portion of the ELM area to perform a more vigorous test of the model. Water Conservation Area 2A was chosen because of the additional data available for that region, primarily on the elevations and habitat types [Jensen, Rutchey et al., 1995]. Since the habitat data had a better resolution than the 1 km² used in ELM, a finer grid of 0.25 km² was developed for this subregion. This also served the goal of testing some of the scaling issues. This truncated ELM model was called CALM - Conservation Area Landscape Model.

The initial plan was to run the model unchanged just switching to the other set of data files. Running the straight Manning's method would require at least a four times smaller time step to compensate for the finer spatial resolution in CALM. But this would take away almost all the gain we got by diminishing the model area. Therefore we were primarily looking at the equilibrating method.

The problem now is that because equation (4) effectively slows down the water fluxing, we are not getting the water redistributed fast enough. This turns out to be especially crucial in our case, where due to the artificially controlled flow through the system of canals and structures, on certain occasions the landscape is receiving tremendous amounts of water, that need to be quickly propagated throughout the model area. In the case of CALM the situation is further aggravated because over the same number of iterations the water is equilibrated over smaller distances, since the cell size is reduced. Pumping station 7 is bringing millions of cubic meters of water daily directly into the area in addition to similar quantities received from Water Conservation Area 1 through a series of structures. Unless the water is moved spatially quickly enough very significant water heads build up. The equilibrium condition (4) works well for purposes of stabilizing the model, but it prevents fast propagation of water.

As an alternative to the two step surface water calculation applied above, we now allow water to be fluxed immediately as the flux is calculated according to (4). Instead of first calculating all the potential fluxes between the cells, storing the results and only then, in a separate loop over all the cells, actually recalculating the water heads, in this new scheme as soon as the flux F is defined, the water heads in the corresponding cells are recalculated:

D_i(t+dt) = D_i(t) - F_i dt / S
D_i+1(t+dt) = D_i+1(t) + F_i dt / S

To do that we also have to consider the NS and EW fluxing separately. That is, first all the landscape cells are updated for the fluxes in the East-West direction, and then over a separate loop all the North-South fluxes are calculated and water heads are updated once again.

In this way, if water rises in ith cell, then the water head in the (i+1)st cell is updated and also rises prior to calculating the flux F_i+1 to the cell (i+2), therefore increasing this flux, and so on. This allows for faster water transfer, but also makes the process dependent upon the order of cell sequencing. In fact, if we were moving in the opposite direction through the cells, we would be getting fairly different results.

Difference generated by the two sequencing algorithms. When running through the cells from left to right, and equilibrating water levels in them, we end up with a distribution different from the one we get, if cells are treated in reverse order.

To mitigate this difference we use a reciprocal algorithm. Over each time step the loop runs first from East to West, then backwards starting from the Westernmost cell and going back to the East. Similarly first the NS loop was run from North to South and then back from South to North. This doubles the number of computations needed, but makes the procedure symmetric with respect to the order that the cells are considered.

In this way we get a better control over the fluxes in the landscape and can distribute water somewhat more rapidly over the area. Another advantage of this method is that condition (5) is no longer needed and therefore there is no longer any potential for mass disbalance in the system. However, because of the condition (4) still in place, occasionally when pulses of pumping are especially high, we may develop unrealistically high water heads. The difference with the proper Manning's method (Method N1) run at dt = 0.001 is still quite substantial, especially when the head difference is large. The equilibrating algorithm (Method N2) is quite efficient in dealing with the small variation of water heads, when there are no significant pulses and associated rises in water heads. In this case the results turn out to be quite close to the base line proper Manning's simulations, and the gain in model efficiency is quite substantial because similar results can be generated with 500 times larger time steps.

Comparison of the water distribution generated by the straight Manning's algorithm with dt = 0.001 and the equilibrating algorithm with dt = 0.5. By day 10 there is still a significant difference in the way the two algorithms accommodate a 1 meter pulse of water. This difference is significantly smaller when the variations of water are small (< 0.5 m).

The method was further modified to accommodate the large pulses, that we had to deal with in CALM and ELM. Once the Manning flux in (5) is so large that it tends to disequilibrate the two adjacent cells, instead of curtailing the flux and limiting it to the amount that would equilibrate the two adjacent cells as in (4), it is assumed that water can be fluxed over more than one cell in one time step. In this way, over one time step, the high water head can be redistributed across as many cells in the landscape as needed to accommodate the excess of water, as long as the heads in those cells are initially lower than the equilibrium thus attained. The propagation of water stops once a cell with a higher water head, or the boundary is reached. This assumption seems to be quite appropriate since over a large time step water in fact will probably travel further than over one adjacent cell.

To maintain some control over the process an additional model parameter is introduced. This parameter is to be specified as the allowed head difference between any two cells. It is also assumed that if the head difference between two cells is initially less than , then this condition is not applied and water is fluxed across these cells, while it is in excess.

Thus, for the ith donor cell we first identify the number of cells N to interact with as recipients. Let T_N = _j=0,N H_i+j and F_i be the Manning's flux calculated according to (1). Then N is incremented while H_i+N <= T_N-1/N - N/2
(to make sure that the water head in the next cell in the row/column was less than the resulting equilibrium, taking into account the allowed head gradient ) and
H_i - T_N/(N+1) - (N+1)/2 <= min (F_i, (H_i - H_i+1))
(to ensure that there is still water to flux according to the calculated Manning's flow and that this flow is not larger than the head difference between the ith cell and its next neighbor).

For these N cells the water heads are then recalculated:

H_i+j = T_N/N + (N-2j-1)/2,
j=0,...,N.

We use essentially as a calibration parameter to achieve distribution patterns close to those generated by the proper Manning's method run with very small time steps on a trial basis. The smaller the allowed head difference , the smoother the propagation of water, the more water is fluxed in one iteration after a pulse has occurred and the slower the water equilibrates afterwards. With = 0 the equilibrium may be attained over one step. In contrast, with sufficiently large the attained distributions are almost step-wise, but eventually the equilibrium is reached faster, because more water is distributed over the cells that are not limited by the condition.

As a result even under conditions of extremely high water deliveries into the system, the model quickly converges to the same water head values as we were to get in the straight Manning's method with very small time steps.

A.

B.
Comparison of the straight Manning's algorithm (A) with the distributing one (B) with dt =0.5. Even though the transfer process is different, by day 4 both methods equilibrate to the same distribution pattern.
The transfer stages do look somewhat different, but in a matter of several days we are getting almost the same distribution. We have tested this method for several scenarios of pumping and even for time steps as large as several days it produced fairly stable results, which were within a 10% error compared with the available observed data. It is impossible to present these results within the format of this publication, since the output is generated as color animations. We therefore refer the reader to our WWW page, were the model output can be viewed.

Conclusion

The methods suggested can be considered as an empirical approach to surface water routing. It is certainly very much based on some empirical assumptions and common sense. However we should remember that the Manning's equation is also empirical and in fact in many applications (1) is generalized to

F = D (6)

where and are offered as empirical constants [Novotny and Olem, 1994], [Beven and Wood, 1993]. Actually this means that in surface hydrology, especially in models of this scale and resolution, there is still quite a lot of empiricism to allow the kind of manipulations that we have applied, especially if the results turn out to be adequate.

The major problem with the proper Manning's equation is that = 1/2. With such and small values of D (D <<1) the change in the flux rate, calculated in (1), becomes relatively much larger than the variations in D. As a result, as we saw above, even with very small time steps there are certain instabilities in the region of small water head differences. In those areas where the really significant water fluxes occur, we see no instabilities even with larger time steps. But it is the area of quasi-equilibrium that causes most of the trouble and requires smaller time steps. It is interesting to note that with = 1 or 1.5 we get a much smoother distribution at significantly larger time steps, which actually very much resembles the picture observed in the cell-equilibrating algorithm (Method N2). Varying we can choose between faster distribution over the landscape and more stability at larger time steps.

Distribution of water in the generalized Manning's method with = 1.5 and dt = 0.01. Even with the larger time step the pattern is smooth (compare), but the rate of propagation has decreased.

The = 1/2 easily justified by Newtonian mechanics, seems to contradict with the kinematic wave approximation for the St. Venant's equation. Neglecting the gravity wave terms in this equation, makes it nearly impossible to reach equilibrium, eventually resulting in the "flip-flop" dynamics, that can continue indefinitely with no trend towards stabilization. On the contrary it has a possible destabilizing effect.

It has been observed that Manning's equation can be easily misused. Kadlec [1990] argues that modeling sheet flow in wetlands is one such case , where the flow is not turbulent enough to use the Manning's equation properly, though it is not quite laminar as well. He suggests to use = 1 in those cases when we find the flow to be more laminar than turbulent. Interestingly, we have arrived at quite similar conclusions based on purely theoretical observations about stability of various modeling methods.

The requirements of a landscape modeling framework, where hydrology is only a part of a much more complex and sophisticated model structure, do not allow the time step to be reduced in order to accommodate the instabilities in the Manning's equation. For that same reason, plus the complicated spatial pattern with multiple internal borders and canals, we cannot employ the stable implicit scheme. The methods suggested certainly sacrifice some of the precision, especially in the transfer processes, but they represent the quasi-equilibrium state well and substantially gain in model efficiency in terms of the CPU time required. In this case we have to rely more on the comparison of the model output with the data available, and be ready to switch from the more process based description to a more empirical one as in (6). As can be seen on the Web page, the results that we have been obtaining for the Everglades seem to match fairly well with the experimental data available on water heads. This speaks in favor of the methods applied.

Certainly by choosing this approach, by diverging from the process-based approach and by allowing more parameter and formulae calibration, we decrease the generality of the model, thus requiring additional testing and calibration when switching to other scales or areas. The question is what is a truly process-based model as against an empirical, regression one. In any process-based model there is a certain level of abstraction at which we actually utilize certain empirical generalizations, rather than true process description. For example, there is hardly an adequate detailed description of the photosynthesis process to be found among the models of vegetation growth, instead some variations of Michaelis-Menten kinetics are applied, which are already empirical generalizations of the process. Nevertheless these models claim to be process-based. As we go to larger systems, such as landscapes, we will need to employ even more generalized formalizations, as was done above for modeling hydrology in the Everglades.

References

• Abbott, M. B., J. C. Bathurst, et al. (1986b). "An introduction to the European Hydrological System -- Systeme Hydrologique Europeen, "SHE" 2: Structure of a physically-based, distributed modelling system." Journal of Hydrology 87: 61-77.
• Beasley, D. B. and L. F. Huggins (1980). ANSWERS (Areal Nonpoint Source Watershed Environment Response SImulation), User's Manual. West Lafayette, Purdue University.
• Beven, K. and E. F. Wood (1993). Flow routing and the hydrological response of channel networks. In: K. Beven, Kirkby, M.J. (Eds.) Chichester,Channel Network Hydrology. John Wiley & Sons.
• Beven, K. J., Kirkby, M. J., Schofield, N. and Tagg, A. F. (1984). " Testing a physically-based flood forecasting model (TOPMODEL) for three U.K. catchments." J. Hyrology 69: 119-143.
• Beven, K. J. and M. J. Kirkby (1979). "A physically-based, variable contributing area model of basin hydrology." Hydrological Sciences Bulletin 24(1): 43-69.
• Boyer, M. C. (1964). Streamflow Measurements. Handbook of Applied Hydrology. V. T. Chow. New York, McGraw-Hill Book Company: 15.1 - 15.42.
• Chow, V. T., Maidment, D. R. and Mays, L. W. (1988). Applied hydrology. New York, McGraw-Hill.
• Fitz, H. C., E. DeBellevue, et al. (1995). "Development of a general ecosystem model for a range of scales and ecosystems." Ecological Modelling: accepted w/revisions.
• Grayson, R. B., I. D. Moore, et al. (1992a). "Physically based hydrologic modeling: 1. A Terrain-based model for investigative purposes." Water Resources Research 28(10): 2639-2658.
• Greco, F. and L. Panattoni (1977). Numerical solution methods of the St. Venant equations. In: T. A. Ciriani, Maione, U., Wallis, J. R. Mathematical Models for Surface Water Hydrology. London, New York, Sydney, Toronto: 181-194.
• Jensen, J. R., K. Rutchey, et al. (1995). "Inland wetland change detection in the Everglades Water Conservation Area 2A using a time series of normalized remotely sensed data." Photogrammetric Engineering & Remote Sensing 61(2): 199-209.
• Kadlec, R. (1990). "Overland Flow in Wetlands: Vegetation Resistance." Journal of Hydraulic Engineering 116(5): 691-706.
• Nalluri, C. and N. D. Judy (1989). Factors affecting roughness in coefficient in vegetated channels. Channel flow and catchment runoff: centennial of Manning's formula and kuichling's rational formula, University of Virginia.
• Novotny, V. and H. Olem (1994). Water Quality. Prevention, Identification, and Management of Diffuse Pollution. New York, Van Nostrand Reinhold.
• Petryk, S., A. M. Asce, et al. (1975). "Analysis of flow through vegetation. " Journal of Hydraulics division HY7: 871-883.
• Voinov, A., C. Fitz, et al. (submitted). "Coupling Raster and Vector Based Models." .