We would like to thank the reviewer for taking the time and effort to provide detailed comments and suggestions. Please see our response below. Line numbers refer to line numbers in the revised manuscript.

This paper by Tang et al. describes a data assimilation system for integrated surface-subsurface hydrologic models, that is capable of assimilating multivariate observations and performing dual state-parameter estimates. The data assimilation system is demonstrated with an ensemble Kalman filter, at a small domain, in a synthetic data experiment, with hydraulic heads and volumetric soil moisture assimilated. The paper is very well written, with the system clearly described and well demonstrated. However, I do have some concerns about the paper:

1. The assimilation interval seems to be one day (it is not very clear from the manuscript: the authors mentioned obtaining daily synthetic observations but did not explicitly mention the assimilation interval), and the errors are evaluated every day. Given the frequent assimilation of observations, it is difficult to evaluate whether the updates of states and parameters truly improved prediction skills. For example, if the observations are assimilated every three days, can the DA runs outperform the open loop runs?

The assimilation interval in this illustrative example is one day. We have added a sentence in the manuscript to clarify this:

Line 460: “The assimilation interval is one day.”

Since the focus of this paper is to show the structure of the HGS-PDAF framework, the synthetic experiment shown in Section 4 is purely an illustrative exercise [LN1] that demonstrates how DA can be achieved via HGS-PDAF, nothing can be generalised from such a synthetic model. We agree that the assimilation frequency can have an influence on the assimilation performance, but since this is entirely illustrative, conducting an analysis of the effect of the assimilation frequency on updating is beyond the scope of this paper.

1. Many hydrologic models are designed to improve flood/drought predictions, which means that stream discharge is the most important prediction. I feel the manuscript could be strengthened by a demonstration of either assimilating discharge observations, or improving predictions of discharge.

As outlined in the reply to the previous comment, the illustrative example is essentially a purely synthetic case tailored towards demonstrating the modular capabilities of HGS-PDAF, not a real world or DA experiment analysing the effects of DA of different observation types on different hydrological predictions. Assimilation/updating of other variables/types of observations such as stream discharge is of course possible with HGS-PDAF and was mentioned in the original manuscript in the conclusions section. As an illustrative case in the paper, we selected two types of observations that are important observations for the hydrogeological modeling and which allow demonstrating how DA can be achieved for HGS with HGS-PDAF. Extending this to more variables/observations is beyond the scope of this methods-oriented paper.

1. It is not clear to me how hydraulic heads and soil water contents are updated separately. Hydraulic heads and soil water contents are connected by water retention curves. If they are updated simultaneously by EnKF, what is being used as initial conditions for the next prediction cycle? This is not explained in the manuscript.

We thank the reviewer for pointing out to us that this was not stated clearly enough in the manuscript. When hydraulic head and soil water content (in terms of saturation as saturation is the directly used variable in HGS) are updated, they are both combined in the state vector and updated simultaneously using the covariance matrix. In the example shown in the paper, when these two variables are updated together, the initial condition for the next prediction cycle was only based on hydraulic head. This is now explained in the manuscript:

Lines 460-461: “When hydraulic heads and soil water saturation are updated together, the initial condition for the next prediction cycle is only hydraulic head.”

We would like to state that the functional relationship between saturation and hydraulic head suggested by the reviewer is only applicable if unsaturated conditions are present. If the groundwater level rises, the head can still change yet the degree of saturation will be at 100%. As we are jointly simulating saturated/unsaturated conditions it is important to consider both saturation and head. Note also that the functional relationships are often associated with large uncertainties and processes such a hysteresis, which is not considered in our models. The consideration of these two variables is therefore not necessarily redundant. Given that our case is a purely illustrative example to demonstrate the modularity of HGS-PDAF, it is therefore out of scope of the paper to analyse the effects of different DA strategies when assimilating both hydraulic heads and soil water saturation simultaneously.

1. I am very curious about why assimilating soil moisture content does not seem to improve the estimates. Have the authors checked the spread of hydraulic heads and soil water saturation of the ensemble, and compared with the errors of hydraulic heads and soil water saturation observations? I feel that my last concern could be also related to this problem.

