This paper developed an automatic quadrilateral meshing system to resolve the numerical simulation problems in the fields of geological mechanics, hydraulics, hydrology, and water resources. An improved combination method was proposed to automatically convert Delaunay triangular meshes into quadrilateral meshes by combining adjacent pairs of triangles through searching advancing-front edges. In order to implement the full conversion, seven combination patterns were established to select appropriate pairs of adjacent triangles. To handle the residual triangles, each nearest pair of residual triangles were moved to each other along a proper path and recombined into a quadrilateral. The movement of residual triangles was implemented through the destruction and reconstruction of neighboring quadrilaterals. A transformation template was established to convert the only residual triangle in each domain into quadrilaterals. To ensure the geometric topology of the resulting mesh, nine topological optimization modes were proposed to improve the topological connections of the degenerated quadrilaterals especially those on the boundaries of the mesh. An area-weighted Laplacian method was used to smooth the interior nodes, and an objective function method was employed to adjust the positions of special nodes. Finally, practical examples of two treatment strategies for river and overland flows were provided to demonstrate the accuracy and reliability of the meshing and optimization algorithms proposed in this paper.