We are grateful to Anonymous Referee #1 for this very positive comment about our manuscript.

We are very grateful to Anonymous Referee #2 (R2) for this very positive comment about our manuscript.

R2 is asking excellent questions about the KCC method. In our original submission, we did not provide much details about the method, since it is described in another published paper, and the focus of this study is applicative. But we are very much willing to provide additional details in the revised manuscript, and fully agree that this will help readers to better understand what is done. We provide some explanation below. Further details will be added to the revised manuscript. 

•  Dimension of x and Sigma_x:
  
  x as defined in Eq (2) is of dimension 5x251+171 = 1426. Indeed, each vector "T" in Eq (2) is a time-series of length 251 years (1850-2100), except T^{ghg}, which only covers 171 years (1850-2020, due to data availability). Consequently, Sigma_x is a matrix of size 1426 x 1426. 
• Estimation of Sigma_x:
  
  The original manuscript said that “mu_x and Sigma_x are estimated as the sample mean and covariance of the CMIP6 model forced responses.” In practice, this is done in three steps. 
  
  i/ For each CMIP6 model considered, we estimate the forced response in each of the "T" vectors shown in Eq (2). So, we estimate the forced response in GSAT, in annual and seasonal mean temperature over France, and also the response to specific forcings (i.e., NAT-only or GHG-only). We use all available members to make this calculation. As a result, we have a sample of 27 estimates of x -- 1 for each CMIP6 model considered.
  
  ii/ We compute the sample mean and variance over this sample of 27 vectors. These are our estimates for mu_x and Sigma_x. This deserves a few remarks. Estimated that way, Sigma_x is not diagonal nor block diagonal -- and covariances play a key role in KCC. The resulting estimate of Sigma_x is not invertible: the rank of our Sigma_x is 26, which is much smaller than its dimension 1426. But, one key point is that Sigma_x does not need being inverted. 
  
  iii/ In computing the posterior, the only matrix which needs being inverted is 
  
          S = (H Sigma_x H' + Sigma_y).
  
  In our implementation, this matrix S is invertible, as Sigma_y is invertible itself (see below).
• Vector y and estimation of Sigma_y:
  
  Vector y is defined in Eq (3). As of 2020, we have 171 observed years for GSAT (1850-2020), and 122 years for the temperature over France (1899-2020). So the dimension of y is 171 + 122 = 293.
  
  Regarding Sigma_y, there are in principle two sources of uncertainty at play: (i) measurement uncertainty, and (ii) internal variability (IV), i.e., variations of the temperature related to intrinsic climate variability (as opposed to forced variability). In our implementation, we assume that measurement uncertainty is null at the regional scale only (i.e., over France). But the second term, related to IV, is relatively large (i.e., we do not assume that "Sigma_y=0"), and it alone guarantees that Sigma_y is invertible. In practice, we assume IV over France to behave like an AR1 process, resulting in a multi-diagonal matrix Sigma_y (this applies only to the block of Sigma_y describing France temperature). 

Please see our responses in the attached file.