Analytical Modeling Tool for Design of Hydrocarbon Sensitive Optical Fibers

Pipelines are the main transportation means for oil and gas products across large distances. Due to the severe conditions they operate in, they are regularly inspected using conventional Pipeline Inspection Gages (PIGs) for corrosion damage. The motivation for researching a real-time distributed monitoring solution arose to mitigate costs and provide a proactive indication of potential failures. Fiber optic sensors with polymer claddings provide a means of detecting contact with hydrocarbons. By coating the fibers with a layer of metal similar in composition to that of the parent pipeline, corrosion of this coating may be detected when the polymer cladding underneath is exposed to the surrounding hydrocarbons contained within the pipeline. A Refractive Index (RI) change occurs in the polymer cladding causing a loss in intensity of a traveling light pulse due to a reduction in the fiber’s modal capacity. Intensity losses may be detected using Optical Time Domain Reflectometry (OTDR) while pinpointing the spatial location of the contact via time delay calculations of the back-scattered pulses. This work presents a theoretical model for the above sensing solution to provide a design tool for the fiber optic cable in the context of hydrocarbon sensing following corrosion of an external metal coating. Results are verified against the experimental data published in the literature.


Introduction
Oil transfer pipelines are used as a means for the transportation of crude oil mixtures to designated processing facilities. Such pipelines are extensively used in the UAE and span large distances [1] while being subject to adverse weather conditions. These conditions are characterized by high temperatures and high levels of ambient humidity during the summer season. As a result, these structures are constantly attacked by the various mechanisms of corrosion both internally and externally, providing motivation to monitor these large structures for corrosion damage.
Conventionally, several standardized monitoring tools are available to pipeline operators for the purpose of inspection [2]. Most internal monitoring tools involve intruding the pipeline. A Pipeline Inspection Gage (PIG) is routinely used for cleaning pipelines. Modern PIGs are instrumented for wall thickness and surface defect detection. The sensors used employ Magnetic Flux Leakage (MFL) as a tool for inferring wall thickness along the pipelines' axial direction [3]. "Smart" PIGs conduct their inspection runs on off-line pipelines, which disrupts operation. Additionally, visual inspection probes may be administered into the pipeline once it is off-line. Closed Circuit TV (CCTV) technology is used to map the internal surface and collect data about its condition. Similarly, the pipeline must be off-line so that the probe may be administered into it for inspection [4].
In order to mitigate the need for disturbing pipeline operation, disposable corrosion coupons have been employed so that they may continuously monitor the corrosion rate as opposed to the periodic inspection intervals offered by the previous tools. Corrosion coupons are inserted into the pipeline and retrieved using a retrieval valve and tool. Once in position, the corrosion current generated by the oxidation of iron may be used to continuously infer the corrosion rate on the internal side of the pipeline [5]. Although this technique offers continuous monitoring of the pipeline conditions without disturbing operation, the required retrieval valves involve extensive retrofitting costs for existing pipelines, making them economically impractical as a monitoring solution. In addition to cost, corrosion coupons have other disadvantages, and they are: • The corrosion coupons may not be experiencing the same flow regime as the pipe surface, thus indicating different corrosion rates than the pipelines.

•
Build-up of scale, asphaltene and wax on the corrosion coupons can impact the true corrosion rate accuracy, particularly if the buildup occurs on the coupons, but not on the interior surface of the pipeline.

