## 1 Introduction  
The construction of the underground transportation system has been gradually increasing in major cities. As the most advanced tunnel construction method at present, the shield method has become the first choice for urban tunnel construction in various regions of the world due to its advantages in low environmental impact, high construction efficiency, and high degree of mechanization. Therefore, a large number of existing shield tunnels have been constructed in urban underground space. Shield tunnels are generally assembled from reinforced concrete segments connected with steel bolts. Owing to their structural characteristics, existing shield tunnels are vulnerable to new tunnel construction underneath [1; 2; 3]. Because of the new tunnel undercrossing, longitudinal tunnel deformations may develop and threaten the safety and serviceability of the tunnel during long-term operation [4; 5]. Therefore, the evaluation of tunneling-induced longitudinal deformations of existing shield tunnels has become a growing concern. The longitudinal performance of a shield tunnel due to new tunneling can be studied by numerical models implemented in such commercial software as ABAQUS, FLAC, or PLAXIS, which are effective tools used to consider complex interaction between the tunnels and the ground [6; 7; 8]. However, the selection of constitutive parameters in the numerical simulation is difficult, and the modeling and computation generally cost substantial time. The complexity and inconvenience of the numerical model limit its application in engineering practice. By contrast, the one-dimensional two-stage analysis method is efficient and convenient for rapidly evaluating the effects of new tunnel excavation on longitudinal behaviors of shield tunnels in the primary design stage [9; 10; 11; 12]. The existing tunnel in these analysis methods is generally assumed as a one-dimensional continuous beam lying on various foundation models. The greenfield soil settlement induced by new tunneling at the existing tunnel location is first determined through empirical equations or analytical solutions [13; 14; 15; 16]. Then, the soil settlement can be used as an input and imposed to the tunnel for obtaining tunnel responses. [PERSON] and [PERSON] [17] introduced a simplified two-stage analysis method based on the Euler-Bernoulli beam theory and the Winkler foundation model to analyze the responses of existing shield tunnels to undercrossing construction. The shearing deformation of shield tunnels is considered in further research by adopting the Timoshenko beam for better simulation of longitudinal tunnel behavior. [PERSON] et al. [18] derived a Timoshenko beam solution to calculate the deformations of existing tunnels owing to new tunneling underneath. [PERSON] et al. [19] proposed an analytical method based on Timoshenko beam theory and a Kerr foundation model for new tunnel construction underneath and perpendicular to the existing shield tunnels. [PERSON] et al. [20] used a Timoshenko beam-based elastic-continuum model to assess the longitudinal responses of existing tunnels to single and twin tunneling. The long-term geological deposition and weathering process result in apparent spatial variability of soil properties [8; 21; 22; 23; 24; 25], which may lead to significant influences on the longitudinal tunnel responses. For the one-dimensional analysis method, the spatial soil variability is reflected by the variations in the soil parameters along the longitudinal direction of the tunnel [26; 27; 28]. However, there are limited studies available at present using the one-dimensional two-stage analysis method to investigate the new tunneling-existing tunnel interaction in spatially variable soil. Compared with other mechanical parameters of the soil, the soil elastic modulus has a relatively large spatial variability [29; 30].