Natural fault damage zones (FDZs) are characterized by complex geometries and topologies, and by strongly-contrasting material properties. Accurate simulation of fluid flow in such systems is dependent on the method of discretization and the mathematical representation of the flow. In this paper, we focus on the conceptual and methodological issues that link a model of a heterogeneous system and its flow response. We study FDZs as our example, where each thin fault strand is a barrier to flow. We examine two contrasting discretization schemes and apply them to 2D FDZ models that contain a realistic array of linear fault traces. Both schemes produce results that are generally in good agreement, and agree with the results calculated by a more accurate (but computationally less efficient) reference scheme. However, differences occur when the discretization approach fails to maintain the fault connectivity (topology) of the input model. It is important to guide the modelling by identifying any continuous flow pathways in the matrix linking fluid inlets and outlets (Through-Going Regions, TGRs). We illustrate a new scheme that identifies all TGRs and determines a grid that is just fine enough to resolve them.