Natural fault damage zones (FDZs) are often composed of complex arrays containing many thin faults. To account for FDZs in large-scale models, it is necessary to represent the flow effects of these many small faults by some form of upscaling. The explicit discretization of the geometry and topology of these thin faults and of the rock matrix typically leads to meshes comprising large numbers of elements, making it difficult to solve the flow equations efficiently. This paper proposes a practical technique based on the mixed finite element method (MFEM), permitting the permeability upscaling in two-dimensional fault damage zones that contain thin and isotropic low-permeability faults. This technique utilizes an implicit discretization scheme that is similar in principle to those used for discrete fracture networks. The scheme described here treats flow across faults accurately, but, unlike the schemes used for fractures, it neglects the flow along the faults because intrafault volumetric flow in thin faults is very small. This approach allows for a simpler numerical discretization and a reduction in numerical calculation effort. The scheme can be readily implemented in a standard mixed finite element code. Compared with the same MFEM based on an explicit discretization scheme, the two techniques are shown to produce results that agree closely. The implicit technique is shown to be more efficient than the explicit technique for a typical fault damage zone model. This technique is readily extendable to three dimensions.