## 1 Introduction  
Although pioneering clinical studies of focused ultrasound were carried out as early as the 1940s [1, 2], it did not attract intensive research interest until the end of the 20 th century and the beginning of 21 st century, during which several theoretical models were developed, improved and then broadly accepted [3, 4, 5, 6, 7, 8, 9]. In the past decades, high intensity focused ultrasound (HIFU) has played an increasingly significant role in the study of non-invasive therapies by demonstrating unique advantages in safety, effectiveness and high efficiency [10, 11, 12, 13, 14]. However, the applications of HIFU are still limited, and clinical treatments are only available for limited sites [11, 14, 15, 16]. One of the challenges confronting HIFU treatments is the spatial precision of tissue ablation. Several techniques such as ultrasound B-Scan and Magnetic Resonance Imaging (MRI) have been combined with HIFU to achieve real-time monitoring of focal areas [17, 18, 19, 20]. With these methods, the actual focal profiles were usually found to deviate from those predicted through theoretical models [21, 22]. As was indicated by [PERSON] et al., the shift of the focal point away from the prescribedposition caused by acoustic aberrations and non-linear wave propagating effects make it mandatory to evaluate the spatial accuracy of HIFU ablation [21]. [PERSON] et al. observed a focal shift of 1-2 mm and ascribed it to the layered distribution of tissues. Several different mechanisms might contribute to these deviations, such as the thermos-lensing effect [23; 24], bubble formation [25; 26; 27], acoustic radiation force [27] and the nonlinear nature of acoustic waves [21; 22]. [PERSON] et al. pointed out that the positioning error of the focal spots could be mainly related to the thermos-lensing effect and nonlinear propagation of ultrasonic waves [23]. When accounting for the complexity of wave propagation, e.g., ribs, abdomen tissues, blood vessels and other celiac organs between the transducers and the targeted area, these could constitute sources of acoustic scattering, diffraction, attenuation and dispersion etc. [28; 29], and negatively affect the precision of HIFU treatments. Therefore, it is very important to take into account the complexities of the acoustic paths by developing theoretical models for higher accuracy. Using the state of the art methods, the spatial distribution of HIFU field can be simulated with the well-known Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation [5; 6], the Westervelt model [22] or the Spherical Beam Equation (SBE) [9]. Because of its reasonable parabolic approximation, the KZK equation has been widely used for describing the propagation of finite amplitude acoustic beams emitted by focused transducers, with the only restriction being that the half angle of divergence of the transducer does not exceed \(16^{\circ}\)[7]. Existing algorithms to solve the KZK equation are usually balanced between accuracy and simplicity, making it possible to calculate and adjust HIFU fields in real time during treatments. In contrast, the heavy calculation burden of the Westervelt equation and the SBE limits their application. Therefore, it is not unexpected that the KZK equation has been preferred against the other two in both clinical situations and industry. In the KZK equation, the effects of diffraction, nonlinearity and attenuation have been taken into consideration and are described by separate terms. Therefore, it has been recognized to be more effective than the Rayleigh integral in simulating HIFU propagation [8]. Meanwhile, the parabolic approximation is also accurate enough in practice, by which the near-axis HIFU field can be calculated relatively accurately [7; 30; 31]. However as [PERSON] et al. have pointed out, it is hard to explain the shifts of focal positions with the traditional KZK equation [24]. Therefore, several questions are still open to discussion.