IMPACT OF OFFSHORE SHIP EMISSIONS 

ON PUERTO RICO AND THE 

U.S. VIRGIN ISLANDS 

FINAL REPORT 

STI-909201-3729-FR 

By 

Neil J. M. Wheeler 

Stephen B. Reid 

Erin K. Pollard 

Kenneth R. Craig 

Yuan Du 

Sonoma Technology, Inc. 

1455 N. McDowell Blvd., Suite D 

Petaluma, CA 94954-6503 

Prepared for 

U.S. Environmental Protection Agency 

Office of Transportation and Air Quality 

2000 Traverwood Drive 

Ann Arbor, MI 48105 

October 27, 2009

iii

TABLE OF CONTENTS 

Section Page

LIST OF FIGURES ....................................................................................................................... iv 

LIST OF TABLES......................................................................................................................... iv 

1. INTRODUCTION.............................................................................................................. 1-1 

1.1 Background............................................................................................................... 1-1 

1.2 Overview of Technical Approach............................................................................. 1-2 

1.3 Summary of Findings ............................................................................................... 1-3 

2. TECHNICAL APPROACH............................................................................................... 2-1 

2.1 Data Selection and Processing.................................................................................. 2-1 

2.2 Creation of Back Trajectories................................................................................... 2-2 

2.3 Analysis of Back Trajectories................................................................................... 2-3 

3. FINDINGS AND RECOMMENDATIONS...................................................................... 3-1 

3.1 Transport Potential.................................................................................................... 3-1 

3.2 Discussion................................................................................................................. 3-9 

3.3 Conclusions and Recommendations....................................................................... 3-10 

4. REFERENCES................................................................................................................... 4-1 

APPENDIX A: FREQUENCY PLOTS BY RECEPTOR SITE ............................................... A-1 
iv

LIST OF FIGURES 

Figure Page

1-1. Ship activity patterns based on the ICOADS/AMVER data sets. ................................... 1-1

1-2. Map of receptor sites for back trajectory analyses........................................................... 1-2 

1-3. Overall spatial probability density map for all receptor sites. ......................................... 1-3 

3-1. Spatial probability density map showing maximum densities for all receptor sites, 

with a 50-m start height and meteorological data for the full year.................................. 3-2 

3-2. Spatial probability density map showing maximum densities for all receptor sites, 

with a 500-m start height and meteorological data for the full year................................ 3-3

3-3. Spatial probability density map showing maximum densities for all receptor sites, 

with a 50-m start height and meteorological data for August through November 

excluded. .......................................................................................................................... 3-4 

3-4. Spatial probability density map showing maximum densities for all receptor sites, 

with a 500-m start height and meteorological data for August through November 

excluded. .......................................................................................................................... 3-5 

3-5. Spatial probability density map for the Aguadilla site, with a 50-m start height and 

meteorological data for the full year................................................................................ 3-6 

3-6. Spatial probability density map for the Aguadilla site, with a 500-m start height and 

meteorological data for the full year................................................................................ 3-7 

3-7. Spatial probability density map for the Aguadilla site, with a 50-m start height and 

meteorological data for August through November excluded......................................... 3-8 

3-8. Spatial probability density map for the Aguadilla site, with a 500-m start height and 

meteorological data for August through November excluded......................................... 3-9 

LIST OF TABLES 

Table Page

2-1. HYSPLIT model inputs. .................................................................................................. 2-3 
1-1

1. INTRODUCTION 

1.1 BACKGROUND 

The purpose of this project is to evaluate the probability that offshore ship emissions 

impact selected onshore sites in Puerto Rico and the U.S. Virgin Islands. To conduct this 

evaluation, Sonoma Technology, Inc. (STI) performed back trajectory analyses, which use 

interpolated meteorological fields to track a parcel of air backward in time from a specified 

location. STI calculated back trajectories from seven locations in Puerto Rico and the Virgin 

Islands and aggregated that data to evaluate the probability that emissions from certain offshore 

locations have an impact on selected onshore sites. STI also generated plots showing 

frequencies of potential impact onshore from ship emissions source regions in the area. As 

shown in Figure 1-1, there is shipping activity, and therefore emission sources, surrounding these 

