Comments by Luis Matias, University of Lisbon

Major comments

The authors apply an earthquake physical model to generated 1 Myr catalogue of earthquakes along the Carboneras Fault (CBF) and adjacent faults. The model is constrained by the assumed fault slip rates and some parameters are tuned to fit an ad-hoc Gutenberg-Richter law. The physical model generates ruptures with variable slip distribution and this feature is explored on a second part of the paper where tsunamis are generated. To our knowledge it is the first time that such physical models for earthquake generation are used in Iberia and surrounding seismically active domains. Such an effort deserves publication on itself, but additional details, and discussion must be provided to encourage the application of the model in other domains. The additional information to be provided may lead to a growth of the paper that could imply splitting it into two parts, Part 1 dedicated to the generation of the earthquake catalogue and Part 2 dedicated to the deterministic evaluation of tsunami hazard in the area. My comments will also be split according to this suggestion.

We agree with the reviewer in the fact that the first part of the paper, devoted to the physcis-based modelling, lacks of several details on the modelling and discussions on the implications and shortcomings on the asumptions made. But as is stated in the text, the presented model is already discussed and detailed in the work by Herrero-Barbero et al., 2021; which has been published in the Journal of Geophysical Research. If, as suggested by the reviewer, we publish again the details and discussions on the physics-based model, we could be incurring in self-plagiarism.

We have included additional information about the physics-based model following suggestions of both reviewers in order to be more concise on some aspects, but trying not to be excessively prolix.

Part 1:The physical model for earthquake generation

My major concern regarding this subject is the lack of relationship between observed seismicity and the Carboneras Fault. This can be inferred from Figure 1 in the paper but is made clear on figure REV-01.

As the reviewer knows, in zones of low or very low tectonic activity, the correlation between instrumental seismicity, of moderate and low magnitude, and the main faults is not direct. On the one hand, location uncertainties can be of several kilometres, and on the other, the epicenters of historical events suffer from a lack of direct observations of shallow fault ruptures. If we also take into account that the seismic cycles of these faults last thousands or tens of thousands of years, it is logical to expect that the instrumental seismicity of a few decades will not reflect the seismogenic behavior of large structures, hence the interest of physics-based models.

The paper mentions that some model parameters are tuned so that the final Gutenberg-Richter (GR) law has a b value equal to 1.0. The paper fails to give the support for this assumption and no information is provided on the a value that also characterizes the GR law. Assessing the ISC catalogue and selecting a generous area surrounding the CBF we obtain the GR law shown in figure REV-02, where the number of earthquakes is scaled to 1 Myr as in the paper.

We obtain a very high b-value, not common for convergent or transcurrent domains, showing that large magnitude events are much less frequent than found on average on the earth. This may be a feature due to the small number of events, but it deserves discussion. The thickness of the brittle layer assumed for the physical model deserves additional discussion in the light of information provided by the earthquake catalogues and deep structure studies in the area.

As the reviewer knows, in zones of low or very low tectonic activity, the correlation between instrumental seismicity, of moderate and low magnitude, and the main faults is not direct. On the one hand, location uncertainties can be of several kilometres, and on the other, the epicenters of historical events suffer from a lack of direct observations of shallow fault ruptures. If we also take into account that the seismic cycles of these faults last thousands or tens of thousands of years, it is logical to expect that the instrumental seismicity of a few decades will not reflect the seismogenic behavior of large structures, hence the interest of physics-based models.The calculation of a Gutenberg-Richter fit requires that the magnitudes of the events used be homogenized in order to be comparable, in addition, a completeness analysis must be done to filter the events by date and the fit should preferably be done with a maximum likelihood adjustment. Nor is it possible to extrapolate the seismicity of a few decades in a seismic cycle of thousands of years, to a behavior of hundreds of thousands or millions of years, for this reason the value of “a” of the Gutenberg Richter law, which depends on the seismic productivity is not used, but the “b” value is compared so that the distribution of the size of the events is similar to the real one.

As explained in Herrero-Barbero et al. (2021), one of the criteria for choosing the best-fit model parameters is that the b-value be close to 1, considering always the same completitude magnitude between several synthetic catalogs. This b-value is justified by the estimations in the same seismogenic zone in previous works based on instrumental seismicity (García-Mayordomo, 2005; IGN-UPM, 2013; Villamor, 2002), and is also a reference value as assumption in numerous papers of synthetic seismicity modeling (e.g., Console et al., 2017; Shaw et al., 2018). These references have been included in the text (lines 151-154).

