Decline of the Black Sea Oxygen Inventory

We show that from 1955 to 2015, the inventory of oxygen in the Black Sea has decreased by 44% and the basin-averaged oxygen penetration depth has decreased from 140m in 1955 to 90m in 2015, which is the shallowest annual value recorded during that period. The oxygenated Black Sea surface layer separates the world’s largest reservoir of toxic hydrogen 5 sulphide from the atmosphere. The threat of chemocline excursion events led to hot debates in the past decades arguing on the vertical stability of the Black Sea oxic/suboxic interface. In the 1970s and 1980s, when the Black Sea faced severe eutrophication, enhanced respiration rates reduced the thickness of the oxygenated layer. Re-increasing oxygen inventory in 1985–1995 supported arguments in favor of the stability of the oxic layer. Concomitant with a reduction of nutrient loads, 10 it also supported the perception of a Black Sea recovering from eutrophication. More recently, atmospheric warming was shown to reduce the ventilation of the lower oxic layer by lowering Cold Intermediate Layer (CIL) formation rates. The debate on the vertical migration of the oxic interface also addressed the natural spatial variability affecting Black Sea properties when expressed in terms of depth. Here we show that using 15 isopycnal coordinates does not free from a significant spatial variability of oxygen penetration depth. By considering this spatial variability, the analysis of a composite historical set of oxygen profiles evidenced a significant shoaling of the oxic layer, and showed that the transient “recovery” of the 1990s was mainly a result of increased CIL formation rates during that period. As both atmospheric warming and eutrophication are expected to increase in the near future, 20 monitoring the dynamics of the Black Sea oxic layer is urgently required to assess the threat of further shoaling.


