VERIFICATION OF A MARINE POLLUTANT SURFACE PLUME MODEL FOR USE IN THE DEVELOPMENT OF AUTONOMOUS VEHICLE TRACKING SYSTEMS

Robotic systems can efficiently and precisely obtain environmental measurements, while operating in locations that are too dangerous or inaccessible for humans (Dunbabin & Marques, 2012)(National Oceanic Service, 2015)(Shukla & Karki, 2016). The next and latest step is to enable these robotic systems to operate autonomously. Autonomous robotic systems are more time and cost effective than traditional survey methods involving a human operator. In some cases, they can be operated at one third of the daily cost and do not require a highly skilled operator such as needed for a remotely operated vehicle (Shukla & Karki, 2016). Additionally, an autonomous robotic system can operate 24 hours a day. With an autonomous robotic system, the operator will have the opportunity to analyze the data in real time instead of operating the vehicle (Shukla & Karki, 2016). To efficiently develop an autonomous robot for tracking marine surface plumes, it is necessary to test the control system in a simulated environment that accurately reflects a spill spatially and temporally.

Previous research on maritime pollutant plume tracing has shown, and the Rhodamine Dye experiment presented here corroborate, that the observed instantaneous concentration values will be significantly different than the time-averaged concentration(Moore & Atema, 1991)(Crimaldi, Wiley, & Koseff, 2002)(Jones, 1983). The fine scale concentration characteristics of the plume complicate the autonomous robotic tracking of the plume movement. The autonomous controller designed for tracking marine pollutant spills will need to be robust enough to handle the concentration variations found in a pollutant plume. It is essential that the model creating the simulated environment accurately predict the fine scale concentration variations the robot will observe.

The particle trajectory model is a commonly used marine transport spill model. Examples include the well-known General NOAA Operational Modeling Environment (GNOME) and MEDSLIKII models. Both of these models have been determined to be accurate at predicting the movement of the pollutant, but have not been evaluated for accuracy in predicting the concentration of the pollutant at a specific location (De Dominicis, Pinardi, Zodiatis, & Archetti, 2013)(Zelenke, O’Connor, Barker, Beegle-Krause, & Eclipse, 2012) or the accuracy of producing the fine scale variations found in the previous studies.

In this study, a simplified particle trajectory model for simulating a neutrally-buoyant pollutant was used to predict the pollutant concentration at a specific location. This model was not designed to take into account all of the physical phenomenon that occur during a pollutant spill; as such it includes the transportation processes, but omits the weathering processes. In an attempt to decrease computation time, the model neglects weathering processes, which are expected to have less effect on surface concentrations over the time scales of the experiment (ITOPF Ltd., 2011)

The model generated simulated pollutant spill was compared with an experimental spill using statistical parameters to determine if the model provides similar general spill characteristics and fine scale concentration variations. A similar study was conducted to compare an atmospheric model to the fine scale structure of an atmospheric pollutant plume (Farrell, Murlis, Long, Li, & Carde, 2002). The statistical parameters used in that study were chosen to verify that the simulated plume resulted in the same general characteristics as the experimental plume in amplitude and temporal statistics. Those plume characteristics identified are the characteristics which complicate the problem of tracking the plume by creating nonconstant fine scale plume concentration variations(Farrell, Murlis, Long, Li, & Carde, 2002). The study presented in this paper uses similar statistical parameters to evaluate the effectiveness of the marine pollutant plume model.

The particle trajectory model uses a combination of Lagrangian and Eulerian specifications to construct the path of the pollutant plume (Al-Rabeh, Cekirge, & Gunay, 1989). The model breaks the total mass of the pollutant pollutant into particles and tracks the movement of each particle. The particles are released into the water evenly in accordance with the predetermined spill rate, and then transported by advection and diffusion. The Lagrangian specification of the horizontal advection-diffusion equation in discrete time can be written as

where P_t is the particle position at the beginning of the time step, P_t+1 is the particle position at the end of the time step, Adv_x/y is the advection in the x and y direction respectively and Dif_x/y is the diffusion in the x and y direction respectively (Chao, Shankar, & Wang, 2003). The model calculates the transportation in the x and y direction independently.

Advection is the change in particle position due to the velocity of external flow field(Reed, et al., 1999). The advection equation can be written as

where V_FF is the velocity of the surface flow field and Δt is the length of each time step. The particles are considered to be point masses during the transportation phase of the model; therefore, drag is not considered to be a factor. The velocity of surface flow field (V_FF) is a combination of the drift due to wind and water currents. Wind drift is approximately 3% of the wind velocity (National Oceanic and Atmospheric Administration, 2002). In this model, the water drift was set as 97% of the water current velocity. The equation for calculating the flow field can be written as

where V_water is the water current velocity and V_wind is the wind velocity (Wang, Shen, & Zheng, 2005).

The diffusion term accounts for both molecular and turbulent diffusion. The diffusion equation is written as,

