Modeling long-term changes of the Black Sea ecosystem characteristics

Abstract. A three dimensional coupled physical-biological model is provided for the Black Sea to investigate its long-term changes under the synergistic impacts of eutrophication, climatic changes and population outbreak of the gelatinous invader Mnemiopsis leidyi. The model circulation field is simulated using the high frequency ERA40 atmospheric forcing as well as assimilation of the available hydrographic and altimeter sea level anomaly data for the 30 yr period of 1971–2001. The circulation dynamics are shown to resolve well the different temporal and spatial scales from mesoscale to sub-basin scale and from seasonal peaks to decadal scale trend-like changes. The biogeochemical model includes the main vertical biological and chemical interactions and processes up to the anoxic interface zone. Its food web structure is represented by two phytoplankton and zooplankton size groups, bacterioplankton, gelatinous carnivores Mnemiopsis and Aurelia, opportunistic species Noctiluca scientillans. The nitrogen cycling is accommodated by the particulate and dissolved organic nitrogen compartments and the dissolved inorganic nitrogen in the forms of ammonium, nitrite and nitrate. The ecosystem model is able to simulate successfully main observed features and trends of the intense eutrophication phase (from the early 1970s to the early 1990s), but points to its modification to simulate better the ecosystem conditions of the post-eutrophication phase.


General remarks:
We are very grateful to editor and referees for their important and useful comments.We have intention to follow major part of their advices.At the same time we would like to point out that the MS is devoted to modeling long-term evolution of the Black Sea ecosystem.Therefore the material presented in the MS about the physical processes reproduced in model simulations is selected to convince the reader that the subsequent reproduction of the ecosystem evolution is based on close to the reality and reasonable dynamics of the Black Sea.We tried to consider the referees and editor comments from this point of view and some part of them to use in our future work.We hope that such approach is acceptable and make possible to avoid consideration in the text of MS C860 insufficient details.
There is an additional fortune for such approach related with publication in English the translation of our paper "Seasonal and Interannual Variability of Black Sea Hydrophys-ical Fields Reconstructed from 1971-1993 Reanalysis Data" by V. V. Knysh, G. K. Ko-rotaev, V. A. Moiseenko, A. I. Kubryakov, V. N. Belokopytov, and N. V. Inyushina in the central peer-reviewed Russian journal Izvestiya, Atmospheric and Oceanic Physics.That paper contains detailed description the preliminary procedures of data prepara-tion, all steps of data assimilation and broader analysis of the reanalysis results than we are able to present in our MS.We have included the reference to the edited version of MS on this new appeared publication: "Seasonal and Interannual Variability of Black Sea Hydrophysical Fields Reconstructed from 1971-1993 Reanalysis Data" by V. V. Knysh, G. K. Korotaev, V. A. Moiseenko, A. I. Kubryakov, V. N. Belokopytov, and N. V. Inyushina ISSN 0001_4338, Izvestiya, Atmospheric and Oceanic Physics, 2011, Vol. 47, No. 3, pp. 399-411. © Pleiades Publishing, Ltd., 2011.First referee comments: General comments.'It seems that assimilation of remotely sensed SST and SLA data is crucial for the model performance.For the period 1971-1993 only the occasional in-situ data were assimilated.I have got an impression that during this period the model was capable to represent seasonal/interannual variations but not the mesoscale variability despite the model was eddy-resolving (7x8 km).The model resolves eddies in the 1994-onwards simulations, when high-res satellite data was assimilated.It suggests that the model skill is heavily dependent on data assim-ilation.For this reason it will be beneficial if the authors presented results for 1994-onwards without data assimilation for comparison.'We insert the next text after line 23: "It should be noted that the model also resolved the mesoscale variability for all period without data assimilation.Two examples of calculation results for 1994-onwards period without data assimilation are presented on Fig. 1, where typical mesoscale eddies structures of the Black Sea are seen very well.And altimetry data assimilation made possible to simulate sea dynamics closer to the C861 reality." Specific comments.'Abstract, last phrase.'The ecosystem model is able to simulate successfully main observed features and trends of the intense eutrophication phase (from the early 1970s to the early 1990s), but points to its modification to simulate better the ecosystem conditions of the post-eutrophication phase.'It should be stated more clearly what the model can and what cannot simulate.The last statement probably means:' the current model is not capable of reproducing the recovery stage after 1990s'.Please clarify.' We meant that the ecological model was calibrated with the data collected during time period from the early 1970s to the early 1990s.Then amount of data for the deep part of the Black Sea was negligible and from late 1990s satellite measurements became the main sources of the data.
1.Introduction.Page 2041, line 16-18 Check the grammar: 'The period 1971-1994 involved relatively dense hydrographic surveys and was replaced by the availability. ..' What was replaced?-the period?
We changed this text: The temperature and salinity data were assimilated for the 1971-1993 because of relatively dense hydrographic surveys took place during this period.
Then the number of surveys was extremely small, but satellite remote sensing data became available.