The seismogenic crust thickness of the model is based on previous seismotectonic studies at Southeastern Spain (García-Mayordomo, 2005, Fernández-Ibañez & Soto, 2008; Mancilla et al., 2013, Grevemeyer et al., 2015). These references have been included in the text (lines 126-128).

If we compare the compute GR law with the earthquake distribution published, as shown in figure REV-03, we remark that, as suspected, the observed seismicity is much less than the modelled catalogue. This feature deserves to be discussed in the paper.

As has been answered previously, the amount of seismicity generated cannot be directly compared to the existent instrumental seismic catalogue; nevertheless details on the fit of the physics-based model, and discussions on the representativeness of the seismicity generated can be found on the Herero-Barbero et al. (2021) work. We have included this aspect in the text in lines 159-162.

Another constrain on the physical model is the average slip rate on the CBF and neighbouring faults. The authors used for the CBF the value 1.3 ± 0.2 mm/yr proposed by Echeverria et al. (2015). We quote here the Echeverria et al. (2015) sentences (CFZ = Carboneras Fault Zone): “The analysis of GPS data in the SE Betics confirmand quantify the ongoing tectonic activity of the onshore segment of the CFZ as a left-lateral strike–slip fault. For the first time, we were able to provide a quantitative measure of the present-day horizontal geodetic slip-rate of the CFZ, suggesting a maximum left-lateral motion of 1.3 ± 0.2 mm/yr. The coincidence of geologic and geodetic strike–slip rates along the CFZ, illustrates how during Quaternary its northern segment has been tectonically active and has been slipping at a rate of 1.1 to 1.5 mm/yr”.

It is clearly suggested that 1.3 is a value valid only for the onshore segment of the CBF and it represents the maximum value. The use of this value deserves further discussion as well as the consequences for its variation along the fault, particularly on the ocean segment. Furthermore, our interpretation of Echeverria et al. (2015) figure 5 is that the CBF slip rate, as measured by GPS, lies between 1.0 and 2.5 mm/yr.

The only way to measure the slip-rate of faults is by means of geodetic data (GNSS or InSAR) or by means of geologic data (paleoseismology, tectonic geomorphology). From the geological point of view Moreno et al. (2015) obtained a minimum strike-slip rate of 1.31 mm/yr for a section in the onshore segment of the Carboneras fault, and a minimum dip-slip rate of 0.18 mm/yr. The offshore segment of the Carboneras fault has been extensively studied also by Moreno (2011); obtaining a minimum strike-slip rate of 1.3 mm/yr and dip-slip rates between 0.1 and 0.3 mm/yr. The geodetic data obtained to date is from the onshore segment by means of high resolution GNSS campaign measurement and permanent stations, with the preferred value of strike-slip rate of 1.3 mm/yr (Echeverria et al., 2015). Taking into account that the onshore and offshore obtained slip rates are the same, and that the onshore geodetic data coincides with the geological data, we don’t have any reason to propose any different slip-rate. We have included this short discussion on the manuscript, in lines 132-142.

All the model parameters used, as well as the associated bibliographical references, will be presented in a Table in the Supplementary Material. They correspond to the compilation made in Herrero-Barbero et al. (2021), shown in the Table 3 of the cited paper.

It may be relevant here for the authors to mention other sources of information on the fault slip of oceanic faults as provided by Neotectonic modelling. While Jiménez-Munt & Negredo (2003) and Cunha et al. (2012) provide slip rates estimates smaller than 0.5 mm/yr, Neres et al. (2016) show a maximum value of 1.7 mm/yr for the CBF.

Thanks to the reviewer for pointing this out. We have included the references in the text in linex 140-142.

We understand that a perfect coupling is assumed for the CBF and neighbouring faults between the kinematic constrain (slip rate) and the earthquake generation. The typical seismic coupling of major plate boundary types has been discussed by Bird and Kagan (2004). They showed that for continental convergent boundaries it lies between 0.51 and 1.00 (1.00 preferred value) and for continent transform faults it lies between 0.38 and 1.00 (0.72 preferred value). The author’s choice of 1.0 deserves some discussion and the consequences of using a different value should be addressed.

It is a good observation that we have missed to mention in the text. It has been included in lines 170-171.

