Oil and gas, along with agriculture, form the backbone of the Alberta economy. Owing to the wide range of geography in the 2 industry sectors, there is an overlap of activities virtually everywhere in the agricultural regions of the province. Sour gas is natural gas containing excessive quantities of hydrogen sulfide (H₂S) and other corrosive compounds. This gas must be processed (cleaned) before it can be delivered to end-users. Following processing, extracted and unwanted components are mixed with gases of high heating value and combusted within, or on top of, a stack. The emitted plume is designed to disperse atmospheric pollution in concentrations unlikely to exceed the legislated maximum ground level concentrations (1,2). As of July 1995, there were 231 gas plants licensed in Alberta to process sour gas and emit combusted by-products into the atmosphere.

For more than 30 y, public concerns have been expressed over perceived negative health effects relating to the Alberta sour gas industry (3,4). These concerns have referred to all spectrums of emissions; from acute, uncontrolled releases of raw sour gas (for example, well blowouts, pipeline and facility failures), through to controlled releases of low-level combustion emissions from sour gas processing plants.

Investigations have been commissioned in response to human health concerns over acute events, such as well blowouts (5). Public advisory boards have made recommendations with respect to set-back distances from sour gas facilities to human residences. Emergency response protocols have been developed to deal with both industrial and residential health hazard risks from acute events (6). Spitzer et al (7) conducted a community health study comparing health outcomes in several communities near and far from sour gas processing facilities. Recently, government regulators have required baseline surveys of both human (8) and animal (9) health in communities located near very large, newly approved facilities.

In terms of animal health studies, several investigations have been conducted into acute blowout events (10,11). Church (10) followed an acute exposure event (a sour gas well blowout), in west-central Alberta. He found several common health complaints (for example, respiratory) that he argued could just as easily be attributed to other confounding factors; such as, weather, management practices, and genetic differences. Coppock et al (11) likewise reported complaints of upper (ocular, nasal) and lower (pneumonia, exercise intolerance) respiratory problems in food animals and horses located downwind of a sour gas well blowout in east-central Alberta. These authors concluded that increased rates of pneumonia and decreased feeding performance could reasonably be assumed to have some relation to the blowout emissions. Small-area longitudinal health surveys have been conducted on swine (12) and beef cattle (13,14). Lore et al (12) conducted a field trial with several generations of swine housed upwind and downwind from a sour gas processing facility. These authors found no difference in fertility, piglet survival, or infectious disease rates between the 2 exposure groups. The 2 studies by Waldner et al (13,14) represent some of the first prospective veterinary epidemiological studies to be published on the subject. While representing data from only 7 beef cattle operations in a small area of the province, the studies proceeded in a methodical and longitudinal manner to provide some of the first insights into health and productivity (particularly reproductive) outcomes associated with oil and gas industry activities. These authors found increased risks (in several years of a multi-year study) for stillbirth associated with elevated sulfation levels on pasture (attributable to a wide range of sources of reduced sulfur compounds, including sour gas processing facilities) (14). In addition, Waldner et al (13) identified an increased risk for non-pregnancy associated (in some years) with proximity to certain oil field facilities. These 2 papers examined a wide range of potential oil and gas industry pollution exposure sources, resulting in a cumulative index of exposure based on sulfation, H₂S deposition, or both, (14) and proximity and density of the sources (13).

A province-wide investigation into the effects of chronic exposure to sour gas processing plant emissions on cattle health and productivity has never been conducted. In response to public concerns, and in an attempt to address unanswered questions, a province-wide longitudinal study of cattle health and productivity in relation to chronic levels of sour gas exposure was undertaken by this research team. As part of a large-scale retrospective field investigation into both dairy and beef cow-calf health and productivity (15), we assessed historical exposures to sour gas processing plant emissions on 1372 dairy farms. The objective of the study was to determine whether or not there is evidence of an impact of sour gas processing plant emissions on the retainment, survival, or both, of cows in dairy herds in Alberta. For the purposes of this study, retainment refers to annual herd-level risks for culling, mortality, or both, while survival refers to productive longevity of individual cows in the herd until these cows are either culled or die.

All individual dairy cows whose birthplace could be traced to the study farm in question, and for whom a lifetime sour gas exposure estimate profile could be established, were included in this study. This data set was combined with a data set listing the details of the last lactation for sold and dead cows from 1985 to 1994. The outcomes considered were either time (months) to culling or time (months) to mortality. A mortality event is of the same consequence to the inventory of the herd milking string as is a culling event. In addition, the only censored events remaining with this approach to the analysis included right-censored end-of-study cows and cows sold for dairy purposes. Time was measured from the month of the first calving for both culling and mortality.

