Modeling photosynthesis of Spartina alterniflora (smooth cordgrass) impacted by the Deepwater Horizon oil spill using Bayesian inference

To study the impact of the Deepwater Horizon oil spill on photosynthesis of coastal salt marsh plants in Mississippi, we developed a hierarchical Bayesian (HB) model based on field measurements collected from July 2010 to November 2011. We sampled three locations in Davis Bayou, Mississippi (30.375°N, 88.790°W) representative of a range of oil spill impacts. Measured photosynthesis was negative (respiration only) at the heavily oiled location in July 2010 only, and rates started to increase by August 2010. Photosynthesis at the medium oiling location was lower than at the control location in July 2010 and it continued to decrease in September 2010. During winter 2010–2011, the contrast between the control and the two impacted locations was not as obvious as in the growing season of 2010. Photosynthesis increased through spring 2011 at the three locations and decreased starting with October at the control location and a month earlier (September) at the impacted locations. Using the field data, we developed an HB model. The model simulations agreed well with the measured photosynthesis, capturing most of the variability of the measured data. On the basis of the posteriors of the parameters, we found that air temperature and photosynthetic active radiation positively influenced photosynthesis whereas the leaf stress level negatively affected photosynthesis. The photosynthesis rates at the heavily impacted location had recovered to the status of the control location about 140 days after the initial impact, while the impact at the medium impact location was never severe enough to make photosynthesis significantly lower than that at the control location over the study period. The uncertainty in modeling photosynthesis rates mainly came from the individual and micro-site scales, and to a lesser extent from the leaf scale.