C862
We changed lines 19-21: "Reanalysis fields of temperature, salinity, and currents of the Black Sea made it possible to reveal different tendencies in their interannual variability which coincide with the variability of the external parameters of atmospheric forcing (Knysh et al., 2011).We would like to change the reference of "Moiseenko et al ( 2009)" to the reference on new, just published "Knysh et al., (2011)", which was mentioned above.The former devoted to the preliminary results of reanalysis and the latest one contains more details.Moreover it is translated in English.The reanalysis presented in the MS is quite the same as in Knysh et al., (2011).That is why we have included in the MS only short description of the procedure of reanalysis and analysis of simulations just to convince the reader that the subsequent reproduction of the ecosystem evolution is based on close to the reality and reasonable dynamics of the Black Sea.[V.V. Knysh, G. K. Korotaev, V. A. Moiseenko, A. I. Kubryakov, V. N. Belokopytov, and N. V. Inyushina.ISSN 0001_4338, Izvestiya, Atmospheric and Oceanic Physics, 2011, Vol. 47, No. 3, pp. 399-411. © Pleiades Publishing, Ltd., 2011.] 2. Methodology of simulations Page 2042 last para ' ERA-40 . . .for the period 1958-2002 was used' .This section is about the period 1971-1993, why did the authors 'use' the data for 1958-1970 and 1994-2002?More careful editing is required.
We have changed: "ERA-40 .Based on the questions of referees, it seems that the issues of preparation of data for reanalysis and the reanalysis itself do not presented clear enough in the MS.We have introduced an additional subsection to section 2, which, in our opinion, will make clearer the distinction between assimilation of climatic seasonal arrays and real observations.15 years is time to reach periodic seasonal cycle when climatic data are assimilated.Simulated climate has a periodic seasonal cycle, so, we stop simulations when periodic oscillations with annual period are achieved.We intent making this item more evident in the text of the MS.Thus, section 2 will be like this (English editing will be done later on): 2. Methodology of simulations 2.1.Preparation of data for assimilation.
The procedure for preparation of the data included the construction of climatic annual cycle of temperature, salinity and current velocity through assimilation of monthly climate arrays of temperature and salinity in the upper 300-meters layer of the Black Sea based on observations data.Assimilation of climatic fields of temperature salinity performed cyclically for 15 years.By the end of the 15-year periodic regime of sea dynamics was resulted, which was identified with the climatic annual cycle.Thus constructed climatic annual cycle was used as a base for the construction of monthly temperature and salinity for the 1971-1993 period by the optimal interpolation (Knysh et al., 2008).Autocorrelation functions of temperature and salinity climatic fields were calculated for different directions every 100 for winter (February), spring (May), summer (August) and autumn (November) on the horizons 10, 50, 105, 200, 500, 1000 and 1500m using the simulated climatic array (Moiseenko and Belokopytov, 2008).
The spatial correlation of temperature and salinity fields is well approximated by an ellipse on all horizons and for all seasons.On all horizons and in all seasons (except C864 for the horizon 50m in winter), the orientation of temperature field isocorrelates is close to zonal.The large ellipse semi-axis has zonal direction and is equal ∼ 330km, the small semi-axis is equal ∼ 160km.Two-dimensional correlation functions of salinity field are stretched in 25 zonal direction and large ellipse semi-axis is equal ∼ 260km, the small semi-axis is equal ∼ 75km.Two-dimensional correlation functions were approximated by analytical expressions.We used the space-time optimal interpolation with a strobe pulse ±45 days to build the monthly arrays of temperature and salinity by means of optimal interpolation.The contribution of measurements with a distance of ∆t from a given month was determined with the help of a time correlation function of the exponential form exp(-α|∆t|).This method for deriving monthly measurement data of temperature and salinity allows each individual measurement to be repeatedly engaged in future data assimilation.This also makes it possible to fill the gaps when measurements were not available.In this way we manage to significantly reduce the negative impact of the inhomogeneity of the spatial and temporal distribution of measurement data on the calculation results.The monthly fields of temperature and salinity measurements were reconstructed at 36 levels: 0, 5, 10, 16, . .., 52, 60, 70, . . ., 100, 110, 120, 135, 150, 175, 200, 250, 300, 400, . . ., 600, 750, . . ., 1150, 1350, 1550, . . ., and 2150.Climatic temperature and salinity data are used deeper 300 m.It is important to note that the optimal interpolation method permits to evaluate a standard deviation of every interpolated value of temperature and salinity.Monthly arrays of temperature and salinity as well as monthly fields of dispersion errors are interpolated for each time step model caculation.Then they are used for correction of temperature and salinity fields simulated by POM using the Kalman filter algorithm formalism (Gandin and Kagan, 1976).Such approach when the error statistics is evaluated a priori is usually named as the data assimilation by optimal interpolation.Let us mention that the weight of the POM simulation correction depends on the standard deviation of optimally interpolated temperature and salinity fields/ 2.2.Reanalysis of marine dynamics of 1971-1993 Reanalysis of the Black Sea dynamics for 1971-1993 was performed by assimilating the temperature and salinity profiles C865 into the circulation model (Moiseenko 5 et al., 2009) based on the POM (Princeton Ocean Model) model (Blumberg and Mellor, 1987).The POM regional model includes the Mellor and Yamada level 2.5 turbulence module (Mellor and Yamada, 1982).The river discharges to the sea are taken into account in the model, as well as the water exchange with the Azov Sea through the Kerch strait and with the Marmara Sea through the Bosporus, where water flows from the Black Sea in the upper layer (upper Bosporus current) and into the Black Sea in the lower layer (lower Bosporus current).Considering the lack of knowledge on the interannual variability of the lower Bosporus flow, it is estimated by assuming that the Black Sea water volume does not changed from one year to another, i.e. annual rivers inflows, water discharges through the Bosporus and the Kerch strait and also precipitations and evaporation altogether are equal to zero.The simulations were carried out on a horizontal rectangular grid of 8.1 km along the zonal and 6.95 km along the meridian directions.Climatic temperature and salinity arrays are based on the all available in-situ measurements.We will add appropriate remark in the text.Unfortunately the Med Atlas data for the Black Sea are very poor quality.We assimilated all in situ data measurements, and comparison model data with in situ data, which are assimilated in model, is not quite correct.Instead of comparison with data, which are assimilated in model, we use indirect method of estimation of quality of model calculations by comparing variability of sea temperature in the upper layers with variability of heat flux from atmosphere, and variability of kinetic energy of currents with variability of wind stress vorticity (Knysh et al.,2011).We will correct on "modeling" through all text.

