The authors developed the process-based model that coupling carbon, erosion, transport and deposition processes. They brought in the lateral movement of carbon into process-based land carbon modelling. The community has been expecting this advancement, especially in large scale studies, for long, but with limited progresses partly due to the high computational cost. The authors presented a detailed and exciting case. The well-illustrated their model formulations, how they tackled on the computational bottleneck, the calibration and validation, with adequate discussion of limitations and future improvements. The study is well-designed and informative. The writing is generally clear despite some part might be a little lengthy. My criticises are between major and minor. Please check below. 

It is not clear to me how the matrix-relevant techniques helped the current study. The equations for the lateral carbon fluxes are presented mostly in carbon balance equations (i.e., no need for the matrix form). Are the matrix techniques only used for constructing the ORCHIDEE emulator? If the ORCHIDEE output only an input to CE-DYNAM, or any parameter changes that require a re-do of model spin-up that requires computation resources? 

Parameter values for soil discretization are optimized.  How about parameter values for other parts of the model? I might miss some part, but is there any table or supplementary information that documents values of relevant parameters used in this study?   

Minor comments: 

Line 1-11. Is the background for abstract a little lengthy? 

Line 24. I would suggest DayCent. For people without background, they don’t know what DayCent refers to. DayCent is not mentioned in texts other than here. Besides, the CENTURY vs. DayCent is another layer that needs background. So DayCent here is not necessary to add more information to this already complex manuscript.  

Line 88 add in before Table 1

Line 92. Is it better to write ORCHIDEE CENTURY as ORCHIDEE (CENTURY-Carbon type) or other better rewordings. ORCHIDEE CENTURY is confusing. 

Line 93. Is the first order kinetics necessary? Current application is with linear-models due to computational cost. In theory, for example, if our studying region is small, the coupling with an nonlinear carbon model is possible, right?   

Line 125. Could we write it as “ induced by the terrain slope (S[x,y]) and the flow accumulation (w[(x,y)])” to reduce confusing? 

Table 1. Could you use other symbols for ω[(x,y)] vs w[(x,y)] (?)? They looks the same 

Equation 6. Lines 200-205, and across the manuscript. By erosion b – t, you mean from lower soil layer to the upper soil layer and the flux is there between adjacent soil layers? Please clarify, by b – t, it refers specifically from the third to the first layer, in your context. So it is not clear to me why “losses from the layers below must be added to the layers above”

Line 525. Is the breakdown of aggregates a transport process? Please clarify. Aggregates breakdown could happen without the transport process. 

Line 525-530. If the “halfway between….”a cause of the difficulty in finding the optimal resolution, or the computation cost and applications? 

The authors provide a spatially upscaled version of an eriosion, transport and deposition model (CE-DYNAM) to European scale, while sticking to a high spatial resolution. This is an important scientific contribution because lateral C transport has largely been ignored so far in process-based models because usually models consider different grid cells as independent. Such lateral transports can have a significant effect and should be investigated because nonlinear effects could potentially lead to drastically new insights and better understanding compared to existing models. 

The spatial upscaling was computationally feasible only because the authors emulated the original model by a matrix approach. This allows the application of sparse matrix models as well as an improved application of parallel computation methods. Another advantage of matrix models is that they allow a rigorous mathematical analysis, something the authors did not do in this manuscript (it was not their goal) but can be done in the future based on the matrix reformulation. This was not possible with the original implementation. Furthermore, I appreciate all the effort the authors showed in reimplementing an existing model with a matrix approach, emulating the original simulation results very vell. In this regard I like also the explanation of the matrix shape in Section 3.3. Nevertheless, I wished I could have "seen" a matrix, at least as a schematic block matrix after the authors speaking so much about matrices.

The presentation and the writing is clear, explaining the model calibration and simulation results as well as model limitiations and future opportunities very well. Sometimes though the text appears too lengthy in my opinion. This starts with the abstract and continues with quite some overlap in the sections about calibration, results, simulations and limitations. Furthermore, in particular in later sections I was overwhelmed by an extensive use of potentially unnecessary numbers.