islands. 

Figure 1-1. Ship activity patterns based on the ICOADS/AMVER data sets.1 

(Puerto Rico is highlighted in red.) 

 

1The ship traffic patterns shown in Figure 1-1 are based on two combined data sets: the International 

Comprehensive Ocean-Atmosphere Data Set (ICOADS) managed by the National Oceanic and Atmospheric 

Administration (NOAA) and the Automated Mutual-Assistance Vessel Rescue System (AMVER) data set managed 

by the U.S. Coast Guard. 
1-2

1.2 OVERVIEW OF TECHNICAL APPROACH 

STI performed back trajectory analyses for four populated areas in Puerto Rico (San 

Juan, Pole Ojea, Ponce, and Aguadilla), two populated areas in the Virgin Islands (Charlotte 

Amalie on the island of St. Thomas and Christiansted on the island of St. Croix), and an 

additional location on the island of St. John (see Figure 1-2). 

Figure 1-2. Map of receptor sites for back trajectory analyses. 

The tasks required to complete the back trajectory analyses are listed below and 

described in detail in Section 2. 

1. Data selection and processing – STI reviewed available meteorological data for the 

region of interest, acquired the best available data for the years 2001-2005, and 

reformatted the data for use with the Hybrid-Single Particle Lagrangian Integrated 

Trajectory (HYSPLIT) model (Air Resources Laboratory, 2008). 

2. Creation of back trajectories – STI used the HYSPLIT model to calculate back 

trajectories for the seven receptor sites in Puerto Rico and the Virgin Islands. 

3. Analyses of back trajectories – STI aggregated the back trajectories to provide data and 

plots detailing the probability that trajectories originated in various source regions around 

each receptor site. STI also compared these plots to ship activity and emissions data 

provided by the EPA. 
1-3

1.3 SUMMARY OF FINDINGS 

The analysis methods used for aggregating an ensemble of trajectories are sufficient to 

provide an overall picture of transport in the region around Puerto Rico and the U.S. Virgin 

Islands and to inform an initial assessment of the potential impact of offshore ship emissions on 

populated areas in the region. Winds are predominantly from the east in this region, and the 

analyses confirm this observation, with the highest percentage of trajectories originating from 

locations east of Puerto Rico and the Virgin Islands. Figure 1-3, which displays the maximum 

probability that trajectories originating in a given location would impact any one of seven 

receptor sites, shows that the highest probability values are clustered around each of the receptor 

sites. Analyses indicate that transport from areas west of Puerto Rico is limited compared to 

transport from the east, north, and south. Dominant easterly flows may transport emissions from 

shipping hundreds of kilometers to populated areas of Puerto Rico and the Virgin Islands. This 

is especially true for shipping to the east, northeast, and southeast of these islands. 

Figure 1-3. Overall spatial probability density map for all receptor sites. 

2-1

2. TECHNICAL APPROACH 

The objective of this analysis is to evaluate the probability that offshore ship emissions 

impact populated areas in Puerto Rico and the Virgin Islands. To perform this evaluation, STI 

identified the sources of meteorological data for the region of interest, selected the most 

appropriate data sets, validated and processed the data, and used the HYSPLIT model to create 

back trajectories that could be used to generate frequency plots showing potential emissions 

impacts on selected sites. Descriptions of data sources and analysis steps for these tasks are 

provided in the sections below. 

2.1 DATA SELECTION AND PROCESSING 

After a review of available meteorological data for the region of interest, the National 

Center for Environmental Prediction’s (NCEP) North American Regional Reanalysis (NARR) 

data set (Mesinger et al., 2006) was chosen as the source of three-dimensional meteorological 

data for the HYSPLIT trajectory calculations. The NARR is a regional long-term, dynamically 

consistent, high-resolution weather and climate data set for North America. The data set was 

created by running NCEP’s Eta regional prediction model at 32-km resolution with a land 

surface physics parameterization package and the Eta three-dimensional variational (3D-VAR) 

Data Assimilation System (EDAS). EDAS assimilates a full complement of quality-controlled 

surface and upper air meteorological observations, plus additional precipitation analyses, satellite 