•
Removal of the corrosion coupons from the pipelines for inspection can lead to deaths due to the release of toxic gases (such as H 2 S, etc.) from pipelines.
Due to the above shortcomings of conventional inspection tools and due to the volatile operating environment within pipelines, a sensing element and communication medium that are intrinsically safe for live pipelines have been targeted as a likely solution. Optical Fiber Sensors (OFS) have been used for Structural Health Monitoring (SHM) in the past and would likely overcome the above shortcomings. Optical fibers are the most natural and logical sensors to be used in the oil and gas industry since data are transmitted through the fiber using pulses of light, the data transfer is extremely fast and light cannot cause sparks and fires. Optical fibers are also immune from electromagnetic interference.
OFS have been categorized into long gage, point, quasi distributed and distributed sensors [6]. Long Period Grating (LPG) sensors are used in civil structures to monitor displacement of a taught sensing fiber through phase modulation against a reference fiber using interferometry. Fiber Bragg Grating (FBG) sensors are used to monitor local strains and temperatures in civil structures using wavelength modulation of a broadband light pulse in relation to distortion of the grating spacing. Other point sensors such as micro and macro bend sensors have also been reported [6,7]. Some of the above-mentioned sensor topologies have found applications specifically in corrosion sensing. Optical fibers have been used as pH sensors via immobilized dyes that have an absorption spectrum sensitive to hydrogen anions. Spectroscopic methods may be used to detect the hydrogen ion concentration and hence the pH level that is strongly related to corrosion. Other approaches utilize LPGs whose cladding is spectrally sensitive to sodium chloride concentrations. Stripped cladding sensors utilizing thermally-deposited metal coatings may be used to detect its removal as transmitted power is restored to previous levels before applying the metal cladding [8][9][10]. Point reflectivity probes may be used to scan the surfaces of pipes and walls to detect surface damage or bacterial colonies [6,11].
Finally, truly distributed fiber optic sensing networks have been realized using backscatter reflectometry in which the time delay between firing a laser pulse and receiving the backscattered pulse along various locations on the fiber is resolved to provide spatial information. The spectral profile of the backscattered light pulse may be used to infer temperature and strain effects by virtue of Rayleigh, Brillouin or Raman scattering [6,12].
Distributed fiber optic sensors have found several applications in pipeline leakage detection. One technique performs distributed temperature measurements via Raman spectroscopy in order to detect leakage from a nearby pipeline transporting fluids at a higher temperature relative to the surroundings [13]. Another technique performs distributed acoustic measurements via the Refractive Index (RI) modulation due to strains. The RI change and its location can be inferred using distributed Sagnac interferometry by interpreting the generated phase signal [14].
Further investigation into distributed fiber optic sensors yielded an approach used for hydrocarbon leakage detection via Optical Time Domain Reflectometry (OTDR). A polymer-clad fiber is sensitive to the hydrocarbons due to the swelling of the cladding and the changes in its RI as a result of hydrocarbon absorption. OTDR is used to measure the loss in light intensity due to chemical perturbations in a polymer-clad fiber that alter its RI [15].
Due to the abundance of hydrocarbons in oil pipelines and the need for corrosivity detection in the internal environment of the pipeline, a hybrid method combining metal-coated fibers with the polymer-clad fibers to detect excessive corrosivity levels is theoretically investigated. The technique employs a high power laser diode that fires a laser pulse through the fiber under inspection [16]. OTDR is used to resolve the location of metal loss from the coating following hydrocarbon absorption into the now exposed cladding. A change in the attenuation coefficient characterizes the losses of a traveling light pulse.
A theoretical investigation of the above phenomenon was explored by the authors in [17]. The results presented are limited in terms of validation and do not provide sensing predictions in the context of internal corrosion sensing of pipelines by virtue of hydrocarbon contact. This work extends [17] by providing a more thorough theoretical tool that is verified against experimental results from the literature and is used to provide a calibration curve that is useful for designing the proposed sensing solution for maximum sensitivity in crude oil pipelines.