It seems that the generation of “characteristic earthquake” recurrence models is a feature of the physical model used. This model has not been used of PSHA and PTHA in Europe and additional discussion should be provided. Is it a model feature or is it explained by some characteristics of the CBF domain?

The pseudo-characteristic earthquake behavior arises from the model configuration. Similar results have been observed with other physics-based models (Console et al., 2021; Rafiei et al., 2022; Shaw et al., 2022). The references have been included in the text in lines 335-336.

Another feature of the physical model for earthquake recurrence applied to the CBF system is that the maximum magnitude exceeds the estimations made by several authors. It is argued that the maximum magnitude value lies at the extreme boundary of some estimates. Is it a feature of the model? Why does it happen? Some additional discussion is needed here.

The maximum magnitude exceeds that of some authors, but coincides with the estimation of others. Additional comments are included in the text in lines 182-185.

The physics-based earthquake generation model, besides the definition of the geometry (that itself deserves additional support), requires many parameters, some to be defined and others used as constrains. A list of the main parameters and information on its choice must be well presented, which is not the case on the current version of the paper. If possible, several runs of the model could be used, first to assess the random uncertainty and then to assess the epistemic uncertainty due to the choice done for some selected critical parameters. One set of the parameters, or features, of the physics-based earthquake generation model that is not explained at all, and is relevant for the higher magnitude events, is the one that rules the multi-segment propagation.

As has been mentioned earlier, all the information required by the reviewer is already published in Herrero-Barbero et al. (2021). We consider it not necessary to repeat the same discussions and details here because we could be incurring in self-plagiarism. Additionally we have included additional information as supplementary material.

Part 2:Tsunami

Here the major comment regards the absence of any reference to previous works on tsunami hazard assessment in the area. The major study requires some reference and comparison is the NEAMTH18 (Basili et al., 2019, 2021). The results of this study can be assessed online 1 . The area investigated by the paper is shown in figure REV-04 as well as the probability results for the Almeria forecast point.

As the reviewer comments, the approximation of NEAMTH18 is probabilistic, incorporating all the potential sources estimated in the area, as such, the results are not directly comparable with our deterministic approach. Nevertheless we have mentioned NEAMTH18 in the text in lines 259-263.

Additional comments

On the pdf provided most of the figures are small and difficult to read.

We think that this is a problem related to the LaTeX manuscript template and in the final version the figures will be bigger.

Lines 21-22: this ability depends on their mode of seismic rupture

In fact, this is not the major parameter defining the potential for tsunami generation by an earthquake, as can be verified by the decision matrix adopted for the ICG/NEAMTWS (2011). The main parameters are:

• Top of the fault depth (focal depth if that information is not available)

• Location in relation to the coastline

• magnitude

Given the high directivity of tsunami propagation we may add also the strike of the structure. We suggest the authors to frame better the above-mentioned sentence.

We have to disagree with the reviewer. It is true that in the used decision matrix the rupture characteristics are not taken into account, but basically because the obtention of such rupture characteristics for an earthquake is a process that lasts several minutes in the best case, and consequently cannot be used in an emergency decision matrix. The above assertion is backed by several studies referenced in the text. It is true that the depth of the fault and the magnitude are of prime importance and we omitted it from the text as is obvious, but following the suggestion of the reviewer we have modified the text accordingly (lines 21-22).

Line 27: Although the lower capacity of strike-slip faults to generate tsunamis is a proven fact,

Mention here the strike-slip generated tsunamis in the Gloria Fault, one domain close to the one investigated in the paper and both belonging to the Nubia-Eurasia plate boundary, as presented by Baptista and Miranda (2009).

We have included the reference in the text (line 29).

Line 33: based on the simulation of tsunamis generated by ruptures of simple, rectangular

This is one occasion to mention the NEAMTHM18 model (Basili et al., 2019, 2021) that covers the investigated area and didn’t use simple rectangular sources.

Ok.

Line 43:

Add to the list the NEAMTHM18 model (Basili et al., 2019, 2021)

In NEAMTHM18 no dynamic rupture models are used and consequently cannot be included as that in line 43.

Figure 1

The location of Almeria is missing, and it is needed.

Location of Almeria has been included.

The geographic projection is “Plate Carrée” which is very unusual. Renders the comparison with other maps difficult. Why this projection was used?

The map shown in figure 1 has been generated with QGIS using a standard Mercator projection. Maybe it has been slightly modified by the vector drawing program used (Inkscape) but we think that for the purpose of the map (a location map) is precise enough.