data, and other special observational data, to minimize the error between the observed and 

modeled atmospheric state. The NARR is an improvement in both resolution and accuracy 

(Mesinger et al., 2006) over the NCEP/National Center for Atmospheric Research (NCAR) 

Global Reanalysis data (Kalnay et al., 1996). 

The NARR was selected for this analysis because it is currently the only multi-year, 

medium-resolution, three-dimensional data set with sufficient coverage of the eastern Caribbean 

Sea region. Global-scale data sets, such as the NCEP/NCAR Global Reanalysis, were not used 

in this analysis because those data sets were developed at a spatial and temporal resolution 

coarser than that of the NARR, and they do not incorporate recent advances in modeling and data 

assimilation systems. 

NARR data products are distributed in gridded binary (GRIB) format on the NCEP 

Advanced Weather Interactive Processing System (AWIPS) 221 Lambert Conformal Conic grid 

at 32-km horizontal resolution on 29 vertical pressure levels at 3-hr intervals. The GRIB data are 

available through the NOAA National Operational Model Archive and Distribution System 

(NOMADS) and through the NCAR mass data store. STI had obtained 27 years of NARR data 

during previous work; thus, no additional data acquisition effort was required for this project. 

The NARR domain covers all of North America and extends into the Atlantic Ocean 

several hundred kilometers east of Puerto Rico and the U.S. Virgin Islands. A 2144-km by 

2335-km (67 by 73 NARR grid cell) analysis subgrid, roughly centered on Puerto Rico, was 

chosen for the HYSPLIT analysis. The eastern boundary of the analysis subgrid was 

intentionally set several grid cells from the eastern boundary of the NARR domain to avoid 
2-2

missing data near the NARR domain edge due to mismatch between the native Eta model grid 

and the AWIPS 221 grid. 

NARR data from the analysis domain were extracted from the GRIB data files and 

translated into the NOAA Air Resources Laboratory (ARL) packed binary file format required 

by HYSPLIT. The code developed to perform this task was based on codes developed by ARL 

to perform similar translations of other model data sets. This code also performed all unit 

conversions necessary to conform to HYSPLIT input requirements. Vertical velocity was 

converted from Pascals per second (Pa/s) to millibars per second (mb/s), while mean sea level 

pressure was converted from hectoPascals (hPa) to millibars (mb). Also, NARR wind fields 

were rotated to grid-relative coordinates. Unlike most modeled data sets, wind fields in the 

NARR have been rotated relative to true geographic north and must be re-rotated prior to use by 

a modeling system. 

Several quality assurance steps were taken to ensure that the NARR data were properly 

translated to ARL format. First, ARL-formatted data were processed through the chk_arl_data 

file check program developed by ARL to verify the data encoding. Vertical profiles at selected 

geographic locations from the ARL-formatted NARR data files were then examined to ensure 

the reasonableness of the vertical temperature, wind, and moisture profiles and to verify unit 

conversions. These vertical profiles were also qualitatively compared for consistency against 

vertical profiles extracted at the same locations from the operational 40-km North American 

Model (NAM) Data Assimilation System, a gridded data set produced by NCEP and translated to 

ARL format by ARL for use in HYSPLIT. Finally, selected HYSPLIT trajectories were 

calculated, using the ARL-formatted NARR and NAM data sets, and compared for consistency. 

2.2 CREATION OF BACK TRAJECTORIES 

The HYSPLIT model was developed by the NOAA ARL (Air Resources Laboratory, 

2008). HYSPLIT uses archived three-dimensional meteorological fields to compute air parcel 

trajectories, and the model can either be run interactively on the web or be run using executable 

code downloaded to a personal computer. 

On the basis of discussions with the EPA Work Assignment Manager (WAM), STI ran 

HYSPLIT using a back trajectory length of 72 hours and trajectory start heights of 50 and 500 m 

above ground level. The 50-m start height was selected to match the effective stack height for 

ocean-going marine vessels (Mason et al., 2008), while the 500-m start height represents the 

middle of the typical mixed layer, which generally has a top that is from 800 to 1,000 m high in 

the region of interest. Running HYSPLIT with these two start heights provides some insight into 