This individual outcome variable was transformed to reflect the actual life period of interest. That is, a survival time (failing to be culled or die, during a given month) of 1 was assigned to the first month following the animal's first calving date. Each month of subsequent survival within that herd (unless censored) added an extra month (+ 1) to the survival time. The distribution of these modified survival data was assessed using software (SAS, version 6.12 Proclifetest; SAS Incorporated, Cary, North Carolina, USA) (18). The natural logarithm of the survival time was then regressed against time (months) and the R² (coefficient of determination) value was assessed. For R² values near 1, an exponential distribution of the survival data could be assumed and a parametric exponential survival distribution could then be utilized in this analysis (19).

Monthly test day summary records of the number of animals on each farm were provided by DHI. This data set was used to compute the annual average number of animals at risk for culling or mortality within each study herd. The variable representing the number of cows in the herd reflected the number of animals in the herd on the previous test day, minus the number of animals that were sold or died, plus new animals entering the herd (either freshening primiparous heifers or purchased stock). Thus, the annual average number of animals at risk reflected the sum of the number of animals in the herd on each test day, divided by the number of test days within a calendar year.

The test day summary files also provided data on the number of animals entering and leaving each herd. The exit value was used to test the validity of the number of animals leaving each herd, from voluntary sale, culling, or mortality during a calendar year, computed from the individual cow previous lactation data set. Where the number of animals recorded as having left the herd (for example, the exit value) did not equal the sum of animals culled, plus animals that died, plus animals sold for other purposes, then the herd-year record was deleted if either: 1) the sum of the 3 categories of departure from the herd exceeded the exit value, or, 2) the sum of the 3 categories fell short of the exit value by more than 5% of the total. For the purposes of this study, only animals that were culled or died were considered as having experienced an event. Those animals that were sold for export, dairy purposes, or were present in herds at the time that the owner had unsubscribed to DHI services, were not considered as having been culled or died. However, for all months during which they were present in the herd, these sold animals were considered to have been at risk for either culling or mortality. The new locations for dairy cows sold to other dairies could not be determined. Annual culling and mortality risks were computed for each farm (criteria specified above) from the number of cows culled (or died) divided by the average number at risk during the study year, respectively.

Sour gas processing facilities are not distributed evenly across Alberta. The highest concentrations of sour gas are found along the foothills, encompassing regions of montane and boreal forest. Because the soil type, vegetative covering, and dominant climatic conditions of these regions were expected to impact on the health and productivity of dairy herds, the potential confounding effects of agro-ecological region were incorporated into all analyses. A 4-level indexed map of Alberta was provided by Agriculture Canada (20). The regions represented by this index included a) boreal forest, b) montane, c) grassland (prairie), and d) parkland. The regional index was appended to each dairy herd data set using the GIS.

The geographical location, emissions designs (for example, stack heights, diameters, exit speeds, and temperatures), and monthly emissions (SO₂ in tonnes) for each of the 231 sour gas plants in Alberta from January 1, 1985 through to December 31, 1994, were provided by Alberta Environmental Protection. The data were cross-checked with equivalent sources from the Alberta Energy and Utilities Board (AEUB) to ensure maximum validity of the working data sets.

Meteorological information was accessed specifically for the purposes entailed in operating the USEPA ISCLT3 dispersion model (see below) (21). The dispersion model file input requirements dictated specific weather data file formats which were provided by Environment Canada, through a network of 14 surface (usually airports) and 1 upper air weather stations in Alberta. Joint frequency distribution data were utilized from 1985 to 1994 for wind speed, stability, and direction. Surface air temperature data and upper air mean afternoon mixing heights were derived from 25-year averages also compiled by Environment Canada (22,23). The mixing heights are used to estimate the height of the ‘boundary’ layer under which vertical mixing of pollutants with atmospheric gases most readily occurs. The height tends to be lower during the winter months (23).

Polar grid plot files from each monthly dispersion model run were converted into latitude (from Y) and longitude (from X) coordinates within the GIS (Spans-Gis version 5.3; Intera Tydac Technologies, Nepean, Ontario). The GIS was subsequently used to interpolate exposure isopleths for each point source. Map levels were coded to represent 0.1 μg/m³ of SO₂ per classification level. All individual gas plant exposure estimate maps for a given month and year were combined using an additive function in the GIS; likewise, annually averaged maps were created. Each of the monthly and annual maps were then assigned to a herd file for each of the dairy farm locations.

