Arny L Blanchard
email: arnyb@ims.uaf.edu.

Howard M Feder

Max K Hoberg

P O Box

Temporal Variability of Benthic Communities in an Alaskan Glacial Fjord, 1971-2007

Keywords: Benthos, Earthquake, Port Valdez, Prince William Sound, Marine Oil Terminal, Multiple Stressors, Disturbance, Long-Term Monitoring, Climate Change

Temporal trends of deep subtidal macrofauna in Port Valdez, Alaska, were assessed with respect to multiple environmental stressors. Effects from a magnitude 9.2 earthquake in Prince William Sound, Alaska, 1964, were reflected in recolonization of the basin of the fjord, increased abundance and number of taxa over time, and moderately increased variability in abundance through 1990, stabilizing 26 years after the earthquake. Long-term climatic variability and local physical processes were important sources of spatial and temporal variability. Correlative evidence suggests that indirect effects of juvenile salmon from a shoreline salmon hatchery and deposition of adult salmon carcasses moderately enhanced deep basin benthic communities. Effects on the deep benthos from a marine oil terminal were negligible. Overall, faunal trends deviated from the stability expected for benthic communities in other fjords. Physical characteristics of the fjord were important in mediating the effects of stressors and in delaying the readjustment process.

Introduction

The coastal fringe is a multiple-use region where frequent environmental perturbations from natural and anthropogenic stressors can have large, long-lasting effects (Lindegarth and Hoskin, 2001). Human activities can result in temporally persistent changes and stress in coastal environments through disturbance as varied as nonindigenous species introductions, alteration of drainage patterns by industrial and urban development, dredging, and waste disposal [START_REF] Nichols | Abundance fluctuations among benthic invertebrates in two pacific estuaries[END_REF][START_REF] Nichols | Persistence of an introduced mudflat community in South San Francisco Bay, California[END_REF][START_REF] Pearson | Long-term changes in benthic populations in some west European coastal areas[END_REF][START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF][START_REF] Perus | Long-term changes in macrozoobenthos in the Åland archipelago, northern Baltic Sea[END_REF]. Changes in nearshore environments often reflect increasing stress from urban growth and development with the direction of change towards greater habitat loss and long-term disturbance to marine communities [START_REF] Schiff | Southern California[END_REF]. In view of large-scale changes in the marine environment resulting from climate change, human activities, and landscape-scale changes, understanding how communities vary in response to stressors becomes critical.