potential differences between near-surface and elevated winds and helps to bound the potential 

impacts of off-shore ship emissions. In addition, HYSPLIT was run using the most recent five 

years of NARR data (2001–2005). Meteorological data from a five-year period was used in 

order to minimize the influence of any invalid data or unusual meteorological events such as 

hurricanes. Table 2-1 provides a summary of HYSPLIT model inputs. 
2-3

Table 2-1. HYSPLIT model inputs. 

Parameter Inputs 

Meteorological Data NARR; 32-km resolution; 2001-2005 

Trajectory Length 72 hours (backward) 

Receptor Locations 7 sites: San Juan, Puerto Rico; Pole Ojea, Puerto Rico; Ponce, Puerto 

Rico; Aguadilla, Puerto Rico; Charlotte Amalie, St. Thomas; 

Christiansted, St. Croix; and St. John Island 

Trajectory Start Heights 50 m and 500 m above ground level 

Vertical Motion Constant sigmaa

aWith the constant sigma vertical motion option in HYSPLIT, a trajectory remains on its original terrain-following 

sigma surface. STI also tested the isobaric option, which derived vertical motion from isobaric data, but there was 

no impact on the resulting aggregated trajectories.

2.3 ANALYSIS OF BACK TRAJECTORIES 

After generating back trajectories for 2001-2005 for each of the seven locations of 

interest, STI used geographical information system (GIS) software to aggregate the trajectories 

and calculate the percentage of time that the trajectories originated in various offshore regions 

around Puerto Rico and the Virgin Islands. STI’s analyses extended at least 500 km in all 

directions from the receptors, and the analyses were performed for two periods: (1) the entire 

year and (2) the year with the rainy season (August-November) excluded. Once these analyses 

were complete, STI prepared frequency plots showing the percentage of time that trajectories 

originated in various regions around the site. Plots were prepared for each receptor site 

individually and for the seven sites overall. These overall plots represent the maximum 

probability that trajectories originating in a given location would impact any of the receptor sites. 

These plots were developed by combining spatial probability values for all seven sites in a GIS 

environment and using GIS to select the maximum probability value among the seven data sets 

for a given location. In addition, plots were prepared for trajectory start heights of 50 m and 

500 m. 

3-1 

3. FINDINGS 

3.1 TRANSPORT POTENTIAL 

Figures 3-1 through 3-4 show overall spatial probability density plots for combined 

trajectories from all seven receptor sites based on 2001-2005 meteorological data (note that the 

black outline around each of the frequency plots shows the extent of the NARR data used for 

these analyses2). These plots show the maximum probability that trajectories originating in a 

given location would impact any of the receptor locations. Note that these plots do not consider 

additive impacts. For instance, if emissions from a given location have a 50% chance of 

reaching Aguadilla and a 20% chance of reaching San Juan, a value of 50% is recorded on the 

plot, rather than a higher value approaching 70%. In addition, emissions reaching shore between 

receptor site areas are not considered (receptor site areas are 32-km grids). As such, these 

overall spatial probability density plots may understate the true impact on land. 

The overall spatial probability density plots are shown with trajectory start heights of 50 

and 500 m for all meteorological data (Figures 3-1 and 3-2) and for meteorological data 

excluding data from the rainy season (Figures 3-3 and 3-4). Multiple trajectory start heights 

were used to represent the effective stack height for ocean-going marine vessels (50 m) and the 

middle of the typical mixed layer in the region (500 m). Running HYSPLIT with these two start 

heights provides some insight into potential differences between near-surface and elevated winds 

and helps to bound the potential impacts of off-shore ship emissions. Data from the rainy season 

were excluded because the meteorological characteristics of these periods may bias long-term 

averages, and the increased mixing and deposition that is likely to occur during stormy weather 

would alter the impact of off-shore emissions on receptor locations. However, the similarity of 

Figures 3-1 through 3-4 indicates that the selection of trajectory start heights and the exclusion of 

rainy season (August-November) meteorological data do not significantly alter transport 

patterns. 

In general, these figures show that easterly flows dominate the region, as the highest 

percentage of trajectories originates in locations east of Puerto Rico and the Virgin Islands. As 

