Finite difference models such as MODFLOW are very flexible in representing geological detail because hydrogeological properties can be assigned on a cell by cell basis. However, large cells in regional models prevent accurate modeling of surface water and groundwater interactions, particularly the effects of wells near water bodies. Analytic element models can more realistically represent these detailed interactions (no grid), but are less flexible in representing the detailed three-dimensional geology. Using the advantages inherent in both methods, we coupled the analytic element model GFLOW to MODFLOW. GFLOW simulates groundwater flow in the shallowest part of the flow system where the surface water and groundwater interactions occur, while MODFLOW simulates flow in the rest of the system. For coupling, GFLOW uses a grid of leakage elements, coinciding with the MODFLOW cells, to model the leakage from the GFLOW layer into the upper MODFLOW layer. However, depending on the hydrogeology, the MODFLOW cell size used to discretize the GFLOW leakage may be too large to accurately represent this leakage. Some simple one-dimensional flow problems served to develop leakage grid size criteria. For cases where the MODFLOW cells are too large, GFLOW uses a finer leakage grid by subdividing the MODFLOW cells, while the average leakage in the sub-grid is assigned to the corresponding MODFLOW cell and passed on to MODFLOW.