A numerical model for simulating variably-saturated flow in fractured porous media is established using the composite element method. The model has the following features. The fractures do not need to be discretized into specific elements

but are inserted into the porous matrix elements according to their geometric positions to form the composite elements. Within the composite element

the calculation equations for the fracture flow and the porous matrix flow are built respectively

which can be linked by the fracture-matrix interface and integrated into the composite element equation. The composite element equation has the same format as the conventional finite element equation

so the solving techniques for the conventional finite element equation can be applied to the composite element one. Using the techniques of under-relaxation iteration

mass matrix lumping and adaptive time-stepping

a calculation program for variably-saturated flow in fractured porous media is developed. The rationality and applicability of the model are verified by simulating the one-dimensional infiltration into dry soils and the flow in complex fractured aquifers. The simulation results provide a theoretical basis for further understanding groundwater flow characteristics in unsaturated fractured aquifers.