expected, the highest probability values for emissions reaching any one of the seven receptor 

sites are clustered around each of the receptor sites. Analyses indicate that transport from areas 

west of Puerto Rico is limited compared to transport from the east, north, and south. Dominant 

easterly flows may transport emissions from shipping hundreds of kilometers to populated areas 

of Puerto Rico and the Virgin Islands. This is especially true for shipping to the east, northeast, 

and southeast of these islands. 

Similar patterns are observed when back trajectories for individual receptor sites are 

considered. For instance, Figures 3-5 through 3-8 show spatial probability density values 

calculated from back trajectories from the Aguadilla site. Plots for all seven receptor sites are 

provided in Appendix A. 

 

2

 The polygon around the analysis domain appears skewed because the NARR data are provided in a Lambert

Conformal projection, which becomes distorted in regions distant from the projection’s two reference parallels. 
3-2 

Figure 3-1. Spatial probability density map showing maximum densities for all 

receptor sites, with a 50-m start height and meteorological data for the full year. 
3-3 

Figure 3-2. Spatial probability density map showing maximum densities for all 

receptor sites, with a 500-m start height and meteorological data for the full year. 
3-4 

Figure 3-3. Spatial probability density map showing maximum densities for all 

receptor sites, with a 50-m start height and meteorological data for August 

through November excluded. 
3-5 

Figure 3-4. Spatial probability density map showing maximum densities for all 

receptor sites, with a 500-m start height and meteorological data for August 

through November excluded. 
3-6 

Figure 3-5. Spatial probability density map for the Aguadilla site, with a 50-m 

start height and meteorological data for the full year. 
3-7 

Figure 3-6. Spatial probability density map for the Aguadilla site, with a 500-m 

start height and meteorological data for the full year. 
3-8 

Figure 3-7. Spatial probability density map for the Aguadilla site, with a 50-m 

start height and meteorological data for August through November excluded. 
3-9 

Figure 3-8. Spatial probability density map for the Aguadilla site, with a 500-m 

start height and meteorological data for August through November excluded. 

3.2 DISCUSSION 

Overall, the spatial probability plots shown in Section 3.1 and in Appendix A offer a 

reasonable picture, obtained by using the best available meteorological data, of transport in the 

region of interest. However, some limitations are associated with these analyses. First, the 

32-km resolution of the NARR data may not capture local circulation patterns (e.g., land-sea 

breezes that may have a more westerly component). In addition, caution should be exercised 

when interpreting probability values within about 100 km of the eastern boundary of the analysis 

domain. This is because, for each data point, a search radius of 50 km was used to calculate the 

fraction of trajectories that passed through a given grid cell; near the eastern boundary of the 

analysis domain, however, the 50-km search radius extends into regions where no meteorological 

data are available. 
3-10 

3.3 CONCLUSIONS 

The analysis methods used for aggregating an ensemble of trajectories are sufficient to 

provide an overall picture of transport in the region around Puerto Rico and the U.S. Virgin 

Islands and to inform an initial assessment of the potential impact of offshore ship emissions on 

populated areas in the region. Analyses indicate that transport from areas west of Puerto Rico is 

limited compared to transport from the east, north, and south. Dominant easterly flows may 

transport emissions from shipping hundreds of kilometers to populated areas of Puerto Rico and 

the Virgin Islands. This is especially true for shipping to the east, northeast, and southeast of 

these islands. 
4-1

4. REFERENCES 

Air Resources Laboratory (2008) HYSPLIT model. Web page of the National Oceanic and 

Atmospheric Association, Air Resources Laboratory, Silver Spring, MD. Available on 

the Internet at http://www.arl.noaa.gov/ready/hysplit4.html. 

Kalnay E., Kanamitsu M., Kistler R., Collins W., Deaven D., Gandin L., Iredell M., Saha S., 

White G., Woollen J., Zhu Y., Leetmaa A., Reynolds B., Chelliah M., Ebisuzaki W., 

Higgins W., Janowiak J., Mo K.C., Ropelewski C., Wang J., Jenne R., and Joseph D. 

(1996) The NCEP/NCAR 40-year reanalysis project. Bull. Amer. Meteor. Soc. 77 (3), 

