## 2 Materials and Methods
### Vibration Energy Comparison Method  
This section describes the specific details of the method used to identify the formation time of the new free surface. First, the observed wave \(S(t)\) and the predicted wave \(S^{\prime\prime}(t)\) of a multihole delay blasting were collected. The former was measured on-site and single-hole waves superposed the latter. The Hilbert transform found the instantaneous vibration energy of both waves and then the energy value was integrated to obtain the cumulative value at any moment. To compare the two accumulated energy curves and find the divergence point as the formation time of the new free surface, the details are as follows.  
Take the observed wave \(S(t)\) as an example. The VMD method can decompose \(S(t)\) into many intrinsic mode function (IMF) components. Each IMF component is a time-amplitude function \(c_{i}(t)\). The \(S(t)\) can be expressed as the sum of each \(c_{i}(t)\) and the residual \(r(t)\), as shown in Equation (6):  
\[S(t)=\sum c_{i}(t)+r(t), \tag{6}\]  
where \(S(t)\) is the observed wave, \(t\) is the time, \(c_{i}(t)\) is its IMF components, \(i\) is the number of the component, and \(r(t)\) is the residual.  
The Hilbert transform is performed on each IMF component \(c_{i}(t)\) to obtain \(H[c_{i}(t)]\). Then build the analytic function \(z_{i}(t)\) as shown in Equation (7), which is expressed in plural form:  
\[z_{i}(t)=c_{i}(t)+jH[c_{i}(t)]=a_{i}(t)e^{i\[PERSON]_{i}(t)}, \tag{7}\]\[a_{i}(t)=\sqrt{c_{i}^{2}(t)+H^{2}[c_{i}(t)]},\ \[PERSON]_{i}(t)=\tan^{-1}H[c_{i}(t)]/c_{i}(t).\]  
where \(a_{i}(t)\) is an amplitude function, \(\[PERSON]_{i}(t)\) is a phase function, \(d\[PERSON]_{i}(t)/dt=f_{i}(t).f_{i}(t)\) is the instantaneous frequency, e is the Napierian base, and \(j\) is the imaginary unit.  
Equation (8) obtains the Hilbert spectrum \(H(f,t)\) and removes the residual \(r(t)\):  
\[H(f,t)=\text{Re}\sum_{i=1}^{n}z_{i}(t)=\text{Re}\sum_{i=1}^{n}a_{i}(t)e^{i\Phi _{i}(t)}, \tag{8}\]  
where Re denotes taking the real part of the plural.  
By integrating the frequency, Equation (9) defines the Hilbert instantaneous energy \(IE(t)\):  
\[IE(t)=\int H^{2}(f,t)df, \tag{9}\]  
where \(IE(t)\) is the instantaneous energy, \(H(f,t)\) is the Hilbert spectrum, and \(f\) is the frequency.  
Equation (10) computes the cumulative value \(E(t)\) of the instantaneous energy \(IE(t)\) of the wave at each moment:  
\[E(t)=\int_{0}^{t}IE(\tau)d\tau, \tag{10}\]  
where \(E(t)\) is the cumulative value of the instantaneous energy.  
The observed wave \(S(t)\) is affected by the variation of the free surface condition while the predicted wave \(S^{\prime\prime}(t)\) does not consider the variation in the calculation. After the formation of the new free surface, the two accumulated energy curves \(E(t)\) and \(E^{\prime\prime}(t)\) of the instantaneous energy will show a divergence. Given that one of the sides involved in the comparison used a linear superposition model, the difference between the observed and predicted results can be measured using the linear correlation. Equation (11) calculates this linear correlation. Before the formation of the new free surface, the two curves will remain similar with a linear correlation close to one. After that, the blasting process is no longer linear and the similarity between the two will decrease.  
\[\rho(t)=\frac{\text{COV}[E(t),E^{\prime}(t)]}{\text{STD}[E(t)]\times\text{STD }[E^{\prime}(t)]}, \tag{11}\]  
where \(\rho(t)\) is the linear correlation of the two curves before the moment \(t\), COV is the covariance of the two variables, and STD is the standard deviation of the variables.  
Figure 2 shows how to identify the new free-surface formation time.  
Figure 2: Identification flow.