Figure 8 shows a typical reading of the developed fiber optic strain in the shearing process (specifically demonstrated for experiment T1). The fiber strain response is also affected by the bending behavior of the ring, and therefore, a post-processing stage was performed in which the fibers were superimposed to nullify (or at least reduce) the bending contribution before they were introduced to the stress evaluation algorithm (which does not include such modes of deformations). Strain changes also appeared outside the ring, although to a much smaller extent. These relate to the fact that the base of the carriage translates and forces movement upon the "free" fiber inside the box, as well as outside (leading to minute displacements and small strains outside the ring area). Note that no temperature compensation was considered, as the experiments were performed in a temperature controlled room.

**Figure 7.** Conventional output from DSS tests.

**Figure 8.** Example of raw measurements along the optical fiber.

The strain data was accumulated for the six tests, with measurements being taken at translation intervals of about 0.1 mm until overall translation of 3 mm, and 0.5 mm intervals until translation of 8 mm. The analysis and interpretation of the results are provided in the following section.

#### **5. Interpretation of the Test Results and Discussion**

The interpretation here focuses on the stress condition at the time of shearing, assuming that a *K*0 condition prevailed prior to shearing. A *K*0 value of 0.54 was considered using Jaky's formula (*K*0 = 1 − sin *φ* ) [49], based on the measured friction angle of the sand, leading to an initial horizontal stress of 264 kPa.

Figure 9 shows a graphical representation of the horizontal stress distribution applied to the soil sample, as analyzed by the Fourier-based optimization at different stages of testing, expressed by the shear strain, *γ*.

**Figure 9.** Contact stress distributions at different stages of the shear tests.

Specifically, values at *γ* = 0.01, 0.02, 0.04, and 0.08 were plotted based on analyses with 3 and 5 modes. The radial distance from the center circle represents the amplitude of the horizontal normal stress to the soil sample (equal to the stress applied to the confining ring). The *K*0 (initial) condition is marked by a blue line, where the length of the light blue arrows represents its magnitude. A scale is given for reference, and the light red arrows extend to the largest calculated stress. Let us refer to a Cartesian coordinate system in which *z* represents the vertical direction, and *x* the direction of translation (*y* being the other horizontal direction). These are plotted in the figure.

As can be seen, with increasing shear strain, the horizontal stress increases in all directions. Overall, the magnitude of the developed stress profiles and their distribution are similar among the different tests. The magnitude of variation is similar to that observed in the applied shear stress in Figure 7. The analysis with five modes shows greater discrepancy between the tests, specifically in the higher spatial frequencies. This is because higher frequencies are more sensitive

to measurement errors. While the five modes analysis indicates a "square" type of distribution, the stress values in the translation direction (*σx*) and the perpendicular one (*σy*) are rather similar to those based on the three modes analysis. Hereafter, all analyses are performed with the five modes data, although similar results and trends are observed from the three modes analyses.

Figure 10 shows the development of the mean *σx*, taken as 0.5(*f*(*θ* = 0) + *f*(*θ* = *π*)). As can be seen, as the shearing process progresses, the horizontal stress increases, eventually becoming larger than the vertical stress (equal to 485 kPa) by a factor of almost two.

**Figure 10.** Mean *σx* throughout the tests, taken as 0.5(*f*(*θ* = 0) + *f*(*θ* = *π*)).

Assuming that this mean *σx* represents the average normal stress applied to sheared vertical planes, a stress path can be plotted in terms of (*σx*, *τxz*). This stress path is unknown from conventional DSS testing, and is a result of the new suggested approach. It can supplement the conventional (*σz*, *τzx*) stress path to allow constructions of Mohr stress circles throughout the shearing process.

Figure 11a shows the various stress paths established from the analysis of all tests. The Mohr circles associated with the peak shear stress, *τzx*, are presented, as well as that of the initial *K*0 condition (identical for all samples). The (*σz*, *τzx*) stress path is essentially identical for all samples, with a slight variation in the peak *τzx* value (as also demonstrated in stress-strain curves in Figure 7). The interpretation of the different samples is rather similar, and as an example, Figure 11b shows the interpretation for test *T*3, entailing the stress paths, initial and final Mohr circles, their poles, and a tangent Mohr-Coulomb cohesionless failure criterion.

**Figure 11.** (**a**) (*σz*, *τzx*) and (*σx*, *τxz*) stress paths for all tests; (**b**) interpretation illustrated on T3.

It can clearly be seen that the (*σz*, *τzx*) stress path reaches the Mohr-Coulomb failure line at tangent point (although this condition was not forced). This indicates that failure occurred along horizontal planes. The principal stress directions rotated significantly throughout the shearing process of the sample, by about 60 degrees, leading the pole to change its position (marked in the figure as "initial *Op*" and "*Op* at failure").

Repeating the same process on the other samples resulted in the same trends. The obtained friction angles, by the tangent line to the Mohr circle of each test, ranged from 27 to 29 degrees among the test. These values are similar to those obtained with the conventional setup and slightly less than those measured on a range of sands in direct shear [50].

