The widely used groundwater flow model MODFLOW offers a range of software packages to simulate the interaction between streams and groundwater in aquifer systems. However, these existing packages lack a general method to address the chaotic simulation sequences of stream segments and require these segments to be ordered by modelers as input to the code. Therefore, it is challenging to simulate a stream network divided into a large number of segments such as a canal irrigation system. In this study, the Streamflow Automatic Routing (SAR) package was developed, and an effective method is proposed to automatically determine the segment simulation sequence. The stream segment order in the SAR input file is arbitrary, which allows modifications of the stream network by removing segments directly and adding segments at the end of the segment group. This mainly includes two processes: scanning all the outlet channels of the water system and calling the recursive algorithm for each outlet channel of the water system. The SAR package was tested using a hypothetical stream–aquifer system and applied to a complex flow field in Aiding Lake of Turpan Basin, China. In the results, a close fitting between the simulation and observations shows that the SAR package precisely simulated the exchange flux between the steams and aquifer. The SAR package can significantly improve the efficiency of simulations in a complex stream network, and it can be widely used as a subroutine package of MODFLOW in agricultural irrigation areas where rivers and canals are interlaced.