An asymptotic exponential decay function, similar to that present in a Gaussian plume dispersion model ground-level concentration profile, was approximated for each point source, while ignoring the impact of wind speed and direction as well as gas plant engineering specifications. Estimating ground-level concentrations (as opposed to concentrations at 2 m above ground) is appropriate for grazing animals, such as cattle and sheep (as opposed to humans). By ignoring wind speed and direction, and plant engineering specifications, we assessed the joint effect of proximity to, and density of, recorded emission levels from processing facilities within a specified radius of all sources. This is in direct contrast to the PD model (above) which fully accounted for wind speed, direction, and plant emission engineering features. A direct comparison of the methods used, and the exposure estimates obtained from these 2 models, has been published elsewhere (27). The exposure data also have been applied in a parallel study of sour gas emissions and dairy cattle reproduction (28).

The potential-mapping method was applied to the known mass of SO₂ emitted monthly from each gas plant, to a distance of 70 km from each plant. The distance of 70 km was empirically determined to be the distance at which SO₂ concentrations approached 0 under the PD model. Farm-site exposure estimates from the potential-mapping model are hereafter termed DD. The exposure estimate contours (classifications) surrounding each gas plant were dependent upon the monthly average emission rate, plus the influence of any other plants falling within a 70 km radius of the estimated surface points (up to a maximum of 32 other plants). The map classification scheme provided for 200 equi-spaced exposure intervals had a minimum value of 0 and a maximum of 2000 (dimensionless on an interval scale). The units were not directly comparable to the PD model estimates. The DD estimates were assigned to the appropriate dairy farm sites using the GIS.

Each exposure variable and postulated confounding variable was subjected to initial bivariable analyses using the survival analysis framework above. In addition to the exposure variables, covariates included in the survival analysis were the breed of cow (Breed), the season of birth of the cow (Season), and the agro-ecological zone (Ecoregion). No individual cow production parameters were included in these analyses as they were considered to be intervening variables (or effect modifiers) rather than potential confounders. Each of the explanatory variables was assessed by comparing the log likelihood of the intercept model to the log likelihood of the model containing the intercept and variable in question. The likelihood ratio statistic, with an asymptotic chi-square distribution and d.f. equal to the number of parameters associated with the effect, was initially used to assess the significance of associations with the outcome at P < 0.05. These were not initially subjected to generalized estimating equations (GEE) adjustments. Variables that were not significant at this level were not used in further analyses.

Multivariable survival analysis methods incorporating GEE were then utilized to adjust for the dependence of responses within herds (29,30). The cluster correlation was specified as independent for this analysis, providing conservative estimates of the standard error for inferring associations. Failure to account for clustered responses within herds may lead to mistaken inferences of association in cases where insufficient information exists in the data to permit such decisions. Full main-effects models, for each of the exposure estimate variables and all potential confounders (deemed to be significantly associated with the outcome at the initial bivariable assessment), were developed within the GEE framework for time to culling data, mortality data, or both. Backwards elimination was used until only variables significant at P < 0.05, utilizing the T score statistic (for GEE), remained (31). All possible first-order interactions between any single exposure variable deemed significant in main- effects models and the remaining confounding variables were likewise assessed for significance at P < 0.05.

Table I.

Table II.

Figure 1. Frequency distribution of cows culled in lactation lifespan (months from first calving date until culling date) in Alberta dairy herds (n = 70 052 cows culled). The clustered vertical bars represent stratification at the median herd PD exposure (more ...)

Figure 2. Frequency distribution of cow fatality in lactation lifespan (months from first calving date until death) in Alberta dairy herds (n = 8743 cows died). The clustered vertical bars represent stratification at the median herd PD exposure level (more ...)

Figure 3. Frequency distribution of annual culling risk in Alberta dairy herds (n = 7810 annual herd reports). The clustered vertical bars represent stratification at the median herd PD exposure level (1.0 × 10−1 μg/m3 of SO2) (more ...)