Line 73: this fault has been proposed as source of the 1522 Almeria earthquake

This sentence is in contradiction with the location of the event in Figure 1. Clarify.

As is stated in the figure 1 caption, the location of the 1522 event is from the IGN-UPM (2013) catalog; while the proposition of the Carboneras fault as source is from Reicherter and Hubscher (2007). It has been clarified in the text as suggested (lines 79-81).

Line 79: Although the tsunami simulations done to date

Add to the list the NEAMTHM18 model (Basili et al., 2019, 2021) and comment briefly its methodology.

In NEAMTHM18 the modelling approach for faults as the Carboneras fault is done by means of “Gaussian-shape unit sources”, and consequently the surface deformation is approximated based on these predefined seabottom deformations. There is no direct modelling of the Carboneras Fault itself. We have rewritten the sentence (lines 87-89).

Lines 108-109: fault depths between 8 and 12 km.

It is no clear what the authors mention as “fault depth”. Is it the width of the fault? On which data is this information based?

The sentence has been partially rewritten (lines 126-128).

Figure 2

Explain in the caption the shorthand terms used, e.g., ASMF, PF, …

Done

Line 115: Besides the input kinematic data

Our interpretation is that the authors model the ruptures on the CF and also on other faults to the NE. What is the slip rate on these faults? What was the source of information?

This information, as is stated in the text, can be found in the work describing the physics-based model of Herrero-Barbero et al. (2021). We have included the reference in the text.

Line 118: We define reference rate-and-state values based on experimental data

This experimental data was likely sampled at shallow depths. How do they apply to the expected seismogenic depths? If they are assumed to be identical on the whole fault system, what might be the consequences of this simplifying hypothesis?

The experimental data refers to fault rock mechanics laboratory tests to describe the fault zone rocks frictional behavior expected at seismogenic depths. These experimental values are the starting point before testing several synthetic catalogs with variation of the rate- and state-dependent coefficients, until achieving the best-fit values (Herrero-Barbero et al., 2021). We have clarified the sentences in the text (lines 148-151).

Line 122: rate-and-state friction parameters a=0.001 and b=0.010;

There might be a confusion with the a and b parameters of the Gutenberg-Richter law.

Yes we know, but these are the terms used in the literature. We added a caveat (lines 155-156).

Line 125: also suggested possible creeping sections in the Carboneras fault

How would this hypothesis affect the paper results? Some discussion on the assumptions and simplifications made should be provided. See main comment earlier in this document.

We added some discussion in the “Earthquake ruptures simulation” section concerning the interevent time dependence on fault coupling (lines 159-162, 170-171).

Figure 3

Given that the physics-based earthquake generation model is applied to a set of faults, it is not clear if the histogram applies only to CF or to the whole system as depicted in Figure 2. It is assumed that the a priori average slip rate of the CF is respected by the model, but a sentence presenting the a posteriori compute slip rate is needed.

As is stated in the caption the data presented is for the Carboneras fault.

Line 132: excluding the aftershocks

Since it seems that the authors are discussing the 6.5 to 7.0 magnitude interval, what is the definition of aftershocks used?

We have to apologize for the confusion, the data presented here includes aftershocks. The reference has been deleted.

Line 137: therefore the released seismic moment.

What is the shear modulus value used to compute seismic moment? Justification for that value? Is it uniform along the fault and over depth? What are the consequences of using a single average value for the modelling?

The description of the details of the model can be found on Herrero-Barbero et al. (2021). We have included here the information on the shear modulus used in line 164, which is 30 GPa, a value widely used as average shear modulus for the upper continental crust.

Line 138: The epicentres

How does the code obtain the slip initiation? Is it a feature of the code or the epicentre is computed a posteriori from the rupture distribution?

The rupture initiation in the code is governed by the rate-and-state formulation, when a “nucleation” state is reached by an element the code, spontaneously, computes the rupture propagation to the neighboring elements on a pseudo-dynamic approximation. The details can be found on the works describing the RSQsim code (Dieterich, 1995; Dieterich and Richards-Dinger, 2010; Richards‐Dinger and Dieterich, 2012). Included in lines 177-180 and in an brief explanation of the physics behind the model in lines 105-123.

Lines 154-155: As the sea-floor deformation generated by the earthquake is usually transferred instantly to the elevation of the water free surface