Introduction
The Black Sea deep waters constitutes the world's largest reservoir of toxic hydrogen sulphide.100 meters of ventilated surface waters are all that separate this reservoir from the atmosphere.This situation results from the permanent halocline (Öszoy and Ünlüata, 1997) that separates the surface layer (of low salinity due to river inflow) from the deeper layer (of high salinity due to inflowing Mediterranean seawater), restraining ventilation to the upper layer (Fig. 1).
In the lower part of the halocline, a permanent suboxic layer separates the Black Sea surface oxygenated waters ([O 2 ] > 20 µM) from the deep sulphidic waters ([H 2 S] > 20 µM) (Murray et al., 1989;Tugrul et al., 1992;Murray et al., 1995).More precisely, Murray et al. (1989) considered a threshold of 10 µM of oxygen because they analyzed high quality oxygen data.The threshold of 20 µM of oxygen was applied later to analyze historical oxygen data of lower quality.The upper (O 2 disappearance) and lower (H 2 S onset) interfaces of this suboxic layer are controlled by different biogeochemical and physical processes (Konovalov et al., 2006;Stanev et al., 2014), and undergo uncorrelated vertical migrations (Konovalov and Murray, 2001).Sinking organic matter is mainly respirated aerobically within the oxycline, ie. the lower part of the oxygenated layer where oxygen concentration decreases downwards to 20 µM.Increasing flux of organic matter, induced by a period of high nutrient load from the 70s to the late 80s, resulted in higher oxygen consumption above the suboxic layer and a shoaling of the upper suboxic interface (Codispoti et al., 1991;Konovalov and Murray, 2001;Tugrul et al., 2014).
After reduction of nutrient inputs around 1990 (Kroiss et al., 2006), the Black Sea was described as a recovering ecosystem (Mee et al., 2005;Oguz et al., 2006).This perspective was supported by improved eutrophication indices in the open sea (Kideys, 2002) as well as the stabilization of the upper suboxic interface in the 90s (Konovalov and Murray, 2001).However, the time scale of the expected recovery, ie. the time scale associated with the chain of biogeochemical mechanisms relating oxycline penetration depth to riverine nutrient loads, is not quantitatively understood.Several processes cause the oxycline depth to respond with a time lag to the reduction of riverine nutrient inputs.First, nutrients are mainly delivered to the northwestern shelf, where the accumulation of organic matter in the sediments buffers the riverine inputs, with slow diagenetic processes controlling and delaying the nutrient outflow across the seaward boundary (Capet et al., 2013).Second, the intermediate oxidation-reduction cycling of nitrogen, sulfur, manganese, iron and phosphorus that separates oxygen from hydrogen sulphide (Shaffer, 1986;Codispoti et al., 1991;Konovalov et al., 2006;Yakushev et al., 2007) can delay the response of the lower suboxic interface to changing nutrient fluxes by several years (Konovalov et al., 2006).
In addition to these biogeochemical factors, the dynamics of the upper and lower interfaces of the suboxic layer are controlled by physical processes (Konovalov et al., 2006;Stanev et al., 2014).In the Black Sea, dense waters formed by winter cooling and mixing Staneva and Stanev (2002) do not sink to the deepest layer, as in the Mediterranean sea, but accumulate on top of the permanent halocline (Fig. 1).The resulting Cold Intermediate Layer (CIL) is a major feature of the Black Sea vertical structure.Cold intermediate water formation and advection by the cyclonic basin-wide Rim Current (Öszoy and Ünlüata, 1997;Capet et al., 2012) ventilate the oxycline and thereby influence variability in the depth of the upper suboxic interface (Konovalov et al., 2006).Recently, atmospheric warming (Oguz et al., 2006) was shown to reduce the ventilation of the lower oxic layer (Tugrul et al., 2014;Pakhomova et al., 2014).At deeper levels, the dense sinking plume formed by the Mediterranean inflow through the Bosporus, which entrains water from the overlying CIL, injects fingers of oxygenated water directly into the deeper part of the suboxic layer and upper sulphidic layer and thus acts to control the depth of the lower suboxic interface (Konovalov and Murray, 2001;Konovalov et al., 2003;Glazer et al., 2006;Konovalov et al., 2006).
Previous long-term analyses of the vertical migration of the suboxic interfaces either ended (1955-1995(Konovalov and Murray, 2001)) or started (1985( -2015( (Pakhomova et al., 2014))) with the eutrophication period, excluding the large-scale overview required to grasp the interactions of eutrophication and climate factors.Those analyses lacked a comprehensive consideration of the natural spatial and seasonal variability of the vertical distribution of oxygen.
In the presence of large gradients, uneven data distribution may induce artificial signals when inter-annual trends are assessed from direct annual averages.In the stratified Black Sea, properties expressed in terms of depth coordinates (m) present a high spatial variability due to mesoscale features (Kempe et al., 1990) and to the general curvature of Black Sea isopycnals (Öszoy and Ünlüata, 1997;Stanev et al., 2014).As an alternative, using density (isopycnal levels, σ θ ) as vertical coordinate is generally considered a stable solution to assess the vertical migration of the chemocline on a decadal scale (Tugrul et al., 1992;Saydam et al., 1993;Murray et al., 1995).However, the spatial confinement of the lateral oxygen injections associated with the Bosporus plume, as well as the spatial variability of diapycnal ventilating processes (Zatsepin et al., 2007), imposes an horizontal structure to the oxygen penetration depth when expressed in terms of density (Stanev et al., 2004;Glazer et al., 2006).As this spatial gradient might scale with the temporal variations (a range of 0.17 kg m −3 was observed during the Knorr2003 campaign, Glazer et al. (2006)), it has to be considered when deriving interannual trends.
The present study describes the application of the DIVA (Data-Interpolating Variational Analysis) detrending procedure (Troupin et al., 2012;Capet et al., 2014) to untangle the temporal and spatial variability of three indices related to the Black Sea oxygenation status: the depth and density level of oxygen penetration and the oxygen inventory.These values were diagnosed from a composite historical dataset of oxygen vertical profiles.We reviewed the evolution of those indices through the past 60 years and discuss the respective controls of eutrophication and climate factors.

Data
We gathered a composite set of 4385 ship-based vertical profiles (oxygen, temperature and salinity) obtained between 1955 and 2005 in the Black Sea using CTD rosette bottles, continuous pumping profilers (Codispoti et al., 1991) and in situ analyzers (Glazer et al., 2006)  To complement the analysis of ship-based casts, we considered profiles originating from ten Argo autonomous profilers (May 2010-December 2015).Only good quality-checked real-time data were considered (Carval et al., 2014).Two of these floats (Argo ID 7900465 and 7900466) have been presented and discussed by Stanev et al. (2013), where the consistence and comparability of Argo and historical profiles is asserted within a 1µM error range.
Several studies address the error of Argo real-time oxygen data (e.g., Bittig and Körtzinger, 2015;Takeshita et al., 2013;Johnson et al., 2015).Demonstrating that the Black Sea real-time Argo data are precisely (ie, at fine scales) comparable with historical Winkler data, or identifying the relevant correction, is beyond the scope of the present study which addresses monthly to decadal time scales.
Evenly distributed small scales error (eg., difference between ascending and descending profiles due to sensor time response) were thus filtered by the temporal smoothing.However, a systematic error is not strictly excluded which could reach an underestimation of 10 µM (Virginie Thierry, IFRE-MER, personal communication, January 2016).Therefore, we evaluated a "worst-case" scenario in the analysis of Argo data by considering a systematic underestimation of oxygen concentration by 10µM.
Although most of the floats drifted along the basin periphery, some were also advected in the central part (Fig. 3).These trajectories highlight the range of spatial variability for the diagnostics described in Sect.2.2.

Profile analysis
From each profile we derived (1) the depth and (2) the potential density anomaly σ θ where oxygen concentration went below 20 µM and (3) the oxygen inventory, integrated above this limit (Fig. 1).The threshold value of 20 µM used to define the upper interface of the suboxic layer was suggested to compare oxygen observations issued from sensors with different detection limits (Konovalov and Murray, 2001).To evaluate how a 10 µM underestimation by Argo profilers would affect the main conclusions, oxygen penetration depths an density levels for Argo were also computed using a threshold of 10 µM.
where ρ is the density and c the heat capacity and T CIL = 8.35 • C (Stanev et al., 2013).

DIVA analysis
Climatologies for the whole period and interannual trends were identified for the three oxygen diagnostics by applying the DIVA detrending algorithm on the ship-based data set (see details in Appendix A).
In short, the DIVA interpolation software (http://modb.oce.ulg.ac.be/mediawiki/index.php/DIVA Troupin et al. (2012)) computes a gridded climatology obtained by minimizing a cost function which penalize gradients and misfits with observations.The DIVA detrending algorithm (Capet et al., 2014) computes trends for each year, i.e. the average difference between data pertaining to this year and the spatial analysis at these data locations.This procedure allows one to account for the sampling error associated with spatial/temporal variability.

Spatial variability
The spatial distribution of the oxygen penetration depth (Fig. 4a) reflects the general curvature of the Black Sea vertical structure.A range of approximately 70 m was observed between oxygen penetration depth in the periphery (150 m) and in the central part (80 m).
A significant spatial variability remains when expressing oxygen penetration in terms of potential density anomaly, σ θ , (Fig. 4b).The ranges of spatial variability derived from theses spatial analysis agreed with those depicted by the Argo profilers (Fig. 5), bearing in mind the different time scales under consideration.

Temporal variability
Between 1955-2005, the oxygen penetration depth rose by an average rate of 7.9 m per decade (Fig. The temporal signals departed from these linear trends between 1988 and 1996, during which deeper oxygen penetration (both in terms of depth and density) and higher oxygen content were observed.

Oxygen inventory and CIL cold content
Positive relationships between oxygen inventory and CIL cold content were obtained for all periods (Fig. 7).Considering a given level of CIL cold content, the corresponding oxygen inventory decreased significantly from period 1955-1975to period 1986-1998 (Fig. 7b (Fig. 7b).
The relationship between oxygen inventory and CIL cold content for the period 1999-2015 does not differ significantly from that obtained for the period 1986-1998 (Fig. 7b).This comparison should be considered with caution, however, as oxygen profiles for the period 1999-2015 originate mainly from Argo floats whose sampling rate is much higher than ship-based casts.
High CIL cold content are much more frequent during the period 1986-1998, while low CIL cold content are more frequent during 1999-2013.

Discussion
The spatial analysis of oxygen penetration depth showed that the use of density coordinates does not eliminate the sampling error associated with uneven spatial coverage (Fig. 4).Deeper pycnal oxygen penetration in the Bosporus area were expected, in relation with the intermediate lateral injections associated with the Bosporus plume.In addition, deeper pycnal oxygen penetration in the southern and eastern periphery suggests the occurrence of diapycnal ventilation along the steep bathymetry (Zatsepin et al., 2007).The aggregation of the most recent ship-based profiles in the Bosporus area and in the southeastern region (Fig. 3), might have led to an overestimation of the basin-average oxygen penetration depth in the last decade, hence to an underestimation of the shoaling trend of the Black Sea oxic layer.
Considering spatial variability revealed a clear shoaling trend for oxygen penetration depth.This shoaling can be seen on both depth and density scales (Fig. 6a,b).This confirms the hypothesis that the shoaling of oxygen penetration depth is not due to a general shoaling of the main halocline, but is associated with a shifted biogeochemical balance in the oxygen budget (Codispoti et al., 1991;Konovalov and Murray, 2001;Tugrul et al., 2014).
Using σ θ coordinates depicts clearer temporal variations (Fig. 5 and Fig. 6).The shoaling rate varies in time and was more intense during 1970-1985 and from 1996 onwards.Argo diagnostics using different oxygen threshold show a larger discrepancy in the case of pycnal coordinates.The co-occurrence of density and oxygen gradients (Fig. ??) results in a higher sensitivity to the sensor accuracy for the σ θ diagnostic for oxygen penetration.However, even a systematic underestimation by 10 µM of oxygen concentration by Argo profilers does not invalidate our results.
The positive correlations between CIL cold content and oxygen inventory observed for all the periods illustrate the ventilation of intermediate layers by CIL formation and advection (Fig. 7b).
In the early 90s, the transient recovery of the three oxygenation diagnostics (Figs. 6a,b,c,7a) provided arguments supporting the stability of the oxic interface (Tugrul et al., 1992;Buesseler et al., 1994).This stabilization matched the convenient perception of a general recovery of the Black Sea ecosystem after the reduction of nutrient load around 1990 (Kroiss et al., 2006).However, Fig. 7 indicates that the oxygenation diagnostics obtained for the period 1986-1998 were associated with much higher ventilation rates (ie, higher CIL cold content) than during the previous periods.If, in response to nutrient reduction, the biogeochemical oxygen consumption terms had been lower during the period 1986-1998 than previously, the increased ventilation during that period would have resulted in higher oxygen inventories.Instead, oxygen inventories observed during 1986-1998 are lower than those observed in the previous decade for similar levels of CIL cold content.We conclude that high CIL formation rates during this period (Piotukh et al., 2011;Capet et al., 2014) provided enough ventilation to mask ongoing high oxygen consumption.
The fact that the relationship between oxygen inventories and CIL content for the last period 1999-2015 is similar to that of 1986-1998 indicates a stabilization in the biogeochemical oxygen consumption terms.Higher air temperature in this last period (Oguz and Cokacar, 2003;Oguz et al., 2006;Pakhomova et al., 2014), by limiting winter convective ventilation events (Capet et al., 2014) led to the lower oxygen inventories ever recorded for the Black Sea (Fig. 6c).
Fore-casted global warming, without excluding transient high ventilation periods, will limit CIL water formation (Capet et al., 2014) and reduce the oxygenation of the Black Sea intermediate layers.
At the same time, uncertainties remain regarding the capacity of re-flourishing economies of the lower Danube watershed to recover their productivity in a more sustainable, less polluting form.
Economic development in the Danube Basin could reverse the improving situation of eutrophication if nutrients are not managed properly (Kroiss et al., 2006).Under these conditions, there is no reason to expect that the oxycline shoaling observed over the past 60 years will stabilize.
There are reasons to worry about a rising oxycline in the Black Sea.First, biological activity is distributed vertically on the whole oxygenated layer, as indicated by zooplankton dial migration Ostrovskii and Zatsepin (2011).The reduction of the oxygenated volume described in this study could therefore have impacted on Black Sea living stocks by reducing carrying capacity and increasing predation encounter rates.It would come in timely to estimate now the impact that a further shoaling of the oxic interfaces would bear on the Black Sea resources for the fishing industry.
Second, under present conditions, a massive atmospheric release of hydrogen sulphide caused by a sudden outcropping of anoxic waters remains unlikely, due to the stability of the Black Sea pycnal structure.Such outcropping event of sulphidic waters would have dramatic ecological and economical consequences (Mee, 1992).On the 27th October 2005, an anomalous quasi-tropical cyclone was observed over the western Black Sea that led, in a few days, to the outcropping of waters initially located at 30 m depth (Efimov et al., 2008).Two years earlier, sulphide was measured in the same area (western central gyre) around 80 m (Glazer et al., 2006).Because global warming is expected to increase the occurrence of extreme meteorological events (Beniston et al., 2007), every meter of oxycline shoaling would bring the Black Sea chemocline excursion events closer to the realm of possibility.

Conclusions
The present study evidenced the decline of the Black Sea oxygen inventory during the second half of the XXth century and first decade of XXIth and highlighted the threat that further atmospheric warming casts upon the vertical stability of the Black Sea oxygenated layer.
Further works are urgently required to assess how actual nutrient emission policies adequately prevents, in the context of fore-casted warming, the ecological and economical damages that would arise from a further shoaling of the oxic interface.
Spatially-resolved biogeochemical models are needed to integrate explicitly the interacting processes affecting the Black Sea oxycline.here.
It is also essential (1) to determine to which extent the shoaling of the oxygen penetration depth entrains a shoaling of the the sulphidic onset depth; (2) to set up a continuous monitoring of the Black Sea oxygen inventory and the intensity of winter convective ventilation (through CIL cold content); and (3) to clarify and quantify the interplays of diapycnal and isopycnal ventilation mechanisms and, in particular, the role played by the peripheral permanent/semi-permanent mesoscale structures and how this relates to the intensity of the Rim Current (Stanev et al., 2014;Kubryakov and Stanichny, 2015).We propose that these objectives might be answered by maintaining in the Black Sea a minimum population of both moored and drifting autonomous profilers equipped with oxygen and sulphidic sensors.
Appendix A: The DIVA detrending algorithm DIVA (Data-Interpolating Variational Analysis) is a method for spatial interpolation.Its principle is to construct an analyzed field ϕ that satisfies a set of constraints expressed in the form of a cost function over a spatial domain Ω.The cost function is made up of (1) an observation constraint, which penalizes the misfit between data and analysis, and (2) a smoothness constraint, which penalizes the irregularity of the analyzed field (gradients, Laplacian etc).
Let us assume that we work with data anomalies, i.e. a reference (or background) field is subtracted from the data points prior the analysis.For N data anomalies d i at locations (x i , y i ), the cost function reads, in Cartesian coordinates: where µ i , α 0 and α 1 are coefficients related to characteristics of the dataset.∇ is the horizontal gradient operator and ∇∇ϕ : ∇∇ϕ = i j (∂ 2 ϕ/∂x i ∂x j ) (∂ 2 ϕ/∂x i ∂x j ), the generalization of the scalar product of two vectors.
The first term of (A1) measures the spatial variability (curvature, gradient and value) of the analyzed field and is identified as the smoothness constraint.The second term is a weighted sum of data-analysis misfits and is identified as the observation constraint: it tends to pull the analyzed field towards the observations.The coefficients of (A1) can be determined from: (1) the relative weights w i attributed to each observation d i , (2) the correlation length L and (3) the signal-to-noise ratio λ (Troupin et al., 2012).The analyses presented in this study were achieved with equal weights w i = 1, L = 0.8 • and λ = 0.5.The minimization of A1 is solved over Ω with a finite-element technique (Brasseur et al., 1996) which excludes data influence across land points Troupin et al. (2010).
The detrending algorithm, presented in (Capet et al., 2014) with synthetic and real case studies, proceeds as follows.Input data can be classified amongst the different classes C j (e.g.1990, 1991, . . . ) of a given group C (e.g. the year).The observation constraint of the functional (A1) can then be rewritten by including an unknown trend value for each class (d C1 , d C2 , . ..): If the function ϕ(x, y) were known, minimization with respect to each of the unknowns d Cj would yield and similarly for the other classes: the trend for each class is the weighted misfit of the class with respect to the overall analysis.
Using an analysis without detrending as a first guess for ϕ, trends are computed for each classes in each group and subtracted from the original data.Following this, a new analysis is performed, the trends are recalculated, and the iterations continue until a specified convergence criterion is fulfilled.
The procedure can be generalized with several groups of classes.The present study considered years and months.
The DIVA software and up-to-date related informations can be found on http://modb.oce.ulg.ac.be/mediawiki/index.php/DIVA.shift in the oxygen budgets (higher consumption) that are independent of the intensity of CIL formation.
While the central part bears typical values of 15.75 kg m −3 , a deeper anomaly (in terms of density) can be seen in the area of the Bosporus plume (16.1 kg m −3 ), which then decreases along the southern (15.85-15.9kg m −3 ) and eastern periphery (15.85 kg m −3 ).These result in a range of spatial variability of 0.35 kg m −3 .The spatial distribution of the oxygen inventory (Fig. 4c) follows that of the oxygen penetration depth.The range of spatial variability reaches 12 mol O m −2 , ie. between 17 mol O m −2 in the central part and 29 mol O m −2 in the periphery.
6a).The basin average was of 140 m in 1955 (ship-based), 100 m in 2005 (ship-based) and 90 m in 2015 (Argo).Considering a systematic underestimation by 10 µM in the Argo data would result in an oxygen penetration depth around 95 m for 2015 (Argo).This shoaling was also observed on the potential density scale (-0.074 kg m −3 per decade, Fig. 6b).The basin average was of 16.05 kg m −3 in 1955 (ship-based), 15.6 kg m −3 in 2005 (ship-based) and around 15.3 kg m −3 in 2015 (Argo).Considering a systematic underestimation by 10 µM in the Argo data would result in values of 15.5 m for 2015 (Argo).The oxygen inventory, integrated from the surface down to the suboxic upper interface, decreased by 44% during the last 60 years (Fig. 6c), considering the ship-based estimate for 1955 (27 mol O m −2 ) and the Argo estimate for 2015 (15 mol O m −2 ).The few ship-based profiles available after the mid 90s revealed the lowest oxygen inventories recorded during the time frame covered by the present study.

Figure 1 .
Figure 1.Typical profiles of temperature, salinity, Brunt-Väisälä frequency (N ), potential density anomaly (σ θ ) and oxygen concentration in the central Black Sea (May).Note the two peaks in the vertical stratification: the thermocline, which is seasonal and corresponds roughly to the upper limit of the Cold Intermediate layer and the halocline, which is permanent, and correspond roughly to the lower limit of the Cold Intermediate layer and the upper boundary of the suboxic zone.Red dotted lines and shaded areas illustrate the diagnostic values derived from each profiles (Sect.2.2).

Figure 2 .
Figure 2. Temporal distribution of the ship-based oxygen profiles merged from the World Ocean Database, R/V Knorr 2003 and R/V Endeavor 2005 campaigns.Only the profiles containing at least 5 observation depths, one observation above 30 m depth and one record with [O2] < 20 µM were considered.

Figure 3 .
Figure 3. Distribution of the ship-based oxygen profiles (Fig. 2) available for each period (black dots).The last panel displays the trajectories of the ARGO floats.Number of profiles for each period are given in the text.Map data: ©Google 2015.

Figure 4 .Figure 5 .
Figure 4. Annual climatologies of (a) oxygen penetration depth (where [O2] = 20 µM), (b) potential anomaly at oxygen penetration depth and (c) oxygen inventory.These spatial climatologies were constructed from the ship-based dataset (1955-2005), accounting for the temporal variability of these diagnostics and the uneven distribution of data (see Sect. 2.3).16

Figure 6 .
Figure 6.Trends of (a) oxygen penetration depth, (b) oxygen penetration density level (σ θ ) and (c) oxygen inventory deduced from (dots) DIVA analysis of ship-based casts and (blue) ARGO floats.On (a) and (b), the diagnostics from ARGO are also shown for the lower threshold of 10 µM to acknowledge a potential bias between Winkler and Argo data.Red lines: the linear trends assessed from the ship-based data set are -7.9 m decades −1 , -0.074 kg m −3 decades −1 and -1.44 mol O m −2 decades −1 for (a), (b) and (c), respectively.Error bars on DIVA estimated trends indicate the standard error associated with the estimation of the mean misfit for each year (see Appendix A).

Figure 7 .
Figure 7. Impact of convective ventilation on oxygen inventory.Frequency distributions of (a) oxygen inventory and (c) Cold Intermediate Layer (CIL) cold content diagnosed from ship-based and Argo profiles for different periods (color legend).(b) Loess regressions (second degree polynomials, span=0.75,Cleveland et al. (1992)) between oxygen inventory and CIL cold content for the different periods (confidence interval α = 0.99).The positive relationships observed during each period illustrate the ventilating action of CIL formation as a source of oxygen to the intermediate levels.The shift of these relationships towards lower oxygen inventories indicates