Introduction
Coastal wetlands are among the most productive ecosystems in the world, providing benefits such as carbon sequestration, storm protection, flood control, habitat, improved water quality, and recreational and aesthetic opportunities (Engle 2011, Jordan andPeterson 2012). However, coastal wetlands and their ecosystem services are not only threatened by upland stressors (Peterson and Lowe 2009), but they are also at risk from ocean-side stressors, in particular accelerated sea-level rise, erosion and anthropogenic pollutants like occasional oil spills.
The Deepwater Horizon (DWH) oil spill in the northern Gulf of Mexico (GOM), from 20 April to 15 July 2010, was the largest marine oil spill in the history of the petroleum industry (available at http://en.wikipedia.org, accessed on 16 March 2012). This accident released about 4.4 × 10 6 ± 20% barrels (7.0 × 10 6 m 3 ) of crude oil into the ocean from the well at 1500 m depth (Crone and Tolstoy 2010). Crude oil degrades in water via evaporation, dissolution, emulsification, sedimentation, biodegradation (microbial oxidation) and photooxidation (e.g., Hunt 1996, Fingas 1999, Plata et al 2008, Liu et al 2012 before coming ashore. On the Mississippi Gulf Coast, small oil slicks (<100 m 2 ) arrived sporadically from early June 2010 through 2011 but slicks arrived onshore more frequently and as larger patches during active tropical systems (e.g., Hurricane Alex, TS Bonnie and TD #5; http://www.nhc.noaa.gov/2010atlan. shtml, accessed on 16 March 2012).
The response of vegetation in coastal wetlands to weathered crude oil is complex and variable, ranging from short-term reductions in photosynthesis and rapid subsequent recovery, to complete mortality and long-term wetland loss (Pezeshki et al 2000, 2001, Roth and Baltz 2009, Engle 2011. Oil can affect plants directly by coating the leaves and blocking stomata, it's chemical toxicity disrupts plant-water relations or directly impacts the living cells, and it reduces oxygen exchange between the atmosphere and the soil with negative consequences for root health (Baker 1970, Pezeshki et al 2000, Ko and Day 2004. Spilled oil can also increase temperature stress, which combined with reduced photosynthetic gas exchange, will cause acute impacts (Baker 1970, Lin and Mendelssohn 1996, Ko and Day 2004. The acute impacts vary depending on the amount and type of oil, oiling frequency, the weather and hydrologic conditions, the species of plant, the spatial extent of oil coverage, season, soil composition and cleanup activity (Lin and Mendelssohn 1996, Hester and Mendelssohn 2000, Pezeshki et al 2000. Spartina alterniflora (smooth cordgrass) was more sensitive to partial oil coating than Juncus roemerianus (black needlerush; Pezeshki and DeLaune 1993), both plant species are important on the US Gulf Coast. In contrast, moderate oiling from the DWH oil spill had no significant effect on Spartina although it significantly lowered live above-ground biomass and stem density of Juncus in the Bay Jimmy, northern Barataria Bay, Louisiana (Lin and Mendelssohn 2012). Furthermore, S. alterniflora and Sagittaria lancifolia (bulltongue arrowhead) were shown to be more sensitive to oiling during the spring/summer growing season than during the pre-dormancy or dormant season in winter (Pezeshki et al 2000, Mishra et al 2012. This notwithstanding, salt marshes in the northern GOM region can recover rapidly DeLaune 1993, Pezeshki et al 2000).
However, oil spills can cause adverse chronic consequences on coastal wetlands and adjacent aquatic habitats due to oil penetration into the sediment and persistence in it (Krebs and Tanner 1981, Alexander and Webb 1987, Lin and Mendelssohn 1998. The severity of chronic impacts is influenced by the volume and chemical nature of the oil, the oil's amount of contact with and penetration of the soil, the physical nature of the coastline (high or low energy) and the composition of the plant community (Dicks and Hartley 1982, Ko and Day 2004. Once oil penetrates the sediment, the recovery to reference conditions may take 3-4 years (Alexander and Webb 1987, Mendelssohn et al 1993, Hester and Mendelssohn 2000 or longer (Bergen et al 2000, Michel et al 2009. Under extreme circumstances, recovery may never occur due to sediment removal (Baca et al 1987, Gilfillan et al 1995 or accelerated erosion after vegetation morality . Even after plants resume growth, oil may continue to adversely affect soil microbial functions in coastal wetlands (Burns andTeal 1979, Pezeshki et al 2001). Understanding the acute and chronic impacts of oil on coastal wetlands is important to be able to assess their resilience.
Net primary productivity is an aggregate measure of energy available to support ecosystem structure and function (Cardoch et al 2002); therefore, photosynthesis can be used to assess the impact of oil on ecosystem health. However, inferring how photosynthesis is affected by chronic and indirect exposure to oil can be difficult because of spatial and temporal complexity. On the one hand, factors at multiple spatial scales (from cm to km) will influence photosynthesis, from individual-plant level traits (cm), to site-specific variation (m), to the broader landscape level (km). On the other hand, photosynthesis changes over time as a response to seasonal patterns, with more photosynthesis during the growing season and less in the dormant season. The difficulty inherent in relating the change in photosynthesis rates directly to a specific oil spill will increase with the scale and complexity of the ecosystem as well as additional disturbance/stress factors such as sea-level rise (Fabricius and De'Ath 2004) unrelated to the oil spill event. A coherent way to link these various processes occurring simultaneously at multiple spatial and temporal scales is therefore necessary to better elucidate the oil-related impacts on photosynthesis during the chronic stress phase of recovery.
Recent advances in computational statistics have produced new tools for inference and prediction, such as the hierarchical Bayesian (HB) modeling approach. This approach has the capacity to exploit diverse sources of information, can accommodate uncertainties, and has the ability to draw inferences on large numbers of latent (inferred) variables and parameters that describe complex relationships (Clark 2005). The approach accommodates complex systems by decomposing their high-dimensional relationships into levels of conditional distribution within a consistent framework: data level, process model level and parameter level (Clark et al 2001). In particular, the process model is made stochastic, consistent with the system complexity, and our limited understanding of the processes (e.g., photosynthesis in our study). The stochasticity allows us to impose conditional independence at the data level and assimilate the relationships among state variables at the process level (Wu et al 2010). Finally, the simulation results are rich (e.g., posterior distributions instead of point estimates) and the interpretation is simplified because uncertainty can be readily quantified using credible intervals from the resulting posteriors.
In this study, we aim to (1) Develop an HB statistical model to simulate photosynthesis of S. alterniflora plants at three locations representative of a range of oil spill impacts based on field measurements.
(2) Apply the model to infer when photosynthesis of S. alterniflora at the oil impacted locations recovered to the status of the control location with inherent leaf, individual, micro-site and seasonal variability accounted for, and which spatial scale contributed the most to variability in the photosynthesis simulations.  figure 1(A)). One location was not impacted by crude oil (control location, figure 1(B)), the second location experienced medium impact by crude oil (some oil observed on the sediment and plants, figure 1(B)) and the third location was heavily impacted by oil (plants covered by crude oil extensively, figure 1(B)). The spatial extent of the heavy impact area was not larger than about 10 m × 5 m (50 m 2 ). Mississippi Gulf Coast's climate is subtropical maritime (Eleuterius 1998), with mean spring, summer, fall and winter temperatures being 19.4 • C, 27.4 • C, 20.4 • C and 11.1 • C, respectively (Reuscher 1998). The area receives an average of 1488 mm of rainfall per year (Eleuterius 1998) with the Hurricane season being from June through November.

Field monitoring
We randomly chose 10-15 individuals of S. alterniflora at each of the three locations monthly between July 2010 and November 2011 and measured the photosynthesis rate on the middle portion of the second leaf from the tip of each individual using an LI-6400 XT portable photosynthetic Figure 2. The conceptual model to illustrate the hierarchical structure with complexity decomposed into stages of data, process and parameters (vertical direction) and the association of different spatial scales (horizontal direction). (Adapted from Clark and LaDeau (2006) and McMahon and Diez (2007).) gas exchange and chlorophyll fluorescence system. The photosynthetic rates for each of the leaves were measured by CO 2 exchange under five identical photosynthetically active radiation (PAR) intensities: ∼400 µmol photon m −2 s −1 , 800 µmol m −2 s −1 , 1200 µmol m −2 s −1 , 1600 µmol m −2 s −1 and 2000 µmol m −2 s −1 to standardize the photosynthesis rate under saturating light intensities. For each leaf we also measured chlorophyll fluorescence (F v /F m ) after at least 30 min dark adaption and used this as an indicator of leaf and individual stress or health (Naumann et al 2007). Air temperature was recorded by the LI-6400 system at the time of measurement.
No summer photosynthesis data was collected at any locations in 2011 because of an instrument exposure to salt-water that required extensive repairs. Data were not collected at the control location in August and September 2010 since we assume that photosynthesis and its variability are similar to those in July, as the daily average air temperature for the three months were similar (28.3 • C in July, 28.9 • C in August and 26.7 • C in September). In order to assess the recovery of photosynthesis at the impacted locations (medium, heavy), we assume that their close proximity to the control location provided similar environmental conditions.

Developing the spatio-temporal hierarchical Bayesian model
In order to model photosynthesis, we started with a saturation function to simulate the response of photosynthesis to PAR intensity, the photosynthesis-irradiance (P-I) curve. We made the model stochastic to reflect the fact that there is inherent natural variability in the saturation function that describes photosynthesis versus PAR. Furthermore, we acknowledged that the variability among individual plants, micro-site variability and temporal variability can all affect photosynthesis in unpredictable ways.
The recognition of these multiple sources of natural spatial and temporal variability guided the development of the process and parameter levels in our HB model (figure 2 and  table 1). The goal of the model was to estimate photosynthesis rates as a function of (1) the covariate at the leaf scale: PAR, (2) the covariates at the individual scale: air temperature and stress/health level, (3) the covariate at the micro-site scale: the level of crude oil contamination, and (4) the covariate at the temporal scale: the number of days after the beginning of the crude oil impact (8 July 2010).
To represent these functions, let i represent month, j represent location, k represent an individual, l represent PAR intensity (unit: µmol photon m −2 s −1 ), P ijkl represent the photosynthesis rate (µmol C m −2 s −1 ) at PAR intensity l for individual k at location j at month i (as we have photosynthesis measurements for each month). Let P ijkl .µ represent the mean of P ijkl , σ 2 L represent the variance of the photosynthesis rates among the 5 different PAR values (∼400, 800, 1200, 1600, 2000 µmol photon m −2 s −1 ) at the leaf scale. N represent a normal distribution, ∼ represent 'is distributed as', so P ijkl is modeled using equation (1): (Note: photosynthesis can be negative when the leaves were covered by crude oil.) Air temperature (unit:

PAR
Photosynthetically active radiation (unit: µmol photon m −2 s −1 ) We applied a Michaelis-Menten equation to model P ijkl .µ as a function (f L ) of available PAR at the leaf scale (equation (2)): where Pmax ijk (µmol C m −2 s −1 ), PAR0 ijk (µmol photon m −2 s −1 ), and ω (µmol photon m −2 s −1 ) are the three parameters that need estimation. Pmax ijk represents the maximum net photosynthesis rate, PAR0 ijk represents the PAR intensity where net photosynthesis becomes positive (called the compensation irradiance), and ω represents the PAR intensity when the photosynthesis is half of the Pmax ijk (assuming PAR0 ijk is 0; called half-saturation irradiance).
We reparameterized the nonlinear Michaelis-Menten equation to make it linear so we could derive the posteriors more efficiently , Clark 2007. To do this, the new variable that replaces PAR ijkl is z ijkl = PAR ijkl ω+PAR ijkl , so that P ijkl .µ can now be modeled as a linear function of z ijkl (equation (3)): The new parameters β 0ijk and β 1 represent the intercept modeled from the individual scale and the coefficient for the transformed variable z ijkl , respectively, and they are related to the original parameters through equations (4) and (5): Therefore, P (photosynthesis rates) for n 1 PAR intensities per individual, n 2 individuals per location at n 3 locations over n 4 months is modeled as in equation (6): where ∝ denotes 'is proportional to'. At the individual scale, where β 0ijk .µ is the mean of β 0ijk , and σ 2 K is the variance at the individual scale. In equation (7) β 0ijk .µ is a function of air temperature (T ijk ) and stress level (S ijk ) (measured as F v /F m ) for individual k at location j at month i, as defined by equation (8): where α 0ij is the intercept modeled from the micro-site scale (among the 3 locations), α 1 and α 2 are the coefficients for T ijk and S ijk respectively, Therefore, β 0 is modeled as equation (9): At the micro-site scale, α 0ij is modeled using equation (10): where α 0ij .µ is the mean of α 0ij , σ 2 MS is the variance at the micro-site scale. In equation (10), α 0ij .µ is modeled as a function of the days after the initial oil impact (Time i ), and two dummy variables representing the locations that have medium (Loc m ) and heavy impact (Loc h ) from crude oil contamination as defined in equation (11): where θ 0 is the global intercept, θ 1 is a vector of coefficients for the days after the initial impact (Time) at each location (Loc), θ 2 is the coefficient for the dummy variable of medium location Loc m , and θ 3 is the coefficient for the dummy variable of heavy location Loc h . The dummy variables Loc m and Loc h will only take the values of 0 or 1. If Loc m is 1, then α 0ij .µ represents the mean of α 0ij for the medium impacted location; if Loc h is 1, then α 0ij .µ represents the mean of α 0ij for the heavily impacted location; if both Loc m and Loc h are 0s, then α 0ij .µ represents the mean of α 0ij for the control location; at no time can Loc m and Loc m be 1 at the same time. Therefore, α 0 is modeled as equation (12): To complete the Bayesian model, we require prior distributions for unknown parameters (αs, βs, θ s and σ s) described above. We used priors that were conjugate with the likelihood for computation efficiency (Calder et al 2003), meaning that prior and posterior distributions had the same form. For instance, the variance parameters σ 2 L , σ 2 K , σ 2 MS had inverse gamma (IG) priors, which is conjugate for the normal likelihood distribution. For the other parameters, we used normal distributions as priors. The priors are 'non-informative', meaning that the prior distributions were rather flat and only weakly influenced the parameter estimates since we knew little about the unknown parameters (Hartigan 1998, Clark andBjornstad 2004).
By combining the parameter, process, and data models as defined above in equations (1)-(12), we derived the joint posterior (equation (13)): p(β 0 , β 1 , α 0 , α 1 , α 2 , θ 0 , θ 1 , θ 2 , θ 3 , σ 2 L , σ 2 K , σ 2 MS | P, PAR, T, S, Time, Loc m , Loc h , priors) From the model just described, we could infer (1) the impact of the PAR intensity (β 1 , ω), air temperature (α 1 ) and stress level (α 2 ) on photosynthesis rates, (2) the temporal trend of photosynthesis at each location (θ 1 ), (3) how photosynthesis at the impacted locations compared to the control location over time (α 0 ), and (4) which spatial scale contributed most to variability in photosynthesis rates (τ ). In particular, we are interested in using the model to infer when the photosynthesis at the impacted locations started to recover to similar rates as the control location, and which spatial scale contributes the most to variability in the photosynthesis simulations.

Implementing the hierarchical Bayesian model
We implemented the model described in equation (13) using Markov Chain Monte Carlo (MCMC) simulations (Gelfand and Smith 1990) in the software OpenBugs 3.2.1 (Rev. 781). We evaluated convergence of the parameters (αs, βs, θ s and σ s) by simulating three MCMC chains from three different sets of initial values (figure S1 available at stacks. iop.org/ERL/7/045302/mmedia). The chains converged based on Gelman and Rubin's convergence statistics (Brooks and Gelman 1998) and convergence required 90 000 iterations of MCMC. These pre-convergence 'burn-in' iterations were discarded and an additional 90 000 iterations were saved for the subsequent analysis. The upper 95% confidence limit of the convergence statistics for the parameters, latent variables and variances based on the post-burn-in iterations are all less than 1.2, indicating that the chains converged to the stable posterior distribution (Clark 2007).
We used an index of agreement (IoA) to assess the match between simulated medians of photosynthesis and the measured photosynthesis (equation (14): Willmott et al 1985Willmott et al , 2012: where Pr i is the model predictions, O i is the observed (or measured) data,Ō is the average of observed (measured) data, and n is the number of data points. IoA measures the agreement between predictions and observations on a one-by-one basis. This dimensionless index ranges from 0.0 (indicating no agreement) to 1.0 (perfect agreement). It approaches 1.0 slowly as Pr approaches O, and therefore, provides a reliable tool to assess the model's performance.

Modifying the hierarchical Bayesian model
We also explored possible improvements to the model by modifying equation (11) to include Julian days as a covariate to represent seasonality (equation (15)): For the effect of Julian day, we wanted a model that would have a single peak value to represent the highest photosynthesis rates in the late spring or summer (between Julian days 120 and 240), so θ 5 should be positive in equation (15).
In particular, we assessed the significance of the coefficients for the added covariates and how the added covariates affect the coefficient estimates for the other covariates.

Measured photosynthesis rates
Not all field measurements obtained were used in the model; 17.7% of the data were removed from our data analysis and model development. Some examples of data that were exclude are: too small a leaf area in the leaf chamber, negative values which could not be explained by health status of the leaves (as indicated by the field notes and F v /F m measurements), or the response of photosynthetic rates being much smaller under high PAR levels compared to low PAR levels but do not represent photoinhibition. At the control location, measured photosynthesis rates were high (41.29 ± 13.17 (standard deviation) µmol C m −2 s −1 ) compared to medium (26.27 ± 10.71 µmol C m −2 s −1 ) and heavy impact locations (−1.53± 0.783 µmol C m −2 s −1 ) in July 2010 (figure 3). Late fall and winter (Nov 2010-Feb 2011) photosynthesis rates were low (15.84 ± 6.33 µmol C m −2 s −1 ), and began to increase in March 2011 at the onset of the growing season of S. alterniflora ( figure 3).
Photosynthesis at the medium oiled location was lower than at the control location in July 2010 and the mean rate  3). Measured photosynthesis was negative (i.e., respiration only) at the heavily oiled location in July 2010 and photosynthesis rates at this location started to increase in August and September 2010. At that time, new leaves were observed to grow back from the roots unaffected by oil toxicity and replace dead leaves coated by crude oil. During the winter months of 2010-2011, the contrast between the photosynthesis rate at the control location and the two impacted locations was not as obvious as just after the oil impact in the summer of 2010 (figure 3). S. alterniflora undergoes senescence of above-ground tissues in winter as it is not an evergreen (Edwards and Mills 2005). Photosynthesis increased through spring 2011 at all three locations during the period of active new growth, no summer data were available due to instrument malfunction, and photosynthesis decreased again starting in October at the control location, and September at the medium and heavily impacted locations.

Model comparison.
The simulation results showed that the coefficients for the quadratic term of Julian day θ 5 in equation (15) was significant larger than 0 (95% credible interval, 2.5% quantile to 97.5% quantile: 1.7555 × 10 −6 -5.5201 × 10 −4 ), meaning that seasonality represented by Julian day significantly affected simulated photosynthesis rates. Meanwhile the coefficient for air temperature α 1 was not significantly different from 0 (95% credible interval: −0.4977-0.1150), meaning that air temperature did not significantly affect simulated photosynthesis once the seasonality has been accounted for by Julian day. When the original equation (11) was applied, the 95% credible interval of the coefficient for air temperature α 1 was positive (table 2) meaning that air temperature significantly affected simulated photosynthesis. These indicated that air temperature accounted for seasonality in the original model while Julian day accounted seasonality in the revised model. However, the penalty due to the inclusion of more parameters related to Julian day made the revised model have a larger Bayesian deviance information criterion (DIC = 8452) compared to the original model (DIC = 8430), suggesting that Julian day was not a better explanatory covariate for seasonal photosynthesis response than air temperature alone. The posteriors of the coefficient for the quadratic term of air temperature (α 3 ) in equation (16) (95% credible interval: −0.04351-0.01788), and the coefficient for the interaction of air temperature and stress level (α 4 ) in equation (17) (95% credible interval: −1.2685-0.9009), ranged from negative to positive. This indicated there were no significant nonlinear effects of air temperature, or the interaction of air temperature with health/stress level. Therefore, we believe the air temperature in the original model effectively accounted for temperature variability in the simulated photosynthesis rates over the different seasons. Since no modification improved the original model, the following simulation results were based on the original model (equation (13)).

Model performance.
The density distributions of the medians of simulated photosynthesis rates from the HB model matched those of the measured photosynthesis well across the three locations, especially at the heavy impact location (figure 4). The good match between the simulated medians and the measured data was also shown in the IoA, which was 0.835 at the control and medium oiled locations, and 0.842 at the heavily oiled location. The posteriors of the simulated photosynthesis captured most of the variability of measured photosynthesis (figure 5). These all suggested the HB model performed well in simulating the observed photosynthesis.

Simulation posteriors.
The photosynthesis rate at the control location showed a decreasing trend from July 2010 to November 2011 as the 95% credible interval of the parameter (θ 1 [control]) was negative (table 2 and figure S2 (available at stacks.iop.org/ERL/7/045302/mmedia); note, missing measurements August and September 2010 and June-August 2011). At the two impacted locations, there was not an apparent trend of the simulated photosynthesis rates over the same time period, though the medians of the change rates were positive at the heavy and negative at the medium locations. The 95% credible intervals for θ 1 [control] and θ 1 [heavy] did not overlap, showing a significantly different temporal trend (table 2) at the control and heavy locations. The overlapping of the 95% credible intervals between the medium and control locations was small, and there did not show a significant difference in temporal trend between the medium and control locations or the medium and heavy locations.
The medians for photosynthesis rate at the leaf scale (P) and the intercept at the individual scale (β 0 ) were of a similar magnitude (20.70 and 22.62), and the median of the intercept from the micro-site scale (α 0 ) was about twice as large as those at finer scales (39.38). In terms of variability, the individual scale (τ K ) showed lower precision (inverse of variance, see table 1) and precision/median ratio (normalized precision) compared to the leaf scale (τ L ), which indicates that individual variability contributes more to the simulated photosynthesis rates than does leaf scale variance (table 2 and figure S5 available at stacks.iop.org/ERL/7/045302/mmedia). At the micro-scale (τ MS ), which includes both spatial and temporal variability, the precision and the precision/median ratio ranged from low to high (table 2). The 97.5% quantile of precision/median at the micro-site scale (τ MS ) is lower than 2.5% quantile of precision/median at the leaf scale (τ L ) (table 2), showing the normalized precision at the leaf scale is significantly higher than that at the micro-site scale. These indicate the uncertainties in modeling photosynthesis rates mainly came from the variability at the individual and micro-site scales, less from the leaf scale.
3.2.4. The monthly trend of the difference between the impacted locations and the control location. The monthly trend of photosynthesis modeled by the HB at each location can be summarized using the intercept at the micro-scale, α 0 (figure 6 and equations (10)- (12)). The value of α 0 was different from the actual photosynthesis as it was standardized for the covariates of PAR intensity, air temperature and F v /F m . At the beginning of the study (after the acute oil impact in July 2010), the 95% credible intervals of α 0 for the control (−37.893, −4.839) and heavy impact locations (−72.532, −40.283) did not overlap. This indicates that there existed significant differences in photosynthesis between these locations. The medians at the heavy (−55.415) and medium (−32.540) locations were smaller than at the control location (−20.265) in July 2010, indicating that median photosynthesis was lower at the two oil impacted locations than the adjacent control location. Over the course of the study, the amount of overlap of the 95% credible intervals among the three locations increased through time, and the medians of the photosynthesis rates at the two impacted locations approached the control location suggesting recovery by the plants.

Discussion
Studies of the impact of oil spills on salt marsh plants have been focused on Louisiana due to more frequent oil spill events and large marsh areas there (Pezeshki and DeLaune 1993, Pezeshki et al 2000, Hester and Mendelssohn 2000, Ko and Day 2004, Mishra et al 2012. During the DWH oil spill, marsh shorelines and not the marsh interior were exposed to the weathered oil . Limited quantitative data are yet available on the impact of DWH oil spill on salt marsh vegetation. Mishra et al (2012) assessed the ecological impact on the salt marshes along the southeastern Louisiana coast using photosynthetic capacity and physiological status through satellite and ground truth data, and found extensive reduction in photosynthetic activity during the peak of the growing season in 2010. Lin and Mendelssohn (2012) documented variable impacts depending on oiling intensity in the Bay Jimmy, northern Barataria Bay, Louisiana. As of the fall of 2011, many of the most heavily oiled shorelines had minimal to no recovery . In Mississippi, some salt marshes experienced crude oil impact, such as in Davis Bayou, Grand Bay, Garden Pond at Horn Island and in Waveland, but the impact and recovery have not attracted enough attention or been documented well. This paper enhances our understanding of the impacts of the DWH oil spill on photosynthesis of the dominant marsh plant, S. alterniflora, from Mississippi.
To improve our understanding of oil's impact on coastal salt marsh plants, we used an HB modeling approach because it can assimilate data and account for variability at multiple spatial and temporal scales in a coherent framework (e.g., McMahon and Diez 2007, Clark et al 2011, Wilson et al 2011). An HB model has been applied to monitor long-term harbor seal abundance changes in Prince William Sound, Alaska impacted by the 1989 Exxon Valdez oil spill (ver Hoef Figure 7. Photosynthesis difference between the impacted locations ((A) heavy location; (B) medium location) and the control location over time. The thin lines represent 2.5% and 97.5% quantiles, and the think lines represent medians of the difference. The vertical space between 2.5% and 97.5% quantiles shows the 95% predictive intervals. The gray line represents the horizontal line at 0 and vertical line at day 140 when photosynthesis at the heavy oiled location stated to recover to the level at the control location. and Frost 2003). However, HB models have not been applied to assess the impact of oil spill's on coastal salt marsh.
Other types of saturation curves and mechanistic models are available to model photosynthesis (e.g., Jassby and Platt 1976, Chalker 1981, Thornley and Johnson 1990, Farquhar et al 1980, Farquhar and von Gaemmerer 1982, Farquhar et al 2001, Marino et al 2010. The application of the Michaelis-Menten function adequately accounted for the relation between photosynthesis and PAR intensity, and facilitated the inference on the temporal trend of photosynthesis and its contrast between the impacted and the control locations, as the leaf scale contributed the least to the variability in the photosynthesis simulations. During the winter months in 2010, the contrast between the control and the two impacted locations in the field data was not as obvious as during the growing season earlier in 2010. This could be explained either by complete recovery of photosynthesis from oil impacts, or it could be due to seasonally lower air temperature being a limiting factor keeping photosynthesis uniformly low. Due to the uncertainty surrounding the interpretation of the field data, we relied on the HB model to derive a temporal pattern of simulated photosynthesis at the three study locations standardized for seasonality, PAR intensity, and individual variability. The linear coefficient θ 1 in the model may oversimplify the temporal trend at each location. However, it captured the general trend of photosynthesis change, although smoothing over small fluctuations (Ver Hoef and Frost 2003).
In particular, we used the HB model to predict the time course of photosynthetic recovery at the impacted locations to values similar to the control location. In order to derive when photosynthesis had recovered to the status of the control location, a finer temporal resolution than month would be required. We applied the posteriors of θ 0 − θ 3 and τ MS in equations (10)-(12) from our HB model to generate the predictive α 0 at each day after the initial impact at each location. Then, the difference of α 0 between heavy and control locations and between medium and control locations were calculated. The 2.5% and 97.5% quantiles and medians of the differences of each day were derived (figure 7). The generally increasing trend over time in the differences between the impacted and control locations was largely due to the decreasing trend at the control location. If the 95% credible intervals of the difference included 0, then the photosynthesis rates between the impacted and the control locations did not show a significant difference. Based on this, we could derive the photosynthesis rates at the heavy location recovered to the status of the control location about 140 days (4.7 months) after the initial impact, which was in early December 2010 ( figure 7(A)). On the other hand, the oil impact was never severe enough to make the photosynthesis rates at the medium location significant lower than that at the control location (figure 7(B)).
There are two possible reasons photosynthesis recovered so quickly. First, the small concentration and very patchy distribution of crude oil associated with salt marsh plants (Biber et al in review) showed the roots were possibly free of contamination, which led to quick regrowth of leaves at the heavily oiled location. At the beginning of the crude oil impact (July 2010), the concentration of oil range organics (C19-C36) ranged from 37 000 to 48 900 mg kg −1 (ppm) on the plants and 32.4-44.2 ppm in the sediment at the heavy impacted location. The detailed chemical analysis of saturated hydrocarbons on the oil collected confirmed that it most likely came from the BP DWH spill as the GC-FID chromatograms showed a similar distribution pattern of resolved peaks and unresolved complex mixture compared to the MC252 reference oil source from BP (Liu et al 2012). The concentration dropped to 686-1070 ppm on the plants and <50 ppm in the sediment at the impacted location by October 2010. Second, quick natural degradation and possible physical removal by waves and tides at the study area reduced the amount of crude oil that washed ashore in the salt marsh. In particular, areas exposed to the predominant southeasterly winds, common during the summer in our study area, experienced substantial tidal 'cleaning' with oil being removed during high tides and periods of strong wave action corresponding to tropical storms in the GOM. Warm temperatures, often exceeding 35 • C during the day are also common in the region from June through September, promoting rapid volatilization of lighter carbon fractions which correlates to reduced toxicity (Irwin et al 1997).

Synthesis
Our study demonstrated a relatively new approach to assess ecological recovery-the HB modeling approach. It shows a promising tool to assess wetland's resilience to disturbance by accounting for uncertainties from different sources at different spatial and temporal scales in a coherent framework. The application of our HB model facilitated a better understanding of oil impacts to S. alterniflora and generated inference we could not have obtained from the empirical data alone. For example, (1) air temperature and PAR positively influenced photosynthesis, whereas the leaf stress level negatively affected photosynthesis, (2) the overall temporal changes in photosynthesis rates with standardized covariates over 17 months had a negative trend at the control location, and ranging from negative to positive trend at the impacted locations, (3) the photosynthesis rates at the heavily impacted location recovered to the control location about 140 days after the initial impact whereas the impact at the medium location was never large enough to make photosynthesis significantly lower than that at the control location, and (4) the uncertainties in modeling photosynthesis rates mainly came from the variability at the individual and micro-site scales, less from the leaf scale. However, even though photosynthesis had recovered to the status of the control location, we recommend continuous monitoring of photosynthesis rates as (1) longer data will help determine if the decreasing trend at the control location persists and help to explain it, (2) the oil residual buried in the sediment at the fine scale may affect the plant root zone chronically and pose a long-term threat, and (3) accumulating baseline data for assessing the next possible disturbance. Thus, our model approach and results can guide our future sampling efforts in photosynthesis of oil impacted salt marsh plants.