This is not true in general, though it applies to the modelling of far source tsunamis. For locally generated tsunamis there are two effects that are not considered in the paper that deserve a comment: i) the finite compressibility acts as a filter when computing the sea surface deformation (e.g. Lotto & Dunham, 2015); ii) the horizontal movement of the sea bottom, in areas of relief, generate an initial velocity on the water that, in some circumstances, must be considered.

We have reworded the sentence (lines 198-199).

3. Tsunami modelling

There are a few general questions that must be addressed by the authors.

1) What are the boundary conditions used for the water borders?

This information has been added in line 236.

2) What are the elastic parameters used to compute the seafloor deformation? Justification?

This information has been added in line 194.

3) What happens close to the coastline? Is friction used? What are its characteristics?

This information has been added in lines 236-237.

4) Is there inundation?

This information has been added in lines 236-237.

5) How is the tsunami amplitude computed? It is recommended that the tsunami wave amplitude to be computed at cells with water depth no smaller than 50 m. The reason is explained in Kamigaichi (2011): “To represent the tsunami waveform correctly in a shallow sea area, very fine bathymetry data mesh is necessary (in a strict sense, 20 or more grid points are necessary within one wave-length [31]), and a vast time is required for the completion of such detailed calculations. To overcome this difficulty, the numerical simulation with the long-wave approximation is applied only to points which are a few to a few ten kilometers seaward from the coast (“forecast points”) where sea depth is about 50m. Then, tsunami amplitude at the coast is calculated by using Green’s law described in the next section.”

We have modelled inundation at the coast and consequently there is no use of the green’s law.

Line 204: maximum wave elevation

The meaning of this parameter must be well explained. See my previous comment.

I think I don’t fully understand the reviewer concern with the term. It is the widely used term to describe the maximum elevation reached by the free water surface at a point of the calculation grid on a propagation.

Line 211: relevant local inundations

It is not explained how “maximum elevation” is converted to inundation. See previous comments.

The inundation is computed by means of the COMCOT numerical model.

Figure 8

Explain the dashed contours.

Added to the caption.

Line 218: have been taken for the 5 m depth isobath

wThis explanation should have been provided earlier. Given previous comments 5 m seems not appropriate. What happens if there is no cell at 5 m depth? Use the Green law to convert it to 5 m?

As the inundation is computed there is no need to use the Green law.

Lines 224 – 234

Given that a single, randomly generated, catalogue was used, I fail to see the relevance of the discussion on these details of the tsunami amplitude histogram.

Would another catalogue generate the same features?

The range in maximum magnitudes would be probably the same as it depends on the dimensions of the structures. If we vary the slip rate or the seismic coupling, the interevent times would be different, and consequently the statistical distribution of maximum elevations and inundations.

The crucial point here is that, when dealing with deterministic approaches usually a single worst case scenario is used. In the best case a set of worst case scenarios considering the uncertainties is done. With these models we can obtain a physically coherent seismic behavior through hundreds of seismic cycles that produce hundreds or thousands of destructive events that can be modeled to characterize not only the worst-case scenario, but also characterize scenarios linked to probabilities of exceedance or return periods.

We have extended our reasoning in the text in lines 304-312.

Line 239: generate locally damaging tsunamis

Define “local” tsunamis. The term “local” has a very specific meaning in the tsunami warning systems (ICG/NEAMTWS ,2011)

We think that the term “locally damaging tsunamis” accompanied with the maps showing the areas where the maximum waves are expected is clear enough. The terminology used in warning systems are crucial for the development of warning systems, but not for the object of our study.

Line 243: If we compare the results of this work with previous results

Given that the “frequency” of tsunamis was mentioned in line 241, comments on the NEAMTHM18 model (Basili et al., 2019, 2021) are appropriate here.

We refer here to the results of our work as we think is clearly stated in the text.

Line 248: Mw 7.62, which is close to the maximum magnitude proposed by Moreno

(2011).

This value and reference were not mentioned in the introduction and they are relevant.

It has been included in line 84.

Lines 252 – 253: On the other hand, the rake used in our models is 10° while Alvarez-Gomez et al. (2011a) used 15o and Gomez de la Peña et al. (2022) used 0°.

What are the consequences of the uncertainty on the rake to the paper results?

Additional comments have been included in lines 302-303.

Lines 264 – 265: we can obtain the maximum magnitudes in a robust manner from

a statistical point of view.

Given the larger number of simplifications and approximations used in the physics-based earthquake generation model, given that a single catalogue was generated without assessing aleatoric and epistemic uncertainties, I cannot classify the results as “robust”, though deserving to be published.