The fact that the (*σz*, *τzx*) stress path reaches the failure envelope at the same location as the tangent to the Mohr circle infers that a direct evaluation of the strength properties from a conventional analysis assuming *τzx*/*σz* = tan *φ* is appropriate for these specific tests. However, this may not always be the case, especially in undrained (or zero volumetric strain) conditions, where failure occurs along vertical planes and the greatest shear stress experienced does not correspond to the maximum shear stress measured [44]. The ability to evaluate such a stress path development, as well as rotation in principal directions, may be important for analysis of undrained shearing and cyclic loading, where there might be interest in the change of confining mean stress due to compaction and dilation at shear direction reversal [51].

#### **6. Conclusions**

This paper presented a new approach for evaluating the horizontal stress condition in direct simple shear tests using high-resolution fiber optic sensing coupled to a Fourier based elastic solution of the confining rings. An optimization problem was formulated aiming to evaluate the best stress profile (in terms of Fourier coefficients) generating an analytical stress profile that would match the measured strain profile.

The approach was demonstrated with a commercial DSS apparatus, for which a new set of fiber optic-embedded confining rings was produced using 3D printing. Six drained direct simple shear tests were performed on clean sand to demonstrate the applicability of the approach. It was demonstrated how a stress path and Mohr circles can be evaluated with the new information derived from fiber optic sensing. Specifically, it was shown that in these particular (drained) tests, the soil truly fails along horizontal plane when considering a Mohr Coulomb yield condition.

The main focus of the paper was to conceptually develop a suitable approach for utilizing high-resolution fiber optic sensing for laboratory geotechnical testing, and to demonstrate its feasibility. Clearly, further research could be performed to improve the suggested approach. For example, studies could be performed to seek the best position of the optical fiber within the rings. Investigation of different ring materials or shape could be considered, as well as utilization of other high-resolution sensing techniques. These investigations, as well as investigations into more complex soil behavior, are beyond the scope of this paper, which simply attempts to take the first steps of advancing the use of high-resolution fiber optic sensing for laboratory geotechnical testing.

**Author Contributions:** Conceptualization, A.K.; Investigation, A.K., M.R., I.R. and I.P.; Methodology, A.K. and M.R.; Supervision, A.K. and I.R.; Writing—Original draft, A.K., M.R. and I.R.; Writing—Review & editing, A.K., M.R., I.R. and I.P.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Klar, A.; Levenberg, E.; Tur, M.; Zadok, A. Sensing for smart infrastructure: Prospective engineering applications. In Proceedings of the International Conference on Smart Infrastructure and Construction, Cambridge, UK, 27–29 June 2016; pp. 289–295, doi:10.1680/tfitsi.61279.289. [CrossRef]

c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Localization and Discrimination of the Perturbation Signals in Fiber Distributed Acoustic Sensing Systems Using Spatial Average Kurtosis**

#### **Fei Jiang 1,2, Honglang Li 2,\*, Zhenhai Zhang 1,\*, Yixin Zhang 3 and Xuping Zhang 3**

Received: 13 July 2018; Accepted: 21 August 2018; Published: 28 August 2018

**Abstract:** Location error and false alarm are noticeable problems in fiber distributed acoustic sensing systems based on phase-sensitive optical time-domain reflectometry (Φ-OTDR). A novel method based on signal kurtosis is proposed to locate and discriminate perturbations in Φ-OTDR systems. The spatial kurtosis (SK) along the fiber is firstly obtained by calculating the kurtosis of acoustic signals at each position of the fiber in a short time period. After the moving average on the spatial dimension, the spatial average kurtosis (SAK) is then obtained, whose peak can accurately locate the center of the vibration segment. By comparing the SAK value with a certain threshold, we may to some degree discriminate the instantaneous destructive perturbations from the system noise and certain ambient environmental interferences. The experimental results show that, comparing with the average of the previous localization methods, the SAK method improves the pencil-break and digging locating signal-to-noise ratio (SNR) by 16.6 dB and 17.3 dB, respectively; and decreases the location standard deviation by 7.3 m and 9.1 m, respectively. For the instantaneous destructive perturbation (pencil-break and digging) detection, the false alarm rate can be as low as 1.02%, while the detection probability is maintained as high as 95.57%. In addition, the time consumption of the SAK method is adequate for a real-time Φ-OTDR system.

**Keywords:** fiber distributed acoustic sensing; Φ-OTDR; localization; discrimination; kurtosis

#### **1. Introduction**

Fiber distributed acoustic sensing systems based on phase-sensitive optical time domain reflectometry (Φ-OTDR) have been widely used in many fields such as oil and gas pipeline monitoring, structure health monitoring, and perimeter security [1–4] due to their high sensitivity, large dynamic range, fully distributed manner, and simple configuration. The principal of Φ-OTDR is based on the interference of Rayleigh scattering lights returned from the fiber. The phase of the Rayleigh scattering lights can be affected by external vibration events, which will result in the amplitude variation of the Rayleigh scattering traces. In Φ-OTDR systems, locating the vibration events is an essential problem. Ideally, the vibration location can be obtained by subtracting a Rayleigh trace from an earlier stored trace [1]. However, this method may cause false locations due to the amplitude fluctuation in Rayleigh traces and the system noise of Φ-OTDR.

Various signal processing methods have been introduced to solve this problem. Some focus on the time-domain signal denoising. For example, Lu et al. proposed the moving averaging and moving differential method to reduce the noise power and also increase the frequency response range of the system [5]. Qin et al. introduced a wavelet denoising method to reject the random noises induced by varied polarization states in different position and detectors [6]. Shi et al. [7] and Qin et al. [8] proposed using empirical mode decomposition (EMD) to improve the signal-to-noise ratio (SNR) of location information for the vibration events. Qin et al. also proposed a curvelet denoising method to reduce the time domain noise and improve the detection performance of Φ-OTDR systems [9]. Ölçer et al. presented an adaptive temporal matched filtering method to reduce the effect of fading noise without any impact on the frequency response of the detection system [10]. Some locate the perturbations using the 2-D image characters of Rayleigh traces. For instance, two-dimensional edge detection can improve both the spatial resolution and the SNR of location information [11,12]. Adaptive 2-D image processing method of bilateral filtering algorithm can also enhance the SNR of vibration/intrusion location in Φ-OTDR [13]. Some pay attention to the spatial frequency information of the Rayleigh signals while locating, such as [14–16], which can simultaneously extract the vibration location and frequency information with an improved SNR. In addition, correlation dimension method has been proposed to solve the coherent fading problem and improve the multiple location results' standard deviation [17].

However, the methods above are sensitive to all kinds of perturbations. Perturbations caused by environmental interference may be located as well. For instance, the fiber cables buried underground may be influenced by pedestrians, vehicles, and trains; the fiber cables hanging in the air are usually perturbed by the wind and rain; and fibers that are embedded in concrete and stuck to pipelines could be affected by natural vibrations of the structures. All of these harmless perturbations may also cause the same disturbance results, which may obscure the real harmful perturbations and cause high false alarm rates (FARs). A new method based on multi-scale wavelet decomposition and backpropagation neural network [18] has been proposed to solve this problem. It can locate the perturbations by multi-scale wavelet decomposition and determine them by a 3-layer BP neural network. The results show that the FAR is decreased efficiently. In addition to this method, various event recognition methods [19–22] have been proposed to discriminate the signals in Φ-OTDR systems.

In this paper, a new localization and discrimination method based on the spatial kurtosis of Rayleigh signals is introduced for Φ-OTDR systems. The spatial kurtosis (SK) along the fiber is first obtained by calculating the kurtosis at each position of the fiber in a short time period, and then spatial average kurtosis (SAK) can be obtained by computing the moving average of SK on the spatial dimension. As a result, both the SK peaks and the SAK peaks may show the perturbation locations. Comparing with SK method, the SAK method can reduce the location error caused by coherent fading problem. We can also simply evaluate whether the perturbation is instantaneous destructive by comparing the SAK value with a certain threshold to some degree. This is because instantaneous destructive perturbations generally have high positive kurtosis, while stationary environmental noises usually have low negative kurtosis. Experimental results show that, the SAK method has higher SNR and lower location standard deviation compared with the previous localization methods. The time consumption of SAK method is low enough to achieve perturbation localization and discrimination in real-time, which is promising for a practical Φ-OTDR system.

#### **2. Methods**

#### *2.1. Kurtosis of Acoustic Signals in Phase-Sensitive Optical Time-Domain Reflectometry*

The kurtosis is a fourth-order statistic that can be used to describe the shape of a probability distribution of a random variable. It is defined as:

$$Kurt[X] = \frac{E\left[\left(X - \mu\right)^{4}\right]}{\left(E\left[\left(X - \mu\right)^{2}\right]\right)^{2}} - 3,\tag{1}$$

where *μ* is the mean of random variable *X*, and *E* [*t*] represents the expected value of quantity *t*. The kurtosis of any univariate normal distribution is 0. Thus we can evaluate a distribution by comparing its kurtosis with 0. Distributions with a positive kurtosis (called leptokurtic distributions) have infrequent extreme deviations (or outliers), while distributions with a negative kurtosis (called platykurtic distribution) produce no outliers. In order to analyze the kurtosis of Φ-OTDR signals, we simulated these three main types of signals in Φ-OTDR systems:

**Figure 1.** The simulated signals in phase-sensitive optical time-domain reflectometry (Φ-OTDR) systems. (**a**) Gaussian noise; (**b**) Transient signal; (**c**) Signal with the frequency of 5 Hz; (**d**) Signal with the frequency of 200 Hz; (**e**) 5–15 Hz band-limited white noise; (**f**) 45–65 Hz band-limited white noise.

We can see that the kurtosis of random noise is around 0, and the reason why it is not totally equal to 0 is that the length of it is too short (only 300). We can also see that signals caused by certain stationary environmental interferences have negative kurtosis, while transient signal (corresponding to instantaneous destructive perturbations) has high positive kurtosis. Thus we may discriminate such types of signals in Φ-OTDR systems by comparing their kurtosis with zero.

#### *2.2. Spatial Kurtosis of Rayleigh Optical Time-Domain Reflectometry Traces*

In Φ-OTDR systems, the Rayleigh backscattering signals can be composed into a matrix. The data picked up from each Rayleigh backscattering trace with the same time delay will form the time series of the corresponding sensing point. Suppose that each Rayleigh backscattering trace has *L* sampling points, which are determined by the sampling rate of the data acquisition card and the fiber length, the matrix with *N* received Rayleigh backscattering traces can be defined as:

$$\mathbf{R}\_{\mathbf{N}\times\mathbf{L}} = \begin{bmatrix} \mathbf{x}\_1^1 & \mathbf{x}\_2^1 & \cdots & \mathbf{x}\_L^1\\ \mathbf{x}\_1^2 & \mathbf{x}\_2^2 & \cdots & \mathbf{x}\_L^2\\ \vdots & \vdots & \ddots & \vdots\\ \mathbf{x}\_1^N & \mathbf{x}\_2^N & \cdots & \mathbf{x}\_L^N \end{bmatrix} \tag{2}$$

where *xn l* represents the amplitude of Rayleigh signal at position *l* at time *n*.

We accumulated *w* adjacent Rayleigh traces as one processing unit to calculate the kurtosis of the amplitude signal. The kurtosis of the signal at position *l* at time *n* can be defined as:

$$k\_{l}^{n} = \frac{\frac{1}{w} \sum\_{i=n-w+1}^{n} \left(\mathbf{x}\_{l}^{i} - \frac{1}{\mathcal{M}} \sum\_{i=n-w+1}^{n} \mathbf{x}\_{l}^{i}\right)^{4}}{\left(\frac{1}{w} \sum\_{i=n-w+1}^{n} \left(\mathbf{x}\_{l}^{i} - \frac{1}{w} \sum\_{i=n-w+1}^{n} \mathbf{x}\_{l}^{i}\right)^{2}\right)^{2}} - 3,\tag{3}$$

where *w* can be also regarded as the window size of kurtosis computing, which determines the temporal resolution of the localization and discrimination. The spatial kurtosis of the Rayleigh traces at time *n* is defined as:

$$\mathbf{SK}(n) = [k\_1^n, k\_2^n, \cdots, k\_L^n]. \tag{4}$$

By identifying the local maxima and minimum of SK at time *n*, we could estimate whether there are perturbations and locate them accurately on the sensing fiber. Moreover, we introduced a time step *s* to reduce the computation in real-time Φ-OTDR systems. It means that the time interval of calculating SK is *s* × *fs* (where *fs* represents the pulse repeat frequency). Composing the SK vectors at different times into a new matrix

$$\mathbf{SK} = \begin{bmatrix} k\_1^1 & k\_2^1 & \cdots & k\_L^1 \\ k\_1^{1+s} & k\_2^{1+s} & \cdots & k\_L^{1+s} \\ k\_1^{1+2s} & k\_2^{1+2s} & \cdots & k\_L^{1+2s} \\ \vdots & \vdots & \cdots & \vdots \end{bmatrix} \tag{5}$$

a two-dimensional image of SK over the monitoring duration can be obtained. Then the perturbation locations could be identified clearly from this image.

#### *2.3. Moving Average on the Spatial Dimension*

In Φ-OTDR systems, the accurate location of perturbations can be figured out by [7]:

$$z = z\_{\rm cent} - \frac{T\_p \upsilon}{4},\tag{6}$$

where *zcent* is the center of the vibration response segment, *Tp* is the optical pulse width, and *v* is the light velocity in fiber. In order to obtain the accurate location of the perturbation in Φ-OTDR systems, we should figure out the center of the vibration segment. However, some subsections of the fiber vibration segment may be in insensitive states due to the coherent fading problem [23]. Figure 2a shows 100 consecutive averaged Rayleigh traces, which were acquired by our Φ-OTDR system. The optical pulse width used here is 100 ns, which means that the spatial resolution is 10 m. The vibration segment is from 5213 m to 5223 m. We can see that two subsections of the vibration segment are in sensitive state while the middle subsection is in insensitive state. This may make it difficult to locate the accurate vibration segment center and result in a location error and a bad locating repeatability.

For the SK localization method this problem exists as well. As shown in Figure 2b, the SK curve of the original Rayleigh traces splits into two peaks in the vibration segment. The two SK peaks correspond to the two sensitive subsections in Figure 2a, while the SK valley corresponds to the insensitive subsection in Figure 2a. It is hard to figure out the accurate vibration location by the SK curve. Virtually, we can solve this problem by calculating the moving average of SK on the spatial dimension. The result can be called SAK, which is defined as:

$$\mathbf{SAK}(n) = \frac{1}{M} \left[ \sum\_{l=1}^{M} k\_l^n, \sum\_{l=2}^{M+1} k\_l^n, \dots, \sum\_{l=L-M+1}^{L} k\_l^n \right] \tag{7}$$

where *M* represents the moving average number, which depends on *Tp* and the data acquisition card (DAQ) sampling rate (SR) of the Φ-OTDR system. It should be equal to the number of sampling points in an optical pulse width (*M* = *Tp* × *SR*). When the vibration occurs, the location of SAK peak should be the center of the vibration segment. We can see from Figure 2b that the SAK peak virtually locates at the vibration segment center. Then we can obtain the accurate location of the vibration according to Equation (6).

**Figure 2.** The vibration locating results when coherent fading problem exists. (**a**) 100 averaged Rayleigh traces; (**b**) The spatial kurtosis (SK) and special average kurtosis (SAK) curves of the Rayleigh traces.

#### **3. Experimental Results**

#### *3.1. Experimental Setup*

The experimental setup of our Φ-OTDR system is shown in Figure 3. A narrow line-width laser (NLL) is used as the light source. The laser light is divided into two parts by a 10:90 optical coupler (OC1). Then 10% of the light is used as the local light to perform a coherent detection, and 90% of the light is modulated into optical pulses with a 200 MHz frequency shift by an acoustic optical modulator (AOM). After being amplified by an erbium-doped fiber amplifier (EDFA), the optical pulses are launched into an optical circulator. The Rayleigh backscattering light from the sensing fiber then interferes with the local light at a 50:50 optical coupler (OC2) and result in a beat light. The beat light is then detected by a balanced photodetector (BPD) through a differential way, which could improve the SNR by 3 dB [4]. The output electric signal from the BPD is then sampled by a DAQ with a 1 GS/s sampling rate. The Rayleigh trace amplitude is finally demodulated by Hilbert transform and downsampled to 250 MS/s to reduce the computation load. In our experiments, the pulse repeat frequency is 1 kHz and the pulse width is 100 ns, which means that the maximum frequency response and the spatial resolution of our system are 500 Hz and 10 m respectively.

**Figure 3.** The experimental setup. NLL: narrow line-width laser; OC: optical coupler; AOM: acoustic-optic modulator; EDFA: Erbium-doped fiber amplifier; Cir: circulator; FUT: fiber under test; BPD: balanced photo detector; DAQ: data acquisition; PC: personal computer.

Both indoor experiment and outdoor experiment were conducted to evaluate the proposed method. In the indoor experiment, a piezoelectric transducer (PZT) cylinder glued with 10 m long fiber was added after the ~5.074 km long fiber under test (FUT) in our system. The PZT was used to simulate the environmental interferences such as structural free vibration and vehicle engine. Pencil-break vibration was used to simulate the civil structure crack, which was felt bya2m long fiber loops glued to a flame retardant (FR-4) plate. The pencil-break vibration was conducted by an HB pencil lead with 0.5 mm diameter. The fiber loops was put at ~5.164 km. In the outdoor experiment, a ~100 m long fiber cable was placed outdoor. The first ~10 m of it was buried beneath a flower bed with 0.3 m deep, while the rest ~90 m of it was placed in a corner. Digging the soil near the buried fiber cable was inducted as the instantaneous destructive perturbations. The flower bed is near to the road, so the fiber cable buried in the soil might be affected by the pedestrians and vehicles, which could be regard as the harmless perturbations.

#### *3.2. Localization Results of Instantaneous Destructive Perturbations*

A total of 200 tests, including 100 pencil-break tests (50 of which were under 100 ns pulse width and 50 of which were under 500 ns pulse width) and 100 digging tests (50 of which were under 100 ns pulse width and 50 of which were under 500 ns pulse width) were carried out to evaluate the localization performance of proposed method. The pencil-break tests and digging tests were performed at two fixed positions (~5164 m for pencil-break and ~5220 m for digging).

In this part, we compared the proposed method with four previous methods, including moving average and differential (MAD) [5], two-dimensional edge detection (TED) [11], wavelet denoising (WD) [6], and multi-scale wavelet decomposition (MWD) [18]. The comparisons include locating SNR, locating repeatability, and time consumption. The parameters of these methods are tuned to get the relatively high locating SNR. The moving average number of MAD was 30; the Sobel operator size of TED was 3 × 3 for pencil-break and 5 × 5 for digging; the wavelet basis function of WD and MWD was db3, and the decomposition level was 6. The window length used to compute kurtosis was 300, the time step was 150, and the moving average number of SAK was 25.

#### 3.2.1. Signal-to-Noise Ratio Improvement of the Locating Information

Locating results of pencil-break vibration and digging are shown in Figures 4–7. Here, SNR is defined as the ratio between the signal peak intensity (S) and the maximum background noise intensity (N) with the equation of 20 lg(S/N) [18]. Figures 4 and 5 respectively show the locating results of pencil-break vibration and digging under 100 ns pulse width. We can see that the SAK method obtains the highest SNR (23.68 dB for pencil-break and 24.32 dB for digging) among all the methods. Figures 6 and 7 respectively show the locating results of pencil-break vibration and digging under 500 ns pulse width. Likewise, the SAK method has the highest SNR (22.28 dB for pencil-break and 23.36 dB for digging) among all the methods. We can see that when the pulse width is increased from 100 ns to 500 ns, the SNR of MAD, TED, and MWD significantly reduce, whereas the SNR of SK and SAK maintain high. The SNR of SAK method is higher than SK, which means that moving average on the spatial dimension may improve the locating SNR.

**Figure 4.** The locating results of pencil-break vibration under 100 ns pulse width. (**a**) The result of moving average and differential (MAD) method; (**b**) The result of two-dimensional edge detection (TED) method; (**c**) The result of multi-scale wavelet decomposition (MWD) method; (**d**) The result of wavelet denoising (WD) method; (**e**) The result of SK method; (**f**) The result of SAK method. SNR: signal-to-noise ratio.

**Figure 5.** The locating results of digging event under 100 ns pulse width. (**a**) The result of MAD method; (**b**) The result of TED method; (**c**) The result of MWD method; (**d**) The result of WD method; (**e**) The result of SK method; (**f**) The result of SAK method.

**Figure 6.** The locating results of pencil-break under 500 ns pulse width. (**a**) The result of MAD method; (**b**) The result of TED method; (**c**) The result of MWD method; (**d**) The result of WD method; (**e**) The result of SK method; (**f**) The result of SAK method.

**Figure 7.** The locating results of digging event under 500 ns pulse width. (**a**) The result of MAD method; (**b**) The result of TED method; (**c**) The result of MWD method; (**d**) The result of WD method; (**e**) The result of SK method; (**f**) The result of SAK method.

The average locating SNRs of each method on the 200 locating tests are shown in Table 1. We can see that, comparing with the average of the MAD, TED, WD, and MWD methods, the SAK method improves the pencil-break locating SNR by 13.95 dB (under 100 ns pulse width) and 16.6 dB (under 500 ns pulse width); and the digging locating SNR by 10.60 dB (under 100 ns pulse width) and 17.3 dB (under 500 ns pulse width).

**Table 1.** The average locating SNRs of each method on the 200 locating tests.

#### 3.2.2. Locating Accuracy

The box plot of locating results is shown in Figure 8. We can see that for both pencil-break and digging, the range of the locating results gets wider when the pulse width increases from 100 ns to 500 ns. Thus the increasing pulse width may degrade the locating accuracy. We can also see that the SAK method has more consistent results than other methods for both pencil-break vibrations locating and digging events locating. Quantitatively, the averages and standard deviations of pencil-break and digging locating results by the 6 methods under two different pulse widths are shown in Table 2. We can see that the averages of the locating results by these methods show no significant differences, which is because for all of the methods the locating results may fluctuate above and below the identical real locations. However, the standard deviations of the locating results differ. We can see that the SAK method has the lowest standard deviation for both pencil-break and digging locating. When the pulse width is 100 ns, the standard deviation of the pencil-break and digging locating results by SAK method

reduce from 1.7 m and 1.4 m (the average locating standard deviation of the previous methods) to 0.8 m and 0.7 m, respectively. Two-sample F-tests show that this improvement is statistically significant for all of the previous methods (*p* < 0.0056). When the pulse width is 500 ns, the pencil-break and digging locating standard deviation of SAK method are only 1.9 m and 1.8 m respectively, while the average locating standard deviation of the previous methods are as high as 9.2 m and 10.9 m respectively. Two-sample F-tests show that this improvement is also statistically significant for all of the previous methods (*p* < 10−19). Comparing with SK method, the SAK method decreases the localization standard deviation of pencil-break and digging under 500 ns pulse width by 8.4 m and 10.4 m, respectively, which means that moving average on the spatial dimension can effectively improve the locating accuracy.

**Figure 8.** The statistic results of locating results. (**a**) The result of pencil-break vibration under 100 ns pulse width; (**b**) The result of pencil-break vibration under 500 ns pulse width; (**c**) The result of digging under 100 ns pulse width; (**d**) The result of digging under 500 ns pulse width.

#### 3.2.3. Time Consumption

For a real-time Φ-OTDR system, the locating speed should be faster than the data generating speed. Therefore, the time consumption of the localization algorithm should be as less as possible. We compared the time cost of SAK method with the previous locating methods under different fiber length, shown in Figure 9. The methods are performed under 5000 Rayleigh traces (corresponding to 5 s long data in our system). The CPU used in our PC is an Intel Core i7-4930K. The algorithms are executed by MATLAB R2017a (MathWorks, Inc., Natick, MA, USA). We can see that the time consumption of SAK method is close to that of the MAD and TED methods and it outperforms both the WD and MWD methods. When the fiber length is 30 km (the spatial sampling interval is 1 m), the SAK method takes only 1.77 s. Consequently, SAK method can meet the time consumption requirement of the real-time Φ-OTDR system.

**Figure 9.** The time consumption of various locating methods under different fiber length.

#### *3.3. Discrimination of Perturbations*

In Section 2.1, we analyzed the kurtosis of different types of signals in Φ-OTDR. The simulation results indicated that we may discriminate the random noise, instantaneous destructive perturbations, and stationary environmental interferences by kurtosis in simple scenarios. In this section, we performed two experiments to validate the simulation results. In addition to using pencil-break vibration and digging as the instantaneous destructive perturbations, we drove the PZT with various frequencies to simulate the interference signals.

Figures 10 and 11 show the locating results of various tests using SAK method. Figure 10a shows the SAK curves of 5 Hz PZT vibration plus pencil-break vibration. The SAK trough at ~5075 m and the SAK peak at ~5164 m locate the PZT vibration and the pencil-break vibration respectively. The SAK image of PZT plus pencil-break is shown in Figure 10b. One conspicuous dark stripe can be observed at ~5075 m, which indicates the PZT vibration. One bright point is observed at ~5164 m, which indicates the pencil-break vibration. Figure 10c,d show the locating results of 200 Hz PZT vibration plus five digging events. The SAK trough at ~5075 m and the SAK peak at ~5220 m in Figure 10c correspond to the PZT vibration and digging events respectively. The dark stripe at ~5075 m and 5 bright points at ~5220 m shown in Figure 10d, respectively, locate the PZT vibration and digging events as well, which can also be seen clearly in the three-dimensional SAK image (shown in Figure 11). It should be noted that when the PZT starts to work and stops working, the SAK may be high as well. Thus such "onset" and "offset" events could also be regarded as instantaneous destructive perturbations in that moment.

**Figure 10.** The location results of piezoelectric transducer (PZT) and instantaneous destructive perturbations (pencil-break and digging) through proposed method. (**a**) The SAK curves of 5 Hz PZT vibration plus pencil-break vibration; (**b**) The SAK image of 5 Hz PZT vibration plus pencil-break vibration; (**c**) The SAK curves of 200 Hz PZT vibration plus digging; (**d**) The SAK image of 200 Hz PZT vibration plus 5 digging events.

**Figure 11.** The three-dimensional SAK image of 200 Hz PZT vibration plus 5 digging events.

In order to generally evaluate the discrimination performance of proposed method, we computed and analyzed the SAK value of a set of various tests (including 126 pencil-break vibrations, 168 digging events, 294 PZT vibrations with various frequencies, and 302 environmental noises). The statistic SAK result of various Φ-OTDR signals is shown in Figure 12a, where IDP means the instantaneous destructive perturbations (pencil-break vibration and digging event); Nmax and Nmin are the SAK maximum and SAK minimum of the noises; IF means interference (PZT vibration with different frequencies including 5, 50, 100 and 200 Hz). We can see that different types of signals have different SAK value interval.

**Figure 12.** The statistic results of perturbation discrimination by SAK method. (**a**) The SAK distribution of different types of signals; (**b**) The receiver operating characteristic (ROC) curves of perturbation discrimination. IDP: instantaneous destructive perturbations; IF: interference.

We used true positive rate (TPR) and false positive rate (FPR) to evaluate the identification ability of the SAK method, where the TPR and FPR are defined as:

$$\text{TPR} = \frac{\text{TP}}{\text{TP} + \text{FN}'} \text{ FPR} = \frac{\text{FP}}{\text{TN} + \text{FP}}. \tag{8}$$

Here, TP is the count of correct positive prediction, FP is the count of wrong positive prediction, TN is the count of correct negative prediction, and FN is the count of wrong negative prediction. The TPR is also known as the probability of detection, while the FPR is also known as the probability of false alarm. By plotting the TPR against FPR at various threshold settings, the receiver operating characteristic (ROC) curve is obtained, shown in Figure 12b. Four-fold cross-validation was used to evaluate the TPR and FPR of the perturbation discrimination by SAK. The SAK thresholds of discriminating the pencil-break (or digging) and PZT vibration were fitted on the training data of each fold using a logistic regression model. The evaluation results on test data of each fold are show in Table 3.

We can see that the average TPR, FPR and accuracy of discriminating pencil-break or digging from noise are 95.57%, 1.02%, and 97.32%, respectively, while the average TPR, FPR, and accuracy of discriminating PZT vibration from noise are 92.54%, 4.38%, and 94.13%, respectively.

**Table 3.** The true positive rate (TPR), false positive rate (FPR), and accuracy (Acc.) of signal discrimination using SAK method.

#### **4. Discussion**

Compared with the previous locating methods for Φ-OTDR systems, the SAK method has the following advantages:

The SAK method has two main defects. One is that the localization time interval is directly related to the time step *s* in Equation (5). When the time step is 150 and the sampling rate is 1 kHz, the time interval of localization will be 0.15 s. However, for an Φ-OTDR threat alarm system this is totally acceptable. The other is that the SAK method only captures the short-time characteristics of the perturbations. Instantaneous destructive perturbations and stationary interferences might be discriminated by SAK, whereas long-term perturbations still require more features to determine whether they are a threat.

#### **5. Conclusions**

In this paper, we proposed a new perturbation localization and discrimination method for Φ-OTDR systems. Pencil-break vibration and digging events are used as the instantaneous destructive perturbations while PZT vibration is used to simulate interference. Experimental results show that, comparing with previous methods (MAD, TED, WD, MWD), the proposed method can locate the instantaneous destructive perturbations with a higher SNR. We also introduced moving average on the spatial dimension to improve the locating accuracy when coherent fading problem occur, thus leading

to the smaller standard deviation of locating results. In addition, the SAK method can discriminate the instantaneous destructive perturbations from system noise and stationary environmental interference to some extent. Due to the little time consumption, SAK can be a promising real-time locating and perturbation discrimination method for Φ-OTDR systems.

**Author Contributions:** Conceptualization, F.J. and H.L.; Formal analysis, F.J.; Funding acquisition, H.L., Z.Z. and X.Z.; Methodology, F.J.; Resources, Z.Z. and Y.Z.; Software, F.J.; Supervision, H.L.; Validation, F.J.; Writing-original draft, F.J.; Writing-review & editing, F.J. and Y.Z.

**Funding:** This work is funded by National Natural Science Foundation of China under Grant No. 61627816 and No. 61773059.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Distributed Temperature Sensing Monitoring of Well Completion Processes in a CO2 Geological Storage Demonstration Site**

#### **Dasom Sharon Lee 1,2, Kwon Gyu Park 1,\*, Changhyun Lee 1,\* and Sang-Jin Choi 1**

Received: 8 November 2018; Accepted: 2 December 2018; Published: 3 December 2018

**Abstract:** The Distributed Temperature Sensing (DTS) profiles obtained during well completion of a CO2 monitoring well were analyzed to characterize each well completion process in terms of temperature anomalies. Before analysis, we corrected the depth by redistributing the discrepancy, and then explored three temperature calibration methods. Consequently, we confirmed the depth discrepancy could be well corrected with conventional error redistribution techniques. Among three temperature calibration methods, the conventional method shows the best results. However, pointwise methods using heat coil or in-well divers also showed reliable accuracy, which allows them to be alternatives when the conventional method is not affordable. The DTS data revealed that each well completion processes can be characterized by their own distinctive temperature anomaly patterns. During gravel packing, the sand progression was monitorable with clear step-like temperature change due to the thermal bridge effect of sand. The DTS data during the cementing operation, also, clearly showed the progression up of the cement slurry and the exothermic reaction associated with curing of cement. During gas lift operations, we could observe the effect of casing transition as well as typical highly oscillating thermal response to gas lifting.

**Keywords:** distributed temperature sensing; temperature calibration; depth correction; fiber optic sensing; well completion monitoring

#### **1. Introduction**

The most distinctive advantage of fiber optic Distributed Temperature Sensing (DTS) technology over other point sensing and wireline logging technologies is its 'distributed sensing' feature. The distributed sensing, as the fiber itself is the sensing medium, enables obtaining a temperature profile for the entire fiber length with much higher spatial density than point sensors in a period. It also enables continuous real-time temperature profiling at any time throughout the well, once it is permanently installed inside the well [1]. Unlike the production logging run, which obtains temperatures after the event, the DTS can identify in-well performance changes as the events occur during production and shut-in [2]. In addition, fiber optic distributed sensing has many more advantages compared to conventional point sensors or wireline logging as it can be deployed in any harsh or unusual environment and it is small, light, corrosion–resistant, and has a long sensing range, good sensitivity, and electromagnetic immunity.

In May of 2017, as a part of an onshore pilot-scale CO2 sequestration demonstration project in Korea, multifunctional monitoring wells with discrete sensors and fiber optic cable for DTS and Distributed Acoustic Sensing (DAS) have been completed. A key consideration factor in developing injection and/or monitoring wells for CO2 geological storage is the integrity of the well that prevents the wells from becoming a conduit for CO2 leakage. Generally, many sensors and pieces of equipment are installed in CO2 monitoring well depending on the aim of the monitoring. Therefore, more care is needed for well design, installation, and completion process in order to ensure that the many sensors and the equipment perform at their best as well as to maintain integrity. Since we applied a permanent method using gravel packing and cementing for the completion, which is irreversible unlike a non-permanent method that uses a packer, monitoring each well completion process is critical to determine whether each process has gone as planned.

For the purpose of monitoring each well completion process in terms of temperature anomaly patterns, we monitored DTS data during whole completion processes without any intervention. However, on-site analysis was limited to qualitative interpretation because temperature pre-calibration and depth correction were not possible due to the field schedule and conditions.

In this paper, we applied a depth correction method and three temperature post-calibration methods for more quantitative and accurate analysis. Then, we characterized the thermal characteristics of each well completion process in order to explore the feasibility of DTS as a well completion monitoring and control tool.

The uncertainty of depth accuracy is inevitable in practical application of DTS to boreholes because it is difficult to detect any depth errors due to stretching and wrapping of the cable during installation. In this paper, conventional depth correction is applied by considering both stretching and wrapping of cable [3]. Then, we confirmed the accuracy of depth correction with our gas lift temperature profile using the depth of casing transition as a priori information.

The conventional way of temperature calibration is submerging both ends of the fiber in ice or hot water to calibrate the temperature offset to an absolute value and its rate of change along the fiber length [3]. It is commonly recommended to apply temperature calibration before installation. In many cases, however, calibration before installation is not possible or limited due to practical factors like site schedule and conditions. This was our case, therefore, we corrected temperature after installation. In addition to a conventional temperature calibration method of submerging the long length in water, we also explored the possibility of two types of pointwise post-calibration methods that may overcome the practical limitation of the submerging method on the site.

During well completion, real-time DTS data were analyzed to monitor the entire well completion process in terms of temperature changes and to control or optimize the process. However, it was limited to a qualitative analysis because the depth correction and temperature calibration were not done. After the depth correction and temperature calibration, thermal events of each well completion process could be characterized more precisely, and each well completion process could be identified with distinctive temperature anomaly pattern. The progress up the well of sand and cement slurry were monitored with clear step-like temperature change. Thermal events due to exothermal reaction during cement curing process was also clearly observed. Consequently, we could conclude that observing the entire well completion process with pre-calibrated DTS data can be a promising tool of real-time control and optimization of well completion process.

#### **2. Materials and Methods**

#### *2.1. Site Description*

The thermal data were collected from the monitoring well, JG-M, which was completed at the end of May 2017 as a part of onshore pilot-scale CO2 geological storage demonstration program at Janggi field. The Janggi field site is located at Pohang city which is the southeastern part of Korean Peninsula as illustrated in Figure 1. The total depth (TD) of the monitoring wells is 1092.44 m, and the well is slanted by 2◦ to the southeast direction.

**Figure 1.** Location of Janggi field site including relative position of auxiliary monitoring well (JG-7-1), monitoring well (JG-M) with its wellhead, and injection well (JG-I).

#### *2.2. Well Schematic and Completion Process*

The monitoring sensors and equipment deployed in JG-M include 32-level electrode array, 24-level 3C geophone strings, a U-tube [4] fluid sampler, two piezo pressure/temperature sensors (P/T), an experimental pH sensor, and a fiber cable for DAS, DTS, and heat-pulse monitoring [5]. The schematic of the Janggi JG-M monitoring well and the completion assembly employed is presented in Figure 2 to help understand the thermal events associated with completion, equipment and materials.

**Figure 2.** Schematic of the monitoring well, JG-M and the completion assembly employed with gravel packing and cementing.

The 32-level electrodes were installed for surface-to-borehole and cross-well Electrical Resistivity Tomography (ERT) to image the distribution and migration of the CO2 plume. Due to fact the 32 electrodes cover a 93 m interval around the reservoir formation (JGF = Janggi conglomerate formation), almost all materials below 803 m should be non-metallic or electrically insulated. Therefore, instead of steel casing, 3-1/2" Fiberglass Reinforced Plastic (FRP) casing was used at 300.54 m interval below 791.90 m, where the slotted (depicted by dashed line in Figure 2) and non-slotted casing alternated to guarantee fluid circulation and to maintain strength against pressure. The 3C geophones were installed for the purpose of acquiring Vertical Seismic Profiling (VSP) and passive micro-seismic monitoring data. The P/T sensors were installed at the top and bottom boundary of JGF reservoir, which were also used for pressure and temperature during the well completion. The fiber cable turnaround which protect spliced part of fiber core for dual-ended mode measurement was installed at 1082.56 m. The gas-lift manifold for fluid production and cementing stage tool for cementing were installed at 605.55 m and 909.94 m, respectively.

All sensors and instruments are installed in the tubing by the convey method by attaching sensors and control line to 3-1/2" steel and FRP casing. A series of centralizers and clamps were used to protect the cables during installation and cement jobs. For the FRP casing section, we used nylon decentralizers which were attached to the casing with 3/8-inch set screws. The mid joint clamp was used to hold the stainless-steel tubes. For the steel casing section, cannon cross coupling protectors with fins provide centralization during cementing.

The two different patterns in the well schematic shown in Figure 2 depict the well completion schemes behind the casing. The bottom section below 910 m indicates the gravel packing section. The basic requirement of the reservoir section is to guarantee the circulation of pore fluid and CO2. The ERT and VSP sensors should be well coupled with the surrounding formations for sensing, which is the reason for choosing the gravel packing over a packer. Above 910 m of the well was completed with cementing to ensure the pressure isolation of the reservoir section and the mechanical integrity of well that prevents the well from becoming a conduit for CO2 migration.

After landing the tubing hanger, DTS system and permanent pressure temperature gauge was used to monitor the process of gravel packing and cementing operation in May 2017. We directly injected 12/20 sand into the annulus and recirculated the fluid up with the 1.9-inch IJ tubing. After the gravel packing procedure, cementing is conducted by circulating slurry into the well with a stage cementing tool at 909.94 m.

#### *2.3. Fiber Optic and DTS Interrogator Configuration*

The fiber optic control line used in JG-M consists of four optical fiber lines for DTS and DAS and two copper conduct lines for heat-pulse monitoring [5], which are protected by stainless steel sheaths and polypropylene jackets to prevent them from being damaged during installation. Among the four optical fiber lines two single-mode fibers are for DAS and the other two are multi-mode fibers for DTS. Temperature profiles are obtained in dual-ended mode by splicing two multi-mode fibers at the end of the control line. The schematic of optical fiber control line and specification, and turnaround sub are as shown in Figure 3.

Silixa's XT-DTS DTS interrogation system was used to monitor the temperature changes for five days during completion and for a month during curing. The data was obtained using a passive mode where no active heating was done to the cable itself. Our interrogation unit can measure temperature profiles up to 5 km, with minimum spatial and temporal resolution of 0.25 m and 5 s, respectively. The temperature profiles were collected with a spatial sampling interval of 0.25 m, with a temporal interval of 600 s during well completion monitoring from May to June of 2017 and 60 s during background monitoring and temperature calibration tests from November, 2017 to April, 2018.

**Figure 3.** Schematic of fiber optic control line and specifications (**above**), and photo of turnaround sub for dual-ended configuration (**below**).

#### **3. Depth and Temperature Correction**

There are many potential factors that can cause errors when using the DTS system. The errors which arise from the instrumentation, the fiber, the nature of the installation, and other sources show up as depth discrepancies and temperature errors [3]. In order to specify the exact location of thermal events and to analyze the events more quantitatively, we applied depth correction and temperature calibration before interpreting thermal events for the well completion processes.

#### *3.1. Depth Correction*

Depth mismatch is an unavoidable issue in the field. Unlike temperature correction, depth discrepancies cannot be pre-calibrated before installation. Even though the cables were carefully handled during the installation, many factors such as cable coiling up around the tubing, stretching of fiber optic cable itself, the usage of a decentralizer, and inner fiber optic twisting may inevitably introduce depth discrepancies. Although it may depend on the application, Smolen and van der Spek [3] stated that generally a depth mismatch of ± 0.3 m is realistic for most applications. With the depth issue comes the issue of determining the end of the fiber depth. The length measured by DTS should be mapped for depth correction to a real physical profile length of interest. In dual-ended or partially returned configurations in borehole applications, the end of the fiber can be easily identified with the peak temperature due to ambient geothermal gradient and it can be confirmed in the raw data by a sudden intensity drop around the splicing point.

In our case, the peak temperature was observed mostly around 1197.114 m in fiber length. On the other hand, the maximum intensity drop occurred between 1196.606 m and 1196.860 m with the difference of 21.17 A.U. (relative decrease of 10.4%) and the second maximum drop occurred between 119.860 m to 1197.114 m with the difference of 18.62 A.U. (relative decrease of 9.3%). Although there is a difference of 0.5 m between the peak temperature point and the maximum intensity drop point, we used the peak temperature point as the fiber-end to determine total fiber length because such a difference is negligible considering our spatial sampling interval of 0.25 m and the width of the probe pulse. With respect to this fiber end, the fiber cable length (*Df*) that corresponds to total casing length of 1082.56 m is 1095.599 m. Therefore, in our case, a depth mismatch (ΔD) of 13.039 m was identified, which is 1.2% longer than the total casing length.

A simple way to match depth discrepancies is to redistribute the amount of depth mismatch to the total length. Thus, we redistributed total depth error throughout the total length under the assumption that the rate of depth mismatch is constant throughout the total length.

For the confirmation of depth correction accuracy, we compared the temperature profiles before and after correction for gas lift process in May and November, 2017 as shown in Figure 4. Since the casing material is changed from steel to FRP casing at 791.90 m as shown in Figure 2, any temperature anomaly associated with the heat conductivity difference of two casing material should be identified around this depth. The effect of casing transition is clearly identifiable with the temperature increase pattern from data in both May and November, regardless of depth correction.

**Figure 4.** Comparison of depth correction effect with gas lift data from May (**a**) and November (**b**).

After correcting the depth, the depth of thermal changes due to casing transition moved up 9.494 m from 801.430 m to 791.936 m. The depth after correction, 791.936 m matches well with the true depth of the casing transition point (791.90 m) within one sampling-point difference, which is 0.25 m in terms of the DTS system's sampling interval. In both May and November, the difference to the true depth was 0.036 m.

Considering the true depth of casing transition point is known and the position of the thermal anomaly due to casing transition can be clearly identifiable, we may redistribute total depth error to each section above and below the transition point section by section under the assumption that the rate of depth mismatch is constant but different in each section. However, the recalculated sampling interval turns out to be the same to the third decimal place regardless of the section division. Thus, no major improvement in depth calibration using sections was seen.

The new sampling interval was applied to the cement injection data for the verification of the depth correction effect again as shown in Figure 5. Although there is no known absolute position of the top of the gravel pack and the bottom of cementing, the top of the gravel pack should be near 911 m and bottom of cement injection near 910 m according to the completion schematic shown in Figure 2.

**Figure 5.** Temperature profiles before and after applying depth correction to data at the end of cement injection in May.

Once the cement slurry meets the gravel pack, it will permeate slowly into the gravel pack. This process can be delineated by a temperature transition between the low temperature of the cooler cement slurry and the high temperature caused by the warmer gravel pack as shown in Figure 5, and this thermal transition remains until the temperature goes back to the normal geothermal gradient.

Considering this, we can interpret the start and the end of the thermal transition zone to the top of gravel pack and the bottom of cement slurry, respectively. The top numbers and the bottom numbers in the figure indicate the depth of the bottom of cement and the depth of the top of gravel pack, respectively. After correcting the depth, the top of gravel pack moved up 10.934 m from 922.937 m to 912.003 m that is 2.003 m difference with expected depth of 910 m. The bottom of cement moved up 10.949 m from 924.208 m to 913.259 m that is 2.259 m with expected depth of 911 m.

#### *3.2. Temperature Correction*

Many DTS acquisition systems have an internal calibration system, however higher accuracy is needed for geoscience operating environments where the fiber optic cable may be under conditions of rapid change and large temperature fluctuations. Under static conditions, the accuracy of calibrated temperatures is generally within ± 0.3 ◦C, however under rapid operating conditions, the accuracy range becomes widened to ± 2 ◦C [6]. The temperature measurement accuracy also declines as the length of the fiber is increased due to the power loss of the laser pulse [3]. In our experiment, pre-calibration could not be done before installation, therefore we suggest a post-calibration method that is within the acceptable range of error, ± 0.3 ◦C.

Three calibration methods were investigated in the JG-M borehole in April 2018 and the experimental setting is displayed in Figure 6. In a dual-ended configuration, only one absolute temperature point is needed because knowing the temperature at a point results in having two points in this configuration. The "water tub" is the conventional method of correction by submerging about one or more pulse width lengths of fiber cable into water to find the stable absolute temperature. However, submerging enough length of fiber longer than a pulse width length into water may not be possible once the cable is installed. Therefore, two pointwise temperature correction methods are investigated. The 'heat coil' heats a portion of fiber cable near the well with a heat coil. The "in-well" is another pointwise method of correction by using the position of a heat coil to apply a slope correction and using the absolute temperature value inside the well to apply the offset correction. For all cases, the slope correction was applied to temperature profile with no events in the well, and background temperature data was collected for 15 min before the calibration method experiments.

**Figure 6.** Schematic diagram and setting of three temperature calibration experiments.

#### 3.2.1. Long Point Calibration Method with Mater Tub

In the water tub experiment, 22 m of fiber optic cable, which is longer than a pulse width length of a laser, was submerged in water. Figure 7 shows the process of the temperature being stabilized in the calibration bath. The temperature fluctuation became smaller as time went by, and about 1 h later, finally stabilized to that of water inside the calibration bath.

**Figure 7.** Temperature change with time at the fiber optic cable submerged in water.

The position of the cable submerged in water could be clearly seen both on the left and right side of the DTS data, due to the dual-ended installation of the fiber. The stable range on the left side of the data ranges from 48.896 m to 64.148 m, and the stable range on the right side of the data, though it is not displayed in the figure, ranges from 2329.826 m to 2345.078 m. Since the submerged part of the fiber optic cable must have the same temperature, the inner 60 points were averaged to eliminate fluctuation of the temperature value on each side. The midpoint on the left side was at 56.522 m, and midpoint on the right side was at 2337.706 m. To make sure the right midpoint matched the left midpoint, same number of 4887 points to max peak were used.

In order to find a better slope correction factor, the temperature change of water within an hour was considered. Even though the specific heat of water is lower than air, the minimum difference of average temperature between the left side and the right side was considered for better estimation of the reference temperature. The ten-minute range which gave the minimum difference between the left and the right side is from 17:53 to 18:02. From the best ten-minute range, the minimum, maximum, and average slope value was used to find the slope equation. The maximum slope gave the minimum difference between the left and the right value. The minimum slope gave the maximum difference between the left and the right value. The slope curve of minimum, maximum and average slope is shown in Figure 8a. In the figure, "all" represents using all the data points without considering the temperature change of water, "DTS" represents the best ten-minute data in terms of minimum difference between the left and the right endpoints. The "Diver" represents the most stable temperature measured with the diver in the water tub for ten minutes.

**Figure 8.** Comparison of long point water tub methods in terms of temperature difference (**a**) and slope range (**b**).

The "DTS" method gave the best result when 60 points were averaged to eliminate fluctuations of the reference temperature as shown in Figure 9a. This was confirmed again by comparing all three methods as displayed in Figure 8. The averaging of 60 points was the most effective in all three cases as it had smaller slope and difference range.

**Figure 9.** (**a**) Slope curve of long point method with best 10-minute range of DTS 60-point average from 17:53 to 18:02 and (**b**) heat profile after long point water tub temperature correction.

For all three cases, the slope equation was applied to the background data from April 12 before any experiment because it is best to apply to the data without any events. The temperature accuracy improved after temperature correction, as the difference to true temperature decreased from 2.6 ◦C before correction to 0.02 ◦C, as illustrated in Figure 9b.

#### 3.2.2. Pointwise Calibration Method with Heat Coil

The heat-coil was used to locate the position on the DTS data for slope correction. A heat coil of 0.5 m in length was wrapped around the fiber optic cable near the well and insulated with a sponge tube of 1.5 m and wrapped in tape to prevent the temperature from fluctuating due to the wind as shown in Figure 6. The distance from the middle of heat coil to the well was 70 cm. In the heat coil experiment, obtaining a stabilized absolute temperature was difficult because the heat coil was shorter than the pulse width of interrogation system and air temperature varied too quickly. Therefore, the ambient temperature value at the heat coil position was used for slope correction. On April 12, the fiber optic cable was heated with heat coil to about 50 ◦C for about 10 minutes and the cooling observed for 27 min.

The slope curves are displayed in Figure 10a, which shows very different slope difference from the long point water tub method. In this case, only slope correction was applied because the absolute temperature was not determined.

**Figure 10.** (**a**) Slope curve of pointwise method with heat coil and (**b**) heat profile after pointwise heat coil slope correction.

Even with only slope correction, we could evaluate if the calibration is done correctly by comparing the two endpoints because two endpoints on the left and the right should have the same temperature value after correction. However, we can identify in Figure 10b, that temperature at the right endpoint, even with an increase of 1.1 ◦C after correction, still has a 1.5 ◦C difference with that of the left endpoint.

#### 3.2.3. In-Well Pointwise Calibration Method

To compensate the difficulty of obtaining a stable absolute temperature with the heat coil method, three divers were lowered into the well at 10 m, 20 m, 30 m below groundwater which has lower specific heat constant than the air. The diver position on the DTS could be located by counting from the known heat-coil position. The position could also be confirmed by temperature change when the fiber optic meets groundwater, which was only 1 m to in-well water from the heat-coil position.

The slopes found with the pointwise in-well method are presented in Figure 11. Compared to the long point water tub method, the variation ranges are wider than the long point method but the average slopes are similar to that of the long point method. The slope and its range did not vary greatly with the depth of diver. The temperature after correction with pointwise in-well diver at 10 m is presented in Figure 12. The temperature accuracy improved as the difference to true temperature is only 0.06 ◦C compared to the difference of 2.8 ◦C before correction.

**Figure 11.** Slope curves of pointwise calibration method by in-well diver at (**a**) 10 m, (**b**) 20 m, (**c**) 30 m.

**Figure 12.** Temperature profile after correction by pointwise method with 10 m in-well diver.

#### 3.2.4. Discussion

The slope and temperature difference curves of three calibration methods are summarized in Figure 13. The long point calibration method with water tub gave the best result. Although the pointwise calibration method using the diver showed a larger range than the long point method, it gave an average slope and the average temperature difference comparable with that of the long point method. Comparing the variation range of the heat coil method and the in-well method, we can conclude that it is crucial to obtain accurate and stable absolute temperatures.

**Figure 13.** Slope curves (**a**) and temperature difference (**b**) of post-calibration method.

The stability analysis was carried by dividing the total round-trip length of the fiber into 10 sections in terms of the same number of sample points per section. The downward trip sections, from wellhead to well bottom, are marked as L5–L1 and the upward trip sections are marked as R1–R5. The number of samples per each section is 897 points for the long point method, 855 points for the pointwise method with in-well diver, and 862 points the pointwise method with heat coil.

The stability curve of all three calibration methods is as shown in Figure 14, which displays the stabilities at 13 points including two top points (Endpt), a bottom point (Maxpeak), and ten midpoints of each section. The slope stability curve of long point calibration method has a range of +0.06 ◦C and −0.03 ◦C. The slope stability curve of pointwise calibration method with in-well of 10 m has a range of +0.35 ◦C and −0.25 ◦C. The slope stability curve of pointwise calibration method with heat coil has a range of +0.7 ◦C and −0.6 ◦C. The temperature difference increases as the length of fiber increase in all three calibration methods as expected due to the power loss that increases as the traveling distance become longer.

**Figure 14.** Slope stability of (**a**) long point calibration with water tub (**b**) pointwise calibration with in-well and (**c**) pointwise calibration with heat coil.

The appropriate temperature error range is within 0.3 ◦C, and both water tub and in-well data fall within the range. However, the error range with water tub slope correction is within 0.3 ◦C with stability of 0.05 ◦C and the temperature error range with diver slope correction is within 0.1 ◦C with a stability of 0.3 ◦C. Therefore, even though the long point calibration works best, the pointwise calibration method can be a useful alternative where long point calibration cannot be done.

#### **4. DTS Monitoring of Well Completion Processes**

The uniqueness of the Janggi field DTS data is the continuous real-time monitoring of the whole well completion process from gravel packing to complete curing without well intervention. The heat map obtained from May 11th to June 4th is shown in Figure 15. Anomalies numbered from 1 to 4 correspond to gravel packing, cement injection, gas lift, and cement curing, respectively.

**Figure 15.** Heat map obtained during well completion where the gravel packing is denoted by 1, the cement injection by 2, the gas lift operation by 3 and the cement curing process by 4.

#### *4.1. Gravel Packing*

There has been monitoring of production in gravel packed completions but monitoring the gravel packing process itself was not yet done [7]. At the Janggi field site, fiber optic cable was connected directly after casing installation to the DTS unit, which allowed real-time monitoring of the gravel packing process throughout the entire length of the well. The detailed temperature profile during gravel packing can be seen in Figure 16 which outlines the major events during the gravel packing process.

The reference (black line) refers to the temperature profile at the end of well completion. The gravel packing process can be characterized by the decrease of temperature below the reference line. The temperature decreases because of the cooling effect due to the injection of sand mixed with fluid that is cooler than the borehole fluid inside the well. In the deeper parts of the well, high temperature anomalies with oscillating pattern can be seen due to water circulating but leaving sand to be accumulated. The sand deposited in the annulus plays a role of thermal bridge between the surrounding formations and fiber cable on tubing. Because sand has higher heat conductivity than water, such a thermal bridging effect explains the high temperature anomalies. The section where sand is not accumulated does not have this connectivity between formations and cable still shows a low

temperature trend. The point where the temperature clearly drops indicate the top level of the gravel pack. The progression up of this anomaly along the well shown in the figure (May 11 19:39 and May 12 14:51) clearly shows that the top of gravel pack can be monitored well by DTS. After gravel packing, the top of gravel pack estimated from DTS data (May 12) matches closely with the planned depth of 910 m as mentioned in Section 3.2.1.

**Figure 16.** Schematic of gravel packing operation (**a**) and temperature profile during the gravel packing operation (**b**).

In the low temperature anomaly section, small peaks that seem like noise are identified. However, they are anomalies due to the clamps being spaced at approximately every 9 m. The clamps provide more direct connection between inside and outside of tubing wall. The effect of clamps appears as high temperature anomalies in the unpacked section while it appears as low temperature anomalies in the gravel-packing zone. In the case of the high temperature anomalies, inside the casing is warmer than the annulus because the sand was directly injected into the annulus. Consequently, the clamped point of the DTS cable is directly connected to the warmer environment while the other parts are not. In the case of low temperature anomalies, inside the casing is cooler in the packed region. Thus, the clamping point has the reverse effect compared to the previous case. The anomalies associated with clamps indirectly shows high resolution of DTS and its usefulness in dynamic condition.

From the heat map in Figure 15, we can identify that two gravel packing processes that were carried out over the two days have similar patterns to each other. The temperature is increasing above the gravel-packed zone and decreasing in the unpacked zone during the operation. The short burst of cooling at the top may be explained due to either flow change or colder fluid injection or breaks during gravel packing operation.

The ability to monitor the progression of sand volume allows the possibility to control the amount of sand volume to reach the top of the intended depth and the capability to check for possible plugging as the sand is poured. All these factors help optimize the gravel packing process at the field level as the change in temperature is displayed without interrupting the process. The gravel packing process also showed potential for detecting sand production and the response to sand production can be optimized to diminish the risk of reduced well productivity.

#### *4.2. Cement Injection and Curing Process*

After the gravel packing operation, the cement slurry was injected to the casing to complete the well. The temperature profile of the cement injection process is displayed in Figure 17. The temperature before cement injection has a lower temperature than the reference, which indicates the effect from gravel packing operation still remains and needs more time to recover to an undisturbed normal geothermal gradient. As the slurry is injected, the temperature is being gradually lowered until the end of the injection because the cement slurry is cooler than the borehole fluid. The top of the gravel pack and the bottom of cement can be determined as mentioned in Section 3.2.1.

**Figure 17.** Schematics of cementing injection process (**a**) and temperature profiles during the cement injection process (**b**).

The major events during curing process from right after injection to a month after injection are shown in Figure 18.

**Figure 18.** Temperature profiles during the curing process in detail.

In the yellow line denoted as type 1, though the high temperature anomaly due to gas lifting remains at depths shallower than 850 m, we can identify the peak near the stage cement tool which can be interpreted as the start of the exothermic process of the curing. The effect of the gas lift anomaly is gone in type 2 blue line; however, another peak appears near 450 m. The peak indicates the start of curing from the upper part. The type 3 green line displays two peaks from the top and bottom boundary. These boundaries move toward each other, indicating the curing process propagates from the boundary to the intermediate. The type 4 red line displays the two peaks finally meet and start settling down.

From the heat map in Figure 15, the major trends of the whole cementing operation can be observed. We can identify when and where the curing starts, when the curing finally done, and how long it took for complete curing from the heat map provide by constant DTS monitoring.

The cementing is the final critical steps to ensure the integrity of the well. Therefore, monitoring the entire process from cement slurry injection to curing in real-time can provide key information in order to control and optimizes the process. For example, the progression of injection outlines both injection rate and volume of cement, which offers the ability to check for possible plugging during injection. The constant monitoring of the entire well also allows correlating unusual anomaly to exothermic reaction due to the settling of the cement, which can provide better insights to the different stages of exothermic reactions. As the time of complete curing could be accurately determined by the homogeneity in temperature, monitoring curing process help move on to the next stage with confidence.

#### *4.3. Gas Lift*

The gas lift process was carried out after cement injection but during cement curing. During the gas lift, the N2 gas was released, and the fluid was pulled upwards due to pressure. The major events of one cycle of gas lift operation is displayed in Figure 19. The temperature increased during gas lift operation in the data below 791.90 m. The increase of temperature is due to pulling up of the deeper and warmer reservoir fluid. The high temperature anomaly associated with gas lift operation has a highly oscillating pattern. In our case, however, its amplitude is different at two regions above and below 791.90 m. The difference in amplitude of the oscillation is due to heat conductivity difference of the casing materials. The 791.90 m is the transition point of the two casing materials from steel to FRP. Since steel has higher heat conductivity than FRP, we can observe a higher amplitude of oscillation in the steel casing region than the FRP casing region.

**Figure 19.** Schematics of gas lift process (**a**) and temperature profiles during the gas lift in May (**b**).

The heat map during the gas lift cycles which shows the global trend during the gas lift is displayed in Figure 20. The number of repeated patterns during gas lift matches the number of gas lift cycles, which was eight cycles. The width of each repeated pattern is the duration of each cycle. Even after the gas lift operation was finished, its effect lasted a day to return to temperature before the gas lift operation began. The temperature did not completely recover due to the starting of curing process near 900 m during the middle of gas lift operation.

**Figure 20.** Heat map of entire gas lift operation in May

As mentioned in Section 3.1, our gas lift data shows not only a typical highly oscillating thermal anomaly pattern of gas lifting event but also the effect of casing transition, which supports the confirmation of the accuracy of depth correction. The DTS data has the potential to detect multiple points of injection and determine gas rate through each injection point. All these are possible as the event is happening and require no shut-in when measuring DTS temperature data [8]. Therefore, DTS helps capture the temperature response throughout the entire length of the well simultaneously as the event is happening.

#### **5. Conclusions**

We have presented continuous real-time DTS monitoring data of the whole well completion process from gravel packing to complete curing without well intervention, and discussed the observed thermal anomaly patterns during each well completion process.

The real-time DTS monitoring data revealed that each well completion process can be characterized by its own distinctive temperature anomaly patterns. In addition, it is remarkable that the upward progression of sand and cement slurry and the thermal anomaly due to curing of cement can be monitored. This clearly shows that pre-calibrated DTS, due to its continuous real-time and high spatial resolution nature compared to conventional point sensors and wireline tools, can be a powerful tool not only for monitoring, but also for controlling and optimizing well completion processes without intervention in each process. However, because the thermal event characteristics associated with each well completion process are site- and condition-dependent, further studies to quantify the temperature anomaly depending on factors like injection rate, hydraulic and thermal properties of the formation and fluid etc. are required for quantitative control and optimization of well completion process by DTS.

For depth calibration, the conventional method distributing the error to the total length is enough. The accuracy could be enhanced more if more points where the depth are known are available in addition to the top and the bottom of well. In the case of temperature calibration, pre-calibration using the conventional method is best. However, our investigation shows that the pointwise method can be an alternative with acceptable accuracy of less than 0.3 ◦C and stability less than 0.3 ◦C.

**Author Contributions:** Data curation, D.S.L.; Formal analysis, D.S.L., S.-J.C., K.G.P. and C.L.; Funding acquisition, K.G.P.; Investigation, D.S.L., S.-J.C. and K.G.P.; Methodology, K.G.P.; Supervision, K.G.P.; Validation, C.L. and K.G.P.; Writing – original draft, D.S.L.; Writing – review & editing, K.G.P.

**Funding:** This work was supported by the Korea CCS R&D Center (Korea CCS 2020 Project) grant funded by the Korea government (Ministry of Science and ICT) in 2017 (KCRC-2014M1A8A1049287).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Research on the Earth Pressure and Internal Force of a High-Fill Open-Cut Tunnel Using a Bilayer Lining Design: A Field Test Using an FBG Automatic Data Acquisition System**

#### **Tianyuan Xu \*, Mingnian Wang, Li Yu \*, Cheng Lv, Yucang Dong and Yuan Tian**

Key Laboratory of Transportation Tunnel Engineering, Ministry of Education, Southwest Jiaotong University, Chengdu 610031, Sichuan, China<EMAIL_ADDRESS>(M.W<EMAIL_ADDRESS>(C.L.);

<EMAIL_ADDRESS>(Y.D<EMAIL_ADDRESS>(Y.T.)

**\*** Correspondence<EMAIL_ADDRESS>(T.X<EMAIL_ADDRESS>(L.Y.)

Received: 30 January 2019; Accepted: 20 March 2019; Published: 27 March 2019

**Abstract:** When there are railway tunnels on both sides of a valley, a bridge is usually built to let trains pass. However, if the valley is very close to an urban area, building an open-cut tunnel at the portal and then backfilling it to create available land resources for the city and to prevent excavation slag from polluting the environment would be a wise choice. This has led to the emergence of a new type of structure, namely, the high-fill open-cut tunnel. In this paper, by performing an automatic long-term field test on the first high-fill open-cut tunnel using a bilayer design in China, the variations of earth pressure and structural internal force during the backfilling process were obtained, and different tunnel foundation types were studied. The results showed that the earth pressure significantly exceeded the soil column weight, with a maximum earth pressure coefficient between 1.341 and 2.278. During the backfilling process, the earth pressure coefficient increased at first and then decreased slowly to a relatively stable value, and a stiffer foundation would make the structure bear higher earth pressure (1.69 times the normal one observed during monitoring). The change of internal force had two stages during backfilling: before the backfill soil reached the arch crown, the internal force of the lining changed slowly and then grew linearly as the backfill process continued. Moreover, the axial force ratio of the inner and outer linings was close to their thickness proportion, and the interaction mode between the two layers was very similar to the composite beam.

**Keywords:** high-fill open-cut tunnel; bilayer lining; FBG sensors; automatic monitoring system; earth pressure of open-cut tunnel

#### **1. Introduction**

In China, the high-speed railway network is rapidly extending into the mountainous area of the southwest. Due to the geological conditions of the southwest region, there has been a substantial increase in the number of tunnels needed, some of which are near cities. There is a special situation for some of these cities: sometimes, there is a valley between two tunnels that is very close to urban areas. In a normal situation, a bridge must be built above the valley to let trains pass. However, if an open-cut tunnel were built at the portal and then backfilled with the engineering spoils from the tunnel excavation, available land resources can be created for cities, which would prevent waste from polluting the environment and lead to the emergence of high-fill open-cut tunnels, as shown in Figure 1. Compared with the traditional open-cut tunnels, a high-fill open-cut tunnel has the following two characteristics:

(1) The height of backfill layer can reach up to 30–40 m, which is almost five times larger than a common one;

(2) In order to bear such high earth pressure, the thickness of the tunnel lining can reach up to 2 m or more, which is 4–5 times thicker than those with shallow overburden. Details of existing high-fill open-cut tunnels in China are shown in Table 1.

**Figure 1.** Open-cut tunnel built in a valley close to a city.

**Table 1.** Details of three existing open-cut tunnels in China.

These two characteristics present two research difficulties:

(1) A large amount of hydration heat can be generated after construction if a monolithic design is adopted for these superthick concrete linings, causing a significant temperature difference between the internal and external surfaces of the lining, which would result in structural shrinkage cracks and safety issues. Therefore, a more reasonable structure design must be adopted.

In order for the lining structure to meet the bearing capacity requirement for such a high backfill layer and to minimize the shrinkage cracks at the same time, a bilayer design is used for high-fill open-cut tunnels which divides the overthick lining into two parts—outer and inner linings. Two layers are concreted separately at 14-day intervals, resulting in a significant decrease of hydration heat generated by each layer of lining. However, the bilayer lining design initially appeared in shield tunnels, which consist of segmental and secondary linings, and has never been used in open-cut tunnels.

Yang et al. [1] carried out a numerical study on the performance of segmental and secondary linings as well as the stress transmission between them. They reported that linings with rebar had a combined bearing capacity, while linings with membranes had separate bearing capacities. Vogel et al. [2] examined double-shell linings in respect to direct shear stress capacity. Their results revealed that a spray-applied waterproofing membrane was able to transfer stresses between both concrete linings but could not transfer any shear stress. Su and Bloodworth [3] presented laboratory

tests on beam samples cut from bilayer linings. They compared the influence of different membrane thicknesses and substrate roughness and provided parameters based on test results for further research and design. Their subsequent research [4] developed a composite mechanical behavior quantification method. Chuan et al. [5] presented several load tests on bilayer shield tunnel models, the results of which revealed that the secondary lining bore most of the bending moment. Their following study focused on a calculation model for bilayer linings [6].

The abovementioned studies provide a better understanding of the mechanical characteristics of the bilayer lining tunnel design, as well as the calculation model and load shearing ratio needed to design the structure. However, the construction sequence of open-cut tunnels is completely the opposite of that of shield tunnels, in which the lining is built at last. Moreover, the size and shape of open-cut tunnels are also quite different from shield tunnels (circular design), as shown in Figure 2.

**Figure 2.** Comparison of size and shape of open-cut and shield tunnels using the bilayer design.

(2) Under such a high backfill layer, it still remains unclear whether the earth pressure on the tunnel crown equals the overburden pressure. Moreover, a deep foundation has been adopted for high-fill structures in some situations, which can also affect the earth pressure of the structure.

A high-fill situation usually occurs in culverts and pipelines in mountainous areas and has never been adopted in railway open-cut tunnels. Qiang et al. [7] analyzed several calculation methods of earth pressure on a culvert crown. The results showed that under a high embankment fill condition, the calculated results by the "neutral point" method were closer to the values measured in situ. Li et al. [8] carried out a field monitoring and numerical simulation to investigate the state and distribution of stress on the exterior surface of slab-culverts under high embankments. The results showed that the earth pressure on culverts was larger than the self-weight of filling, and the pressure distribution was uneven. Chen et al. [9] presented a new formula to calculate the vertical earth pressure on a culvert, the calculation of which results were compared with field tests. Zhang et al. [10] carried out a centrifugal experiment on the earth pressure distribution of a culvert top. The results showed that a box culvert with a pile foundation bore higher earth pressure than a normal one, and they suggested that soil could be backfilled first on both sides of the culvert before it was constructed. Some researchers developed a load calculation method for deeply buried culverts by summarizing the displacement distribution characteristics on the culvert crown [11,12]. Additionally, controlling displacement was considered to reduce the vertical pressure on culverts. Moreover, backfill materials and some engineering measures also affect the earth pressure of high-fill structures. Meguid [13] presented an experimental investigation to measure the earth pressure distribution on a rigid pipe backfilled with tire-derived aggregate. The average measured earth pressure above the crown of the pipe was found to be as low as 30% of the overburden pressure. An optimum soft zone geometry for imperfect trench installation was proposed to maximize the reduction of the earth pressure on buried corrugated steel pipes [14], and the maximum wall stress was reduced by 69%.

According to previous research, the earth pressure of an open-cut structure is closely related to the side slope angle and foundation type. As a more rigid foundation usually results in a higher earth pressure, it is suggested to avoid using a deep foundation in the design of culverts. However, some culverts still use a strengthened foundation to achieve better performance in particular situations [15]. For railway open-cut tunnels, the track surface subsidence must be strictly controlled to meet the operation requirements of high-speed trains, making it necessary to use deep foundations for railway open-cut tunnels in poor geological conditions. In conclusion, a high-fill open-cut tunnel using a bilayer lining design is very different from the current culvert or shield tunnel in terms of construction method, shape, size, and foundation type, and its mechanical characteristics are not clear yet. Besides, current research rarely contains long-term field measurement data, especially measurements conducted when the structure is built on different types of foundations.

In geotechnical and structural domains, fiber Bragg grating (FBG) sensors are applied to measure structural strain [16–19], seepage pressure [20,21], temperature [22], and vibration [23]. However, the monitoring work for tunnels does not continue after the lining has stabilized. Therefore, vibrating wire sensors are more widely used to reduce unnecessary costs, and these sensors typically require manual data acquisition. However, for dynamic or long-term data acquisition from tunnels, such as the real-time safety monitoring of tunnels in special geological conditions and measuring the vibrations generated by train travel [24,25], it is more sensible to use FBG sensors, as they have strong anti-interference and long-term stability advantages. In our research, FBG pressure cells and strain gauges were adopted to acquire the soil pressure and structural internal force for further analysis.

Moreover, for an open-cut tunnel, once the backfilling process begins, personnel have to enter the tunnel from the undercut portal, which sometimes is more than 10 km from the test section. Therefore, an automatic data acquisition system was adopted to provide more timely and accurate data. These test results can be a useful complement to both high-fill open-cut and bilayer structures to help provide a deeper understanding of them.

#### **2. Field Test Procedure**

#### *2.1. Tunnel Description*

The field test was carried out at the first high-fill open-cut tunnel in Fengdu, China, where the bilayer lining design was firstly adopted, with a 0.5-m-thick C35 concrete inner lining covered by a 1.4–2.8-m-thick C35 concrete outer lining. Further, a waterproof layer consisting of nonwovens and a polyethylene sheet was applied between the inner and outer linings as well as on the surface of the outer lining. The outer lining was constructed first. When the outer lining was stabilized, the inner lining was constructed, and finally, the backfilling process was carried out.

This open-cut tunnel is 373-m long with a backfill height up to 22–28 m (from tunnel arch crown to the backfill surface) and a 100-m horizontal backfill range on both sides of the tunnel axis. Both sides of the tunnel were filled symmetrically with 0.5-m-thick earth in each cycle, which took 6–7 days, and the whole backfill process lasted about 9 months. Also, a C30 concrete dam was adopted as the foundation for some part of the open-cut tunnel to control the subsidence of the tunnel bottom. The remaining parts of the open-cut tunnel were set on bedrock with weatherproof protection. More design details are shown in Figures 3–5.

#### *2.2. Automatic Data Acquisition System*

The application of intelligent test systems and methods in tunnels has been a growing trend in recent years [26–29]. In our field test, an automatic data acquisition system was developed and adopted, using FBG sensors to measure the structural internal force and soil pressure. The system was developed by the research team in collaboration with Sensorlead Technology Co., Ltd, Shanghai, China. It can capture the wavelength data of the sensors at 100 Hz and allows users to input different formulas to convert wavelength data into stress, pressure, displacement, etc. Moreover, it has a remote operating subsystem to allow users to view the data anytime by personal computers or even smartphones.

**Figure 3.** Design details of the tunnel. (**a**) Cross section of the tunnel and the concrete dam foundation; (**b**) Profile view of the whole project.

**Figure 4.** Backfilling range (half side).

**Figure 5.** Open-cut tunnel in the backfill process.

Step 1: Adjacent sensors were connected by optical fibers at both ends and were finally connected to a fiber optic cable closure, as shown in Figure 6.

Step 2: All the closures were connected in a series and were finally connected to one main fiber cable.

Step 3: The main fiber cable was connected to FBG demodulators (Figure 7), which transmitted the data to a remote-control computer, as shown in Figures 8 and 9.

**Figure 6.** Fiber optic cable closure.

**Figure 7.** FBG demodulators.

**Figure 8.** Remote-control computer.

Figure 9 shows the flow chart of the whole data acquisition system. In order to reduce the test error caused by construction disturbance, the data acquired at 4:00 a.m. was used for subsequent analysis.

**Figure 9.** Automatic data collecting and delivering system.

#### *2.3. Sensor Layout and Basic Structure*

Two sections were chosen to test the earth pressure and structural internal force, using bedrock and the concrete dam as their foundation, respectively, as shown in Table 2.

In order to obtain the variation law of earth pressure and structural internal force during the backfilling process, earth pressure cells were installed at the arch crown and arch rib on the outer lining surface, and strain gauges were installed in the inner and outer linings, respectively, as shown in Figure 10. For those pressure cells installed on the outer lining surface, metal cable boxes were used to protect the fiber cables from being damaged by the falling backfill soil, as shown in Figure 11.

**Figure 10.** The layout of pressure cells and strain gauges.

**Figure 11.** Pressure cells on the outer lining surface with metal cable boxes to protect the optical cables.

The gauges and cells adopted in this research were made by Sensorlead Technology Co., Ltd. Shanghai, China. Their model and basic structure are shown in Figures 12 and 13. Changes in stress and pressure can cause changes in the wavelength of the sensor, which could be measured by the gratings. Further, real-time ambient temperature could be read by the temperature compensation gratings in the sensors, correcting errors due to seasonal temperature differences and providing more accurate data.

**Figure 12.** Model and basic structure of pressure cells.

**Figure 13.** Model and basic structure of gauges.

#### **3. Numerical Model**

For further study and design improvement of a high-fill open-cut tunnel using bilayer linings, a numerical model was built in this study, as shown in Figure 14. The elastic model was adopted for the lining and concrete foundation of the open-cut tunnel, and the Drucker–Prager elastic–plastic model was used for the bedrock and backfill soil. The contact elements were added between the open-cut tunnel and the backfill to simulate the mutual squeezing and slip between them (Figure 15). According to the site situation, a waterproof layer consisting of nonwovens and a polyethylene sheet was set between the outer lining and backfill soil as well as between the outer and inner linings. The friction coefficient between the nonwovens and the polyethylene sheet was 0.23. Moreover, the backfilling process of the open-cut tunnel was simulated by the method of layered activation, that is, finite element units for every single backfilling layer (0.5 m) were activated sequentially. The calculation parameters of the model are shown in Table 3.

**Figure 14.** Three-dimensional finite element model for a numerical simulation from ANSYS 12.0 software.

**Figure 15.** Contact elements between different structures.

In the 3-D numerical model, the node stress in the original Cartesian coordinate system was transformed into the tangential stress of the lining by the coordinate transformation of elastic mechanics. As a result, the bending moment and axial force of the lining were obtained and used for further analysis. The transformation method is described below (Figure 16).

**Figure 16.** Coordinate transformation for node stress of lining unit.

Taking nodes 1 and 2 on the outer and inner sides of a section, and assuming that the angle between the line of the two nodes above and the vertical plane is *θ*, the tangential stress of each node on the lining section can be calculated from Equation (1) [30]:

$$
\sigma\_{\text{ll}} = \sigma\_{\text{x}} \cos^2 \theta + \sigma\_{\text{y}} \sin^2 \theta + \tau\_{\text{xy}} \sin 2\theta \tag{1}
$$

where *σx*, *σy*, and *τxy* are node stress components in the original coordinate system, respectively, and *θ* is the angle between the outer normal of the section and the *y* axis.

Therefore, the tangential stress of nodes 1 and 2 was obtained, and then the axial force *N* and bending moment *M* of the lining section could be deduced according to Equations (2) and (3) [30]:

$$M = bh^2 \frac{\sigma\_1 - \sigma\_2}{12} \tag{2}$$

$$N = bh \frac{\sigma\_1 + \sigma\_2}{2} \tag{3}$$

#### **4. Earth Pressure Test Results and Discussion**

#### *4.1. Earth Pressure Test Results*

The in situ earth pressure was obtained in the backfilling process of the open-cut tunnel. The measured results are represented by the earth pressure coefficient (*c*), which was defined as earth pressure divided by bulk density, as described in Equation (4). The average earth pressure of the left and right arch rib acted to reduce the testing error. The testing results are shown in Figure 17:

$$
\mathfrak{c} = \frac{p}{\rho g h} \tag{4}
$$

where *p* is the earth pressure, *ρ* is the soil density, and *h* is the backfilling height.

**Figure 17.** Variation of earth pressure coefficient during the backfilling process. (**a**) Arch crown; (**b**) Arch rib.

As shown in Figure 17:

(1) The earth pressure coefficient firstly increased and then decreased slowly with the increasing backfilling height. The earth pressure coefficient of sample A1 reached its peak value when the backfilling height was around 8 m and then decreased slowly until stabilizing. However, the maximum coefficient of sample A2 appeared at the stage when the backfilling height was about 14 m and then rapidly decreased until stabilizing. The different change trends of the earth pressure coefficient between samples A1 and A2 means that the interaction process in the soil was gradual during the backfilling phase.

(2) At same backfilling height, the earth pressure coefficient at the arch crown and arch rib of sample A2 was obviously higher than that of sample A1, respectively, as listed in Table 4. The earth pressure coefficient at the arch crown of sample A2 was 1.53 times that of sample A1, and sample A2 was 1.69 times that of sample A1 at the arch rib. The root cause of this phenomenon was the additional shear stress difference between samples A1 and A2. The additional shear stress was caused by the settlement difference between the inside and outside soil columns, as shown in Figure 18. The inside soil column consisted of lining and the above soil. Its stiffness was much higher than that of the outside backfilled soil, which led to the settlement difference between the outside and inside soil columns. So, the outside soil column pulled the inside soil column down, resulting in additional shear stress, as shown in Figure 18. Therefore, compared with the stiffness of sample A1, that of sample A2 was much higher, resulting in greater differential settlement and, finally, a significant increase in earth pressure.

Arch rib 1.341 2.278 1.69

**Table 4.** The earth pressure coefficients of different samples.

In summary, at the same backfilled height, the increasing stiffness and thickness of the base caused there to be higher earth pressure acting on the open-cut tunnel, which increased the bearing capacity requirement for the lining. Nevertheless, rail surface settlement must be strictly controlled to satisfy the operation requirements of railway trains. Thus, a deep foundation is always applied in open-cut railway tunnels which cross through weak and soft ground. So, to avoid excessive loading of the lining, load shedding measures can be used which have already been adopted in high-fill culverts and pipes [31]. For instance, the load shedding layer was installed above the structure to decrease the settlement difference between the inside and outside soil columns, which reduced the loading of the open-cut structure [32].

**Figure 18.** Additional shear stress caused by the settlement difference between inside and outside soil columns.

#### *4.2. The Comparative Analysis of In Situ Testing and Numerical Simulation*

In order to validate the abovementioned conclusion, a comparative analysis of the in situ test and numerical simulation is shown in Figure 19.

**Figure 19.** Earth pressure coefficient (comparison of numerical simulations and in situ testing). (**a**) Arch crown; (**b**) Arch rib.

As shown in Figure 19, the earth pressure coefficient change trend of the numerical simulation was similar to that of the in situ test, which firstly increased and then decreased until stabilizing. Moreover, the vertical earth pressure coefficient of sample A2 was obviously higher than that of sample A1, as listed in Table 5. The earth pressure coefficient at the arch crown of sample A2 was 1.63 times that of sample A1, and sample A2 was 1.59 times that of sample A1 at the arch rib.

**Table 5.** The earth pressure coefficients of different samples in numerical simulations.

The in situ testing results were in good agreement with the numerical simulation results, and the maximum error was 21.7% at the beginning phase of the backfilling process. The reason for this could be that the backfill had not been completely compacted and stabilized during the actual construction.

In summary, both the in situ testing and the numerical simulation demonstrated that the total loading that acted on the high-fill open-cut tunnel consisted of two parts: the additional loading and the soil column loading. Moreover, the value of additional loading was closely related to the stiffness of the base. In order to more accurately determine the earth pressure of a high-fill open-cut tunnel, the lining and foundation below must be considered as a whole structure. The results also verified the accuracy of the numerical model, which can be used to analyze the settlement difference between the inside and outside soil columns, as shown in Section 4.3.

#### *4.3. The Settlement Difference Analysis between the Inside and Outside Soil Columns*

The key factor for the additional loading was the settlement difference between the inside and outside soil columns. However, the change trend or results of this difference could not be obtained exactly through in situ measurement, making it difficult to explain the gradual decreasing phenomenon of the earth pressure coefficient. Therefore, the relationship between earth pressure and settlement difference was analyzed via numerical simulation. The settlement difference between the inside and outside soil columns of two base types is illustrated in Figure 20.

**Figure 20.** Settlement difference between the inside and outside soil columns. (**a**) Settlement difference for sample A1; (**b**) Settlement difference for sample A2.

As presented in Figure 20, the settlement difference of sample A2 was obviously higher than that of sample A1. The settlement difference results from the arch crown in the 2–28-m range are shown in Figure 21, where the settlement difference on the horizontal plane near the arch crown for sample A1 was 18 mm, while that of sample A2 was 107 mm, about six times that of sample A1, which means the stiffness of the base had a significant influence on the settlement difference. Moreover, the settlement difference gradually decreased to a stable value with the increasing backfill height, which clearly explains the change trend of the earth pressure coefficient.

**Figure 21.** Settlement difference variation law as the backfill height increased.

#### **5. Internal Force Test Results and Discussion**

The structural strain data were firstly transformed into tangential stress by multiplying the elastic modulus of the material. The test data of the measuring points in the symmetrical position were averaged to reduce the error. The dynamic change of axial force and bending moment of the bilayer lining during the backfilling process is shown in Figures 22 and 23.

**Figure 22.** Variation of internal force during the backfilling process for sample A1. (**a**) Axial force of outer lining (sample A1); (**b**) Bending moment of outer lining (sample A1); (**c**) Axial force of inner lining (sample A1); (**d**) Bending moment of inner lining (sample A1).

**Figure 23.** Variation of internal force during the backfilling process for sample A2. (**a**) Axial force of outer lining (sample A2); (**b**) Bending moment of outer lining (sample A2); (**c**) Axial force of inner lining (sample A2); (**d**) Bending moment of inner lining (sample A2).

The change trend of axial force and bending moment can be represented as two stages:

Stage 1: Before the backfilling soil reached the arch crown, the structure only bore lateral pressure, resulting in a very small increase of the axial force and bending moment.

Stage 2: When the backfilling soil exceeded the arch crown, the lining bore both lateral and vertical pressure. The axial force showed linear growth, while the bending moment linearly decreased to a negative value. Moreover, the internal force change trend of samples A1 and A2 during the backfilling process was quite similar.

However, the absolute value of the internal force of the tunnel on the concrete dam foundation was much higher than the one without the concrete dam. The axial force of sample A1 was 1.2–1.6 times that of the other, and the bending moment was 1.6–2.4 times that of the other, as seen in Table 6. This further proves that an open-cut tunnel using a concrete dam foundation will bear greater earth pressure at the same backfilling height, and the internal force of the structure will also increase significantly.

In addition, the axial force of the outer lining was about 2.6–3.2 times that of the inner lining, which was very close to their thickness ratio (3:1), as shown in Table 6. This was because the waterproof layer transition between the two linings reduced the friction between the outer and inner linings. The transmission of shear force between the outer and inner linings was blocked, while only the axial force could be transmitted. That is, the mechanical behavior of the double lining could be analyzed according to the combined beam model.

**Table 6.** Relationship between the internal force and the tunnel foundation.

#### **6. Conclusions**

A field test study was carried out at the first high-fill open-cut tunnel using the bilayer design in China. An automatic data acquisition system using FBG sensors was adopted to obtain the earth pressure and internal force of the tunnel, and numerical models were created for further analysis as well. The following major conclusions can be drawn:

**Author Contributions:** Conceptualization, M.W.; Methodology, T.X. and L.Y.; Software, T.X. and C.L.; Validation, M.W. and T.X.; Investigation, T.X. and Y.D.; Data Curation, T.X. and C.L.; Writing-Original Draft Preparation, T.X.; Writing-Review & Editing, T.X., L.Y. and M.W.; Visualization, Y.D., T.Y. and T.X.; Supervision, M.W.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51878568, and the Chongqing Development Plan Project, grant number cstc2014yykfB30003.

**Acknowledgments:** The authors would like to thank Rui Xu, Xin Zheng, Hao Cai and Hanbo Chen for their help in the field test. The constructive comments from the reviewers are greatly appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Field Performance of Open-Ended Prestressed High-Strength Concrete Pipe Piles Jacked into Clay**

#### **Hai-Lei Kou 1,2,\*, Wen-Zhou Diao 1, Tao Liu 3, Dan-Liang Yang 1 and Suksun Horpibulsuk 4**

Received: 31 October 2018; Accepted: 29 November 2018; Published: 1 December 2018

**Abstract:** The behavior of open-ended pipe piles is different from that of closed-ended pipe piles due to the soil plugging effect. In this study, a series of field tests were conducted to investigate the behavior of open-ended prestressed high-strength concrete (PHC) pipe piles installed into clay. Two open-ended PHC pipe piles were instrumented with Fiber Bragg Grating (FBG) sensors and jacked into clay for subsequent static loading tests. Soil plug length of the test piles was continuously measured during installation, allowing for calculation of the incremental filling ratio. The recorded data in static loading test were reported and analyzed. The distribution of residual forces after installation and the effect on the bearing capacity were also discussed in detail. The test piles were observed to be in partially plugged condition during installation. The measured ultimate shaft resistance and total resistance of the test piles were 639 and 1180 kN, respectively. The residual forces locked in the test piles after installation significantly affected the evaluation of the axial forces, and thus the shaft and end resistances. It tended to underestimate the end resistances and overestimate the shaft resistances if the residual forces were not considered in the loading test. However, the residual forces did not affect the total bearing capacity of open-ended PHC pipe piles in this study.

**Keywords:** open-ended pipe piles; field test; static loading test; soil plugging; residual forces

#### **1. Introduction**

Prestressed high-strength concrete (PHC) pipe piles are commonly adopted as deep foundations in the coastal areas of Southeast China. The installation methods usually involve jacking and driving. A vibration and predrilling method is also utilized in the field. Compared with other installation methods, the jacking method is preferred as it avoids noise, vibration and slurry pollution. Hence, it is more suitable to install piles in urban areas. The performance of jacked pipe piles during and after installation has becomes a focus of attention in recent years [1,2], but many issues remain to be solved, in particular for open-ended pipe piles, such as the termination criteria, the bearing capacity, the design methods and so on [3,4]. This is mainly because the physical properties involved in open-ended pipe piles are extremely complicated [5,6].

Many references have reported on the behavior of jacked steel H-piles, including the effect of excess pore water pressure, the set-up, the estimation of shaft and end resistances and so on [7]. As full displacement piles, the steel H-piles have no soil plug effect during installation [8,9]. However,

for non-displacement or partial displacement piles, soil plugging was identified as an important factor in pile behavior [10,11]. A number of field and model tests were conducted to study the effect of soil plugging on the behaviour of open-ended piles jacked into sand [12]. Kou et al. [13,14] revealed that the loading-settlement behavior was closely related to the development of soil plugs. Yu and Yang [15] compared major open-ended pipe pile design methods and proposed the Hong Kong University (HKU) method to estimate the end bearing capacity of open-ended steel piles in sand based on the Imperial College Pile (ICP) and the University of Western Australia (UWA) methods [16,17]. The proposed HKU method was capable of producing satisfactory predictions over a wide range of pile lengths, diameters and slenderness. Moormann et al. [18] adopted Coupled Eulerian-Lagrangian (CEL) method to simulate the soil plug behavior in open-ended steel piles in loose sand to estimate the current design methods. For silty clay, soil plug behavior research focuses on open-ended steel piles with model tests [19] or field tests [20,21]. Few studies were conducted to investigate the soil plugging effect to the performance of jacked open-ended concrete pipe piles. Liu et al. [22] investigated the formation of soil plug in open-ended concrete pipe piles in the field, however, the effect of soil plugging on pile behavior was not discussed.

Further, the conventional strain gauges used in the monitoring of axial strains has some disadvantages [23]. Fiber optic sensors are becoming a well-established technology for a variety of geophysical and civil engineering applications [24]. Schmidt-Hattenberger et al. [25] conducted strain measurements by fiber Bragg grating sensors for in situ pile loading tests. The monitored strain results of fiber Bragg grating sensors agree fairly well with conventional concrete strain gages. Authors have also examined the skin friction between piles and subsoils as well as the settlement behavior of the pile based on the measured strain data. Klar et al. [26] conducted a static loading test on bored piles based on fiber optic technology and discussed the differences between fiber measurements and conventional discrete measurements. Doherty et al. [27] used a fiber sensing system to monitor the axial loads of soil-cement injected precase piles and steel pipe piles in the field. The advantages of fiber optic sensors compared to conventional strain gauges are long-term performance, resistance to corrosive environments, immunity against electromagnetic interferences, array-capability by wavelength-demultiplexing, and miniaturized size. Hence, there is a practical demand to investigate the behaviour of open-ended concrete pipe piles jacked into silty clay considering soil plugging using novel fiber optic sensors [28].

In this study, Fiber Bragg Grating (FBG) sensors were adopted to instrument two open-ended PHC pipe piles to study the soil plugging effect on pile behaviour when jacked into clay. The performance of test piles was continually monitored during the installation and the loading tests. The main objectives of this study were to investigate: (1) the behaviour of open-ended PHC pipe piles jacked into clay; (2) the soil plugging during and after jacking; (3) the residual forces after installation and the effect to the evaluation of total bearing capacity; (4) the loading transfer and settlement mechanisms of the test piles in loading test.

#### **2. Site Conditions and Test Program**

#### *2.1. Site Description*

Full-scale tests were conducted at a site in Hangzhou, China. Before the test pile installation, two auger borings *B*1 and *B*2, and one cone penetration test (CPT) were conducted to investigate the soil information, as shown in Figure 1. Figure 2a,b show the soil profiles from boreholes with the cone resistance *qc*, sleeve friction *fs* and friction ratio *Rf*. It is indicated that the test site is in sequence by 2.1 to 2.4 m of fill layer, 2.1 m of silty clay, 3.5 m of sandy clay, 2.0 m of muddy-silty clay, 2.5 m of silty clay, and 9.2 to 9.5 m thick alluvium layer. The soil properties in embedded length of test piles are summarized in Table 1.

**Figure 2.** CPT data and Su values in test site. (**a**) Soils profiles and CPT values along depth; (**b**) Friction ratio with depth; (**c**) Undrain shear strength with depth.

**Table 1.** Subsoil Properties.

Note: *c* and *ϕ* were determined using the quick shear tests; *α*1-2 was determined using one-dimensional compressibility test.

The coefficient of compressibility *α*1-2, defined as the decrease in volume per unit volume produced by a unit change of pressure, given in the table was determined through one-dimensional compressibility test. The cohesion *c* and friction angle *ϕ* were determined using direct shear test. The undrained shear strength *Su* is calculated from CPT results according to Equation (1) proposed by Schmertmann [29]:

$$S\_{\rm II} = (q\_t - \sigma\_{\rm v0}) / N\_{\rm kt} \tag{1}$$

where *qt* is the modified cone resistance; *σv0* is the total overburden pressure; *Nkt* is a cone bearing capacity factor. A *Nkt* value of 10 was used here.

#### *2.2. Test Piles Details and Instrumentation*

In order to separate all resistance components for partially plugged open-ended piles, the instrumented double-walled pile system is usually used in model tests [30,31] or full-scale field tests [32]. In this paper, the FBG sensors are installed along test piles using embedded methods. However, the survival rate of embedded FBG system in pile instrumentation is lower than that of conventional sensors because the exposed grating points are very fragile and require further cable protection [33]. Hence, an epoxy resin was used in this test to seal the FBG sensors and cables.

The test piles used in this study are two open-ended PHC pipe piles with 400 mm outer diameter, 75 mm wall thickness, and 13.0 m length. Details are given in Table 2.

#### **Table 2.** Details of Test Piles

Note: In size, 400 mm is the outer diameter of the test piles, 75 mm is the wall thickness; 36.0 GPa is the elastic modulus of the test piles material.

Fourteen FBG sensors were installed in each test pile at seven levels to monitor the axial strain along the test piles and then separate the end and shaft resistances. The detailed configuration of FBG sensors instrumentation is shown in Figure 3a. The sensors were installed about 2.5 m intervals and closer together near the pile toe. After two shallow grooves with 50 mm width and 15 mm depth were created, all FBG sensors were attached on the groove bottom according to pre-designed intervals. Then the FBG sensors were sealed with epoxy resin to prevent any damage caused by groundwater and the surrounding soils during jacking, as shown in Figure 3b,c. The axial strain recorded by FBG sensors could be transformed into axial forces. The Young's modulus of the test piles material was taken as 3.6 × 104 N/mm2 in this study without considering pre-stressing. The end resistance was estimated with linear extrapolation method using the measured data. As discussed above, the shaft resistances between soil plug and inner surface were not measured in this program. In other words, the estimated end resistances include annulus and plug resistances.

(**c**) **Figure 3.** Strain gauges installation on test piles. (**a**) Schematic arrangement for strain gauges (Unit: m); (**b**) Installation of strain gauges along test piles; (**c**) Strain gauges sealed with epoxy resin in field.

#### *2.3. Test Program*

Test piles were installed using a ZYJ680A jacking machine (Figure 4), which has a loading capacity of 6.8 tons. The jacking machine includes two components of clamps and hydraulic jacks. In jacking, the clamp with four arc pieces was used to stabilize against the pile vertically. Then a hydraulic jacking force was applied to push the piles. After the hydraulic clamp reached the maximum travel range, the hydraulic clamp was released and moved up to the top position. The step was repeated to push the piles down until the test piles were installed into the ground. Two test piles were both jacked into 13.0 m depth and founded in the layer of marine deposits. The ultimate bearing capacity of the test piles was estimated as 1200 kN from the calculation before loading test.

**Figure 4.** Jacking machine used in the test.

The jacking forces were recorded to assess the drivability of the test piles. The soil plug length during jacking was continuously measured using two different weights, as shown in Figure 5.

**Figure 5.** Measurement of soil plug length during pile jacking.

The loading test was conducted on test piles of T2 at 17 days after jacking, as illustrated in Figure 6. The total applied loads to the piles head were controlled by a hydraulic pump. The vertical displacement was monitored using two dial gauges mounted in the reference beams. The values of all FBG sensors attached to T2 were set to be zero before loading in order to measure the axial forces induced by applied loads. The static loading test was slowly maintained loading test and conducted in accordance with the Chinese technical code for testing of building foundation piles [34], with no unload-load loops.

**Figure 6.** Static loading test of T2. (**a**) Schematic view; (**b**) Static loading test in field.

#### **3. Test Resultis and Discussion**

#### *3.1. End and Shaft Resistance during Jacking*

The jacking forces of test piles during installation are plotted in Figure 7. It can be seen from the figure that the jacking forces increased with penetration depth. At final depth of 13.0 m, the jacking forces were 806 and 815 kN for T1 and T2, respectively. The induced axial forces along test piles are illustrated in Figure 8. It is indicated that the axial forces increased with jacking depth and the distribution shown similar characteristics to each other. It also should be noticed that the induced axial forces in Figure 8 are different from the jacking forces in Figure 7. The induced axial force on the piles head is the same as the jacking forces at the corresponding depth in Figure 7.

**Figure 7.** Jacking forces versus penetration depth.

For a better view of the loading transfer, the end and shaft resistances are separated and shown in Figure 9.

The shaft resistances of T1 and T2 consistently increased with installation, while the end resistances were in a range from 21 to 151 kN. At a penetration depth of 5.5 m, the shaft resistances of test piles reached a level of 457.3 and 471.2 kN, taking on 94.8 and 95.7% of the jacking forces, respectively. At 13.0 m, the shaft resistance for T1 was about 771.6 kN, amounting to 95.8% of jacking forces. For T2, the shaft resistance reached a value of 793.1 kN, while the end resistance was about 22.7 kN. The corresponding ratios of the shaft and end resistances to the jacking forces were 97.2 and 2.7%, respectively. The average shear stresses (*τ*av = shaft load/shaft area) and unit end resistance (*q*b) measured during installation are plotted against depth, as shown in Figure 10.

**Figure 8.** Variation of axial pile forces during jacking: (**a**) T1; and (**b**) T2.

**Figure 9.** Variation of end and shaft resistance during jacking: (**a**) T1; and (**b**) T2.

The *q*b is the force per unit cross section area, which the cross-section area is calculated as a closed-ended pipe piles. In order to compare with the cone resistance *q*c from CPT, the *q*c values at corresponding depths were also marked in the figure. It is apparent that these *q*b values are around 0.8 *q*c, which is in good agreement with those reported in previous database [35,36].

**Figure 10.** Pile installation stresses with depth. (**a**) Average shear stresses profile; (**b**) Average end bearing stress.

#### *3.2. Soil Plug Behavior*