Background information

This page contains information specific the PCB simulation embedded in the MITgcm driven by the ECCOv4 state estimate. 

Code origins

Polychlorinated biphenyls (PCBs), mercury (Hg), and perfluorinated compounds (PFCs) are global contaminants of concern. MITgcm simulations exists for each of these contaminants. Yanxu Zhang, was the first  to implement contaminant cycling in the MITgcm. Yanxu implemented inorganic mercury cycling in the MITgcm ECCOv1 (Y. Zhang et al., 2015) and later implemented organic mercury cycling in the MITgcm ECCOv4 (Y. Zhang et al., 2016). The mercury simulation in MITgcm ECCOv1 was subsequently adapted by Xianming Zhang for PFC cycling. Helen Amos adapted the MITgcm ECCOv4 PFC simulation for PCBs and later updated the PCB simulation to run with MITgcm ECCOv4.

The MITgcm simulates the physical transport of contaminants in the ocean through advection, diffusion, and ocean currents. The ECCOv4 ocean state estimate was created by Gael Forget and is driven by a re-anlysis data set spanning 1992-2012 (Forget et al., 2015a,b). Marine plankton dynamics are simulated in the MITgcm by the DARWIN model (ECCOv1: Dutkiewicz et al., 2012). DARWIN had not been run in ECCOv4. Yanxu Zhang, with the help of Helen Amos and Steph Dutkiewicz, was the first to run DARWIN in ECCOv4, which allowed mercury and PCBs to be updated from ECCOv1 to ECCOv4. The primary advantage of ECCOv4 is that the ocean grid covers the whole globe. Previously, ECCOv1's grid domain terminated at ~65 N and did not include the Arctic Ocean. 

(Helen M. Amos, 26 June 2016)

Migrating from Svante to Odyssey

The PCB simulation in MITgcm ECCOv1 was developed on Harvard's compute cluster "Odyssey" while updating the PCB simulation into MITgcm ECCOv4 happened on MIT's compute cluster "Svante". The PCB MITgcm ECCOv4 code, archived runs, and Matlab packages for visualizing model output were migrated from Svante to Odyssey by GEOS-Chem Support Team members, Bob Yantosca and Lizzie Lundgren. All files copied from svante are archived on Odyssey at the following location:

/n/sunderland_lab/share_root/Lab/svante/

The PCB component of the MITgcm code is incorporated into the MITgcm master branch of the Harvard-BGC MITgcm repository for use on Odyssey. You can access the code by cloning the git repositiory from the Harvard-BGC account on Bitbucket. See instructions for using MITgcm on Odyssey elsewhere on this wiki.

(Lizzie Lundgren, 15 August 2016)

Svante archive on Odyssey

A brief description of directories in /n/sunderland_lab/share_root/Lab/svante/helen/:

• GEOS_Chem/ : Carey Friedman's global atmospheric PCB simulation (Friedman & Selin, 2015). This is my personal copy I had been modifying to work with the MITgcm / GEOS-Chem Coupler. 
• IDL/ : Contains a copy of GAMAP2 visualization software. Used to plot and do quick visualizations with GEOS-Chem model output. 
• MATLAB/ : Contains all Matlab routines to process either GEOS-Chem or MITgcm ECCOv4 model data. 

  • gcmfaces/gcmfaces/gcmfaces_GEOSChem/ : Purpose-built routines to regrid and transform GEOS-Chem data to the MITgcm ECCOv4 grid and vice versa. 

  • gcmfaces/gcmfaces/gcmfaces_utility/: Purppose-built routines to plot observed and modeled PCB ocean concentrations and fluxes.

• MITgcm_ECCOv1.tar.gz : Zipped copy of the PCB simulation implemented in MITgcm ECCOv1. This is an old version of the PCB code and is not recommended for use.