Sensor Deployment
It is envisioned that the optical fiber corrosion sensor, proposed in this paper, will be deployed in the oil pipelines either through the PIG launcher or retrieval stations or through valves. The work in [18] shows a proven method to install fiber optic cable inside gas pipelines using existing valves. The fixtures, proposed in [18], have been designed to provide gas-tight seals around the cable entry and exit points. Just like [18], it is envisioned that our proposed optical fiber corrosion sensor will lay on the bottom of the pipeline very near the surface of the pipeline; therefore, seeing the same flow regime as the pipe surface; thus, experiencing similar corrosion rates as the oil pipelines. If sediments such as scale, wax, or asphaltene form in the interior surface of the pipe, since the optical fiber is also very near the pipe surface, it will experience the same phenomenon.
The transduction principle of the proposed sensor relies on a break in the metal coat due to excessive corrosion. Upon exposure of the polymer cladding underneath, hydrocarbons in the surrounding medium diffuse into the polymer, causing its RI to change. This problem will be modeled by utilizing the established theoretical models for light intensity distribution and a specially-developed numerical model for arbitrary RI profiles as explained in the following sections.
From a practical point of view, it is worth noting that the diffusion of hydrocarbons into the cladding is indeed reversible and does not affect the integrity of the fiber. The work in [15] (Section 3.3) makes use of this phenomena to regenerate the same fiber after testing by evaporating the volatile diffused hydrocarbons in air or acetone. We must note that the model presented in this paper does not account for the effects of swelling.
Another practical aspect of the proposed methodology could be the temperature sensitivity of the cladding RI to the high temperatures encountered within oil transfer lines. A study by [19] on photonic Polydimethylsiloxane (PDMS) crystals shows a very small thermo-optic effect (variation of RI with temperature). PDMS will be the material of choice for the cladding of the fiber for reasons explained in subsequent sections. RI changes with temperature were in the order of 10 −5 , which are unlikely to have an appreciable effect on the attenuation of the signal. The work in [20] mentions in an experimental study and characterization that RI of light crude mixture varies from 1.4530 at 20 • C to 1.4379 at 60 • C. In large transfer lines, the oil will be closer to the surrounding temperature of 40-50 • C, which is within the range of RI values included in subsequent calculations within this paper.
From a fabrication point of view, uniform deposition of metal, similar to the pipeline, can be achieved on optical fibers using the sputtering technique (physical vapor deposition). Many fiber optics companies already provide optical fibers with polyamide, gold, silver, copper, copper nickel alloy, aluminum, platinum and palladium coatings. It should be possible to coat the optical fibers with steel from the parent pipeline (API 5L X-grade material). Coating thicknesses, currently used by the fiber optics companies, range from 15 µm for single-mode fibers to 60 µm for large-core multi-mode fibers. We believe that thicker coatings suitable for this application are feasible because metallic deposition technology is currently widely used for manufacturing of automotive trim, optical components and even bathroom fixtures.
A schematic of the perturbed fiber is shown in Figure 1.