I furthermore do have some issues about an easy reproducibility of the method, because to me it seems that in Section 2 some formulas are incorrect and notation is not precise. Well, either the formulas are incorrect or I understood them wrong, neither option is preferrable. In particular I consider the statement (p. 3, l. 71): "We expect our mathematical development of CE-DYNAM to facilitate its reproduction and incorporation in LSMs such as ORCHIDEE, DayCent, and others." to be pretty bold. Under this point of view formulas and notation should show no flaws.

This starts with Table 1, which is in general very nice, but it is incomplete, some symbols that are used later on do not appear here. This was sometimes annoying for me while reading. For example I* was never properly introduced. On p.5, l. 132 it is stated that inputs are disturbed along D, which I do not understand.

Here my issues with notation and formulas, which should be thoroughly checked:

p. 8. l. 1: Should it be Δ_i? Is it used to compute d₀=0, d₁=d₀+Δ₀ and so forth?

p. 8. l. 183: k_e is missing in Table 1. The "Total mass of soil" under the bracket refers to total mass of soild in hill slopes? Why is it decreasing with h? To me this looks like we have an infinite erosion if we do not have any hill slopes, isn't this going in the wrong direction?

p. 10, l. 202. If there is going to be a new input flux, should there also be a new output flux? The input must come from somewhere. Please also indicate where the new input flux goes, I do not want to guess here.

p. 10, l. 213. The way k_τ (missing in Table 1) is defined two lines later, it seems wrong to me to call it a flux. It is a rate (dimension 1/time). It will only become a flux once it is multiplied with a pool content. This leads to major confusions for me later on.

p. 11, l. 235: Please use other notation than (a, b) for the indexing of the sum, the two letters are already taken. It becomes very confusing this way. Again, I think that if k_τ is a rate rather than a flux, then k_s will be as well because all the other factors in Eq. 11 are dimensionless.

p. 10, l. 241 "all PFTs" should rather be "all PFTs but p₀", right?

p.12, l. 247: "Such input flux". I disagree again, same problem. It is not yet a flux because it is not yet multiplied with a carbon stock, which should be S_{[(a, b), ...]} here?

p. 12, l. 257: The "P" here looks different from the ones introduced on p. 10, l. 221. Furthermore, I am not sure whether it belongs here in the first place, well the absolute value is unecessary in any case. But why multiply by the number of non-zero PFTs in cell (a, b)? Shouldn't this be already included in k_τ already? Could you write it down explictly for yourself without the P but a second sum instead and check whether it is correcct this way?

p. 12, l. 260: In k_s I think that source and target are confused.

p.13, l. 265: What is k_t^*?

p. 13. l. 266. Should k_s with source (a, b) be multiplied with some stock indexed by (a, b) instead (x, y).

I obviously do have some confusions about the firs-order description, where sources and targets seem not to match, at least in my head. So would like to encourgae the authors to carefully check the notatian and the formulas again, along with their implementation.

Small issues:

• In general units are sometimes italic and sometimes not, sometimes with a space between the number and the unit, sometimes without.
• The use of singular and plurar gets mixed up quite often.
• p. 3, l. 86: "approach" --> "an approach"?
• p. 3, l. 88: Table 1 --> (Table 1)?
• Table 1:
•     - Description sometimes ends with a period, sometimes not. 
•     - The adimensional respiration rates of carbon k_r actually do have a dimension: 1/day.
•     - I^{[c_j, c_j]}: One of the j's should be an i. Probably the first one, then please also adapt the description, to make it consistent with k_t.
•     - Then it is a little unfortunate to use ω and w for the depth to bedrock and the flow accumulation, respectively.
• p. 7, l. 145: I am not sure if such a procedure necessarily converges to a pullback attractor (given there is one in the first place). The pullback attractor is reached when starting the simulation earlier and earlier, basically moving toward an infinite simulation history.
• p.7, l. 152: Why now change the notation from S, I, and k to x, B, and A?
• p. 10, l. 243: "soi"
• p.14, l. 280: The seventh-largest what?
• p. 15, l. 308: u_i instead of p_i?
• p. 15, l. 312: "P is Panagos", what does this mean?
• p.15. unnamed formula: what is r(y, m)?
• p. 18, l. 398, 399: What are C-factor and "R factor"? Please also note the different way of writing them.