Assessment of disturbance in marine communities is complicated by multiple stressors acting and interacting over varying spatial and temporal scales (Blanchard and Feder, 2003;[START_REF] Hewitt | Assessing and monitoring ecological community health in marine systems[END_REF]. Preferably, undisturbed community characteristics and information on natural sources and directions of change should be available prior to a disturbance [START_REF] Nichols | Abundance fluctuations among benthic invertebrates in two pacific estuaries[END_REF]. This is difficult since marine communities are continually adjusting to long-term change from varied sources of stress such as nonindigenous species introductions, global climatic change, global spread of contaminants, urbanization, and increasing resource use [START_REF] Steneck | Human alterations of marine communities: students beware![END_REF][START_REF] Dojiri | Changing anthropogenic influence on the Santa Monica Bay watershed[END_REF][START_REF] Hewitt | Assessing and monitoring ecological community health in marine systems[END_REF]. Small, localized disturbances from human activities can often be readily identified [START_REF] Blanchard | Long-term investigation of Benthic Fauna and the Influence of Treated Ballast Water Disposal in Port Valdez, Alaska[END_REF](Blanchard et al., , 2003;;Blanchard and Feder, 2003) whereas determination of community readjustment from perturbations large enough to impact whole marine ecosystems [START_REF] Josefson | Long-term soft-bottom faunal changes in three shallow fjords, west Sweden[END_REF][START_REF] Austen | Comparison of long-term trends in benthic and pelagic communities of the North Sea[END_REF][START_REF] Lay | The great Sumatra-Andaman earthquake of 26 December 2004[END_REF] is more difficult. Although the effects of stressors cannot always be predicted, long-term studies of disturbed systems provide a means to gain insights into the responses to stress and temporal variability of marine communities and long-term trends [START_REF] Hawkins | Recovery of polluted ecosystems: the case for long-term studies[END_REF][START_REF] Whomersley | The use of time-series data in the assessment of macrobenthic community change after the cessation of sewage-sludge disposal in Liverpool Bay (UK)[END_REF].

Physical and ecological characteristics of glacial fjords result in strong environmental sensitivity to natural and anthropogenic stressors. Deep basins in glacial fjords are sheltered environments with restricted exchange with oceanic water and are highly stratified with only seasonal mixing of deep water [START_REF] Muench | Physical oceanography[END_REF][START_REF] Jewett | Effects of the Exxon Valdez oil spill on benthic invertebrates in an oxygen-deficient embayment in Prince William Sound, Alaska[END_REF][START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. The restricted circulation and stratification in silled-fjords limits larval dispersion, nutrient exchange, and deep water renewal as well [START_REF] Jewett | Effects of the Exxon Valdez oil spill on benthic invertebrates in an oxygen-deficient embayment in Prince William Sound, Alaska[END_REF][START_REF] Josefson | Long-term soft-bottom faunal changes in three shallow fjords, west Sweden[END_REF][START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF][START_REF] Molinet | Effects of sill processes on the distribution of epineustonic competent larvae in a stratified system of Southern Chile[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. The fine glacial sediments deposited in glacial fjords form unstable slopes, bury juvenile and adult macrofauna, and dilute particulate organic carbon (POC) [START_REF] Görlich | Suspension settling effect on macrobenthos biomass distribution in the Hornsund fjord[END_REF][START_REF] Feder | The subtidal benthos[END_REF][START_REF] Holte | Some macrofaunal effects of local pollution and glacier-induced sedimentation with indicative chemical analyses, in the sediments of two arctic fjords[END_REF][START_REF] Wlodarska-Kowalczuk | Soft-bottom macrobenthic faunal associations and factors affecting species distributions in an arctic glacial fjord (Kongsfjord, Spitsbergen)[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. Fluxes of organic carbon to deep basins of fjords may be on the order of 1-3 g C m -2 year -1 , as in the Balsfjord in Norway and Port Valdez, Alaska, resulting in carbon-poor benthic systems with low infaunal abundance and biomass [START_REF] Feder | The subtidal benthos[END_REF][START_REF] Wlodarska-Kowalczuk | Spitsbergen glacial bays macrobenthos -a comparative study[END_REF][START_REF] Oug | Soft-bottom macrofauna in the high-latitude ecosystem of Balsford, northern Norway: Species composition, community structure, and temporal variability[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. Spatial gradients in fauna and community structure may show little change on decadal scales but the sensitivity of fauna to long-term climatic change is unknown [START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. As a result of their physical and ecological characteristics, glacial fjords can be very sensitive to anthropogenic disturbance [START_REF] Josefson | Long-term soft-bottom faunal changes in three shallow fjords, west Sweden[END_REF][START_REF] Holte | Some macrofaunal effects of local pollution and glacier-induced sedimentation with indicative chemical analyses, in the sediments of two arctic fjords[END_REF][START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF]Blanchard and Feder, 2003;[START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. The restricted circulation and thus, isolation of deep basins of glacial fjords contributes to their importance as sites for monitoring long-term change [START_REF] Josefson | Long-term soft-bottom faunal changes in three shallow fjords, west Sweden[END_REF][START_REF] Holte | Some macrofaunal effects of local pollution and glacier-induced sedimentation with indicative chemical analyses, in the sediments of two arctic fjords[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF].

A 35-year investigation of the deep-subtidal macrofauna in Port Valdez, Alaska, provides the opportunity to examine the temporal variability of benthic organisms in a glacial fjord during a period with large environmental changes. Environmental perturbations in the study area include a major earthquake in 1964 and increasing anthropogenic stresses [START_REF] Feder | The subtidal benthos[END_REF][START_REF] Mcroy | Natural and anthropogenic disturbances at the ecosystem level[END_REF][START_REF] Wiegers | A regional multiple-stressor rank-based ecological risk assessment for the fjord of Port Valdez, Alaska[END_REF]. Regional climatic variations have been identified as key factors for ecological processes in the Gulf of Alaska and Prince William Sound (PWS) ecosystems [START_REF] Mundy | Introduction[END_REF] contributing to temporal variability of infaunal communities in Port Valdez and PWS. Additionally, limited responses of fauna to short-term increases in sediment hydrocarbons have been observed in shallow sediments adjacent to the marine oil terminal in Port Valdez [START_REF] Blanchard | Long-term investigation of Benthic Fauna and the Influence of Treated Ballast Water Disposal in Port Valdez, Alaska[END_REF]. Objectives of this study were to examine the macrofaunal community in the deep basin of the fjord from 1971 to 2007 to determine if temporal variability represented biologically significant deviations from the stability expected of fjord-basin communities and to gain insights into the influences of climatic variability, natural gradients, and stressors on variability in community structure.

Study area

Sampling occurred in Port Valdez, a glacial outwash fjord in the northeastern corner of PWS (Hood et al., 1973;[START_REF] Colonell | Port Valdez, Alaska: Environmental Studies 1976-1979[END_REF]Shaw and Hameedi, 1988) (Figure 1). Port Valdez has two sills (at 120 and 200 m depth) and a relatively flat bottom in the deep basin varying between 230-250 m depth. In summer, positive estuarine circulation occurs within surface waters to about 15 m [START_REF] Muench | Physical oceanography[END_REF][START_REF] Sharma | Geological oceanography[END_REF]Colonell, 1980a). The counterclockwise movement of water is slow (< 5 cm sec -1 ) with a residence time estimated as 40 days. Deep-water layers below sills (> 120 m) have little exchange with upper water layers in summer but are well mixed in winter (Colonell, 1980a). The sediment is dominated by silt and clay fractions carried to the fjord in glacial melt water from the surrounding mountains. Seasonal variations in sediment flux are large ranging from less than 5 mg cm -2 day -1 in winter up to 48 mg cm -2 day -1 in summer with highest sediment fluxes observed in the eastern end of the fjord [START_REF] Feder | Subtidal benthos[END_REF]Fig. 2). Annual sedimentation rate at the eastern end of the fjord was estimated to be 13.5 cm yr -1 but < 1 cm yr -1 at the western end [START_REF] Naidu | Sedimentation processes[END_REF]. Macrofaunal communities in the fjord's basin demonstrated relatively low numbers overall and moderate gradients in community structure associated with high sediment loads from glacial rivers [START_REF] Feder | Subtidal benthos[END_REF][START_REF] Feder | The subtidal benthos[END_REF]).

An east to west division of faunal communities was related to sedimentation of fine glacial sediments corresponding to the general position of the seasonal sediment plume from glacial runoff in the fjord (Fig. 1) [START_REF] Feder | Subtidal benthos[END_REF][START_REF] Feder | The subtidal benthos[END_REF].

A magnitude 9.2 earthquake in March, 1964, caused catastrophic disturbances to all marine communities in PWS [START_REF] Haven | Effects of land-level changes on intertidal invertebrates, with discussion of postearthquake ecological succession[END_REF][START_REF] Hubbard | Distribution and abundance of intertidal invertebrates at Olsen Bay in Prince William Sound, Alaska, one year after the 1964 earthquake[END_REF][START_REF] Stanley | Effects of the Alaska earthquake of March 27, 1964 on shore processes and beach morphology[END_REF]. The intertidal and subtidal benthic communities of the study area were impacted by scouring from the tsunamis and catastrophic sediment slumping [START_REF] Coulter | Effects at Valdez[END_REF][START_REF] Mcroy | Natural and anthropogenic disturbances at the ecosystem level[END_REF][START_REF] Sharma | Geological oceanography[END_REF][START_REF] Naidu | Sedimentation processes[END_REF]. Displaced sediments were deposited in the basin of Port Valdez, directly impacting benthic organisms through burial.

Anthropogenic stresses on the marine environment in Port Valdez increased greatly over the 35-year investigation [START_REF] Mcroy | Natural and anthropogenic disturbances at the ecosystem level[END_REF]. As the northernmost, ice-free port available for large vessel traffic, Port Valdez was selected to be the site of the marine oil terminal serving as the terminus for the Trans-Alaska pipeline which became operational in 1977. Discharge of treated ballast water from the oil terminal is a continuing stress to some shallow-water biological communities adjacent to the oil terminal [START_REF] Blanchard | Long-term investigation of Benthic Fauna and the Influence of Treated Ballast Water Disposal in Port Valdez, Alaska[END_REF](Blanchard et al., , 2003)). Oily ballast water is unloaded from tankers at the oil terminal, treated onshore to remove hydrocarbons, and discharged into the fjord through a diffuser pipe at 60m to 80 m depth. Discharges are decreasing due to the declining volume of oil transported [START_REF] Blanchard | Long-term investigation of Benthic Fauna and the Influence of Treated Ballast Water Disposal in Port Valdez, Alaska[END_REF][START_REF] Richardson | Economics of human uses and activities in the Northern Gulf of Alaska[END_REF]. A salmon hatchery located on the south side of the fjord, downstream from the oil terminal, began commercial operations in 1982 with first returns in 1983. Over 200 million salmon fry are now raised and released every year with returns of more than 10 million adult fish (predominantly pink salmon Oncorhynchus gorbuscha) [START_REF] White | Alaska salmon enhancement program 2007 annual report[END_REF].

There are too many returning adult salmon to fully harvest, so many salmon die in the fjord and their carcasses are deposited in the intertidal or settle on subtidal sediments.

Growth of the City of Valdez and associated urban development altered nearshore habitats largely through urban development and industrialization such as construction of roads, docks, and a gas refinery [START_REF] Wiegers | A regional multiple-stressor rank-based ecological risk assessment for the fjord of Port Valdez, Alaska[END_REF]Blanchard and Feder 2003).

Ecosystem stress in the fjord increased through various means including urban discharges and uncontrolled urban runoff, habitat alterations, increased marine vessel traffic, and continuing nearshore development [START_REF] Mcroy | Natural and anthropogenic disturbances at the ecosystem level[END_REF][START_REF] Wiegers | A regional multiple-stressor rank-based ecological risk assessment for the fjord of Port Valdez, Alaska[END_REF]Blanchard and Feder, 2003).

Methods

Subtidal fauna were first sampled in Port Valdez in 1971. Six stations sampled over the study period comprise a transect down the center of the fjord with three replicate samples collected at sites on the transect in September 1971 and generally five from late summer to fall in 1976-1977, 1980 to 1982, 1985, 1987 and annually since 1989 to 2007 (not all sites were sampled every year) (Fig. 1). This transect encompasses a gradient in sedimentation rates and environmental processes from the head of the fjord (where glacial rivers enter the fjord and suspended sediments are highest) to the fjord mouth.

Sediment for faunal samples was collected with a van Veen grab (0.1 m 2 ), washed over a 1.0 mm screen, and fauna preserved in 10% buffered formalin. Macrofauna (infaunal invertebrates retained on 1.0 mm mesh sieve) were identified to the lowest taxonomic level practical, mostly to genera or species, counted, and weighed. Faunal identifications were aggregated to the family level or higher for all analyses because of taxonomic uncertainties and changes in identifications from the early years of the study (1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980) to the present. While analysis of more detailed taxonomic information is preferable, the loss of information due to aggregation to family is not too great and ecologically relevant trends will still be apparent [START_REF] Feder | The deep benthos of Prince William Sound, Alaska, sixteen months after the Exxon Valdez oil spill[END_REF][START_REF] Wlodarska-Kowalczuk | Surrogacy in natural patterns of benthic distribution and diversity: selected taxa versus lower taxonomic resolution[END_REF]. Species identifications are noted in the text where possible. In 2002, field sampling also included occupation of many sites from the fjord-wide grid established in 1971 (Fig. 1).

Faunal measures calculated include average abundance (individuals m -2 ), wet weight biomass (g m -2 ), and total number of taxa (includes family or unique higher taxonomic categories). Average abundance and biomass for each year represent the averages for stations on the deep-basin transect for each year. The tabulated statistic total number of taxa was calculated as the number of taxa found within a year (over all stations combined) whereas geostatistical modeling relied on the number of taxa per station per year.

Nonmetric multidimensional scaling (MDS) was used to elucidate communitylevel changes in family abundance estimates over time. To assess the changes in the faunal relationships, surveys from 1971, 1976, 1981, 1985, 1990, 1995, 2000, and 2005 were analyzed. Bray-Curtis similarity coefficients [START_REF] Bray | An ordination of the upland forest communities of Southern Wisconsin[END_REF] were calculated between the individual stations for each year using ln(X+1)-transformed mean family abundance. MDS was applied to the resulting similarity matrices using the PRIMER software package to assess similarities of stations between years [START_REF] Clarke | User Manual/Tutorial[END_REF]. The ANOSIM routine tested for differences between years and the eastern and western stations (as determined by [START_REF] Feder | Subtidal benthos[END_REF]. The SIMPER routine of PRIMER was used to determine contributions of individual taxon categories to year groups.

Inclusion of the fjord-wide sites sampled in 1971 and 2002 allowed for direct comparisons to be made in the abundance and distribution of selected families. Bubble plots and bar charts of the abundance with 95% confidence intervals of the selected fauna were generated to assess changes in the distributions of fauna. Two-way ANOVA was used to assess differences between years and eastern and western fjord stations.

Spatio-temporal modeling was applied to assess changes over space and time in macrofauna abundance, biomass and the number of taxa [START_REF] Stein | Space-time statistics for environmental and agricultural related phenomena[END_REF]. The statistical program R (www.r-project.org) and the spatial tools in the library GeoR [START_REF] Ribeiro | geoR: A package for geostatistical analysis[END_REF]; http://cran.r-project.org/) were used to perform geostatistical analysis of the biotic variables. The data were average abundance, biomass, and number of taxa for each station by year. Only those stations on the transect along the center of the fjord were used for spatio-temporal modeling (Fig. 1). The suite of stations sampled over time represents a one-dimensional measure of distance so distance and time were used to model changes in faunal parameters. The correlation structures of the variables were modeled with the spherical model [START_REF] Cressie | Statistics for spatial data[END_REF]. Abundance values were ln(X+1)-transformed and biomass ln(X)-transformed to meet the assumption of normal errors and all data were detrended by modeling trends as a polynomial functions of distance and year. Contour plots of back-transformed predicted values from the spatial models are presented.

Spearman's nonparametric correlation coefficient was used to determine associations between biotic variables and measures of anthropogenic stress, weather variables, and a climatic index. Biological indices used include average abundance and biomass, and total number of taxa per year. Available measures of anthropogenic stress were salmon fry releases , adult salmon returns , and average total aromatics (PAH: 1976(PAH: , 1977(PAH: (different deep stations for 1976(PAH: and 1977 but also representing a transect down the fjord), 1980-1982, 1985, 1989-2007). The data on salmon fry releases and adult returns were from annual reports to the State of Alaska (Alaska salmon enhancement program annual reports available through Alaska Dept.

Fish and Game, Juneau, AK; e.g., [START_REF] White | Alaska salmon enhancement program 2007 annual report[END_REF] and data for PAH (the sum of 18 polycyclic aromatic hydrocarbons) were from the long-term monitoring study in Port Valdez [START_REF] Blanchard | Long-term investigation of Benthic Fauna and the Influence of Treated Ballast Water Disposal in Port Valdez, Alaska[END_REF] (Table 1). Measures of weather and regional climatic variability include average annual temperature, total precipitation, and annual snowfall from the Valdez weather station (http://climate.gi.alaska.edu/) and the yearly average of the Pacific Decadal Oscillation index (PDO) (www.aoss.org). The PDO reflects anomalies in sea surface temperature (with positive anomalies representing warmer water) and is a measure of climatic variability strongly associated with ecological processes and faunal changes in the northern Pacific ocean [START_REF] Mundy | Introduction[END_REF][START_REF] Pinchuk | Climate-related variability in abundance and reproduction of euphausiids in the northern Gulf of Alaska in 1998-2003[END_REF].

Results

A total of 11 phyla encompassing 77 families were identified. Polychaetes and bivalves comprised 64% and 22% of the total numbers of individuals, respectively, for all years combined. The highest abundance occurred in 1987 with an average of 589 ind. m - 2 , the highest biomass value was in 2005 with 33.79 g m -2 , and the highest total number of taxa (family or greater) across the transect occurred in 1997 with 52 taxon categories followed by 1976 with 51 categories recorded (Table 1). The lowest abundance was in 1980 with 156 ind. m -2 , lowest biomass was in 1985 with 7.02 g m -2 , and the lowest number of taxa was in 1971 with 29 taxon categories identified.

The MDS ordination suggested a high level of similarity among years and sampling locations (Fig. 3). Nevertheless, separations by year were apparent as sampling locations for 1971 and 1976 grouped by year while stations for 1981 and 1985 were mixed together forming a separate group, as was also true for 1990 to 2005. ANOSIM indicated that the year groupings were significant as the overall R = 0.563 (p = 0.001). 1981 and 1985 and 7% in 1990 to 2005 (1990, 1995, 2000, and 2005) (Fig. 3). The decreasing values with increasing distance from the west end of the fjord (Fig. 7). The whisker plot of biomass grouped by year showed a negligible increase in biomass over time and negligible differences in variance from 1971 through 1990 (SD = 7.4 ind. m -2 ) compared to 1991 to 2007 (SD = 5.3 ind. m -2 ). For the number of taxa, the kriging and whisker plots demonstrated increasing number of taxa over time, with a sharp increase in the number of taxa from 1985 through 2007, and decreasing values with distance west to east (Fig. 8). The difference in variability in the number of taxa was negligible from 1971 to 1990 (SD = 8.0 ind. m -2 ), as compared to later years (SD = 6.6 ind. m -2 ).

Correlation analysis indicated moderate relationships between faunal measures and indices of environmental change. Moderate correlations were observed between biomass and salmon fry releases and adult returns (ρ = 0.350 and 0.399, respectively) and the average PDO and abundance and biomass (ρ = 0.377 and -0.311, respectively; Table 2 and Fig. 9). Negligible to small correlations were observed between releases and returns and abundance (ρ = 0.229 and 0.119, respectively) and the number of taxa (ρ = -0.269 and -0.244, respectively). Correlations between PAH and the biological variables were also negligible to small (ρ = 0.049 to 0.292). Correlations between local climatic measures were small for average precipitation and biomass (p = -0.280), average temperature and abundance (p = 0.250), and total snowfall and biomass (ρ = 0.257).

Discussion

Influences of Environmental Gradients in Fjords and Climatic Variability

The macrofauna assemblage in the deep basin of Port Valdez is characteristic of fauna found in glacial fjords worldwide. Deposit-feeding polychaetes are dominant in the study area, as found in glacial fjords elsewhere, including the surface-deposit feeding owenid polychaete Galathowenia oculata, spionid polychaetes (Polydora spp. and Prionospio spp.), and subsurface-deposit feeding polychaetes Heteromastus filiformis (family Capitellidae) and family Maldanidae [START_REF] Gulliksen | The soft bottom fauna in van Mijenfjord and Raudfjord, Svalbard[END_REF][START_REF] Holte | The macrofauna and main functional interactions in the sill basin sediments of the pristine Holandsfjord, northern Norway, with autecological reviews for some keyspecies[END_REF][START_REF] Wlodarska-Kowalczuk | Spitsbergen glacial bays macrobenthos -a comparative study[END_REF]Feder andBlanchard, 1998, Burd et al., 2000;[START_REF] Hoberg | The macrobenthos of sites within Prince William Sound, Alaska, prior to the Exxon Valdez oil spill[END_REF][START_REF] Ríos | Soft-bottom macrobenthic faunal associations in the southern Chilean glacial fjord complex[END_REF]. Deposit-feeding bivalves found in Port Valdez are also common in northern-hemisphere glacial fjords including the families Nuculanidae (Megayoldia/Yoldia spp. and Nuculana sp.) and Thyasiridae (Axinospida spp. and Thyasira spp.) [START_REF] Holte | The macrofauna and main functional interactions in the sill basin sediments of the pristine Holandsfjord, northern Norway, with autecological reviews for some keyspecies[END_REF]. The dominance of deposit-feeding animals is the result of the inability of most surface-dwelling fauna to adjust to the reduced availability of carbon fluxing to the bottom and high sedimentation rates [START_REF] Görlich | Suspension settling effect on macrobenthos biomass distribution in the Hornsund fjord[END_REF][START_REF] Holte | Some macrofaunal effects of local pollution and glacier-induced sedimentation with indicative chemical analyses, in the sediments of two arctic fjords[END_REF][START_REF] Wlodarska-Kowalczuk | Soft-bottom macrobenthic faunal associations and factors affecting species distributions in an arctic glacial fjord (Kongsfjord, Spitsbergen)[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. The combined influence of seasonal stratification, seasonality of production, early consumption of water-column primary production by zooplankton, and dilution of food by glacial sediments contribute to the paucity of carbon in benthic systems of glacial fjords [START_REF] Sargent | Partial characterization of organic material in surface sediments from Balsfjorden, northern Norway, in relation to its origin and nutritional value for sediment-ingesting animals[END_REF][START_REF] Cooney | Water column production[END_REF][START_REF] Oug | Soft-bottom macrofauna in the high-latitude ecosystem of Balsford, northern Norway: Species composition, community structure, and temporal variability[END_REF]. Lack of carbon reaching the benthos in Port Valdez was evident in the low estimate for annual macrofaunal production (0.3 to 1.7 g m -2 yr -1 ) [START_REF] Feder | The subtidal benthos[END_REF] and low abundance and biomass (Table 1).

Macrofauna in glacial fjords are stressed by high sedimentation rates resulting in strong gradients in community structure [START_REF] Holte | Common macrofaunal dominant species in the sediments of some north Norwegian and Svalbard glacial fjords[END_REF][START_REF] Wlodarska-Kowalczuk | Benthic response to chronic natural disturbance by glacial sedimentation in an Arctic fjord[END_REF]. A common trend of increasing abundance and diversity with lower biomass towards the head of a fjord reflects the increased proportions of small disturbance-tolerant organisms and the greater stress found near sources of glacial sediments [START_REF] Holte | Common macrofaunal dominant species in the sediments of some north Norwegian and Svalbard glacial fjords[END_REF][START_REF] Wlodarska-Kowalczuk | Soft-bottom macrobenthic faunal associations and factors affecting species distributions in an arctic glacial fjord (Kongsfjord, Spitsbergen)[END_REF][START_REF] Wlodarska-Kowalczuk | Benthic response to chronic natural disturbance by glacial sedimentation in an Arctic fjord[END_REF]. The general spatial trend associated with sedimentation rates and fine sediments was observed in prior studies of Port Valdez [START_REF] Feder | Preliminary benthos survey[END_REF][START_REF] Feder | The subtidal benthos[END_REF][START_REF] Naidu | Sedimentation processes[END_REF] where, as also indicated by spatio-temporal modeling in the present study, communities in the western end had higher number of taxa and biomass and lower abundance compared to the faunal assemblage of the eastern end (closest to the major sources for glacial sediments). The major faunal trend in van Mijenfjord, a silled fjord in Spitsbergen, reflected disturbance from fine sediments, and gradients in faunal communities were similar between [START_REF] Colonell | Physical oceanography[END_REF]and 2001[START_REF] Wlodarska-Kowalczuk | Benthic response to chronic natural disturbance by glacial sedimentation in an Arctic fjord[END_REF][START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. Like van Mijenfjord, little temporal change in community structure was reported for Balsfjord, Norway although minor increases in abundance were noted, possibly a result of changes in the composition of pelagic fauna and grazing rates [START_REF] Oug | Soft-bottom macrofauna in the high-latitude ecosystem of Balsford, northern Norway: Species composition, community structure, and temporal variability[END_REF]. In general, the physical characteristics and sediment gradients in glacial fjords result in the temporal and spatial persistence of gradients in benthic faunal community structure as well.

The Pacific decadal oscillation (PDO) is an index of climate variability associated with ecologically-important changes in physical processes of the North Pacific Ocean and nearshore coastal waters. The PDO contrasts sea surface temperatures between the Central and Northeast Pacific such that when the PDO is positive, sea surface temperatures are elevated in the coastal Northeast Pacific [START_REF] Royer | Interdecadal variability of Northeast Pacific coastal freshwater and its implications on biological productivity[END_REF]. The warmer waters of a positive PDO leads to increased rainfall and greater coastal freshwater flow along the Northeast Pacific which may have a number of effects on the ecosystem including strengthening of stratification in the water column in summer [START_REF] Neal | Linking the pacific decadal osscilation to seasonal stream discharge patterns in Southeast Alaska[END_REF]Mundy and Cooney, 1995;[START_REF] Weingartner | Freshwater variability and predictability in the Alaska Coastal Current[END_REF]. Increased stratification of the water column will initially result in greater primary production but may prevent mixing of nutrients upward from deep water later [START_REF] Royer | Interdecadal variability of Northeast Pacific coastal freshwater and its implications on biological productivity[END_REF]. Shifts in ecological processes associated with the PDO variations may influence benthic fauna in Prince William Sound (PWS) and Port Valdez through changes in the timing and volume of annual primary production, grazer availability, grazing rates, flux of nutrients to the benthos, and processes influencing larval survival (Mundy and Cooney, 1995;[START_REF] Beuchel | Long-term patterns of rocky bottom macrobenthic community structure in an Arctic fjord (Kongsfjorden, Svalbard) in relation to climate variability (1980-2003)[END_REF]. 2).

The association between abundance and the PDO index were strongest from 1987 to 2007 with less apparent association from 1971 to 1987 when the benthos was recovering from the 1964 earthquake (Fig. 9) (see discussion below). Most primary production of Alaskan glacial fjords reaches deep sediments through undigested matter in zooplankton fecal pellets with most production occurring in the spring before glacial melt waters reduce visibility [START_REF] Cooney | Water column production[END_REF][START_REF] Burrell | Patterns of carbon supply and distribution and oxygen renewal in two Alaskan fjords[END_REF]). This provides one potential pathway for climatic variability to affect benthic communities as greater rainfall associated with positive PDO values may result in a greater snowpack in surrounding mountains. The increased spring runoff may then enhance stratification resulting in greater spring primary production (Mundy and Cooney, 1995;[START_REF] Royer | Interdecadal variability of Northeast Pacific coastal freshwater and its implications on biological productivity[END_REF] and consumption of the increased production by zooplankton and fall of fecal pellets to the bottom would then enhance survival of new recruits in the benthos. Benthic community structure may also be directly influenced by interactions between water temperature and survival of pelagic larvae [START_REF] Barnett | Long-term changes in some benthic species in the Firth of Clyde, with particular reference to Telinna tenuis da Costa[END_REF][START_REF] Pearson | Long term changes in the benthic communities of Loch Linnhe and Loch Eil (Scotland)[END_REF].

The temporal stability of community structure in a protected, silled-fjord in Spitsbergen (with tidewater glaciers) suggested that benthic fauna were largely influenced by interactions between local hydrography, physical characteristics, and faunal life-history traits and may be resistant to decadal-scale environmental variability [START_REF] Renaud | Multidecadal stability of benthic community structure in a high-Arctic glacial fjord (van Mijenfjord, Spitsbergen)[END_REF]. Like the macrofauna of Spitsbergen (with tidewater glaciers), the overall community structure of macrofauna in the basin of Port Valdez has not shown great change since 1990 (see discussion below) since overall faunal composition is controlled by the ecological characteristics of the fjord. However, the present study does suggest that fauna in deep basins can be responsive to regional environmental trends through changes in other facets of community composition (abundance and biomass).

The increased volume of glacial melt waters and runoff in warmer years from the large watershed in the Chugach Mountains surrounding Port Valdez (a glacial outwash fjord) (which would be absent or reduced in fjords with tidewater glaciers) is likely the source for effects on benthic fauna at depth (>150 m). Thus, sensitivity to climatic change in nearshore communities may be dependent on landscape-scale factors of surrounding watersheds rather than just interactions of oceanographic and biological features within a fjord.

Effects of the Great Alaska Earthquake in 1964

The greatest perturbation in Port Valdez in recent history was a magnitude 9.2 earthquake in PWS March, 1964. The earthquake was so great as to have impacted every marine community in PWS by topographic changes, scouring by tsunamis, and sediment slumping and resuspension [START_REF] Coulter | Effects at Valdez[END_REF]NRC, 1971). The deep benthos in Port Valdez was disrupted by the mass sediment slumping triggered by the earthquake and deposition of sediments over the bottom [START_REF] Coulter | Effects at Valdez[END_REF][START_REF] Sharma | Geological oceanography[END_REF][START_REF] Feder | The subtidal benthos[END_REF][START_REF] Naidu | Sedimentation processes[END_REF]. Effects of the catastrophic disturbance in Port Valdez were reflected in differences of community structure among years in the multivariate analyses, the temporal changes in abundance and distributions of fauna, and spatio-temporal modeling (Fig. 3). Additionally, moderately increased variability in abundance (a symptom of disturbance) was apparent for Port Valdez from 1971 to 1989 as also occurred following cessation of sewage-sludge disposal in Liverpool Bay, UK (Fig. 6) [START_REF] Warwick | Increased variability as a symptom of stress in marine communities[END_REF][START_REF] Whomersley | The use of time-series data in the assessment of macrobenthic community change after the cessation of sewage-sludge disposal in Liverpool Bay (UK)[END_REF]. The much weaker association of faunal abundance to the PDO index from 1971 to 1989 is additional evidence for the presence of a strong, ecological process in the benthos related to faunal readjustment from the earthquake rather than climate effects (Fig. 9). Since there are no published studies of long-term trends in benthic fauna following a physical disturbance of a magnitude or scale comparable to the 1964 earthquake in PWS, it is difficult to discern to what extent physical characteristics of Port Valdez contributed to delayed recruitment. However, interactions between fauna and physical factors can be expected to influence temporal trends of faunal assemblages since readjustment is determined to varying degrees by hydrodynamics, circulation patterns or restrictions, and other physical features (e.g., [START_REF] Hall | Physical disturbance and marine benthic communities: life in unconsolidated sediments[END_REF][START_REF] Molinet | Effects of sill processes on the distribution of epineustonic competent larvae in a stratified system of Southern Chile[END_REF].

Faunal readjustment from large disturbances can be a lengthy process in glacial fjords. The benthic community in disturbed sediments in Alice Arm, a silled fjord in British Columbia disturbed by mining activities, was dominated by disturbance-tolerant species until disturbance abated [START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF]. Although considerable faunal readjustment occurred three years after cessation of tailing disposal, up to 18 years was required for biota to return to the ambient faunal community in Alice Arm. Faunal changes included introductions of numerous species into disturbed locations where they became dominant members of the macrofaunal assemblage. Initial colonizers responded quickly to cessation of disturbance but recruitment of stress-intolerant species was delayed [START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF]. The temporal trends in Port Valdez were similar in scope to readjustment patterns of fauna in Alice Arm suggesting that large-scale recolonization occurred, as would be expected following a very large earthquake. In the present study, the absence or low abundances of Maldanidae, Oweniidae, Paraonidae and Thyasiridae in 1971 and subsequent redistributions were striking since these organisms are common in muddy sediments in Prince William Sound (PWS) (Figs. 4 and5) [START_REF] Feder | The deep benthos of Prince William Sound, Alaska, sixteen months after the Exxon Valdez oil spill[END_REF][START_REF] Hoberg | The macrobenthos of sites within Prince William Sound, Alaska, prior to the Exxon Valdez oil spill[END_REF]. Additionally, spatio-temporal modeling indicated that the number of taxa and abundance increased over time, as would be expected following disturbance (Figs. 6 and8) [START_REF] Boesch | Response to Stress in Marine Benthic Communities[END_REF][START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF].

Disturbance to the faunal assemblage was also suggested by the relative dominance of Capitellidae, a disturbance-tolerant, opportunistic polychaete family, in 1971 and 1976 (Fig. 3). Evidence from faunal trends and spatio-temporal modeling indicates that the macrofaunal assemblage in the basin of Port Valdez reached a point of stability in 1990, 26 years after the 1964 earthquake.

Temporal persistence of effects from physical disturbance on benthic systems is dependent on the scale of disturbance, larval and adult recruit supplies, pre and postsettlement processes for recruits, and other factors [START_REF] Hall | Physical disturbance and marine benthic communities: life in unconsolidated sediments[END_REF]. Recovery can be rapid when new recruits and adults are readily available even when disturbance is large [START_REF] Boesch | Response to Stress in Marine Benthic Communities[END_REF][START_REF] Snelgrove | Animal-sediment relationships revisited: cause versus effect[END_REF][START_REF] Reise | Sediment mediated species interactions in coastal waters[END_REF][START_REF] Gimenez | Marine community Ecology: importance of trait-mediated effects propagating through complex life cycles[END_REF]. In fjords, however, seasonal stratification of the water column and restriction of deep-water movement by sills in fjords may contribute to the slow recolonization of benthic organisms. The sills of glacial fjords like Port Valdez interfere with larval transport and recruitment, thereby preventing some larvae from reaching the deep basin while sluggish water movement may slow their dispersal within the basin [START_REF] Muench | Physical oceanography[END_REF]Colonell, 1980a;[START_REF] Burd | Punctuated recovery of sediments and benthic infauna: a 19-year study of tailings deposition in a British Columbia fjord[END_REF][START_REF] Molinet | Effects of sill processes on the distribution of epineustonic competent larvae in a stratified system of Southern Chile[END_REF]. As a result, it can be difficult to predict recovery processes or the time frame required for benthic communities to return to prior environmental conditions and ambient faunal assemblages in fjords. Interactions between fauna, climate variations, and physical characteristics may be even more important when mature fauna are unavailable for larval production and larvae must come from distant sources as might occur with a very large, landscape-scale perturbation. Such a large reduction of mature, reproducing macrofauna populations may have occurred in PWS following the 1964 earthquake, as suggested by the reduced abundance of fauna in 1971 and delayed recruitment of some otherwise ubiquitous polychaetes in Port Valdez.

Anthropogenic Stressors

Migratory salmon are an important link between oceans, estuaries, and watersheds wherever they occur [START_REF] Kline | Recycling of elements transported upstream by runs of Pacific Salmon: I. 15 N and 13 C evidence in Sashin Creek, Southeastern Alaska[END_REF][START_REF] Kline | Recycling of elements transported upstream by runs of Pacific Salmon: II. 15 N and 13 C evidence in Kvichak River watershed, Bristol Bay, Southwestern Alaska[END_REF][START_REF] Finney | Impacts of climatic change and fishing on pacific salmon abundance over the past 300 years[END_REF][START_REF] Naiman | Pacific Salmon, nutrients, and the dynamics of freshwater and riparian ecosystems[END_REF][START_REF] Ray | Connectivities of estuarine fishes to the coastal realm[END_REF]. Evidence for a fjord-wide influence of salmon released from a shoreline hatchery in Port Valdez includes the positive correlations of salmon releases and returns with macrofauna biomass shown in this study (ρ = 0.35 and 0.40, respectively) and the ubiquitous presence of fish bones and otoliths in sediment from the deep benthos (Table 2; A. L. Blanchard, personal observations). The large number of juveniles released into Port Valdez undoubtedly influence benthic systems by feeding on zooplankton thereby altering pathways of energy exchange between ecosystem components with unknown effects on the benthic components of marine ecosystems. Carcasses of adult salmon may directly influence benthic communities by providing an immediate food source for scavenging macrofauna but most carcasses are likely deposited intertidally or in the shallow subtidal with fewer reaching the deep benthos so the pathways through which carcasses may influence fauna in the basin of a fjord are not entirely clear. In freshwater systems, nutrients from salmon carcasses (dissolved or particulate marine-derived nutrients) are dispersed through numerous ecosystem components and multiple pathways [START_REF] Naiman | Pacific Salmon, nutrients, and the dynamics of freshwater and riparian ecosystems[END_REF][START_REF] Ray | Connectivities of estuarine fishes to the coastal realm[END_REF]. Similarly, dissolved nutrients and particulates released by decomposing salmon from the intertidal and shallow subtidal are likely to be dispersed throughout a fjord. Particulate organics and dissolved nutrients indirectly reaching the benthos (e.g., transfers of particulates through water layers as opposed to direct carcass falls) would be available to deposit-feeding organisms in the basin of Port Valdez and the lack of scavenging organisms in the deep benthos suggests that such indirect transfers may be greatest. While the exact pathways for energy and nutrient flow of salmon carcass decomposition products to the deep basin are unknown for Port Valdez, the evidence suggest that a moderate flow (likely through multiple, indirect pathways) is present, as shown from studies of freshwater systems [START_REF] Naiman | Pacific Salmon, nutrients, and the dynamics of freshwater and riparian ecosystems[END_REF][START_REF] Ray | Connectivities of estuarine fishes to the coastal realm[END_REF].

The correlative evidence for associations between faunal measures from the deep basin of Port Valdez and salmon releases and returns suggests small effects at most rather than large changes like the gross organic enrichment caused by fish farming and eutrophication elsewhere (Josefson and Rosenburg, 1988;[START_REF] Muslow | Sediment profile imaging (SPI) and microelectrode technologies in impact assessment studies: example from two fjords in Southern Chile used for fish farming[END_REF]. Dominance by opportunistic species or complete loss of species resulting from organic enrichment was not observed in the deep basin of Port Valdez, although deposition of fish wastes from a processing plant adjacent to the City of Valdez led to such conditions in a limited shallow subtidal area (Blanchard and Feder, 2003). Given the tendency for increasing stress in coastal waters [START_REF] Schiff | Southern California[END_REF], effects from hatchery salmon fry releases (~200 million annually) and returns (~10 million annually) [START_REF] White | Alaska salmon enhancement program 2007 annual report[END_REF]

in Port

Valdez could lead to large interactions in the future with other stressors and potentially stronger, negative effects.

Effects from discharge of treated tanker ballast water from a treatment plant at the marine oil terminal in Port Valdez were not observed in the deep basin. In the present study, average sediment PAH concentrations in the deep basin were low and below concentrations expected to result in alterations of faunal communities [START_REF] Leung | Deriving sediment quality guidelines from field-based species sensitivity distributions[END_REF]. Hydrocarbon accumulations and associated effects on fauna of Port Valdez are limited as most hydrocarbons are diluted within the water column due largely to the steep slope where hydrocarbons are discharged, depth of the fjord, oceanography, and low quantities of organics in sediment to which hydrocarbons may adsorb [START_REF] Shaw | Hydrocarbon accumulations[END_REF].

Spatially and temporally limited effects associated with increased PAH were observed in shallower waters within a few hundred meters of the point of discharge of treated ballast water but similar trends in fauna (e.g., increases or gradients in hydrocarbon-tolerant species near the oil terminal) were not noted for deep basin communities [START_REF] Shaw | Hydrocarbon accumulations[END_REF][START_REF] Blanchard | Long-term investigation of Benthic Fauna and the Influence of Treated Ballast Water Disposal in Port Valdez, Alaska[END_REF]Blanchard et al., , 2003)). As indicated by [START_REF] Wake | Oil refineries: a review of their ecological impacts on the aquatic environment[END_REF], interactions between effluents and the physical environment are important in the fate of hydrocarbons discharged into the marine environment and with adequate treatment, residual hydrocarbons have very limited effects. Effects of hydrocarbons on fauna in Sullom Voe, Shetland Islands, were strongly influenced by interactions with environmental gradients and effects were reduced where accumulations were prevented by hydrodynamics (although hydrocarbon deposition and associated effects increased at sites with greater organics) [START_REF] May | Effects of oil-industry operations on the macrobenthos of Sullom Voe[END_REF]. The lack of effects from hydrocarbons in the deep basin of Port Valdez was likewise due to a combination of physical characteristics, oceanography, and treatment of wastes at the ballast-water treatment facility [START_REF] Shaw | Hydrocarbon accumulations[END_REF]. . #Stat = number of stations per transect, abund. = abundance, PDO = average annual pacific decadal oscillation index, temp = annual temperature (C), snowfall = total annual snowfall (m), precip = total annual precipitation (m), and N/A = not available. Units for the variables are ind. m -2 for abundance, g m -2 for biomass, ng g -1 for PAH, and millions yr -1 for salmon releases and returns. Valdez, 1971Valdez, , 1976Valdez, , 1981Valdez, , 1985Valdez, , 1990Valdez, , 1995Valdez, , 2000Valdez, , and 2005. Data are ln(X+1)transformed average abundance of infaunal families (ind. m -2 ). Circles on the MDS plot represent year groupings. Average abundance and percent of total group abundance of dominant fauna, as suggested by SIMPER rankings, are presented for each year grouping.

  The multiple comparisons supported the observed groupings with R = 0.159 (p = 0.084) for the comparison between 1981 and 1985 suggesting weak differences at most, R 0.169 (p 0.212) for comparisons between the years 1990 to 2005 indicating no differences between years, and R 0.426 (p 0.036) indicating highly significant differences for all other comparisons. The MDS ordination also supported an additional separation of eastern and western fjord stations (as subgroups) as the ANOSIM overall R = 0.583 (p = 0.001) for an overall comparison of east vs. west stations (see Fig. 1 for definition of division boundaries). Support from multiple comparisons was weaker due to sample size limitations. For the east vs. west comparison by year, R = 0.583 for 1971 (small sample with too few permutations for accurate determination of the p-value), R = 0.583 for 1976 (small sample), R = 0.319 (p = 0.011) for 1981 to 1985, and R = 0.248 (p = 0.033) for 1990 to 2005. Numerically dominant families determined by SIMPER for the multivariate year groupings demonstrated large temporal changes. In 1971, dominant species in the fjord included the polychaete family Capitellidae (frequently Heteromastus filiformis), a family known for opportunistic behavior, which comprised 21% and 34% of overall abundance in 1971 and 1976, respectively, with abundance dropping to 6% of total abundance in

  polychaete family Spionidae was moderately abundant in 1971 comprising 12% of overall abundance but was infrequent in later years as abundance was 1% or less of total abundance in 1976 and later. Cumaceans of family Leuconidae (largely Eudorella emarginata) comprised 11% of abundancein 1971, 7% in 1976, and 2% from 1981 to 2005. The abundance of polychaetes of family Lumbrineridae (Lumbrineris spp. and Ninoe gemmea) was moderate in 1971 comprising 9% of total abundance, higher in 1976 at 11% of total abundance but declined to4% in 1981 and 1985 and 6% in 1990 to 2005. Bivalves of family Thyasiridae (Adontorhina cyclia, Axinopsida serricata, and infrequently, Thyasira flexuosa) were not abundant in 1971 comprising only 2% of abundance but increased to7% in 1976, 24% in 1981 and 1985, and 18% in 1990 to 2005. The polychaetes Oweniidae (primarily Galathowenia oculata) and Trichobranchidae (mostly Terebellides stroemi) were absent from the survey transect in 1971. Oweniid polychaetes were also absent in 1976 but present in 1981 and 1985 with 4% of total abundance and 6% in 1990 to 2005. Polychaetes of family Trichobranchidae were 4% of total abundance in 1976, 8% in 1981 and 1985, but only 2% in 1990 to 2005. The polychaete family Nephtyidae (numerically dominated by Nephtys punctata) was present all years in moderate abundance comprising 10% of total abundance in 1971, 4% abundance in 1976, 4% of total abundance in 1981 and 1985, and 6% total abundance in 1990 to 2005. Changes in abundance and distributions of selected macrofauna families demonstrated significant changes over time between the stations sampled over the whole of the fjord in 1971 and 2002. Bubble plots of abundance demonstrate that the polychaetes Maldanidae, Paraonidae, and Trichobranchidae were sparse in the fjord in 1971 but abundant and distributed throughout the fjord in 2002 (Fig. 4). The bivalves Thyasiridae occurred throughout the fjord in 1971 but in low abundance, as compared to 2002. The change in abundance for these taxon groups represents a 15-fold increase in abundance from 1971 to 2002 for Maldanidae and 8-fold, 13-fold, and an 8-fold increases for Paraonidae, Trichobranchidae, and Thyasiridae, respectively. Two-way ANOVA indicated significantly lower abundance values in 1971 than in 2002 for Maldanidae, Trichobranchidae, and Thyasiridae (Fig. 4). Average abundance of Paraonidae was significantly lower for the western stations than for the eastern sites and in 1971 as compared to 2002. The cumacean family Leuconidae (mostly Eudorella emarginata) was abundant in most of the fjord in 1971 with lower abundance in the eastern fjord but in 2002, the cumaceans were more abundant along the margins of the fjord than in the basin (Fig. 4). Leuconidae abundance in 1971 was 2 times that of the abundance found in 2002 with significantly higher abundance at the western stations in 1971 as compared to the other categories. Ophiuroids (mostly Ophiura sarsi) were common in 1971 in the western deep basin but were absent in 2002 (Fig. 4). In both years, the ophiuroid family Amphiuridae was moderately abundant at a few sites and particularly abundant adjacent to a small, outwash delta on the southern edge of the fjord. One polychaete family, Oweniidae (largely G. oculata) demonstrated particularly dramatic temporal changes in distribution. Although oweniid polychaetes were absent from the deep basin of the fjord in 1971, they occurred in very low abundance at one site in Valdez Arm just outside the fjord (Fig. 5). Comparisons of rank-transformed abundance indicated significant differences between 1971 and 2002. Oweniid polychaetes were still relatively uncommon in 1976, more common in 1982, and by 2002, they were a dominant member of the benthic community. Spatio-temporal modeling demonstrated significant trends in abundance, biomass, and the number of taxa from 1971 to 2007. The kriging plot (the graphic summary of the spatio-temporal model) for abundance indicated increasing abundance over time and with distance from the west (station 50 near the mouth of the fjord: 0 km) to the east (station 11 near the head of the fjord: 14 km) with low abundance indicated in the lower left of the plot by orange and yellow and high values in the upper right corner shown in blue and purple (Fig. 6). The whisker plot by year demonstrated more variable abundance through 1990 (SD = 139.1 ind. m -2 based on overall means; Table 1), as compared to 1991 to 2007 (SD = 77.1 ind. m -2 ), with significantly lower values (shown by the non-overlapping confidence intervals) generally occurring from 1971 to 1990 and an overall trend of increasing abundance. The whisker plot by distance indicates significantly lower abundance mid-fjord and significantly higher abundance at distances > 11 km. The kriging plot and the whisker plot by distance for biomass demonstrated

Fauna

  in the deep basin (>150 m) of Port Valdez appear to be sensitive to decadal-scale climatic variability although the pathways through which the PDO may have influenced fauna are unknown. Relationships between the PDO and macrofauna are suggested by the moderate associations of faunal summary measures and the PDO index in Port Valdez (ρ = 0.38 and -0.31, for abundance and biomass, respectively; Table

  Port Valdez from 1971 to 2007 suggest an ecosystem in flux. Although the macrofaunal assemblage in the deep basin was generally similar to that of glacial fjords elsewhere, trends in faunal abundance and community structure deviated significantly from the stability expected for glacial fjords. Natural factors were the strongest drivers of spatial and temporal change. Effects on biotic communities from the catastrophic disturbance of the 9.2 magnitude earthquake in 1964 were large and faunal readjustment was the dominant source of variability in the deep benthos from 1971 to 1989. The readjustment process was demonstrated through introductions and subsequent dispersal of macrofauna families in the basin and trends in the spatio-temporal models. Following readjustment from the earthquake, long-term climatic trends were important sources of variability for benthic fauna from 1989 to 2007. Spatial trends in faunal indices were associated with the sediment gradient from the head to the mouth of the fjord as well. Physical characteristics of the fjord (e.g., the sills, basin depth, limited water movement in the basin, and sedimentation rate) were important structuring factors of faunal communities spatially and temporally during the period of faunal readjustment. Releases of juvenile salmon at the hatchery since 1982 and returns of adult salmon were positively associated with macrofaunal biomass suggesting enhancement of the deep benthic fauna. The Valdez marine oil terminal does not measurably contribute to faunal variations in the deep basin largely due to the fjord's physical and oceanographic characteristics. Overall, benthic communities in the deep basin were in a state of flux up to 1990 while a state of greater stability was reached in 1990 that more closely reflected the stability observed in fjord-basin fauna of undisturbed fjords elsewhere.

Figure 1 .

 1 Figure 1.Sites sampled for benthic fauna in PortValdez, Alaska, 1971-2007. The 1971 and 2002 grid sites are shown with the 1971-2007 transect sites identified by their station numbers. The dotted line represents the division between eastern and western stations as described by[START_REF] Feder | Subtidal benthos[END_REF].

Figure 2 .

 2 Figure 2. Sediment flux in Port Valdez, 1976 To 1978.

Figure 3 .

 3 Figure 3. MDS ordination of abundance data for the benthos transect sites from PortValdez, 1971Valdez, , 1976Valdez, , 1981Valdez, , 1985Valdez, , 1990Valdez, , 1995Valdez, , 2000Valdez, , and 2005. Data are ln(X+1)transformed average abundance of infaunal families (ind. m -2 ). Circles on the MDS plot represent year groupings. Average abundance and percent of total group abundance of dominant fauna, as suggested by SIMPER rankings, are presented for each year grouping.

Figure 4 .

 4 Figure 4.Bubble and bar plots of the abundance (ind.m -2 ) of selected fauna from PortValdez, 1971 and 2002. Bar plots include 95% confidence intervals of average abundance (ind. m -2 ). The dashed lines in the plots represent the division between eastern and western stations. P-values from ANOVA are given following the comparisons for year (Y), location (L), and the interaction (YxL) Figure 5. Bubble plot of the abundance (ind.m -2 ) of the polychaete Galathowenia oculata of the family Oweniidae in Port Valdez, 1971 to 2002. The thin dashed lines represent the boundaries for sampling in each year. Abundance data were ranktransformed prior to ANOVA. See legend of Figure 4 for further details.

Figure 6 .

 6 Figure 6.Spatio-temporal model and whisker plots of abundance (ind. m -2 ) for deep stations fromPort Valdez, 1971-2007. Distances are from station 50 (0 km) in the west near the mouth to station 11 (14 km) in the east near the head of the fjord. Whisker plots are mean abundance with 95% confidence intervals. The dotted lines in the whisker plots are the grand means and the solid lines are estimated curve lines to show the general trends.

Figure 7 .

 7 Figure 7. Spatio-temporal model and whisker plots of biomass (g m -2 ) for deep stations from Port Valdez, 1971-2007. See legend of Figure 6 for details.

Figure 8 .

 8 Figure 8. Spatio-temporal model and whisker plots of the number of taxa for deep stations from Port Valdez, 1971-2007. See legend of Figure 6 for details.

Figure 9 .

 9 Figure 9. Plot of mean abundance (ind. m -2 ) for the deep-basin station transect from Port Valdez and the annual mean pacific decadal oscillation index (PDO), 1971-2007. The dashed line indicates 0 for the PDO index. PDO values above the line indicate warm years and values below the line cold years. Major events occurring in the fjord are listed and indicated on the horizontal axis of the plot by an asterisk.

Table 1 .

 1 Summary measures from Port Valdez, Alaska

Table 1

 1

Table 2 .

 2 Spearman correlations ( ) of benthic indices, measures of anthropogenic stress, and climatic variables, fromPort Valdez, Alaska 1982-2007. Values in bold are significant at = 0.10.

	Table 2
	Variables	PAH Releases Returns PDO Temperature Snowfall Precipitation
	Abundance 0.29	0.23	0.12	0.38	0.25	-0.02	-0.05
	Biomass	0.12	0.35	0.40	-0.31	0.03	0.26	-0.28
	# Taxa	0.05	-0.27	-0.24	0.06	0.07	0.13	0.06


Khurelbaatar Gonchigdanzan

Kamil M Kosiński
email: k.m.kosinski@uva.nl

On the functional limits for partial sums under stable law

Keywords: Almost sure limit theorem, logarithmic average, stable law, product of partial sums

For the partial sums (S n ) of independent random variables we define a stochastic process s n (t

, for some sequence (d n ) and distribution G t . We also prove an almost sure functional limit theorem for the product of partial sums of i.i.d. positive random variables attracted to an α-stable law with α ∈ (1, 2].

Introduction and main result

In the past two decades many interesting extensions of the classical central limit theorem (CLT) have been obtained. One of the extensions is known as almost sure central limit theorem (ASCLT) which is discovered by Brosamler (1988) and Schatte (1988) and has been extensively studied for independent random variables as well as dependent variables. Motivated by ASCLT, almost sure versions of many limit theorems in probability and statistics have been obtained in the past. It is known that for i.i.d. r.v.'s ASCLT holds under the same assumptions as CLT but in general, the existence of the weak limit does not always imply the almost sure limiting result. For more discussions about the early results on ASCLT we refer to [START_REF] Berkes | Results and problems related to the pointwise central limit theorem[END_REF].

In this note we consider the product of partial sums, denoted by S n , of a sequence of random variables attracted to a stable distribution and its limit distributions. [START_REF] Rempała | Asymptotics for products of sums and U-statistics[END_REF] established the limit distribution of the product of partial sums of a sequence of i.i.d. positive r.v.'s with mean µ and variance σ 2 :

n k=1 S k n!µ n µ/(σ √ n) d → e √ 2N(0,1) . (1)
Zhang and Huang (2007) proves a weak invariance principle of (1) for i.i.d. r.v.'s. Recently, [START_REF] Kosiński | On the functional limits of a function of partial sums[END_REF] has shown that the weak invariance principle still holds when the partial sums are attracted to an α-stable law with α ∈ (1, 2] which also generalizes the earlier result by [START_REF] Qi | Limit distributions for products of sums[END_REF].

Throughout this paper, log log x and log x stand for ln ln(max{x, e e }) and ln(max{x, e}) respectively. We also use the notations a n b n for a n = O(b n ) and I(A) for the indicator function on a set A.

Our main result in this note is to establish an almost sure version of the result by [START_REF] Kosiński | On the functional limits of a function of partial sums[END_REF] that can generalize the early results by [START_REF] Gonchigdanzan | A note on the almost sure limit theorem for the product of partial sums[END_REF] and [START_REF] Gonchigdanzan | An almost sure limit theorem for the product of partial sums with stable distribution[END_REF][START_REF] Gonchigdanzan | Almost sure functional limit theorem for the product of partial sums[END_REF].

Recall that a sequence of and i.i.d. r.v.'s {X n : n ≥ 1} is said to be in the domain of attraction of a stable law L if there exist sequences (a n ) and (b n ) such that

S n -b n a n d → L,
where L is one of the stable distributions with index α ∈ (0, 2].

The following theorem is well known (see e.g. [START_REF] Hall | A comedy of error: the canonical form for a stable characteristic function[END_REF].

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Theorem 1.1 (Stability Theorem). The general stable law is given, to within type, by a characteristic function of one of the following forms:

1. φ(t) = exp(-t 2 /2) (normal case, α = 2); 2. φ(t) = exp(-|t| α (1 -iβ sgn(t) tan( 1 2 πα))) (0 < α < 1 or 1 < α < 2, -1 ≤ β ≤ 1); 3. φ(t) = exp(-|t|(1 + iβ sgn(t)2/π log |t|) (α = 1, -1 ≤ β ≤ 1).
It is worth mentioning that in Theorem 1.1, β is the skewness parameter. In our paper, β = 1 since X 1 is a positive random variable.

The first result of this note is the following almost sure functional limit theorem:

Theorem 1.2. Let {X n : n ≥ 1} be a sequence of i.i.d. positive random variables with EX 1 = µ in the domain of attraction of an α-stable law L α with α ∈ (1, 2] and characteristic function as in Theorem 1.1. Define a process

{π n (t) : 0 ≤ t ≤ 1} by π n (t) :=         [nt] k=1 S k µk         µ/a n ,
where (a n ) is a sequence of positive numbers that satisfies (S nµn)/a n d → L α as n → ∞. Then for any real x

1 log N N n=1 1 n I (π n (t) ≤ x) a.s. → F t (x) as N → ∞
where F t is the distribution function of the random variable exp

t 0 L α (s) s ds . Remark 1.3. If X 1 has finite variance equal to σ 2 then α = 2, L α d = N(0, 1) and a n ∼ σ √ n, thus Theorem 1.2 implies
the main result of [START_REF] Gonchigdanzan | Almost sure functional limit theorem for the product of partial sums[END_REF] which in particular yields the result of [START_REF] Gonchigdanzan | A note on the almost sure limit theorem for the product of partial sums[END_REF] Theorem 2 since it is easy to verify that

1 0 L α (s) s ds d = √ 2 N(0, 1).
Moreover, [START_REF] Kosiński | On the functional limits of a function of partial sums[END_REF] showed that for any α ∈ (1, 2]

1 0 L α (s) s ds d = (Γ(α + 1)) 1/α L α ,
hence Theorem 1.2 also yields the result of [START_REF] Gonchigdanzan | An almost sure limit theorem for the product of partial sums with stable distribution[END_REF] Theorem 1.1.

Our next result is the following Berkes-Dehling type of theorem (Theorem 2 [START_REF] Berkes | Some limit theorems in log density[END_REF].

Theorem 1.4. Let {Y n : n ≥ 1} be a sequence of independent random variables and

S n = Y 1 + • • • + Y n . Let (d n ) be a sequence of positive numbers such that d l d k l k γ (l ≥ k ≥ n 0 ) (2)
for some γ > 0 and n 0 ≥ 1 and

E S n -µn d n e γ (log n) 1-ε (3)
for some constant µ and γ ∈ (0, γ). Then for any distribution G t ,

1 log N N n=1 1 n I 1 d n [nt] k=1 S k k -µ ≤ x a.s. → G t (x) as N → ∞ (4)
if and only if

1 log N N n=1 1 n P 1 d n [nt] k=1 S k k -µ ≤ x → G t (x) as N → ∞. (5)
The following three lemmas are needed for the proof of our main result.

Lemma 2.1 (Lemma 2.3, [START_REF] Gonchigdanzan | An almost sure limit theorem for the product of partial sums with stable distribution[END_REF]). Under the assumption of Theorem 1.2 we have

µ a n [nt] k=1 log S k µk - 1 a n [nt] k=1 S k k -µ a.s.
→ 0 as n → ∞.

Lemma 2.2. Under the assumptions of Theorem 1.2 we have

1 a n [nt] k=1 S k k -µ d → t 0 L α (s) s ds in D[0, 1].
Proof This is a particular case of Theorem 2 in [START_REF] Kosiński | On the functional limits of a function of partial sums[END_REF] when f (x) = x.

Lemma 2.3. Under the assumptions of Theorem 1.4 we have

E 1 d n max 1≤k≤n k j=1 log n + 1 j (Y j -µ) log n E 1 d n max 1≤k≤n |S k -kµ| log n e γ (log n) 1-ε .
Proof The first part is Lemma 1 in [START_REF] Gonchigdanzan | Almost sure functional limit theorem for the product of partial sums[END_REF] valid for any sequence of random variables. The second part is Lemma 1 in [START_REF] Berkes | Some limit theorems in log density[END_REF] combined with the assumption (3).

Proofs of the main results

To prove Theorem 1.2 we need the result in Theorem 1.4. Let us prove Theorem 1.4 first, then Theorem 1.2 for convenience.

Proof of Theorem 1.4 According to [START_REF] Berkes | Results and problems related to the pointwise central limit theorem[END_REF]Dehling (1993) (p.1647) it suffices to prove that for any bounded Lipschitz function g on D[0, 1] we have

1 log n n k=1 1 k g s k d k -Eg s k d k a.s. → 0 as n → ∞, (6)
where

s n := s n (t) = k≤[nt] (S k /k -µ).
It turns out that the following estimate is indeed sufficient for (6) (see p.1648 [START_REF] Berkes | Some limit theorems in log density[END_REF], for the proof):

E n k=1 1 k ξ k 2 log 2 n(log log n) -1-ε for some ε > 0, (7)
where

ξ k = g(s k /d k ) -Eg(s k /d k ). Observe that n k=1 (S k /k -µ) = n k=1 b k,n (Y k -µ) where b k,n = n j=k 1/ j. It can be easily seen that s l -s k = b [kt]+1,[lt] (S [kt] -[kt]µ) + b [kt]+1,[lt] (Y [kt]+1 -µ) + • • • + b [lt],[lt] (Y [lt] -µ) for l ≥ k. Obviously s l -s k -b [kt]+1,[lt] (S [kt] -µ[kt]) is independent of s k , so we get Cov g s k d k , g s l -s k -b [kt]+1,[lt] (S [kt] -µ[kt]) d l = 0 for l ≥ k.
Since g is a bounded Lipschitz it follows that

E(ξ k ξ l ) = Cov g s k d k , g s l d l -g s l -s k -b [kt]+1,[lt] (S [kt] -µ[kt]) d l E max 0≤t≤1 |s k + b [kt]+1,[lt] (S [kt] -µ[kt])| d l ≤ E max 0≤t≤1 |s k | d l + E max 0≤t≤1 |b [kt]+1,[lt] (S [kt] -µ[kt])| d l = d k d l E max 0≤t≤1 |s k | d k + E max 0≤t≤1 b [kt]+1,[lt] |S [kt] -µ[kt]| d k .
Moreover, noticing max 0≤t≤1 b [kt]+1,[lt] = log(l/k) and applying Lemma 2.3 we get

E(ξ k ξ l ) d k d l         E max 0≤t≤1 1 d k [kt] i=1 b i,k (Y i -µ) + log(l/k)E max 0≤t≤1 |S [kt] -µ[kt]| d k         = d k d l         E max 0≤ j≤k 1 d k j i=1 b i,k (Y i -µ) + log(l/k)E max 0≤ j≤k |S j -µ j| d k         log l d k d l E max 1≤ j≤k |S j -µ j| d k log l k l γ log ke γ (log k) 1-ε =: c k,l .
On the other hand we also have E(ξ k ξ l ) 1 because ξ k is bounded. Hence we estimate E(ξ k ξ l ) as follows:

E(ξ k ξ l ) 1, if l/k ≤ exp (log n) 1-ε c k,l , if l/k ≥ exp (log n) 1-ε
where ε is any positive number.

Thus we get

1≤k≤l≤n l/k≤exp ((log n) 1-ε ) E(ξ k ξ l ) kl ≤ 1≤k≤n 1 k k≤l≤ke (log n) 1-ε 1 l n k=1 1 k log 1-ε n log 2-ε n (8) and 1≤k≤l≤n l/k≥exp ((log n) 1-ε ) E(ξ k ξ l ) kl ≤ log 2 ne γ (log n) 1-ε 1≤k≤l≤n l/k≥exp ((log n) 1-ε ) 1 kl k l γ ≤ log 2 ne (γ -γ)(log n) 1-ε 1≤k≤l≤n 1 kl log 4 ne (γ -γ)(log n) 1-ε log 2-ε n, (9)
where the last estimation follows because γ ∈ (0, γ).

Since

E n k=1 1 k ξ k 2 1≤k≤l≤n 1 kl |E(ξ k ξ l )|
by ( 8) and (9) it follows (7).

Before proving Theorem 1.2, recall that it is well known that the sequence (a n ) in Theorem 1.2 can be written as

a n = n 1/α L(n)
where L is a slowly varying function.

Proof of Theorem 1.2 We first show the equivalence of (4) and ( 5) under the conditions of Theorem 1.2 setting d n := a n . In fact (3) is a direct consequence of Theorem 6.2 in [START_REF] Deacosta | Convergence of moments and related functionals in the general central limit theorem in Banach spaces[END_REF]. ( 2) can be easily verified using the facts that a n = n 1/α L(n) and L(k)/L(n) (k/n) ε for any ε > 0 where L is a slowly varying function. Thus by

ACCEPTED MANUSCRIPT


N Balakrishnan

E Cramer

K Davies

Some results on order statistics generated by two simulation methods

Keywords: Order statistics, simulation methods, joint distribution, positive quadrant dependence, exceedence probabilities, copula, precedence-type tests

In this paper, we consider joint distributions of order statistics generated by two simulation methods. By using these distributions, we study the nature of dependence and exceedence probabilities between them.

Introduction

Let U 1 , • • • , U n be iid uniform random variables, and let us denote the corresponding r-th order statistic by U r:n and its cumulative distribution function (cdf) by F r:n . Then, for simulating U r:n , we have two options. First is simply to order U 1 , • • • , U n and take the r-th smallest value. Second is to take the first uniform random variable U 1 and apply a quantile transformation to define U * r:n = F -1 r:n (U 1 ); see [START_REF] Nagaraja | Some relations between order statistics generated by different methods[END_REF] and [START_REF] Arnold | A First Course in Order Statistics, Classic Edition[END_REF]. Evidently, the random variables U r:n and U * r:n have the same distribution, but the generation process is clearly different. Moreover, they are dependent since both definitions involve a common random variable U 1 .

In this paper, we are interested in the nature of dependence between U r:n and U * r:n . In Section 2, we derive the joint distribution of U r:n and U * r:n , and show that they are positive quadrant dependent. In Section 3, we derive explicit expressions for exceedence probabilities of the form π * r:n = P (U * r:n ≤ U r:n ) and establish some associated properties. In Section 4, we present some asymptotic approximations to exceedence probabilities.

In the subsequent derivations, we make use of the following relations for cdf's of uniform order statistics, which can be found in [START_REF] David | Robustness of location estimators in the presence of an outlier[END_REF] (see also Arnold et al. 2008, p. 40) :

F r-1:n-1 (t) = F r:n-1 (t) + n -1 r -1 t r-1 (1 -t) n-r (1.1) = F r:n-1 (t) + 1 n f r:n (t), t ∈ [0, 1], F r:n (t) = F r:n-1 (t) + n -1 r -1 t r (1 -t) n-r (1.2) = F r:n-1 (t) + 1 n tf r:n (t), t ∈ [0, 1].
2 Joint distribution of U r:n and U * r:n

We shall denote the order statistics of the sample

U 2 , • • • , U n by V r:n-1 for 1 ≤ r ≤ n -1. Note that the order statistics V 1:n-1 , • • • , V n-1:n-1 are independent of U 1 .
For convenience, let us also define

F n:n-1 (x) = 0 ∀ x ∈ [0, 1].
Theorem 2.1 The joint cdf of U * r:n and U r:n is given by

F * r:n (s, t) = P (U * r:n ≤ s, U r:n ≤ t) = F r-1:n-1 (t)F r:n (s), t ≥ F r:n (s) F r:n-1 (t)F r:n (s) + 1 n tf r:n (t), t < F r:n (s) , 0 ≤ s, t ≤ 1. (2.1)
Proof The r-th order statistic U r:n in the sample U 1 , • • • , U n is related to the order statistics in the subsample U 2 , • • • , U n in the following way:

U r:n =      V r-1:n-1 , U 1 < V r-1:n-1 , U 1 , V r-1:n-1 ≤ U 1 < V r:n-1 , V r:n-1 , V r:n-1 ≤ U 1 .
So, for 0 ≤ s, t ≤ 1, the joint cdf can be expressed as

P (U * r:n ≤ s, U r:n ≤ t) = P (U 1 ≤ F r:n (s), U r:n ≤ t) = P (U 1 ≤ F r:n (s), U r:n ≤ t, U 1 < V r-1:n-1 ) + P (U 1 ≤ F r:n (s), U r:n ≤ t, V r-1:n-1 ≤ U 1 < V r:n-1 ) + P (U 1 ≤ F r:n (s), U r:n ≤ t, V r:n-1 ≤ U 1 ) = P (U 1 ≤ F r:n (s), V r-1:n-1 ≤ t, U 1 < V r-1:n-1 ) + P (U 1 ≤ F r:n (s), U 1 ≤ t, V r-1:n-1 ≤ U 1 < V r:n-1 ) + P (U 1 ≤ F r:n (s), V r:n-1 ≤ t, V r:n-1 ≤ U 1 ) = I 1 + I 2 + I 3 , say .
Now, we obtain the following representations for I 1 , I 2 and I 3 : F r:n-1 (u)du + F r:n-1 (t) F r:n (s) -min{F r:n (s), t} .

I 1 = Fr:n(s) 0 P (V r-1:n-1 ≤ t, u < V r-1:n-1 )du = min{Fr:n(s),t} 0 F r-1:n-1 (t) -F r-1:n-1 (u) du = F r-1:n-1 (t) min{F r:n (s), t} - min{Fr:n(s),t} 0 F r-1:n-1 (u)du, I 2 = min{Fr:n(s),t} 0 P (V r-1:n-1 ≤ u < V r:n-1 )du = min{Fr:n(s),t} 0 n -1 r -1 u r-1 (1 -u) n-r du (1.1)
Summing I 2 and I 3 , we obtain

I 2 + I 3 = min{Fr:n(s),t} 0 F r-1:n-1 (u)du + F r:n-1 (t) F r:n (s) -min{F r:n (s), t}
so that we finally find

I 1 + I 2 + I 3 = F r-1:n-1 (t) min{F r:n (s), t} + F r:n-1 (t)[F r:n (s) -min{F r:n (s), t}].
Thence, we have

P (U * r:n ≤ s, U r:n ≤ t) = F r-1:n-1 (t)F r:n (s), t ≥ F r:n (s) F r-1:n-1 (t)t + F r:n-1 (t)[F r:n (s) -t], t < F r:n (s) = F r-1:n-1 (t)F r:n (s), t ≥ F r:n (s) F r:n-1 (t)F r:n (s) + t[F r-1:n-1 (t) -F r:n-1 (t)], t < F r:n (s) = F r-1:n-1 (t)F r:n (s), t ≥ F r:n (s) F r:n-1 (t)F r:n (s) + n-1 r-1 t r (1 -t) n-r , t < F r:n (s) = F r-1:n-1 (t)F r:n (s), t ≥ F r:n (s) F r:n-1 (t)F r:n (s) + 1 n tf r:n (t), t < F r:n (s)
.

Hence, the theorem.

Remark 2.2 Noting that U * n:n = U 1/n 1

, for r = n, (2.1) yields

P (U 1/n 1 ≤ s, U n:n ≤ t) = t n-1 s n , t ≥ s n t n , t < s n , 0 ≤ s, t ≤ 1,
which is an expression given earlier by [START_REF] Nagaraja | Some relations between order statistics generated by different methods[END_REF]; see also Arnold et al. (2008, p. 106).

In particular, we find from the above representation the copula of U * r:n and U r:n . Since both marginal cdf's are given by F r:n and so, with s = F -1 r:n (u) and t = F -1 r:n (v), we find

C(u, v) = P (U 1 ≤ u, U r:n ≤ F -1 r:n (v)) = F r-1:n-1 (F -1 r:n (v))u, v ≥ F r:n (u) F r:n-1 (F -1 r:n (v))u + 1 n F -1 r:n (v)f r:n (F -1 r:n (v)), v < F r:n (u)
to be the copula of U * r:n and U r:n . Theorem 2.3 U * r:n and U r:n are positive quadrant dependent, and so they are positively correlated, in particular.

Proof U * r:n and U r:n are positive quadrant dependent (PQD) if

P (U * r:n ≤ s, U r:n ≤ t) ≥ P (U * r:n ≤ s)P r(U r:n ≤ t) = F r:n (s)F r:n (t), 0 ≤ s, t ≤ 1.
By using the result (Shaked and Shanthikumar 2007, p. 277)

U i:n ≤ st U j:m ⇐⇒ i ≤ j and n -i ≥ m -j,
we find U r-1:n-1 ≤ st U r:n so that F r-1:n-1 ≥ F r:n . Therefore, we get from (2.1), for t ≥ F r:n (s),

F * r:n (s, t) = F r-1:n-1 (t)F r:n (s) ≥ F r:n (s)F r:n (t). Let t < F r:n (s); then F * r:n (s, t) = F r:n-1 (t)F r:n (s) + 1 n tf r:n (t) ≥ F r:n (s)F r:n (t) iff F r:n (s) F r:n-1 (t) -F r:n (t) + 1 n tf r:n (t) ≥ 0. (2.2) Now, from (1.2), we have F r:n-1 (t)-F r:n (t) = -1 n tf r:n (t), and so (2.2) is equivalent to 1 n tf r:n (t)[1 -F r:n (s)] ≥ 0.
This obviously holds, which implies the positive quadrant dependence of U * r:n and U r:n . The positive correlation then readily follows. Using Hoeffding's expression for the covariance and

1 0 F r:n (s)ds = n-r+1 n+1 , we find 1 12 Cov(U * r:n , U r:n ) = 1 0 Fr:n(s) 0 F r:n-1 (t)F r:n (s) + 1 n tf r:n (t) dtds + 1 0 1 Fr:n(s) F r-1:n-1 (t)F r:n (s)dtds - 1 0 F r:n (s)ds 2 = 1 0 Fr:n(s) 0 {F r:n-1 (t) -F r-1:n-1 (t)} F r:n (s) + 1 n tf r:n (t) dtds + 1 0 1 0 F r-1:n-1 (t)F r:n (s)dtds - n -r + 1 n + 1 2 = 1 n 1 0 Fr:n(s) 0 f r:n (t) {t -F r:n (s)} dtds + (n -r + 1) 2 n(n + 1) - n -r + 1 n + 1 2 = 1 n 1 0 Fr:n(s) 0 f r:n (t) {t -F r:n (s)} dtds + (n -r + 1) 2 n(n + 1) 2 = - 1 n 1 0 Fr:n(s) 0 F r:n (t)dtds + (n -r + 1) 2 n(n + 1) 2 .
Since Cov(U * r:n , U r:n ) ≥ 0 by the PQD-property established in Theorem 2.3, we obtain the following bounds for the covariance:

0 ≤ Cov(U * r:n , U r:n ) ≤ 12 (n -r + 1) 2 n(n + 1) 2 ≤ 12 n ,
which readily reveals that the covariance tends to zero as n → ∞.

In this section, we obtain an expression for the exceedence probabilities π * r:n = P (U * r:n ≤ U r:n ) for r = 1, • • • , n. First, we get the following decomposition of π * r:n :

π * r:n = P (U * r:n ≤ U r:n ) = P (F -1 r:n (U 1 ) ≤ U r:n ) = P (F -1 r:n (U 1 ) ≤ V r-1:n-1 , U 1 < V r-1:n-1 ) + P (F -1 r:n (U 1 ) ≤ U 1 , V r-1:n-1 ≤ U 1 < V r:n-1 ) + P (F -1 r:n (U 1 ) ≤ V r:n-1 , V r:n-1 ≤ U 1 ) = J 1 + J 2 + J 3 , say . (3.1)
It turns out that the solution of the equation

t = F r:n (t), t ∈ (0, 1), (3.2)
is important in the calculation of the probabilities J 1 , J 2 and J 3 . F r:n is an S-shaped function which is unimodal, as seen in Figure 3.1. It is first convex and then concave with F r:n (0) = 0 and F r:n (1) = 1. For 1 < r < n, F r:n has exactly one fixpoint in the interval (0, 1), i.e., there exists exactly one solution y r:n of (3.2). Observe that (3.2) has no solution for r ∈ {1, n}. In these cases, let us define for convenience y 1:n = 0 and y n:n = 1. Geometrically, y r:n is the intersection point where the graph of F r:n intersects the line from (0, 0) to (1, 1). Moreover, (i) y r:n is strictly increasing in r;

(ii) y n-r+1:n = 1y r:n , r = 1, • • • , n;

(iii) For n odd and r = n+1 2 , we have y r:n = 1 2 ; (iv) For 1 ≤ r < n+1 2 , we have y r:n < 1 2 ; (v) For n+1 2 < r ≤ n, we have y r:n > 1 2 . Proof (i) Suppose that y r:n < y r-1:n . Then, we conclude from (3.3) that y r-1:n < F r:n (y r-1:n ). The cdf's F r:n are strictly ordered with respect to r, i.e., for fixed n, F r:n (t) < F r-1:n (t), t ∈ [0, 1] (see (1.1)). Therefore, F r:n (y r-1:n ) < F r-1:n (y r-1:n ) = y r-1:n which results in a contradiction. Thus, y r-1:n ≤ y r:n , 1 ≤ r ≤ n. Moreover, equality is not possible. For r = 1, we have y 1:n = 0 < y 2:n . Assuming that r ≥ 2, and that y r:n = y r-1:n , we deduce from (1.1) that

0 = F r:n (y r:n ) -F r-1:n (y r:n ) = - 1 n f r:n+1 (y r:n ) < 0.
Thus, we obtain a contradiction once again, and so y r:n is strictly increasing in r.

(ii) Using the fact that F r:

n (t) = 1 -F n-r+1:n (1 -t), t ∈ [0, 1], we have F r:n (t) = t ⇐⇒ 1 -t = F n-r+1:n (1 -t).
Therefore, α is a fixpoint of F r:n iff 1α is a fixpoint of F n-r+1:n . (iii) From (ii), we conclude for n odd and r = n+1 2 , that nr + 1 = r. Therefore, y r:n and 1-y r:n are fixpoints. Since we have exactly one fixpoint, we conclude that y r:n = 1 2 . (iv) Assertions (iv) and (v) follow immediately from (i) and (ii).

A C C E P T E D M

These properties of y r:n , 1 ≤ r ≤ n, are illustrated in Figure 3.1. Table 3.1 gives intersection points y r:n for n = 5, 10, 15, 20. The values are given only for r = 2, . . . , ⌊(n + 1)/2⌋ since y n-r+1:n = 1y r:n (see Part (ii) of Lemma 3.1). The computation of y r:n is difficult for large values of n, but we conjecture the following bounds. Notice that r-1 n-1 is the change point of F r:n . Conjecture 3.2 Let 2 ≤ r ≤ ⌊(n + 1)/2⌋ and w r:n be the smallest change point of f r:n given by

w r:n = r -1 n -1 - √ 2r 2 -r 2 n -rn -2r + rn 2 -n 2 + 2n (n -1)(n -2) .
Then,

w r:n ≤ y r:n ≤ r -1 n -1 .
Moreover, for (r n ) n with lim n→∞ rn n = α ∈ (0, 1), we have

lim n→∞ √ n r n n -y rn:n = α(1 -α).
Now, we are ready to establish the following theorem which gives a simple expression for the exceedence probability π * r:n = P (U * r:n ≤ U r:n ). Theorem 3.3 The exceedence probability π * r:

n = P (U * r:n ≤ U r:n ), 1 ≤ r ≤ n, is given by π * r:n = 1 2 + 1 2 n-1 r-1 2 2n-1 2r-1 - 2r-2 r-1 2n-2r n-r 2n-1 n-1 F 2r-1:2n-1 (y r:n ) (3.4) = 1 2 + 2r-2 r-1 2n-2r n-r 2n-1 n-1 2r -1 2n -F 2r-1:2n-1 (y r:n ) . (3.5)
Proof Suppose 1 < r < n. Then, we have from (3.1),

J 1 = 1 0 P (U 1 ≤ F r:n (u), U 1 ≤ u)f r-1:n-1 (u)du = 1 0 min{F r:n (u), u}f r-1:n-1 (u)du = yr:n 0 F r:n (u)f r-1:n-1 (u)du + 1 yr:n uf r-1:n-1 (u)du A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT J 2 = 1 0 P (u ≤ F r:n (u), V r-1:n-1 ≤ u < V r:n-1 )du = 1 yr:n P (V r-1:n-1 ≤ u < V r:n-1 )du = 1 n 1 yr:n f r:n (u)du = 1 n {1 -F r:n (y r:n )} = 1 n (1 -y r:n ), J 3 = 1 0 P (U 1 ≤ F r:n (u), u ≤ U 1 )f r:n-1 (u)du = 1 yr:n P (u ≤ U 1 ≤ F r:n (u))f r:n-1 (u)du = 1 yr:n {F r:n (u) -u} f r:n-1 (u)du.
From (1.1), we have

f r-1:n-1 (t) = f r:n-1 (t) + 1 n f ′ r:n (t), t ∈ (0, 1). (3.6)
So, by partial integration, we obtain

1 yr:n uf r-1:n-1 (u)du - 1 yr:n uf r:n-1 (u)du = 1 n 1 yr:n uf ′ r:n (u)du = - 1 n y r:n f r:n (y r:n ) - 1 n (1 -y r:n ).
Moreover, by (3.6), we have

yr:n 0 F r:n (u)f r-1:n-1 (u)du - yr:n 0 F r:n (u)f r:n-1 (u)du = 1 n yr:n 0 F r:n (u)f ′ r:n (u)du = 1 n F r:n (u)f r:n (u) yr:n 0 - 1 n yr:n 0 f 2 r:n (u)du = 1 n y r:n f r:n (y r:n ) - 1 n yr:n 0 f 2 r:n (u)du.
Thus, summing up J 1 , J 2 and J 3 , we arrive at

P (U * r:n ≤ U r:n ) = 1 0 F r:n (u)f r:n-1 (u)du - 1 n yr:n 0 f 2 r:n (u)du = P (U r:n ≤ V r:n-1 ) - 1 n yr:n 0 f 2 r:n (u)du.
Since, 1

n yr:n 0 f 2 r:n (u)du = 2r-2 r-1 2n-2r n-r 2n-1 n-1 F 2r-1:2n-1 (y r:n ), we get P (U * r:n ≤ U r:n ) = P (U r:n ≤ V r:n-1 ) - 2r-2 r-1 2n-2r n-r 2n-1 n-1 F 2r-1:2n-1 (y r:n ).
(3.7) Proof As above, we can write

P (U * r:n = U r:n ) = P F -1 r:n (U 1 ) = V r-1:n-1 , U 1 < V r-1:n-1 + P F -1 r:n (U 1 ) = U 1 , V r-1:n-1 ≤ U 1 < V r:n-1 + P F -1 r:n (U 1 ) = V r:n-1 , V r:n-1 ≤ U 1 ≤ P F -1 r:n (U 1 ) = V r-1:n-1 + P F -1 r:n (U 1 ) = U 1 + P F -1 r:n (U 1 ) = V r:n-1 = 0.
The second result can be seen as follows:

P (U * r:n , U r:n ) ∈ G = P (F -1 r:n (U 1 ), U r:n ) ∈ G = P U 1 = U r:n = 1 n .
Finally, we find

P (U * r:n , U r:n ) ∈ G > = P U * r:n > F -1 r:n (U r:n ) = P U 1 > U r:n = P U 1 > V r:n-1 = 1 0 F r:n-1 (t)dt = n -r n = 1 - r n .
In an analogous manner, we find P (U * r:n , U r:n ) ∈ G > = r-1 n . By combining the above results, we obtain the following theorem.

Theorem 3.7 For any 1 ≤ r ≤ n, we have

π * r:n = P (U * r:n ≤ U r:n ) = 1 -P (U * n-r+1:n ≤ U n-r+1:n ) = 1 -π * n-r+1:n .
In particular, for n odd and r = n+1 2 , we have so that in the case of the median the probability equals 1 2 .

P (U * r:n ≤ U r:n ) = 1 2 , A C C E P T E D M A N U S C
Proof According to Lemma 3.6, P F -1 r:n (U 1 ) = U r:n = 0 so that we do not have to take care of equality. Then, due to the relations U r:n = 1 -(1 -U ) n-r+1:n and F n-r+1:n (t) = 1 -F r:n (1t) for t ∈ [0, 1], we find

P (U * r:n ≤ U r:n ) = P F -1 r:n (U 1 ) ≤ U r:n = P 1 -F -1 n-r+1:n (1 -U 1 ) ≤ 1 -(1 -U ) n-r+1:n = P (1 -U ) n-r+1:n ≤ F -1 n-r+1:n (1 -U 1 ) = P U n-r+1:n ≤ F -1 n-r+1:n (U 1 ) = 1 -P U * n-r+1:n ≤ U n-r+1:n .
Note that both 1 -U and U have the uniform distribution on [0, 1]. The second assertion follows since, in this case, nr + 1 = n+1 2 = r and so both events are the same.

= P (U * r:n ≤ U r:n ) > 1 2 .
For r = 1, the result follows directly since y 1:n = 0. Let r ≥ 2. Then, from (3.5), it is sufficient to show that 2r

-1 2n ≥ F 2r-1:2n-1 (y r:n ), 2 ≤ r ≤ n + 1 2 .
The latter condition is equivalent to

E (U 2r-1:2n-1 ) ≥ F 2r-1:2n-1 (y r:n ) = P (U 2r-1:2n-1 ≤ y r:n ).
Finally, we conjecture that the sequence π * r:n is increasing-decreasing. Conjecture 3.9 π * r:n is increasing-decreasing in r f or 1 ≤ r ≤ n + 1 2 .

Lemma 4.1 Suppose y r:n is the unique solution of the equation F r:n (t) = t in the interval (0, 1). Moreover, let (r n ) n be a sequence such that lim n→∞ rn n = α ∈ (0, 1). Then, we have lim n→∞ y rn:n = α.

Proof Let t ∈ (0, 1). Suppose S n is a binomial random variable with parameters n and t. Then, with S * n = Sn-nt √ nt(1-t)

, we have

F rn:n (t) = P (S n ≥ r n ) = P S * n ≥ r n -nt nt(1 -t) = P 1 √ n S * n ≥ rn n -t t(1 -t) .
By central limit theorem, S * n is asymptotically standard normal so that Now, for any r and n, F r:n is a continuous function such that the equation F r:n (t) = t has exactly one solution y r:n in the interval (0, 1). Let δ > 0 and ε > 0. From (4.1), we find that there exists n 1 such that F rn:n (z) < δ for all n ≥ n 1 and z ∈ (0, αε).

Analogously, we find that there exists n 2 such that F rn:n (z) > 1δ for all n ≥ n 2 and z ∈ (α + ε, 1). Hence, the solution of F rn:n (t) = t must be contained in the interval (αε, α + ε) for all n ≥ max{n 1 , n 2 }. Since ε and δ are chosen arbitrarily, we find that lim n→∞ y rn:n = α. Since F 2rn-1:2n-1 (y rn:n ) ∈ [0, 1], we obtain the desired result.

F

  r-1:n-1 (u) -F r:n-1 (u) du,

F

  r:n (t) < t for t ∈ (0, y r:n ) and F r:n (t) > t for t ∈ (y r:n , 1). (3.3)For a given n, the solutions y r:n ∈ [0, 1] of Eq. (3.2), for r = 1, • • • , n, have the following properties:Lemma 3.1 Let n be fixed. Then:

Figure 3 . 1 :

 31 Figure 3.1: Plots of cumulative distribution functions of order statistics for n =5, 10, 15, 20.

Figure 3

 3 Figure 3.2: Scatterplot of 2000 simulated random pairs (U 6:15 , U * 6:15 ).

Conjecture 4. 2

 2 Let (r n ) n be a sequence such that lim n→∞ rn n = α ∈ (0, 1). Then, lim n→∞ F 2rn-1:2n-1 (y rn:n ) = v(α) for an appropriate function v : [0, 1] → [0, 1].Theorem 4.3 Let (r n ) n be a sequence such that lim n→∞ rn n = α ∈ (0, 1). Then,

Table 3

 3

	.2 gives exceedence probabilities π * r:n for n = 5, 10, 15, 20. Due to Theo-rem 3.7, we only present the values of π * r:n for r = 1, • • • , ⌊(n + 1)/2⌋.
	Conjecture 3.8 For any 1 ≤ r < n+1 2 , we have
	π * r:n

ACCEPTED MANUSCRIPT

The expression in (3.7) can be evaluated by using (1.2). As

, the theorem follows.

Remark 3.4 For r = n, we deduce

which has been derived by [START_REF] Nagaraja | Some relations between order statistics generated by different methods[END_REF]. For r = 1, with V j = 1 -U j , 1 ≤ j ≤ n, we find