We have reworded the phrase.

Lines 270 – 271: relationship between the size of the earthquake rupture and the slip;

Not shown in the paper. Show as supplement?

The earthquake ruptures are included now as supplementary data.

Line 278: would reflect a GR relation

In fact, one of the most frequent earthquake recurrence laws is the truncated GR relation, not the simple (and open) GR law.

We are referring to both flavors when mentioning the Gutenberg-Richter relation.

Line 285: see supplementary models

These are not available on the documentation provided.

Our mistake. Now is included.

Line 295: the simplified model overestimates

It should be mentioned that nowadays the common procedure is to taper the uniform slip at the borders of the rectangular faults (e.g. Davies and Griffin, 2018). Given this, the comparison between tsunamis generated by irregular and uniform slip faults is unfair, for the tips of the fault as mentioned in the text.

Additional comments are included in lines 359-360.

Figure 13

The labels mentioned in the caption cannot be seen on the figure. Too small?

We think that this is a problem related to the LaTeX manuscript template and in the final version the figures will be bigger.

Line 309: allows a more robust characterization of the scenarios

Given the larger number of simplifications and approximations used in the physics-based earthquake generation model, given that a single catalogue was generatedwithout assessing aleatoric and epistemic uncertainties, I cannot classify the results as “robust”, though deserving to be published.

We have reworded the phrase.

Figure 14

The labels mentioned in the caption cannot be seen on the figure. Too small?

We think that this is a problem related to the LaTeX manuscript template and in the final version the figures will be bigger.

Lines 316, 317: which allows a robust implementation of uncertainty estimation

Given the larger number of simplifications and approximations used in the physics-based earthquake generation model, given that a single catalogue was generated without assessing aleatoric and epistemic uncertainties, I cannot classify the results as “robust”, though deserving to be published.

We have reworded the phrase.

Lines 332, 333: The implementation of these methodologies in the Probabilistic Tsunami Hazard Analyses (PTHA) is a logical and necessary step.

Comment/discussion of the NEAMTHM18 model (Basili et al., 2019, 2021) are needed here since they apply to the same area discussed in the paper.

We do not see the necessity of including in our conclusions this reference again, as has been now included several times throughout the text.

Line 337: GMT (Wessel et aI., 2013) has been used to

Given this information we do not understand the use of the “Plate Carrée” projection in Figure 1.

As is stated in the text: “GMT (Wessel et al., 2013) has been used to perform some calculations and to produce most of the figures.” not all. The map shown in figure 1 has been generated with QGIS using a standard Mercator projection. Maybe it has been slightly modified by the vector drawing program used (Inkscape) but we think that for the purpose of the map (a location map) is precise enough.

Comments by Julián García-Mayordomo, Geological Survey of Spain

Comments:

- A proper cite of the earthquake catalogue used in figure 1 (apart from IGN-UPM (2013)) is missing. You should also cite the paper of Cabañas et al., 2015 in which the homogenization process to Mw is explained. Additionality, it is important to know that the catalogue you are using in fig 1, apart from being homogenized in size to Mw, was also declustered (IGN-UPM, 2013).

The reference and observation about declustering has been included in figure 1 caption.

- Could you extend a bit more on explaining the earthquake ruptures simulation method? I know that a proper understanding would require the reader to go to the original source, but it would be good to show an equation (for example) in which the parameters are shown (a, b,...). If it is possible.

More details on the formulation of the simulation has been added in lines 105-123.

- I assume you are using the rate and state parameters provided in Herrero-Barbero et al(2022) after a testing process that was done in that paper, is that right? If so, could you make it explicit in the manuscript; otherwise it seems that the testing process was performed in the frame of this paper.

We added a clarification in the text in lines 148-149.

- I see that references cited in the body of the manuscript are ordered alphabetically. The authors should check if that is correct according to the journal format; as references are usually ordered in increasing year of publication.

We have used the Copernicus bibliographic style for LaTeX provided by the journal. Nevertheless we will check the consistency of the references.

- I have made few more comments and suggestions for improving clarity in some parts of the text and figures right on the pdf attached. Please, in case you cannot read them (done with the adobe acrobat tool) get me back and I will pass them to you in a different format.

- I found some typos in the manuscript, also marked in the attached pdf.

We have reviewed the minor comments and suggestions from the referee and most have been included in the manuscript.