Figure 4. Frequency distribution of annual mortality risk in Alberta dairy herds (n = 7810 annual herd reports). The clustered vertical bars represent stratification at the median herd PD exposure level (1.0 × 10−1 μg/m3 of SO (more ...)

On the 1382 dairy farms studied, the annual median DD classification ranged from a low of 1.0 units, to a high of 2.0 units (annual MIN = 0; annual MAX = 59). The monthly farm-site classification ranged from a low of 0.0 units to a high of 66.0. There was a similar, though less marked, seasonality for the monthly and annual mean DD-estimated farm-site exposures. This was totally attributable to processing volume variations. The distribution of exposures was highly right-skewed with a peak just slightly above zero (data not shown) (27).

There was a total of 150 210 cows with both a birth date and a first calving date between January 1, 1985 and December 31, 1994. Of these, 70 052 were culled before December 31, 1994. The remaining 80 158 were considered censored; that is, they either: a) died; b) remained in the herds as of December 31, 1994; c) left the herds for dairy or export purposes; or d) the herd withdrew from DHI services. Only 8743 cows died on the farm. The remaining 141 467 cows were treated as censored observations for the mortality analysis.

The average length of productive (lactation) lifespan (first calving date until culling), death, or censoring (end of study or sold for dairying purposes) for the 150 210 dairy cows was 22.8 mo (s = 18.58; MIN = 0; MEDIAN = 18; MAX = 98), ending at an average age of 50.1 mo (s = 18.98; MIN = 22; MEDIAN = 46; MAX = 120). Of these, 70 052 cows were culled at an average of 21.3 mo (s = 16.49; MIN = 0; MEDIAN = 18; MAX = 94) following the delivery of their first calf (average age = 48.7 mo; s = 16.90; MIN = 22; MEDIAN = 45; MAX = 119). Of those deaths on the farms (n = 8743), the average lactation lifespan was 23.4 mo (s = 17.86; MIN = 0; MEDIAN = 19; MAX = 91) and these cows died at an average age of 50.9 mo. Animals right-censored (n = 71 415) from the culling and mortality study appeared to out-live their culled or mortal companions, leaving the herd at an average age of 51.3 mo (s = 20.83) after contributing a lactation lifespan of 25.1 mo (s = 20.40) to the dairy herds. The frequency distributions for the lactation lifespan of cows that were culled or died during the study period are shown in Figures 1 and 2, respectively.

The culling and mortality risk data set was derived from 91 305 monthly herd DHI reports combined into annual summaries of cows that exited the herd (for non-dairy sale purposes), and the average number of cows at risk for culling or mortality during a calendar year. Following the exclusion of herds that did not contribute a sufficient number of monthly records, the final data set included 7810 annual herd estimates of the culling risk. The annual herd culling risks, from 1986 to 1994, averaged 0.28 (s = 0.135; MIN = 0; MEDIAN = 0.27; MAX = 1; 99th percentile = 0.69) (Figure 3). The very highest culling risks (for example, greater than 60%) suggest the erroneous use of DHI disposal codes (for example, coded as culling as opposed to herd-dispersal sales) in a few cases. The mean annual herd mortality risk (adult cows) was 0.04 (s = 0.032; MIN = 0; MEDIAN = 0.03; MAX = 0.32; 99th percentile = 0.08) (Figure 4).

Figure 5. Survival curve of time to culling (months from first calving date until culling date) in Alberta dairy herds (n = 70 052 cows culled). The clustered vertical bars represent stratification at the median herd PD exposure level (0.98 × (more ...)

Figure 6. Survival curve of time to death (months from first calving date until death) in Alberta dairy herds (n = 8743 cows died). The clustered vertical bars represent stratification at the median herd PD exposure level (0.91 × 10−1 (more ...)

With respect to survival (lactation lifespan) until culling, the PD exposure estimate variable was found be marginally significant (β = -0.0018, P = 0.035) under initial (non-GEE) bivariate tests of association. The DD exposure estimate variable, likewise had a parameter estimate in a similar direction, but of lower magnitude (β = 0.0010; P = 0.128), and was not significant. Note that in this analysis a lower hazard ratio (= e^β, then compared to a referent level of exposure) indicates longer survival (or time to culling). The multivariable final model for time to culling did not contain either of the exposure estimate variables (such as, PD or DD) but included both the agro-ecological region (cattle from grassland regions stayed in the herd longer) and season of birth (cattle born during winter months had a shorter lactation lifespan).

For the time to death analysis, initial bivariate associations found to be significant in naïve (non-GEE) statistical tests included: a) the PD exposure estimate variable (β = −0.012; P < 0.0001); and b) the DD exposure estimate variable (β = −0.0079; P < 0.0001). In the multivariable modelling of the time to death, both PD and DD remained significant terms for the time to mortality models. The final parameters for PD (β = −0.0119; P = 0.029) and DD (β = 0.007; P < 0.0001) suggest that cows in herds exposed to higher levels of sour gas emissions died at an older age than cows that were not exposed or received lower exposure. Final models for the time to death analysis did not include any other significant (P > 0.05) covariates.