• MITgcm_ECCOv4/ : The PCB simulation implemented in MITgcm ECCOv4. The directory structure is the same for all ECCOv4 simulations (e.g., PCBs, Hg, PFCs) and is described in MITgcm ECCOv4 Svante Guide in Section 2. Obtain source code. 
  • pkg/pcb/ : Directory containing the PCB Fortran routines making up the basis of the PCB simulation.
  • verification/global_pcb_llc90/build/ : Directory where Fortran routines for PCB simulation are compiled. Note 
• NCL/ : Some templates for NCL plotting routines. 
• net_fs02_d2/ : MITgcm ECCOv4 PCB simulation run directory and archived PCB model output. ECCOv4 model output is very large, so we were asked to run the model on a file server instead of the head node. We were on file server 2 on Svante, hence "fs02". 
  • MITgcm_ECCOv4/verification/global_pcb_llc90/run/ : Working run directory.
  • MITgcm_ECCOv4/verification/global_pcb_llc90/results_pcb/ : Archived PCB model output. 
• public_resources/ : Copies of code and results shared with other MITgcm users. 
• temp/ : Place to hold temporary files. 

(Helen M. Amos, 2 Jul 2016)

 Compiling the PCB simulation

Selecting a congener

The PCB simulation is set up to run one congener at a time. Users must select the congener before compiling the code. The code is set up to simulate the seven IUPAC congeners: PCB-28, 52, 101, 118, 138, 153, and 180. Open the following FORTRAN file to select the congener you want to simulate:

cd /n/sunderland_lab/share_root/Lab/svante/helen/MITgcm_ECCOv4/pkg/pcb/

vi pcb_param.F

A user specifies which congener they want to simulate by uncommenting the physiochemical constants for their desired congeners and commenting out the values for the other congeners. Here's an example for the molar volume for PCB-28:

The user needs to set the following six physiochemical constants for their desired congener:

1. Molar volume, Vb
2. Molecular mass, Mb
3. Octanol-water partition coefficient, logKow
4. Degradation half life, tHalf
5. Enthalpy of air-water exchange, dH
6. Dimensionless gas-over-liquid Henry's law constant, KH0
   
   

(Helen M. Amos, 27 June 2016)

Compiling

The compilation procedure is the same for every congener. After selecting your congener, compile the MITgcm by following these instructions:

• On Svante, MITgcm ECCOv4 Svante Guide, Section 5. Compiling the Code
  
• On Odyssey, 4.3. Compiling MITgcm for the PCB simulation. 

(Helen M. Amos, 2 Jul 2016)

Running the PCB simulation

A copy of my working (on Svante) PCB (in ECCOv4) run directory has been archived on Odyssey here:

/n/sunderland_lab/share_root/Lab/svante/helen/net_fs02_d2/MITgcm_ECCOv4/verification/global_pcb_llc90/run/

The structure of the run/ directory for any contaminant (Hg, PCB, PFC...) is going to be identical, except the input_*/ directory containing atmospheric forcing fields and the data* file pointing to PCB input files. Both will be named for the contaminant of interest. For PCBs, data.pcb contains the file paths pointing to the atmospheric forcing fields in input_PCB/. Resources for running any ECCOv4 simulation are available on these two wiki pages:

• MITgcm ECCOv4 Odyssey Guide, Section 6. Running the Simulation (by Yanxu Zhang)
  
• Running the MITgcm PFOS and Hg Simulations on Odyssey (by Bob Yantosca) 

MITgcm ECCOv4 runs generally have to be submitted to a queue system using a bash or c-shell runscript. My runscript to submit a single PCB job to Svante's queue system is called qsub_itXX.csh. 

Setting up a run for a specific congener

The PCB simulation is set up to run one congener at a time. After selecting the congener to run and compiling the code, you need to make sure the atmospheric forcing fields (deposition and atmospheric concentration) match the congener compiled. The atmospheric forcing files are set in data.pcb. Let's say you want to run PCB-28, the following paths need to point to PCB-28 files:

 # Atmospheric concentrations (mol/m3) 

## 10 yrs of monthly means from GEOS-Chem

 PCBa_concFile = './input_PCB/PCB28_popgconc_1992-2010_eccov4_llc90.bin',

 