437-471. 

Mason R., Dolwick P., Carey P., Kinnee E., and Wilson M. (2008) Emissions processing and 

sensitivity air quality modeling of category 3 commercial marine vessel emissions. Paper 

presented at the 17th International Emissions Inventory Conference, Portland, OR, June 

5, by the U.S. Environmental Protection Agency, Research Triangle Park, NC. 

Mesinger F., DiMego G., Kalnay E., Mitchell K., Shafran P.C., Ebisuzaki W., Jovic D., Woollen 

J., Rogers E., Berbery E.H., Ek M.B., Fan Y., Grumbine R., Higgins W., Li H., Lin Y., 

Manikin G., Parrish D., and Shi W. (2006) North American Regional Reanalysis. Bull. 

Amer. Meteor. Soc. 87 (3), 343-360. 

A-1

APPENDIX A 

FREQUENCY PLOTS BY RECEPTOR SITE 
 
A-3

Spatial probability density map for the Aguadilla site, with a 50-m start height and 

meteorological data for the full year. 
A-4

Spatial probability density map for the Aguadilla site, with a 500-m start height 

and meteorological data for the full year. 
A-5

Spatial probability density map for the Aguadilla site, with a 50-m start height and 

meteorological data for August through November excluded. 
A-6

Spatial probability density map for the Aguadilla site, with a 500-m start height 

and meteorological data for August through November excluded. 
A-7

Spatial probability density map for the Pole Ojea site, with a 50-m start height and 

meteorological data for the full year. 
A-8

Spatial probability density map for the Pole Ojea site, with a 500-m start height 

and meteorological data for the full year. 
A-9

Spatial probability density map for the Pole Ojea site, with a 50-m start height and 

meteorological data for August through November excluded. 
A-10

Spatial probability density map for the Pole Ojea site, with a 500-m start height 

and meteorological data for August through November excluded. 
A-11

Spatial probability density map for the Ponce site, with a 50-m start height and 

meteorological data for the full year. 
A-12

Spatial probability density map for the Ponce site, with a 500-m start height and 

meteorological data for the full year. 
A-13

Spatial probability density map for the Ponce site, with a 50-m start height and 

meteorological data for August through November excluded. 
A-14

Spatial probability density map for the Ponce site, with a 500-m start height and 

meteorological data for August through November excluded. 
A-15

Spatial probability density map for the San Juan site, with a 50-m start height and 

meteorological data for the full year. 
A-16

Spatial probability density map for the San Juan site, with a 500-m start height 

and meteorological data for the full year. 
A-17

Spatial probability density map for the San Juan site, with a 50-m start height and 

meteorological data for August through November excluded. 
A-18

Spatial probability density map for the San Juan site, with a 500-m start height 

and meteorological data for August through November excluded. 
A-19

Spatial probability density map for the Charlotte Amalie site, with a 50-m start 

height and meteorological data for the full year. 
A-20

Spatial probability density map for the Charlotte Amalie site, with a 500-m start 

height and meteorological data for the full year. 
A-21

Spatial probability density map for the Charlotte Amalie site, with a 50-m start 

height and meteorological data for August through November excluded. 
A-22

Spatial probability density map for the Charlotte Amalie site, with a 500-m start 

height and meteorological data for August through November excluded. 
A-23

Spatial probability density map for the Christiansted site, with a 50-m start height 

and meteorological data for the full year. 
A-24

Spatial probability density map for the Christiansted site, with a 500-m start 

height and meteorological data for the full year. 
A-25

Spatial probability density map for the Christiansted site, with a 50-m start height 

and meteorological data for August through November excluded. 
A-26

Spatial probability density map for the Christiansted site, with a 500-m start 

height and meteorological data for August through November excluded. 
A-27

Spatial probability density map for the St. John site, with a 50-m start height and 

meteorological data for the full year. 
A-28

Spatial probability density map for the St. John site, with a 500-m start height and 

meteorological data for the full year. 
A-29

Spatial probability density map for the St. John site, with a 50-m start height and 

meteorological data for August through November excluded. 
A-30

Spatial probability density map for the St. John site, with a 500-m start height and 

meteorological data for August through November excluded.