Culling and mortality risk; herd level — The unconditional (bivariable) associations for each of the exposure variables with both the annual herd culling and mortality risks, were significantly (P < 0.0001) associated with a decreased annual risk. Note that these likelyhood ratio tests were based on naïve significance tests (unadjusted by GEE for within-herd dependence). Neither DD (P = 0.06) or PD (P = 0.26) remained significantly associated with the culling risk in the final GEE model. The agro-ecological zone was the only term in the final model. The final annual multivariable mortality risk model did not include either of the DD or PD exposure estimate variables.

In the present study, 28.1% of the animals across all study herds were culled per year and cattle survived to an average age of 50 mo (maximum = 120 mo). These outcomes are comparable to data reported elsewhere. Sol et al (32) reported an annual culling rate in the Netherlands of 25 to 30 %. Dohoo and Martin (33) reported an average culling rate (per annum) of 20.1% with an average survival to 67 mo of age. Sargeant (34) reported an annual removal risk, similar to the present study, of 28%.

We are aware of no similar scientific investigations on retainment and survival of cows in dairy herds exposed to sour gas processing plant emissions. However, there are several well-documented anecdotal and clinical reports relating to culling of dairy herds in the literature. A practising veterinarian in west-central Alberta has described situations of sub-optimal herd health in the practice area which the veterinarian attributed to several sour gas processing plants in the region (35). In particular, the veterinarian noted dairy cattle problems including unthriftiness, susceptibility to infectious disease, reproductive problems, and downer cows. As a response to these and other factors, culling rates were perceived to be elevated on these farms. In 1993, the Energy Resources Conservation Board (ERCB – the predecessor of the AEUB) ruled against an oil company's application for sour-oil expansion adjacent to a dairy farm in east-central Alberta (36). Specific concerns cited by the complainants included: a) an increased abortion rate; b) slow rate of increase in milk productivity over several years (unrealized genetic potential); and c) subsequent increased culling rate resulting in an inability to meet quota requirements (36). This latter concern illustrates the additional market forces (for example, supply management) that may direct dairy producer decisions pertaining to cow retainment.

In terms of longevity of individual animals within the 1382 study herds (as modelled in the survival analyses), there was no evidence to suggest that sour gas processing plant emissions impacted the time from first-calving to culling. However, both the PD and DD variables indicated a decreased hazard for mortality on the study farms. This suggests that if or when, cows died on high sour gas-exposed farms, they did so at later ages than did cows on low sour gas-exposed farms. It is difficult to explain this finding independent of the knowledge pertaining to reasons for animal retainment in the herds. As mentioned earlier, given such a choice most dairy producers would choose to send a cow to slaughter rather than have that cow die on the farm. In this study, cows that died on farms did so at later ages than cows that were culled. This may suggest a later retention of unhealthy or poor-doing cattle on farms with higher levels of sour gas exposure. The reasons why this may be are unclear, and any such hypotheses would be untestable within the context of the historical data available for the present study. It seems quite plausible that the outcome suggesting a decreased hazard for mortality would likely not be independent of the herd culling risk noted below; even though overall survival (wastage) and time to culling did not appear to be associated with exposure.

This study provided no evidence to suggest that herds in regions of Alberta with higher levels of exposure to air emissions from sour gas processing plants exhibited higher risks for culling or mortality. No exposure variable was associated with herd-level mortality risk. The additional variable found to be associated with culling risk (but not mortality risk) was the agro-ecological region.

The survival to mortality analysis suggests that there may be a very small association manifested between levels of exposure to sour gas processing emissions and retention of animals in the herd. One interpretation could be that herds in high exposure regions exhibit improved retention of animals, in that deaths occur at a later age. However, any such judgement is heavily value-laden and presumes that herd operators are free to make culling and retention decisions as a series of voluntary (as opposed to involuntary) choices. Indeed, target levels (37) for culling reflect the principle idea that rates of culling can also be too low (resulting in slower genetic progress and the retention of low-producing animals), as well as too high. However, without information pertaining to the production or economic (for example, quota) needs of the herd, it is difficult to attribute value-laden terms; such as, improved or decreased, to culling data. At this juncture, we can simply state that there is no evidence of a harmful association between exposure to sour gas processing emissions and the hazard for mortality, in dairy herds in Alberta.