Subsurface drains are important features in many problems involving shallow groundwater flow and have generally been represented in two-dimensional numerical models by specifying hydraulic head or fluid flux at one or more nodes in the computational grid. Here we describe a new approach for representing drains in three-dimensional variably saturated flow models. The method is based on a one-dimensional flow equation for the drain which is analogous to a one-dimensional open-channel flow equation; the numerical contributions to the system of global equations from the one-dimensional discretized equation are added to the three-dimensional, porous medium contributions. The inclusion of drain flow by this technique is readily adaptable to any three-dimensional numerical simulator. The accuracy of the method is demonstrated by comparing results obtained with the new method to steady state and transient results presented by Fipps et al. <1986>. The numerical formulation naturally accommodates flow in heterogeneous media and demonstrates excellent mass balance. Illustrative examples are presented which highlight the usefulness of the algorithm for studying the influence of drains in problems involving subsurface effluent disposal in heterogeneous aquifers and agricultural land drainage. The method will have application to many studies involving water and contaminant movement in drained soils. ¿ American Geophysical Union 1996