where DM is the diffusion magnitude and θ is the angle of diffusion. DM is modeled as a random walk. DM and θ are calculated using

where R_m and R_θ are uniformly distributed random numbers between 0 and 1 and c_d is the effective diffusion coefficient. R_m and R_θ are calculated separately in the x and y directions. The diffusion coefficient is assumed to be constant in the horizontal plane; i.e. in both the x and y directions (Chao, Shankar, & Wang, 2003).

Two particle characteristics were developed during the creation of this model for use in determining the pollutant concentration at a specified location: particle mass and particle distribution extent. The particle mass is the portion of the total mass that each particle has and remains constant throughout the simulation. The particle distribution extent is the standard deviation of the distribution of the particle mass around the center location of the particle. The mass distribution of the particle was found by multiplying a simplified joint probability density function for the bivariate normal distribution (Hogg & Tanis, 2010) with the total particle mass. The formula can be written as

where md is the mass distribution function, m is the particle mass and σ is the standard deviation in both the x and y direction.

The particle characteristics were calibrated for the environment where the experiment was conducted: salt water, near shore, partially protected and with Rhodamine WT dye as the “pollutant.” The calibrated particle characteristics used to compare the experimental and model results are particle mass of 0.01 g and a particle distribution extent of 0.6 m, which simplifies the mass distribution function to

Concentration is the mass of the pollutant per unit volume of a solution. In this simplified model, the pollutant was assumed to stay near the surface of the water in a well-mixed layer. The amount of the particle within the sample area is determined using the particle mass distribution function and integrating it to find the volume of the concentration in the sample area. The function to determine the portion of the particle mass that is within the sample area can be written as

where SM is the portion of the mass that is within the sample area, md is the mass distribution function and m is the total mass of the particle(Hogg & Tanis, 2010). The total concentration is the sum of the portion of mass for all particles that is within the sample volume divided by the sample volume and can be written as

where p is the set of particles, E is the mass of the particles that is within the sample area and v is the volume of the sample area in liters. The volume was calculated using an well-mixed surface layer with an assumed depth of 0.1 m.

Several experiments were conducted over the course of two weeks at the Makai Pier on Oahu, HI to collect data on the pollutant concentration of a marine pollutant pollutant spill. This location was chosen because the breakwater pier creates a small basin, protecting the site from the major ocean swells while still allowing wind and water currents. The pier also provides access from three of the four sides of the observation area. The average water depth of the area is 4 m. Environmental data was collected throughout the experiment. The water currents, up to 2 meters deep, were measured by a Doppler Velocity Logger at multiple locations around the basin, and the wind was measured using an anemometer.

The experiments consisted of two unmanned marine surface vessels (USV) [ Figure 1]; each equipped with two fluorometer concentration sensors. Each fluorometer was located approximately 0.1 m below the surface of the water. The USV’s collected the time, location and concentration of the pollutant at 10 Hz. During the experiments, Rhodamine WT was used to simulate the marine pollutant. It is a nontoxic highly visible dye (Rhodamine WT Liquid Safety Data Sheet, 2015). The dye was released from a known source location at a constant rate of 5 g/min for each experiment [Figure 2]. Four experiments were conducted at four different locations in the pollutant plume. The locations were defined by their distance from the pollutant source: 6.7 m, 22.5 m, 34.9 m and 51.3m.

The recorded data from the experiment provides a time series of dye concentration for the stationary point in the pollutant plume. The data was only examined for the time period where the concentration was at steady state, i.e. after the particle plume front had reached the sensor location. This allows the analysis to find the fine scale concentration structures within the pollutant plume.

A statistical analysis was conducted for both the time series data gathered during the experiment and the data produced by the model. The time series was used to determine the maximum, mean, standard deviation, skewness and kurtosis of the concentration at steady state. The pollutant plume amplitude and temporal structures were characterized using statistical methods, including the coefficient of variation, peak-to-mean ratio, intermittency, burst length and burst return.

The coefficient of variation is a measure of the dispersion of the concentration. It is a non-dimensional number that can be used to compare the plume structures regardless of the magnitude of the concentration. A smaller coefficient of variation implies that the variation in the concentration with respect to time was small with respect to the mean concentration. The equation for the coefficient of variation can be written as

where CoV is the coefficient of variation, C̄ is the mean concentration and σ is the concentration standard deviation.

The peak-to-mean ratio is a measure of the extremeness of the maximum concentration in relation to the mean concentration. The closer the peak-to-mean ratio is to 1.0 the closer the maximum and mean concentrations are together. The peak-to-mean ratio can be written as

where PMR is the peak-to-mean ratio, C_max is the maximum concentration and C̄ is the mean concentration.

Intermittency is the percentage of time the concentration is below a threshold level. The equation for intermittency can be written as

where γ is the intermittency, τ is the threshold value, T is the total time, c_t is the concentration at the time step and I is the indicator function, where I is 1 if the concentration is less than the threshold concentration and 0 if it is greater (Liao & Cowen, 2002). Intermittency is a function of the threshold value; i.e. intermittency increases as the threshold value increases.