Second referee comments
The manuscript presents the results of a 3D numerical model applied to the Black sea for investigating tis long term dynamics.The model assimilates various data sets (temperature, salinity, currents, sea surface height) and the period of simulation has been chosen when data are available.Since the model has the aim to be used for operational and prediction purposes, this is important to know what are the performances of the model without data assimilation.I recommend that the authors assess model errors (e.g.bias, RMS, correlation) when run without assimilation.Indeed, data assimilation has to be used when knowing model deficiencies.It seems in this case that the model performances strongly rely on data assimilation and it is important to know what will be the performances of the model when no such a large amount of data will be available (which will be the case when run in operational).In the simulations presented, we are told that the model assimilates data each day (interpolated data, but the model is corrected daily leaving very little freedom to the model dynamics).Since all the available data are assimilated, the model can not be validated with independent C868 data set and this is not optimal.This prevents testing the adequacy of the assimilation protocol.The manuscript does not investigate which type of data need to be assimilated to improve at the best the results, as well as the spatial and time frequency that is necessary to obtain good performances.This testing is very important because it could give recommendation on the types and frequency of data that are needed in the Black Sea for operational purposes.This manuscript use all the available data sets collected during 1971-1993 and this prevents to validate the approach.Besides, this is sure that in an operational model all these data sets will not be available and this is why I strongly recommend to 1) test the model without assimilation over the period 1971-2001 , 2) to test the assimilation scheme/protocol, 3) to derive recommendations of the frequency/type of data that are needed in operational model.
We will add the next text into the MS concerning to notes: It is known that the sea state, calculated with the help of a mathematical circulation model, describes an inherent trajectory in the phase space that can be quite different from the actual trajectory.
In carrying out lengthy calculations, this difference can reach significant values.That's why circulation models applied to operational forecasting of the real state of the sea, used some procedure of assimilation of observational data, which is just like "attracts" the trajectory of their own models to the real trajectory of the system.It is very helpful to know what kind of contribution make assimilated observations to the dynamics of the system, i.e., how much of the information supplied the data, and how much the model.The closer the model trajectory in phase space to the real trajectory, the smaller the amount of data can be assimilated, and the more successful will be the dynamic interpolation.The quality of the model and the role of used data assimilation procedure can be estimated by comparing the reanalysis data with the data of prognostic calculation without assimilation of any data.We have carried out such comparison for the ten years period of 1985-1994.On the Fig. 2, 3 some results of this comparison are presented.Fig. 2 shows the evolution of the average area of salinity at different depths during the period.Approximately to a depth of 200 m salinity values in the prognostic calculation is higher than in the reanalysis, in the first case clearly shows a positive C869 trend, which is virtually absent in the reanalysis.Salinity with depth at the prognosis is close to the values of the reanalysis, and deeper than 300 meters prognostic salinity is less than the reanalysis one, and its trend is negative.The seasonal signal is clearly visible to a depth of 400 m, and its amplitude at the reanalysis calculation is higher than the prognostic calculation.Deeper than 400 m the seasonal signal at the reanalysis is practically absent, but at the prognostic run this signal can be seen, albeit minor, to a depth of 750 m.The difference between the average area salinity at the prognosis and reanalysis in the upper 100 m reaches a value of approximately 0.4 psu; in the layer of 100 -175 m of about 0.1 psu; in the 175-250 m layer is about 0 psu to the end of the calculation period.Deeper the difference becomes negative: in a layer of 250 -400 m it is around -0.05 psu, and deeper 400 m it reaches values of -0.1 psu.Fig. 3 shows the average for the entire decade vertical profiles of average area temperature, salinity and the kinetic energy of currents.In the upper 25 m layer temperature (Fig. 3a) at the prognostic and reanalysis calculations is the same.In the layer from 25 to 40 meters the temperature is slightly higher at the prognostic run.At the prognostic calculation the Cold Intermediate Layer (CIL) is thinner, the temperature in its core is smaller and the thermocline is sharper than at the reanalysis.Deeper 100 m down to the bottom, the temperature is less at the prognosis.Salinity (Fig. 3b) in the layer from the surface to 150 meters is greater at the prognostic calculation than the reanalysis, and deeper -less.So, the halocline (and, therefore, main pycnocline) is slightly sharper at the reanalysis.The values of kinetic energy at the reanalysis, on the contrary, up to 150 meters great and less in the deeper layers (Fig. 3c).The performed comparison shows that, indeed, the phase trajectories thermohydrodynamic state of the Black Sea in the space of states differs at the prognostic calculation and the calculation with the data assimilation upon at the same external forcing.However, these differences does not lead to catastrophic deviations even during a quite long time calculations.This allows us to hope that the model can be used quite well to carry out hydrodynamic interpolation of observational data in order to restore the true state of the Black Sea.
This is not clear whether data assimilation is also used in the ecosystem model.C870 The ecosystem model does not use data assimilation.We will emphasize this in the text of the manuscript.
The long period of simulation performed with the ecosystem model is really challenging and I would like that the authors comment on the balance of nitrogen in the Black Sea (how the sources and sinks of nitrogen balance after 28 years of integration).This is surprising since the authors do not represent benthic processes which are important for the balance of nitrogen in the Black Sea (see for instance Friedl et al , 1998 andFriedrich et al ., 2002 who show than benthic denitrification was a critical sink of nitrate consuming the huge stock of nitrate brought by the Danube).
The model does not really represent benthic processes.They are important mainly on the north-western shelf of the Black Sea.The main sources of the nitrogen in the model are river discharges and ammonium flux on the lower boundary of the model area in the deep part of the basin.The last parameter varied during the period under review to provide nitrate concentration in the layer of maximum corresponding to that which is known from measurements.These sources of nitrogen are compensated by sinking of nitrogen with diatoms and particulate organic matter.We will include some comments in the text.
General remarks: . . .most of the figures show volume integrated values (volume salinity, phytoplankton, nitrate,..), in the revised version, it would be very relevant to show the long term evolution (1971-1992 period of simulations) of these quantities for the physical model (salinity, temperature, mixing layer) but also for the ecosystem model such as phytoplankton, nitrate, zooplankton and gelatinous, sulfide, oxygen.In all the figures it would be helpful that the authors clarify whether they show model results or interpolated data.At some places in the text, they refer to collected data.
All figures concerning ecological model in the article present results of modeling.We can add diagrams with evolution of the ecosystem quantities as it is suggested above.They look like in figure 4.