Mathematical Modeling of Light Intensity Distributions
The paper will introduce several mathematical models for the diffusion of hydrocarbons into the polymer cladding. A model for the calculation of the resulting modal capacity and light intensity profiles due to the diffusion of the hydrocarbons into the cladding is then developed. Finally, the resulting loss in intensity from a reference fiber configuration (before hydrocarbon intrusion) is computed in each case and compared to experimental measurements found in the literature.
A mathematical model is sought to calculate the light intensity distribution for arbitrary RI profiles in order to compare the resulting power to that carried by a reference case and devise an attenuation coefficient representative of the perturbation. Following the derivations presented in the literature [21][22][23], the following derived models build on the established governing equations to define the ordinary differential equations that describe the light intensity profile while holding the RI parameter, an arbitrary function of radius. Analytical results of the governing equations (Maxwell's equations) in the literature describe intensity distributions for fixed and quadratic core refractive index profiles. However, no general purpose method exists for solving these differential equations in the case of arbitrary RI profiles (especially in the cladding). The resulting RI profiles due to hydrocarbon concentration gradients in the cladding could result in such arbitrary RI profiles. Hence, a numerical tool suited for this case will be developed in Section 2.2.2.

Derivation of General Light Intensity Field Distributions in a Circular Fiber Optic Waveguide
To begin with, the assumed solutions of the light intensity distribution given by Equations (1) and (2) for the core and cladding regions respectively are substituted into the governing equations (please see the Supplementary Information for a review of the derivations of the wave equation from Maxwell's equations).
e t (r, φ) = e cl (r) cos lφx f or r > a core .
(2) e t represents the intensity distribution to be solved for as a function of the radial and azimuthal coordinates used in cylindrical coordinates to define a circular fiber; e co and e cl represent the radial intensity distributions for the core and cladding respectively and are to be solved for. Here, l is the azimuthal mode number, which can be any non-negative integer; 0, 1, 2..... The Helmholtz equation, which is known as the wave equation in cylindrical coordinates, is derived from Maxwell's equations. In terms of the defined field intensity distribution, the wave equation is given as follows: where µ 0 is the magnetic permittivity of free space, ε 0 is the electrical permittivity of free space, ω is the angular velocity of the light wave, which is derived from the wavelength, β is the propagation constant of the fiber, which is a parameter to be derived from the boundary conditions of the problem, and n (r) is the RI profile function. Substituting Equation (1) into Equation (3) yields the first differential equation in the core region: where k 2 = µ 0 ε 0 ω 2 . Similarly, substituting Equation (2) into Equation (3) yields the differential equation in the cladding region: The term k 2 n 2 − β 2 is defined as follows for the core region [21]: where V is the modal capacity of the fiber defined as V = ka core n co 2 − n cl 2 , a core is the fiber core radius and b is the generalized guide index, which is computed from boundary conditions. Equation (6) is used to eliminate β from Equation (4): where n (r) denotes the arbitrary profile of the RI in the core region of the fiber and n co denotes the nominal RI of the core at the central axis of the fiber. A similar manipulation is performed on the cladding differential equation knowing that the term k 2 n 2 − β 2 is defined as follows for the cladding region [21]: Similarly, Equation (5) is left in terms of n (r) to accommodate an arbitrary cladding RI, and Equation (8) is used to eliminate the propagation constant β: where n (r) denotes the arbitrary profile of the RI in the cladding region of the fiber and n cl denotes the nominal RI of the cladding. Equations (7) and (9) represent the general differential equations governing the solution for the radial tangential field, which defines the intensity distribution for any particular mode.

Derivation of State Equations to Be Solved Using Numerical Computation
The development of a numerical model for the arbitrary-index fiber involves solving the governing differential equations for the light intensity. The governing equations after expanding Equations (7) and (9) are given by: However, since the equations are divided on the interval 0 ≤ r ≤ a core and r > a core , each equation needs a change of variables such that both equations share the same parametric variable that ranges from 0-1. This is done by first changing the first equation in the core region using the following change of variables to ψ: The equation in the cladding region is given by the following change of variables: The state variables for the above equations are: While the derivative d q dψ denoted as ∇ q is: From Equation (10) and the above transformations, we have: where Ψ = ψ (b clad − a core ) + a core and 0 ≤ ψ ≤ 1. The above equation may be written in the final form in terms of the elements of the vector q, Equation (11) as: However, the use of a numerical solver would prove difficult in this situation since there are singular terms in the first two elements of ∇ q; − 1 / ψ q [2] and l 2 / a 2 ψ 2 when ψ → 0. This problem is circumnavigated by approximating the value of the function e co and de co dψ near the singular coordinate (zero) by using a Taylor series approximation for the function and its derivative at a point d near the origin as follows: where, e co = de co dψ and e co = d 2 e co dψ 2 . e co (0) = 0 due to the symmetric boundary condition at the center of the fiber core, and e co is obtained from Equation (10) with the coordinates transformed to ψ as follows: to yield the boundary conditions at ψ = 0: where the ψ term has been replaced where appropriate with the infinitesimal variable d to avoid singular terms: p , and p = e co (0) is an unknown parameter that must be solved for. Furthermore, the slope and value of e co and e cl at ψ = 1 and ψ = 0 respectively should be equal: where E l is an arbitrary amplitude coefficient that determines the power input into the fiber. The remaining boundary conditions are set at ψ = 1 for the slope e cl as r → ∞: provided that a sufficiently large value for the cladding diameter b clad is chosen to approximate its end as ∞.
The problem contains two unknown parameters (the generalized guide index b and the parameter p) along with two sorder differential equations requiring a total of six boundary conditions, which are now summarized below: The boundary value problem is also summarized below and written in terms of the elements of the state vector q, the core radius a core , the cladding radius b clad , the core RI n co , the cladding RI n cl , the generalized frequency k = 2π/λ and the modal parameter V = ka core n 2 co − n 2 cl . The model is in terms of normalized radial coordinates ψ = r/a core and Ψ = ψ(b clad − a core ) + a core . Here, l is the azimuthal mode number, which determines the mode configuration being solved for, and E l is an arbitrary amplitude constant that determines the power carried by an individual mode.
Now that the numerical model is established, it is solved using MATLAB's built-in boundary value solver, which utilized three-stage Lobatto methods to obtain the solution. The results of the analytical solution will be used as a guess for the numerical model in order to initialize the solution. The numerical approach is verified against the known step index solution [21][22][23] for the case of step index profiles to check the accuracy of the numerical model.

Analytical Solutions to Step Index Profile Problems
Step index fibers feature constant core and cladding RIs with a step drop to the cladding RI. An analytical solution exists for this type of fiber, which is obtained by solving Equations (7) and (9) directly by making the substitution of n (r) = n co for the core region equation and n (r) = n cl for the cladding region producing the following well-known results [21]: where J l is the Bessel function of the first kind, K l is the modified Bessel function of the second kind and E l and E l are amplitude constants to be determined through boundary conditions. The result of applying the continuity boundary condition at the core/cladding boundary yields the following results for the amplitude coefficients and the generalized guide index b: where B = V/ a core (1 − b). Next, the solutions to Maxwell's equations are verified for the developed numerical model in Section 2.2.2 against the theoretical model in Section 2.2.3.

Verification of Developed Numerical Model against Step Index Solutions
Both models will be solved for the fiber presented below as it is similar to the application described in [15]. The parameters used for the case study are summarized in Table 1. First a step index fiber is assumed and is compared against the numerical arbitrary index profile solver where n(r) is fed into the solver as a constant RI in the core and cladding regions for a few selected modes that can propagate in the fiber.
First, the characteristic equation, Equation (24), is solved for b. It cannot be solved analytically and must therefore be solved graphically by plotting the right-hand side versus the left-hand side to obtain all the possible intersections and, hence, obtain a value for the parameter B from which b may be inferred. The results for the analytical dispersion relation for the fiber above for l = 1, l = 2 and l = 5 is solved for and used as an initial guess for the generalized guide index b to be input into the solver. An example of the solution is shown in Figure 2 for l = 1.  (Table 1), l = 1.
Next, both the analytical solver, Equations (21) and (22), and the numerical solver, Equations (13)- (20), are run simultaneously for the following linearly-polarized modes; LP 1,5 ,LP 2,2 and LP 5,4 where LP l,m is the designation for a linearly-polarized mode with azimuthal number l and radial mode number m. A radial mode number of m = 1 would correspond to the first mode solution beginning from the left of the abscissa. Furthermore, a value of E l = 1 is used in both models in order to normalize the result at the cladding/core interface. The amplitude constant E l is set based on how much power is coupled into a specific mode. This parameter will later on be used to distribute the input power to a fiber among the various possible modes that can be accommodated.
The results of the these comparisons are shown in Figure 3. The value for the propagation constant b is obtained by using the analytical dispersion solution as an initial guess. Since for this case both methods are being used to solve the same case for a step index fiber, the same value of b is expected from the numerical search approach. Table 2 below summarizes and compares the results of both cases for each generalized guide index.
The maximum percentage difference with reference to the analytical results for b is in the order of 10 −3 , which is expected since the numerical solver is already at the convergence point.
After gaining confidence in the numerical model, the solutions to Maxwell's equations developed in Section 2.2 will be used to calculate the light intensity profiles in the fiber before and after contact with hydrocarbons. The results will be compared with those presented in [15], which conducted OTDR experiments on different hydrocarbons characterized by different RIs. Table 2. Comparison between the propagation constant estimates for both the analytic and numerical approaches to step index problems.

Characterization of OTDR Traces Using Established Theoretical and Numerical Models
In order to verify the developed theoretical tool, the attenuation results from [15] will be used. The developed theoretical tools for calculating the intensity distribution will be used to calculate the net intensity profile before and, after contact with hydrocarbons, calculate the net power and the corresponding attenuation coefficient, which will be compared against the attenuation results reported in the literature. Table 3 summarizes the properties of the fiber optic cable used for the experimental setup in [15]. A test was carried out assuming that the chemical in contact with the fiber was trichloromethane (CHCl 3 ) with an RI of 1.4459. Two possible scenarios will be explored for the possible RI distributions due to soaking in CHCl 3 .

Refractive Index Perturbation Mechanisms
The first scenario is considered by assuming a first order diffusion model of the hydrocarbon in question through the cladding layers. According to [24], the mass concentration of the chemical agent in question will grow in time until reaching a steady-state concentration profile that varies linearly between the initial concentration of 0 mol/L at the core/cladding interface and the ambient concentration of the hydrocarbon in the medium surrounding the cladding. This results in the diffusion modeled RI profile depicted by Figure 4.
The second scenario is considered by assuming a mixed RI as a result of combining the optical properties of the hydrocarbon in question (CHCl 3 ) and the Polydimethylsiloxane (PDMS) polymer that forms the cladding. The mixed RI is predicted by calculating the volume fraction of each constituent under the assumption of complete soaking of CHCl 3 into the PDMS cladding. Table 4 below outlines the details of the calculations used to compute the volume fraction of each constituent.  In order to calculate the RI of a mixture of PDMS and CHCl 3 under the assumption that all the present free volume is occupied by the hydrocarbon, the generalized Lorentz-Lorenz relation [27] is adapted for this problem: where n 1→m is the RI of m constituents, n i is the RI of each constituent and φ i is the volume fraction of each constituent. Using the results from Table 3 and the established equation for the free volume fraction in PDMS, the effective cladding RI is calculated as n cl = 1.44435. The results of the previous scenarios are summarized by Figure 4 in order to depict the cases for which the light distribution is to be solved.

Mode Configuration in Reference and Perturbed Fibers
Next, the corresponding permissible mode sets are computed in order to determine the number of modes propagating in the fiber. Each mode is represented by a combination of an azimuthal and radial mode number denoted as l and m, respectively. For each unique solution of the generalized guide index b, a non-negative integer value for m is assigned, with m = 1 corresponding to the largest unique solution for b. Here, l designates the azimuthal mode number, which controls the cylindrical distribution of the radial mode solution obtained as follows: e t (r, φ) = e co (r) cos lφx f or r ≤ a core (26) e t (r, φ) = e cl (r) cos lφx f or r > a core (27) where e co and e cl are radial mode solutions in the core and cladding regions obtained from Equations (19) and (20) or the theoretical solutions [21][22][23]. The parameters l and m are enumerated and fed to the solver that is used to obtain the intensity distribution for each individual mode. The possible enumerations of the mode sets are characterized by the limiting equations for guided and leaky modes [22,28] as depicted graphically in Figure 5b.     Table 3. The same fiber under a hydrocarbon perturbation in the RI as a result of the mixed RI assumption is shown in Figure 5b. It is assumed that light traveling through the fiber is fully developed before encountering the perturbation. As a result, when the modal capacity is reduced, some modes are immediately attenuated, while others travel in the radiation domain denoted as "leaky" modes. These leaky modes are attenuated considerably faster than guided modes, but due to the short length of the perturbation, they are returned to the set of guided modes as soon as the modal capacity of the fiber is restored following the perturbation. The remainder of the modes are permanently decoupled from the fiber. This effect is represented by the black hatched area in Figure 5b. In other words, this area represents the lost energy from the fiber due to the perturbation.

Cumulative Light Intensity Distribution of All Possible Mode Sets
The permissible mode configurations explained in Section 3.2 can now be solved for using the developed solver in Section 2.2.2 and the analytical solution in Section 2.2.3. The solver is run for each mode while setting the amplitude constant E l to unity. E l determines the amount of power coupled into each mode relative to the other supported modes. By setting E l to unity for all modes, an equal power distribution among all permissible modes is enforced. In a practical test setup, this condition is met by the use of a mode scrambler to equally excite all light propagation modes [15]. The intensity distribution e t is recorded, and the corresponding power for each mode is calculated as: [21] All the modes are normalized to carry unity power by dividing the square of each intensity distribution e t 2 by the total power carried as calculated in Equation (28) yielding the square of the normalized intensity distributionê 2 t . The total power in the fiber is therefore given by the sum of all the normalized intensity distributions.
However, the number of modes is not simply given by the number of permissible combinations of l and m as shown in Figure 5. Examining Equations (26) and (27), for l = 0, it can be seen that the cosine terms become unity. However, for l > 0, the cosine terms may assume two identical states in terms of magnitude, cos(lθ) and cos(lθ + π/2) which results in each mode being counted twice due to the presence of two polarization states [28].
All supported modes for each assumption are normalized, summed and integrated by Equation (28). The ratio of the power of the perturbed fiber to the reference fiber is used to calculate the step drop per unit length (i.e., the attenuation coefficient α), which characterizes the step drop: P perturbed P re f erence = 10 −α∆z/10 , α = −10 log 10 P perturbed P re f erence /∆z , where ∆z is the interval of the step drop usually set to a unit meter. With the above procedure, the attenuation coefficients and respective OTDR traces can be generated for both scenarios explained in Section 3.1.

Power Attenuation Results for Diffusion Modeled Fibers
Utilizing the pre-established mode sets in Figure 5a, the numerical model for diffusion modeled clad fibers returned the intensity distributions for individual modes. Figure 6 displays the intensity results for a few selected high order modes who have strongly affected evanescent waves. The drop in intensity is not very significant as shown in Table 5, where the maximum change in the generalized guide index was 2.98%. As a result the cumulative intensity distribution is not likely to show a significant drop in power as shown in Figure 7 where the net intensity distribution slightly drifted from that of the reference step index fiber. The resulting change in power will be translated into an attenuation coefficient and compared with the results of Section 3.5.

Power Attenuation Results for Mixed Refractive Index Fibers
Utilizing the pre-established mode sets in Figure 5b, the analytical solution for step index fibers is used to obtain the light intensity solution for each mode and uses the summation algorithm from Section 3.3 to obtain the cumulative light intensity distribution. The total results are reported in Figure 8. The attenuation coefficient is then reported relative to the reference fiber and summarized in Section 3.6.

Comparison and Discussion of Diffusion Modeled and Mixed RI Fibers as Describing Models of CHCl 3 Perturbation
A comparison between diffusion modeled fibers and mixed RI fibers is performed to determine which of the two models best describes the reference case obtained through experiments in the literature. As noted in Section 3.4 earlier, the diffusion model is likely to yield unrealistic results due to the assumption that the RI at the core/cladding boundary is fixed. Table 6 displays this result showing a very small comparative step drop in power levels. Mixed RI modeling exhibits more promising results as the step drop experienced due to the modeled RI profile closely agrees with the value reported in the literature of 1.46 dBm −1 resulting in more credibility being awarded to the mixed RI model. Due to the close agreement in the results between the reference experimental case and the mixed RI model, a calibration curve for the proposed sensor and the variation of the step drop with various liquids with various RIs will be studied by using the mixed RI model to calculate the step drop for different hydrocarbons.

Sensor Calibration Data and Discussion of Sensor Operation for Various Hydrocarbon Stimulants
Using the pre-established mixed RI model, a calibration curve describing the proposed sensor was sought by solving the model for selected RIs at 5% intervals. The resulting step drops computed from the cumulative light intensity distributions for each cladding interval are depicted in Figure 10, while the intensity profiles at each increment are shown in Figure 9. To assess the level of agreement with experimental data [15], the instrument used for obtaining the OTDR traces (TFS 3031, Tektronix) [15] was investigated in order to obtain a measure of uncertainty in the results.

Design Stage Uncertainty of the Tektronix OTDR Unit
The manufacturer specifications and the uncertainty characteristics of the (TFS 3031, Tektronix) instrument are given in the manufacturer datasheets. Relevant parameters are summarized in Table 7. As a result, the uncertainty in the measurement is given as follows: For the case of a step drop of 1.46 dB, the total design stage uncertainty is given by: The resulting uncertainty is incorporated into the experimental measurements to verify whether the discrepancy between theoretical calculations and experimental measurements is simply due to hardware uncertainty. It can be seen from Figure 10 that the discrepancy is not within hardware uncertainty limits, which justifies a physical explanation in the assumptions made.
A two-term exponential fit can be used for the simulated points to obtain a theoretical calibration curve for various RIs of the form: α(n) = ae b·n clad + ce d·n clad (30) where coefficients are given as: a = 2.109; b = 0.8399; c = −0.2702; d = −0.4631.
Note that this calibration is only valid in the range 1.470 ≤ n clad ≤ 1.4550 where n clad is normalized by the mean (µ) = 1.446 and the standard deviation (σ) = 0.005627 as follows: n i = n clad + µ σ .
Below n = 1.4370, the attenuation coefficient is equal 0 dB, and above 1.4550, complete power loss occurs or, for simplicity, losses in excess of 20 dB are expected. Please note that n clad is calculated from Equation 25 to account for the mixing of the hydrocarbon with the PDMS cladding as demonstrated in Table 4. Figure 10 displays a close agreement between the simulated data and the experimental data points obtained by [15]. The work in [15] tested two hydrocarbon compounds; CHCl 3 and CCl 4 , as shown in the figure. Table 8 shows the experimental discrepancy between the simulated calibration curve and the experimental points. Table 8. Comparison between experimental data points in [15] and the simulated calibration curve.

Effective RI of HC Expected Attenuation Attenuation Difference and Cladding Drop (Simulated) dB
Drop dB (%) * (Equation (25)) (Equation (30) There is an increasing level of discrepancy as the RI approaches the core RI. The reasons for this discrepancy could be attributed to the idealization of the assumed RI profile. However, the increasing level in discrepancy can be attributed to the increasing count of leaky modes as the gap between the core RI and the cladding/hydrocarbon mixed RI diminishes. Figure 11 demonstrates how the number of leaky modes is increased when the RI of the hydrocarbon increases. Figure 11 shows this affect as the number of leaky modes shown in green increases compared to the number of guided modes when increasing the RI from 1.444-1.447.
The attenuation coefficient calculations in this paper assume that the permitted leaky modes are loss-less and behave just like guided modes. However, that is not necessarily valid. According to [29], a leaky mode is defined as a mode that is attenuated to 1% of its original power value. This corresponds to a drop of 20 dB per meter according to Equation (29). If the perturbation length is assumed to be 1 cm, then the drop would only be 0.2 dB, which is small enough, but comparable to the calculated attenuation coefficient for CHCl 3 of around 1.46 dB from Figure 10. In other words, the assumption that the perturbation is infinitesimal is not a valid one and must be accommodated for by studying the attenuation coefficient associated with leaky modes and incorporating them into the established mathematical model for a closer agreement.
The mathematical model provides a powerful design tool for distributed sensors that use optically-sensitive claddings to transduce an attenuation in the light pulse. Based on the optical properties of alkanes present in crude oil, 1.3571 < n < 1.3886 [30], the sensor needs to be redesigned to operate in this region. As a result, the model will be used to design the fiber to be used for future experimental verifications. Data obtained through material handbooks about the optical properties of several of these compounds was obtained as shown in Table 9.
Propagating modes leaky modes permitted leaky modes (b) Possible modes for n clad = 1.447. Figure 11. Comparison of permissible guided and leaky mode sets for (a) 40% of the total RI range (b) 60% of the total RI range. Figure 12 demonstrates the sensitivity of this sensor to crude oil and how the fiber should be redesigned for the best performance in the expected environment known to exist in oil transfer lines. A best fit line of the simulated attenuation coefficients was projected on a coordinate space showing the RI band for certain hydrocarbon species (namely those most abundant in oil and gas pipelines). The sensing fiber in its current configuration is beyond the bandwidths of any of these species (especially alkanes). As a result, when considering field deployment, the chemical composition of the pipeline environment should be known such that the core and cladding RIs are tuned to provide a detectable signal attenuation.  For example, the core RI should be set at approximately 1.38, while the cladding RI should be set at around 1.355 such that the fiber detects contact with alkanes. However, caution should be exercised when making these determinations such that the hydrocarbon compound in question does not possess an RI close to, but lower than that of the core RI. Such a situation would cause a huge attenuation in the signal rendering the remainder of the fiber insensitive as the entire signal has been stripped from the fiber.
However, if the sensing band is chosen to give a reasonable attenuation in the signal, the fiber can detect hydrocarbon intrusion at multiple points along its length allowing for spatial monitoring of the fiber condition without the need for replacement of the already exposed sections of the fiber.

Conclusions and Future Work
An analytical modeling tool was developed to simulate the performance of a hydrocarbon OTDR sensor operating as a corrosion indicator in oil pipelines. The data obtained by the tool agree with several data points from the literature with a certain measure of discrepancy due to a number of idealizations made. The modeling tool describes the proposed sensor's performance for a variety of hydrocarbon stimulants allowing the sensor to be redesigned for operation in crude oil mixtures. It will help to provide insight into the design of the experimental setup needed to accommodate and test the sensor's operation.
Furthermore, mathematical modeling of all the guided and leaky modes provides insight into the attenuation coefficient of each mode, which can be used to infer the length of the perturbation to be made in the sensing fiber. The step drop exhibited following a hydrocarbon perturbation should contain two components: attenuation due to loss of higher order modes and the relatively large attenuation associated with leaky modes. By being able to quantify the latter, the length traveled by the leaky modes throughout the perturbation may be inferred.

Potential Applications of This Technology
In this paper, an optical fiber corrosion sensor with a hydrocarbon sensitive cladding, coated with a metal layer similar to the pipeline material, was introduced to detect internal corrosion of oil pipelines. In this paper, the cladding was sensitive to hydrocarbons, but the cladding can be made sensitive to other chemicals or conditions such as corrosion-causing bacteria, scale, asphaltene, wax or all the above. The potential applications of our proposed fiber optics sensor in the oil and gas industry are enormous.