The subsurface in the Western part of the Netherlands is generally characterized as a regional aquifer, in which most of the lateral groundwater flow takes place, covered by a semi-permeable layer. The semi-permeable layer is often split into a phreatic layer on top of an aquitard (which has an even lower permeability). The phreatic layer and aquitard together are referred to as the top-system. The phreatic layer contains ditches and drains that transport excess water out of the area to, for example, a belt canal or a river. In regional groundwater models upscaling methods are applied to replace the top-system by a linear head-flux relationship (a Cauchy boundary condition) to take into account the interaction between many small surface water features and the groundwater in a lumped fashion. The upscaled system contains two parameters: the effective water level (p*) and the effective resistance (c*). The geology of the phreatic layer has a large influence on the value of these parameters, but existing upscaling methods are limited to homogeneous and isotropic phreatic layers. New formulas are derived for the effective parameters that take into account layered heterogeneity and vertical anisotropy. This multi-layer upscaling method is compared to an existing method, derived by de Lange, “A Cauchy boundary condition for the lumped interaction between an arbitrary number of surface waters and a regional aquifer”, Journal of Hydrology 266 (1999) (p. 250-261), which does not incorporate these aspects directly. Two upscaled models are created using the two upscaling methods and are compared to an analytic element model containing all features explicitly. The comparison is carried out for a cross-section with an extraction in the regional aquifer. In homogeneous isotropic phreatic layers both upscaling methods perform well. The drawdown in the regional aquifer differs by less than 3% as compared to the drawdown in the explicit model. When heterogeneity and anisotropy are introduced to the phreatic layer, the multi-layer method performs better with differences in drawdown remaining below 1%. The method by de Lange shows differences of up to 15%. The values of the effective parameters calculated by both models are similar in heterogeneous top-systems characterized by larger phreatic transmissivities, although the multi-layer method is almost always more accurate. For phreatic layers with smaller transmissivities, the multi-layer method yields significantly better estimates for the effective parameters. The new method does have limitations. It is less accurate in situations where surface water features are wide relative to the distance between them or when the bed resistance of the surface water features is large relative to the resistance of the aquitard. De Lange's approach is, in theory, applicable under those conditions, but does not take into account heterogeneity or vertical anisotropy. For most practical purposes, the multi-layer upscaling method is preferred and yields similar or more accurate results than de Lange's method.