C871
Other Comments: Introduction : please add some references Please provide in the appendix a description of the ecosystem model.For instance: Please specify how hydrogen sulfide is represented in the model (In Oguz work, Iron and Manganese were represented which is not the case here).How sulfide is generated?Please give some details about the representation of gelatinous (are they explicitly simulated or are they forced from the data?)Bottom boundary conditions of the northwestern shelf where sediments are important It is specified that the model only extends to 200m but at page 6 we are told that the model imposes a huge pool of hydrogen sulfide and ammonium, this is contradictory.
We will add in text of the manuscript more expanded description of the ecological model describing these topics.As far as sediments on NWS are concerned, we run a few numerical experiments with resedimentation in the form of nitrate fluxes upward from the bottom.(The model includes sediments of diatoms part of phytoplankton and particulate organic matter).These experiments have shown that this process is important only on the shelf, but it is not as important for the whole basin.The hydrogen sulfide and ammonium zone in the Black Sea starts from the depth of approximately 120m, so on the 200m depth concentrations of hydrogen sulfide and ammonium are high enough.
More important, how the model can reach a steady state conditions?Spin up time?Balance of nitrogen after 28 years of integration?Role of the sediments on the shelf?Do you conclude that sediments degradation of the shelf are not important if one wants to integrate the model over decades?
The spin up time of the model is approximately 5 yerars.To reach a steady state conditions (i.e.normal seasonal sycle) we ran the model with the same hydrodynamical fields, corresponding the first year of our study and then used these biogeochemical fields as initial conditions.Concerning the role of sediments see the text above.
Page 6, section 3: "the reanalysis data collected every five days are used. .."You C872 mean simulated?
We saved simulated data for every 5 days on hard disk and then we used them for analysis of the results.Now we have changed the line 22-24 to: "the reanalysis data stored every five days are used. .." Figure 1 is model results?Or data?I would suggest to make the average over the upper layer where changes are important.This is the model results.We have changed the cutline on: Fig. 1.Interannual and seasonal variability of the volume averaged model temperature.
We have added on line 2, page 2047: "The average over the some upper layers are presented in (Knysh et al.,2011)." Page 6: The authors conclude that "Such characteristics of the thermohaline structure as the autumn-winter cooling and spring-summer heating of surface waters, formation of the upper mixed layer (UML) and the seasonal thermocline, renewal of the cold intermediate layer (CIL) and formation of a new CIL, decrease of cold reserve in CIL by autumn are resolved by model simulations." Please show some comparisons.For instance show the cold content of CIL as defined by Stanev (1995) and compared with data (even if we can expect a good agreement since the model assimilates all available data).Nevertheless this comparison can give information on how the assimilation of temperature data allow to represent the dynamics of CIL which is one of the key processes in the Black Sea dynamics.Same idea with the mixing layer depth, stratification index (potential energy anomaly for instance. ..).Please add this comparison with data.
We will add to MS the reference"(Knysh et al., 2011)" after line 7, page 2047.This paper is translated in English, and extended analysis of CIL is presented there.We think that it is not reasonable to provide here extensive discussion of physical fields and parameters because of the MS is devoted to the ecosystem description (see the C873 general remarks).
Figure 5: Please give a similar plot for the surface salinity.Please explain why the volume average salinity increase from 1971 to 1992 although the Bosphorus inflow is estimated in order to keep the system in balance (as said by the authors above).Why the salinity plot stop at 1992 although simulations where performed until 2009?
Because the model assimilates data, there is a source term in the right hand side of the equation for salinity.So, the system cannot be in balance.We have two periods of reanalysis which differ from each other of the type of assimilated data -hydrographic data and satellite data.Firstly, after 1992 there were not the observations of salinity, and secondly we do not know exact input through Bosporus Strait.Taking into account the general remarks made at the beginning, we have replaced the Russian references on these pages: (Hydrometeorology and hydrochemistry, 1991), Lipchenko et al. (2006), (Repetin et al., 2006) with the reference (Knysh et al., 2011), where detailed analysis and pictures are presented.We are very thankful for these notes.We are going to continue analysis of physical fields.But, the MS is devoted to the ecosystem description (see the general remarks) and we think that enlarging the text by adding detailed description of physical part of work does not appear to be sufficient reasons.
Please compare results shown in Figures 9 and 10 with available observations.Since the model is 3D it would be very relevant to show spatial maps for different years of phytoplankton (ditaoms/flagellates) and zooplankton in order to appraise the spatial structure.to represent these trends identified in Konovalov and Murray( 2001).
The long-term evolution of the nitrate concentration in the layer of maximum (basinaveraged) is presented on figure 12 (blue dots).In the first year of the modeling the nitrate concentration was approximately 3 mmolN/m3, in the late 1980s -early 1990s it reached the value of about 8.5 mmolN/m3 and then decreased to about 4.5 mmolN/m3 in late 1990s -early 2000s in accodance with Konovalov's data mentioned in the text of the manuscript.
Also, I suggest showing the evolution of the vertical nitracline during 1971-1992 as well as for phytoplankton, zooplankton and gelatinous as already mentioned in my general comments.
We are going to add the evolution of the vertical nitracline during 1971-2001 as well as for phytoplankton, zooplankton and gelatinous.These diagrams are presented above.
Vertical plots over 1971-1992 are required.I guess that Figure 12 are data and not model results.Is it correct?Please specify.
Figure 12 presents the model results.We have included this explanation in the text of the manuscript.
Figure 16: Please specify the units.I repeat that in the revised version, the authors may provide similar plots but continuously from 1972 to 2002 in order to see the transition from periods because this is very challenging to have a model that is able to simulate this transition especially without data assimilation for the ecosystem.
The units are mmolN/m3 .We have corrected the text.As we mentioned above we will include additional time diagram to demonstrate time evolution of the simulated parameter in water column.
Page 19, the authors refer to a climatic change (warming tendency).Do you identify a climatic trend with your simulations?

