In this paper we study the diffusion-free liquid phase residence time distribution (RID) of laminar Taylor flow in square mini-channels numerically and analytically. We evaluate the RTD caused by the non-uniform laminar velocity field for co-current upward and downward Taylor flow from detailed numerical simulations obtained with a volume-of-fluid method by diffusion-free Lagrangian tracking of a set of virtual particles. The numerical RID curves are used to develop a compartment model for the RTD of a Taylor flow unit cell (which consists of one gas bubble and one liquid slug) in the fixed frame of reference. While the new model is conceptually similar to models from literature, it is refined with respect to (i) the delay time (i.e. the residence time of the fastest liquid elements, modeled by a plug flow reactor) and (ii) the difference in slope of the RID at small and large residence times (modeled by two parallel continuous stirred tank reactors with different mean residence time). It is shown that the refined model fits the numerical unit cell RTD reasonably well for different flow conditions. From the unit cell RTD model, the RTD for two and three unit cells in series is computed by a convolution procedure. The agreement between the numerical and the convolution-based RID for multiple unit cells is not satisfactory and indicates the inappropriateness of the convolution procedure for computing the diffusion-free RID of multiple unit cells in Taylor flow. The failure of the convolution procedure is attributed to the dividing streamline, which separates in Taylor flow the liquid phase in two regions, one with recirculating and one with bypass flow. In the absence of diffusion, the diving streamline is never crossed. As a consequence, locations of tracer particles in neighboring unit cells are not independent which invalidates the unit cell convolution approach. (C) 2012 Elsevier B.V. All rights reserved.