# Atmospheric deposition (mol/m2/s)

## 10 yrs of monthly means from GEOS-Chem

PCBp_depoFile  = './input_PCB/PCB28_partdep_1992-2010_eccov4_llc90.bin',

PCBa_depoFile  = './input_PCB/PCB28_gasdep_1992-2010_eccov4_llc90.bin',

Binary files for PCBa_concFile, PCBp_depoFile, and PCBa_depoFile for congeners PCB=218, 101, 153, and 180 for the years 1930 through 2010 are archived in run/input_PCB/. Atmospheric deposition and concentration fields are all monthly means and were generated from archived GEOS-Chem atmospheric PCB simulations. The GEOS-Chem simulations were performed by Carey Friedman for Friedman & Selin (2016). Additional details about how the files for PCBa_concFile, PCBp_depoFile, and PCBa_depoFile were generated can be found in the README in run/input_PCB/.

If you are initializing the simulation for zero PCBs in the ocean, there is nothing else you need to do to run for a specific congener. If you are running a "pickup" simulation (i.e., starting your current run from the PCB concentrations at the end of the previous run), you need to make sure your pickup.0000* file in run/ matches the congener you wish to simulate. Copies of pickup.0000* files for congeners PCB-28, 101, 153, and 180 are archived in:

/n/sunderland_lab/share_root/Lab/svante/helen/net_fs02_d2/MITgcm_ECCOv4/verification/global_pcb_llc90/results_pcb/

Submitting consecutive jobs to run an 80-yr PCB simulation 

The PCB simulation can be run for 80 consecutive years (1930 - 2010).  The simulation period (1930-2010) was broken into three stages based:

• Stage 1: 1930 - 1931. This is your "initialization". Simulation is initialized with zero PCB concentrations.
• Stage 2: 1931 - 1991. This is your "spin up". Each new job starts from where previous job ended. Jobs are broken into four-year chunks and repeat over meteorology from 1992, 1993, 1994, and 1995. Four years was chosen as the job length to account for leap years. 
• Stage 3: 1992-2010. This is your "science run". It's run as a single 20-yr job driven by meteorology for years 1992-2010.

The 80-yr simulation was broken into these three stages based on commonalities in the runscripts. In other words, you need one runscript to submit a job that is initialized from zero, you need another script to submit a 4-yr pickup job, and another script to submit a 20-yr job that picks up from where the spinup finished. 

You only need one data file and one data.pcb file to submit a single job to the queue. In the case of the 80-yr simulation, you need a data file for each stage and a data.pcb file for each stage and job:

data.stage1                          1930-31

data.stage2                          1931-91

data.stage3                          1992-2010

--

data.pcb.stage1                   1930-31

data.pcb.stage2_001           1932-35

data.pcb.stage2_002           1936-39

data.pcb.stage2_003           1940-43

data.pcb.stage2_004           1944-47

data.pcb.stage2_005           1948-51

data.pcb.stage2_006           1952-55

data.pcb.stage2_007           1956-59

data.pcb.stage2_008           1960-63

data.pcb.stage2_009           1964-67

data.pcb.stage2_010           1968-71

data.pcb.stage2_011           1972-75

data.pcb.stage2_012           1976-79

data.pcb.stage2_013           1980-83

data.pcb.stage2_014           1984-87

data.pcb.stage2_015           1988-91

data.pcb.stage3                   1992-2010 

In the data.stage* files, the user is specifying (1) how long the job is and (2) whether it's a pickup run. In the data.pcb.stage* files, the user is making sure PCBa_concFile, PCBp_depoFile, and PCBa_depoFile are pointing to the file paths to the correct years. For example, data.pcb.stage2_14 runs calendar years 1984-1987 and should look like:

 # Atmospheric concentrations (mol/m3) 

## 10 yrs of monthly means from GEOS-Chem

PCBa_concFile = './input_PCB/PCB28_popgconc_1984-1987_eccov4_llc90.bin',

# Atmospheric deposition (mol/m2/s)

## 10 yrs of monthly means from GEOS-Chem