C876
Figure 17a demonstrates evolution of the annual-and winter-mean surface temperature, derived from the reanalysis (i.e.results of the modeling).
Figures 19and 20 Please give some quantitative estimations of the model errors in agreement with the quality required for MyOcean products (for instance is annual average comparison enough for a MyOcean product ?) We will add the revised version with some quantitative estimations of the model errors, i.e rms deviation of surface chlorophyll derived from modeling from satellite data.
Conclusions: They authors may provide some conclusions about the ability of the model to be sued for operational purposes in the Black Sea investigating the type f data that are needed to reach good performances.
Now the ecological model is used in the Black Sea Forecast system.It assimilates satellite chlorophyll data.Other parameters which are needed for good performance are inorganic nitrogen in river estuaries and values of ammonium and hydrogen sulfide on the lower boundary of the model area.Unfortunately these parameters are hardly known.We will add these comments in the text of MS.

M. Bell (Editor)
The two main difficulties with the paper are somewhat linked.1) I believe that figures 1-10 and 12-18 only show results from the model integrations.The comparisons with satellite data in figures 19-21 are useful but not particularly reassuring.It is difficult for the reader to assess the reliability of the results as they are currently presented.
We are going to add this part of the article with a plot showing RMS deviation of the model data from SeaWiFS data.
2) Some of the more important references on the in situ observations are to papers written in Russian.These are relatively inaccessible to western scientists.
We have changed the references in Russian to the reference "(Knysh et al., 2011)" in C877 English the MS.
Section 2.1 Para 4: The Knysh and Moiseenko references are in Russian.It would be helpful to give a better picture of the data from the hydrographic surveys.Some indication of the number of profiles taken per survey and the months surveyed.Perhaps one or two figures showing the distribution of observations.
We have such figures but taking into account the general remarks made at the beginning we insert the reference to (Knysh et al., 2011), where this information is available.We are thankful for your comments which concern to the preparation of data for reanalysis.We will insert the additional subsection "2.1.Preparation of data for assimilation" devoted to this subject (see answer to #1 referee and general remarks).
Section 2.3 Para 1: Most of this paragraph is redundant (covered in 2.1).Para 2: Typo "ne-dimensional" in first line.
We have changed the text of MS.
Last para: River discharges are evidently crucial.Is the Ludwig reference readily accessible ?Last sentence: Is all nitrate upwelling at 200 metres set to zero ?Is nitrate C878 upwelling confined to a shallower layer (e.g. the UML).
The Ludwig's data were accessible for the 'SESAME' project participants.200 metres level is a lower boundary of the model area.Nitrate concentration on it is set to zero.Nitrate upwelling is important in shallower layer (upper nitroclyne).
Section 3.1 Figure 1 is it possible to show a corresponding time-series from your monthly arrays (built solely from observed values) ?Clearly there would be gaps in that time-series so the results would have to be presented as values for individual months.In the winter months when the water cools and convention forms new water of CIL, it is assumed that the upper bound of the CIL coincides with the sea surface.Thank you.We will correct.
Para 6: change "enormously warm" to "very warm" or "exceptionally warm".Last 2 C879 sentences: the 50 metre level is right at the base of a very sharp thermocline.Would it not be better to focus on the volume of water colder than 8oC in summer which does show the variations you intend to describe.
We will change "enormously warm" to "exceptionally warm" and "cold reserve" to "cold storage".
Para 7: The trend in layer 500 -1000 m is really small isn't it ?(0.0002 K over 20 years).I don't think it is worth highlighting.
Yes, the variability of temperature on these levels is very small.We would like to show our model result.
Figure 4: Please label the plots a) to e).The scale on d) needs correction (8.6 for all values at the moment).
We will correct.
Para 8: "The temperature increase on 200 m horizon": presumably means to refer to 100 -300 m depth range and figure d).These changes again appear to be very small and not worth highlighting.
We will demonstrate our model result.It seems important.
Section 3.2 Figure 5: Same comment as figure 1.
We will demonstrate our model result.It seems that it is possible conclude because of evaporation decreased, precipitations (in annual values) grew, the trend of rivers' inflow appeared to be not significant, as noted in MS.No, 0-50m is correct.
Next para: "Thimplis" please check spelling (I think it is "Tsimplis" ? ) Yes, you are right: it is "Tsimplis", we will correct.
Figure 7: This figure shows salinity decreasing over the period at the surface and increasing at depth.This is consistent with the water becoming more stratified over the period and less winter-time mixing.Should your discussion capture this better ?
We are very thankful for these notes.We are going to continue analysis of physical fields.But, the MS is devoted to the ecosystem description (see the general remarks) and we think that enlarging the text by adding detailed description of physical part of work does not appear to be sufficient reasons.
Section 3.3 Does Figure 8 show variations due only to basin-scale circulations i.e. are the contributions from mesoscale variations negligible ?The units of the trends might be clearer as m2s-2 per year (note use "s" not "c") We assimilate smoothed values of temperature and salinity.So, there are not the mesoscale variations.1971-1993 and 1994-2001 affect the simulated sea thermodynamics [http://catalogue.myocean.eu.org/static/resources/myocean/quid/MYO-WP10-QUID-V1-BS-MHI-BS-MHIR_DAVER_PHYS-RAN-v1.0.pdf].As example we demonstrate this on temperature fields.Prepared for the first period, the average monthly temperature fields are smooth.There are no mesoscale structure in these fields (except eddies to right of RIM current) because of the surveys in the most cases do not cover the all entire area of the sea, and the average monthly gridded arrays were got by means of from the optimal interpolation deviations in situ data from climate data.At the second period dynamic altimetry sea level was assimilating in model, There are a lot by synoptic eddies in the structure of sea level field.The relevant structures were induced in the fields of temperature as a result of assimilation of these kind observations.Variability of area average in vertical during two periods is presented on Fig. 8.It can be seen the main processes which form the thermohaline structure of Black Sea waters during both periods.But, for example, the CIL is located deeper and its cold storage is great when altimetry sea level is assimilated (Fig. 8b).This is very clearly seen for the 1993 year.Upper mixed layer (UML) depth is ranging from 20 to 64 m.Increased thickness of the UML (up to ∼ 54 m) was in 1976, 1985, 1987, 1989, 1992, 1993 -1995 and 1998 1993-2002. Fig.1b . Fig.1b shows a tendency to reduce of cold storage of CIL and increasing of temperature in the active layer of the sea for this period." Para 1 3 lines above figure 10: "lateral support".More evidence would be needed to substantiate this assertion.
Phytoplankton biomass is in general higher on the NWS of the Black Sea than in deep part of the basin due to Danube supply with nutrients (see for example figures 20, 21 illustrating surface distribution of the chlorophyll-a).As a result, the phytoplankton biomass in deep western part of the Black Sea is higher than in eastern one due to lateral exchange with NWS.Last 2 sentences above figure 17: This explanation is not very convincing to my mind.
The slopes of the blue and red curves in figure 17b are nearly equal and opposite after 1992 but the red slope is larger than that of the blue before 1992.
It seems that the text in these two sentences is misleading without a reference to figure 18.We will add this reference in the text.Section 5 General question: Are there major sources of nutrients from other rivers than the Danube ?
We used sources of nutrients from 12 rivers, but the Danube river makes the greatest contribution.
We will correct the text.

Aurelia
26 _-surfaces were used along the vertical: 0, −0.003, −0.006, −0.009, −0.012, −0.015, −0.020, −0.025, −0.030, −0.035, −0.040, −0.045, −0.050, −0.055, −0.060, −0.067, −0.075, −0.090, −0.140, −0.200, 20 −0.330, −0.500, −0.670, −0.830, −0.910, −1.000.The coefficients of horizontal turbulent exchange of momentum, heat and salt were assumed to be: AM = 300m−2 s−1, AH = 60m2 s−1, respectively.The coefficients of vertical turbulent viscosity KM and diffusion KH are represented by: KM = SMql ; KH = SHql, where q2 is twice the kinetic energy of turbulence (per unit mass), and l is a master scale of turbulence.SM and SH are non-dimensional characteristics depending on the Richardson number, Ri:Ri = −(Nl/q)2 ,where N is the local buoyancy frequency.For calculation of SM and SH, the empirical formulas are used(Mellor, 2001).[Mellor, G. L., One-dimensional, ocean surface modeling, a problem and a solution.J.Phys.Oceanogr., 31, 790-809, 2001.]C866The products of Mediterranean array of global reanalysis ERA-40 created by the European Center of Mid-term Weather Forecasts ECMWF(Uppala et al., 2005)  with spatial resolution 1.125_×1.125_and temporal resolution 6 h for the period 1971-2002 was used for the atmospheric forcing of the model.They include wind direction and speed at 10m height, temperature on the sea surface, fluxes of short-wave and long-wave radiation, and sensible and latent heats, precipitations (continuous and convective ones) and evaporation.The wind stress vector is defined using the bulk formula.P2042 line 22-24.'The coefficients of turbulent exchange of momentum...' ->' The coefficients of horizontal turbulent exchange of momentum...' 'The coefficients of vertical turbulent viscosity KM and diffusion KH are expressed through the turbulence kinetic energy and the stability parameters which are the functions of the Richardson number.'As the parameterization of the vertical diffusion/viscosity is critical to the accuracy of simulations, the equation and parameters used or a reference should be given.We will add corresponding text -see the answer above.Page 2044, line 20-22 '. ..the climatic sea surface topography obtained from the assimilation of the climatic temperature and salinity arrays in POM. ..'. Are the 'climatic temperature and salinity arrays' based on observations or these are the simulated climatologies as in(Moiseenko and Belokopytov, 2008)?Please clarify.
Page 2046, line 21- 25  This section discusses intricate variations sometimes as small as 0.05 degrees of temperature.It would greatly benefit if the uncertainties of simulations are estimated ( e.g. by error bars on the graphs) and validation against observations is provided.This in contrast to the biochemical component, which is validated against SeaWiFS data.All information in this section is based on the model output.At least some comparisons with in situ data measurements are required to assess the accuracy of C867 the model results.Temperature and salinity data for the Black Sea over the period 1971-2002 are available from the Med Atlas and other databases.
Minor issues: Page 2070, Fig 11.Harmonize the spelling (UK/US) :'Modelling' in figure 11 but 'modeling' in the caption and further in the text.

Figure 7 :
Figure 7: The variability is much lower than in Figure 5. Did you average the salinity field?Fig.7 presents the annual mean values.We have changed on the next: "Figure 7. Interannual variability of the annual mean salinity averaged in the layer -20 m (a), 20 -150 m (b), 150 -300 m(c), 500 -1000 m (d)."Pages 10 and 11: the authors analyze the trends present in the salinity fields shown in Figures5 and 6.They enumerate a lot of possible factors that can explain these trends but do not identify which is/are the factor(s) actually responsible for these trends.Since the authors have a model, they can use the model in order to select the cause.This is the advantage of using a model compare to only data analysis.Moreover, this is surprising that the authors refer toRepetin et al (2006), Lipchenko et al (2006)  although they are not using these data sets (they are using ERA-40 and Ludwig river data).Also, the authors may analyze the trends in the river data they use and atmospheric precipitation they use to force the model in order to investigate the existence of trends in these forcings to explain the long time trends they simulated.Section 3.3 (pages 12-13): The authors analyze the Black Sea circulation at different depth.For clarity it would be very helpful to have plots of this circulation at different depth and for different C874

Figures 9
Figures 9 and 10 demonstrate evolution of area-averaged parameters.Unfortunately, we do not have sufficient in-situ date to compare with results of modeling.Some comparison of the modeling results with data is presented in the section "Calibration of the ecological model" in the paper 'G.K. Korotaev, T. Oguz, V.L. Dorofeev, S.G.Demyshev, A.I. Kubryakov, Yu.B. Ratner, Development of Black Sea nowcasting and forecasting system, Ocean Sci., 7, 1-21, 2011".(We will add the reference to this article in the text).We will also present in revised version of the manuscript the maps of space distribution of phytoplankton and zooplankton for different years, corresponding early eutrophication, strong eutrophication and post eutrophication phases of the Black Sea ecosystem evolution.Page 15, 1st paragraph: The authors enumerate the conclusions of theKonovalov and  Murray (2001)  paper on the long term evolution of the nitrate content.Since the authors have realized long term simulations, it is important to check whether the mdoel is able C875

Para 4 :
The description of the optimal interpolation is concise.Could you add a brief description of the following details: (a) How are the covariances specified for the calculation of the monthly climate arrays ?(b) You describe how climate values are calculated for each day.How are these values assimilated (e.g.what value do you use for the nudging coefficients).Para 5: Could you briefly describe how the variances are specified.Do you interpolate in the vertical ?Are the temperatures and salinities interpolated independently ?Para 6: By a strobe pulse of +/-45 days do you mean a 90-day time window centred on the middle of the month ?Could you describe in more detail how the fields are assimilated (what weight is given to the observed field relative to the model background).Section 2.2: I presume that Korotaev et al 2003 describes the assimilation for this period in detail.If so please refer to it in the last 2 sentences.

Fig. 1
Fig.1 presents the model result.Monthly arrays based on observations were prepared for assimilation only in upper 300m layer.Deeper we used climatic temperature and salinity data (see added subsection 2.1).So, it is possible to construct the pictures only for upper 300m layer.We have changed to: Fig. 1.Interannual and seasonal variability of the volume averaged model temperature.Para 3: Are the values for UML quoted basin averages ?If not what are they ?Yes, these are basin averaged.Para 4: In winter Figs 2a and 2b appear to show that the temperature increases monotonically with depth in January.The definition of CIL given then only applies for months when there is an intermediate temperature minimum.Because of my lack of familiarity with the Black Sea I was initially rather confused by this paragraph.

Para 5 :
Fig 3.2 should say Fig 2b.Near the end of the para the ref to "(Fig 3)" could helpfully say "(Fig 3a and 3b)".

Figure 6 :
Figure 6: Why are the years 1976 and 1984 chosen here when figures 2 and 3 used 1993 and 1981.Typo in caption "1076".A consistent choice would increase the value of the figures.We will correct the text.The noted figures correspond to different weather conditions: Fig.2,3 -to normal, cold, warm and very warm conditions; Fig.6 -to cold and warm winters.C880

Figure 11 :
Figure 11: Why are the model scales (on lhs of figure) 1000 times smaller than the measurements ?Does the figure show data from the Purcell model or the integrations described in this paper ?

Figure 16 :
Figure 16: Please define the "central part of the Black Sea" more precisely.Diagrams presented on this figure correspond to the point in the central part of the Black Sea with coordinates: 34E, 43N.

Figures 20 and 21 .
Figures 20 and 21.Does the chlorophyll edge in the SeaWIFS data follow bathymetric contours ?Why is the chlorophyll in the model so much more diffuse ?Is the bathymetric slope in the model greatly smoothed (to avoid pressure gradient errors) ?Yes, the chlorophyll edge in the SeaWIFS data follows bathymetric contours along the western boundary.Bathymetry in the model was smoothed and resolution of the model is much poorer than in the SeaWIFS data.

Fig. 1
Fig.1 Simulated sea level calculated without data assimilation.
That is why we assimilate altimetry data from the Topex/Poseidon and ERS missions for the 1994-2001 years.
Ibid, lines 19-21.'Reanalysis of the Black Sea dynamics for 1971-1993 performed byMoiseenko et al. (2009)indicates that the seasonal and interannual variability of temperature, salinity and current fields are well resolved'.What does it mean 'well resolved'?What was the metrics?As Moiseenko et al (2009) was published in a difficult to find journal, it is worth to give these details in the MS.What is the difference between reanalysis performed byMoiseenko et al (2009)and this study?
. .for the period 1971-2002 was used" Page 2043, lines 4-15.The text does not provide a clear description of what and how assimilation was done.As it follows from the paper, assimilation is a crucial component of the model and should be better described in the MS.authors state that they use ECMWF 6-hourly reanalysis data, here it is monthly climatic forcing.Please clarify what meteo forcing was actually used.'The model fields demonstrated periodic oscillations at the end of integration and we have considered them as the climatic one'.Not clear what the authors mean-please clarify.
'The monthly climatic arrays of temperature and salinity ...were first interpolated on the model grid and then temporally for each day of a year by means of harmonic functions of time.'Whatdoes it mean?How many harmonics?An equation will be helpful here.'Theywerethen assimilated in the model'.What was assimilated-the climatic data or ' ...three to ten monthly hydrographic surveys per year...' ( see line 4)?'The simulations were carried out on time period of 15 yr using climatic atmospheric forcing(Staneva and Stanev,  1998).'The length of the 1971-1993 period is 23 years, not 15.Please clarify.On the C863 previous page the And to illustrate this we will add the next text with figure by the end of section 3.1: "The difference in the types of assimilated data for two period Section 4 Previous figures (e.g. 4 and 7) have only covered1971  -1993 .Figures  in this section cover 1971  -2002.I have some concerns that the change in the data available for assimilation will affect the trends previously presented and the ecosystem dynamics.I think you should present the earlier figures for the full period (to 2002) if you have the data available.Of course, the change in the data available for assimilation affects the simu-C881 lated sea dynamics.
. The maximum depth of UML (∼ 64 m) was in 1993.The upper and lower boundaries of the CIL identified by the location of the isotherm 8 • C. Seasonal thermocline formed because of the spring-summer warming of water.It manifests itself mainly at depths of 10-40 m.It can be seen that CIL thickness decreases in the spring and summer, mainly due to penetration of its upper boundary by surface waters warming, vertical advection and turbulent diffusion of heat.Temperature in the CIL increases mainly by the fall.Behavior of the average temperature on deeper of CIL likes the seasonal and interannual variability of temperature near the bottom of the CIL.The 1987 and 1989 years are characterized by medium, and 1976, 1985, 1992  and 1993 -the cold winter thermal conditions (Knysh et al., 2011).