It is important to know how often the concentration in the plume could be observed as zero and thus the robot could observe a zero concentration and still be within the plume. The control system needs to be able to account for that. Thus, the threshold for calculating intermittency was set at 1.5 ppb, which is the sum of the sensor’s lower detection limit and error. This will give the percentage of time that the sensors could see a zero concentration value. A high intermittency value indicates that a large portion of the time the concentration is below the sensor reading threshold.

The burst length is the time that the concentration is continuously above a threshold value and the burst return is the time the concentration is continuously below the amount. As these two values approach each other, the concentration is spending an equal time above and below the threshold value. There is no standard threshold value to use in comparison of the burst length and return. To analyze the effect of the particle mass and particle distribution on the burst length and return, the threshold was set as the mean pollutant pollutant concentration. This value was different for each experiment and model observation. As the burst length and burst return approach each other, the average frequency around the mean concentration is equal to the burst length and return values.

Table 1 outlines the results of the statistical analysis performed on the experimental results. The results show a general decrease in parameter value as the distance from the pollutant source increased. The temporal parameters (intermittency, burst length, and burst return) do not follow that pattern and show a wide variety of values in sensors located at the same distance from the pollutant source.

Table 2 outlines the statistical parameters for the model concentration time series. The results show a similar pattern to the experimental results with a general decrease in parameter value as the distance from the pollutant source increased; however, the amplitude parameters (coefficient of variation and peak-to-mean ratio) in the model statistical results show an increase in value as the distance from the source increases. This is opposite of what was seen in the experimental results.

The results show that the model was able to predict some fluctuations in the concentration; however, there were some significant differences in the values calculated from the model results and the experimental observations. The percent error was used to compare the statistical parameters [Table 3]. The statistical results of maximum, mean, standard deviation, skewness and kurtosis for the concentration predicted from the model were significantly different from concentration time series observed during the experiment, over 100% error. The structural statistics were more similar with less than 90% error observed.

The model preformed the best in mirroring the plume amplitude structural statistics: coefficient of variation and peak-to-mean ratio (average 51% error). Even though the maximum, mean and standard deviations of the time series were significantly different (>300%), the statistical parameters that capture variation and temporal structural statistics, CoV through BR, have the lowest amount of relative error and were similar enough that the model reproduces to some extent the variable nature of the concentration in the plume. If the autonomous robotic control system is robust enough to handle the model, it will be sufficient to operate in a marine environment that is similar to the model environment, near shore in an enclosed basin. More testing is required to determine if this model is sufficient for modeling in the open ocean.

When modeling the plume temporal structural statistics, the model performed worst in intermittency, error of 88%. This occurred because there was always some intermittency in the experimental results; however, the model rarely produced intermittency, which is a potential important factor of control system testing. This can be seen in the Figure 3, which shows the concentration time series for both the experimental and model data.

Figure 3 also shows the model produces a high frequency variation that is not present in the experimental data. This high frequency variation occurred because the number of particles inside the sample area was highly variable. This model does not take into account spreading, which can be an important factor during the first several hours of a pollutant spill (ITOPF Ltd., 2011). It is important to note that spreading is primarily associated with hydrophobic pollutants, such as oil. Here, spreading is being considered a change in the model’s particle mass distribution extent as a function of time. Increasing the particle mass distribution extent over time will decreases the variability in the number of particles within the sample window and reduce the high frequency variations in the model.

A simplified particle trajectory model was used to predict the concentration of a pollutant at a point location. The simulated data produced by the particle trajectory model was tested against experimental data to determine is accuracy and reliability for use in the development of autonomous robotic control systems. The experimental location was a near shore basin, partially protected by a pier. Rhodamine WT was used in place of a toxic pollutant pollutant. Sensors were placed at a different location in the plume during each experiment. Those locations were at different distances from the dye source: 6.7 m, 2.5 m, 34.9 m, and 51.3 m. A time series of the dye concentration was gathered at each location when the dye plume had reach steady state.

The scenarios were analyzed using statistical parameters. The amplitude comparison statistics used were the peak-to-mean ratio and the coefficient of variation; the temporal comparison statistics were intermittency, burst length and burst return. The comparison of the simulated concentration data with the experimental data showed that the amplitude statistics had lower normalized errors than the temporal statistics. Although the model only loosely matches experimental data, after tuning, this simplified model was able to relatively quickly produce simulations of real-time, fine scale pollutant plume structure that can be used for testing control algorithms in simulation.

For future work, the effects of each of the model components (or future added components e.g. spreading), from a computational speed standpoint, need to be quantified. Then the model should be improved to better match the temporal statistics of the experimental data; most importantly by determining the minimum number of additional model considerations, again from a computational speed standpoint, that need to be added to the model to reduce the high frequency variation in concentration and improve the intermittency.