PCBp_depoFile  = './input_PCB/PCB28_partdep_1984-1987_eccov4_llc90.bin',

PCBa_depoFile  = './input_PCB/PCB28_gasdep_1984-1987_eccov4_llc90.bin',

Please see simulations archived in

/n/sunderland_lab/share_root/Lab/svante/helen/net_fs02_d2/MITgcm_ECCOv4/verification/global_pcb_llc90/results_pcb/ for examples of data files from 80-yr three-stage runs. 

You are ready to submit your 80-yr run after modifying all the data.stage* and data.pcb.stage* files. There are four runscripts:

• runscript_stage1.csh : c-shell runscript to run Stage 1 
• runscript_stage2.csh : c-shell runscript to run Stage 2 
• runscript_stage3.csh : c-shell runscript to run Stage 3 
• run_stages123.csh    : wrapper to submit Stages 1, 2, and 3 to the queue all at the same time

Each stage can be run independently by submitting its own runscript to the queue (e.g., qsub runscript_stage1.csh to run Stage 1). If you want to submit all three stages to the queue at once, execute run_stages123.csh on the command line with:

./run_stages.csh

Disclaimer: These shell scripts were developed on Svante and have not been tested on Odyssey at the time of writing this post. 

(Helen M. Amos, 2 Jul 2016)

Analyzing and visualizing output with the gcmfaces Matlab package

The gcmfaces Matlab package forms the backbone of all Matlab scripts used to visualize and analyze PCB results from ECCOv4 simulations. Basic descriptions of gcmfaces and a users manual can be found on MITgcm ECCOv4 Odyssey Guide, Section 6.9 Final Remarks. Credit for the original gcmfaces modules goes to Gael Forget (MIT).

Purpose-built routines to plot observed and modeled PCB ocean concentrations and fluxes are in:

/n/sunderland_lab/share_root/Lab/svante/helen/MATLAB/gcmfaces/gcmfaces/gcmfaces_utility/

Included is this directory is a routine gcmfaces_benchmark_pcb.m, which is the driver for the PCB benchmark. The benchmark creates figures comparing modeled and observed seawater PCB concentrations to help assess model performance. 

Lizzie Lundgren has been expanding the plotting and analysis capabilities for PCB model output. One of the major efforts has been developing a code to calculate the PCB budget for each ocean basin. An archived copy of Lizzie's PCB routines is here: 

/n/sunderland_lab/share_root/Lab/svante/lizzie/pcb_budget

The pcb_budget repo is under Git version control and hosted on BitBucket. Contact Colin Thackray (thackray@g.harvard.edu) for an invitation to Bitbucket to get an update on the status of this work. 

(Helen M. Amos, 2 Jul 2016)

Additional documentation

I kept a work log in Evernote as a I developed and tested the PCB simulation in MITgcm ECCOv1 on Odyssey and then in MITgcm ECCO v4 on Svante. The work log is completely digital and word searchable. You can find a public copy here. This is a great resource of you are troubleshooting. 

(Helen M. Amos, 27 June 2016)

References

Dutkiewicz et al., 2012, Interconnection of nitrogen fixers and iron in the Pacific Ocean: theory and numerical simulations, Global Biogeochemical Cycles, doi: 10.1029/2011GB004039

Forget et al., 2015, On the observability of turbulent transport rates by argo: supporting evidence from an inversion experiment (2015), Ocean Science Discussions (under review)

 Forget,et al.,  (2015), Ecco version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation. Geoscientific Model Development Discussions (under review)

Friedman, C. & N. Selin (2016), PCBs in the Arctic atmosphere: determining important driving forces using a global atmospheric transport model, Atmospheric Chemistry & Physics, 16:3433-3448.

Y. Zhang et al., 2015, Biogeochemical drivers of the fate of riverine mercury discharged to the global and Arctic oceans, Global Biogeochemical Cycles, doi:10.1002/2015GB005124

Y. Zhang et al., 2016, Modeling methylmercury formation and trophic transfer in the ocean. in preparation.