Existing fracture network models are used almost exclusively to model flow and transport in saturated rock masses. When modelling shallow groundwater systems, the existence of the water table will have a significant effect on the flow behaviour and hence a fully saturated modelling approach is inappropriate. Thus, extension of an existing fracture network code to include the influence of the water table in unconfined conditions has been undertaken. The baseline code used is NAPSAC, a powerful finite-element code written for radioactive waste disposal research. Variably saturated flow module is characterised by the use of a single-variable capillary curve unsaturated flow model. These non-linear relationships are solved by the inclusion of a Newton-Raphson linearisation scheme coupled with direct-frontal or iterative PCCG matrix solution procedures. Additional features useful for unconfined flow simulation such as a Generalised Flux boundary condition are also included. Successful verification of the new code is demonstrated via a series of progressively more complex test cases, examining one, two and three-dimensional fluid flow. Fracture data from the southern UK chalk are used to demonstrate flow and transport capability in field geometries. The new code was found to operate successfully and to be suitable for future field and laboratory application.