In this specific example, we did not explicitly simulate evapotranspiration, and the thin unsaturated zone thus only exists when the groundwater level decreases. Therefore, assimilating the soil water saturation only affects the unsaturated zone while assimilating the heads will change the head in both the saturated and unsaturated zones. In this very specific example, therefore, water saturation is not as informative as hydraulic heads. However, again we would like to stress that the focus of this paper was to show the structure of the HGS-PDAF framework. By considering saturation observations we intended to show the possibility of assimilating multiple observation types simultaneously (i.e. multi-variate assimilation), in this case hydraulic heads and soil water saturation. The point of the example is not to investigate the suitability of DA strategies, or the effect (and many difficult choices to be made) by jointly updating of hydraulic heads and soil water saturation. Doing this would be out of the scope of this paper.

Specific comments

1.   L83: “the coupling was neither modular nor user-friendly for…” “For” is redundant.

Deleted as suggested by the reviewer.

2.   L89: “PDAF makes it very easy to switch between different assimilation methods without the need for additional coding.” Does the observation array needs to be re-coded if the assimilation method has changed?

No, the observation array is independent of the assimilation method. Changing the assimilation method within the PDAF doesn't affect the observation array.

3.   In the manuscript, hydraulic conductivity is chosen to be modified. Based on my past experience, the parameters controlling the water retention curves can be even more important. Have the authors considered this?

We agree that the parameters controlling the water retention curves, such as alpha and n in the van Genuchten model, are important when considering variably saturated flow in the aquifer. In our case we have predefined the pressure-saturation relationship table and therefore the water retention curve doesn't change during the assimilation. These parameters are taken into account and can be flexibly added to the HGS-PDAF framework in future applications. 

1. L129: “a dual dual-aquifer configuration.”

The first dual is deleted. 

1. L159: “typical states that are considered for updating are hydraulic heads, surface water discharge, soil moisture, evapotranspiration or solute concentrations.” Discharge and evapotranspiration are not states, but fluxes.

Corrected as suggested by the reviewer. 

1. Equation (2) “the observations are perturbed by a reasonably chosen representative observation error.” This is interesting. I don’t think the classic EnKF requires the perturbation of assimilated observations though.

The classical EnKF requires perturbed observations. This was clarified by Burgers et al. (1998) but it was missing in Evensen (1994). 

1. L181: The authors may need to explain the term filter divergence. People outside of the DA community may not know what it stands for.

We have added one sentence to explain filter divergence:

Lines 187-190: “Filter divergence refers to the situation where the estimated state of the system becomes increasingly inaccurate or divergent from the true state over time. This divergence occurs when the filtering algorithm fails to effectively incorporate new observations or when the model's dynamics do not properly represent the underlying system.”

1. L431: “These observation time series were subsequently stochastically perturbed by a normally distributed error with a standard deviation of 5cm for hydraulic heads and 1% for soil water saturation.” How were the errors determined?

These observation errors are based on the prior knowledge and the tuning experiments, e.g. 5 % and 10 % have also been tested as the saturation error. As the example shown in the paper is based on a synthetic model setup and the observations are also generated synthetically, we use a relatively small observation error to better illustrate how HGS-PDAF works. We appreciate the reviewer's suggestion and this has now been clarified in the manuscript:

Lines 445-447: “The values of the observation errors are determined by our prior knowledge and tuning experiments. Different percentages such as 5% and 10% were tested and subsequently defined to provide a most illustrative use case.”

1. Figure 9: I must admit I got lost when looking at Figure 9. I am not sure what the x- and y-axes are. They look like spatial maps and I assume they are the spatial x and y directions but I am not sure.

Yes, the original Figure 9 (which in the revised manuscript is now Figure 10) is a spatial map of the model domain. We have added the x- and y- axis legend as well as the flow direction for this figure. We hope it is now clear. 

References:

Burgers, G., Jan van Leeuwen, P., and Evensen, G., 1998, Analysis Scheme in the Ensemble Kalman Filter: Monthly Weather Review, v. 126, no. 6, p. 1719-1724.

Evensen, G., 1994, Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics: Journal of Geophysical Research: Oceans, v. 99, no. C5, p. 10143-10162.

We would like to thank the reviewer for taking the time and effort to provide detailed comments and suggestions. Please see our response attached as one table is added into the manuscript and we think it would be better to display the table in a pdf file. Line numbers refer to line numbers in the revised manuscript.

We would like to thank the reviewer for taking the time and effort to provide detailed comments and suggestions. Please see our response in the attached pdf file as some equations are included. The line numbers refer to the line numbers in the revised manuscript.