Obtaining an accurate estimation of surface settlement during tunnel excavation is challenging due to the many factors that influence settlement. This study developed a mathematical model using multi-gene genetic programming for predicting the maximum surface settlement associated with earth pressure balance shield tunneling. Based on the field-monitored data collected from a metro construction project, five input parameters, namely, buried depth, face pressure at the top and center of the shield, advance rate, and grouting pressure, were used for the model development. Three statistical metrics for performance evaluation indicated the potential of the proposed model to estimate the maximum surface settlement. A parametric study was performed to explore the variation trend of the maximum surface settlement induced by different parameters, which validated the reasonability of the proposed mathematical model. A global sensitivity analysis revealed that the cover depth and grouting pressure were the two more influential parameters in developing the maximum surface settlement during shield tunneling than other selected features. The uncertainty analysis proved the robustness and reliability of the model in predicting the maximum surface settlement during tunneling. Comparative analysis demonstrated that the proposed model is easy to implement in routine design with acceptable error and could supersede simple empirical models.