This paper presents a nice example of combining theory based models and machine learning to efficiently identify parameters of an ecosystem model, exploiting observation data recorded at multiple sites. The approach is valid and the results are interesting. However, the documentation of data and methods is currently deficient on a level that makes it hard to grasp the main messages and interpret the results. Section 2 of the paper does in my yes require a thorough revision, including new explanatory figures, restructuring and replacement of text blocks. For this reason I recommend a major revision or rejection with an invitation to resubmit. 

Major comments

1. I assume a key point of the developed framework is that it enables to directly backpropagate from the outputs through the model equations to the neural networks. This is not clear from the paper at all. Much of the framework description seems like you feed NN predictions of parameters through a black box physics based model, which is a standard approach. I suggest a dedicated subsection, possibly including a figure, to clarify this detail.

2. The datasets used for training and testing are not properly documented. We don't know how many datapoints are included over which time periods. The random holdout suddenly appears in the results, and in general we don't know how training/validation/testing splits are defined. CLM4.5 standard parameters play a central role in the results, but we know nothing about where they come from / how they are defined and if, for example, all or a subset of values are used for comparison. 

3. The explanation of the ecosystem model suffers from a clear struggle between trying not to include the entire set of equations in the paper, while providing sufficient detail. For me the level of detail provided in the paper was actually confusing, because it required constant looking up in the appendix to understand the context, distracting from the main messages. I think a way out could be to include a figure that summarizes the main blocks of the model (including what parts correspond to f1 and f2), include only the changed equations in the paper, and otherwise keep the full model description in the appendix. On a sidenote: is f2 not the same as an observation equation, that is commonly used in state space models?

4. Details on hyperparameters (neural network # of layers, activation functions, learning rates etc.) are not provided at all. Some key information should be provided in the paper, and a reference to supporting information or the code should be provided for details.

Detailed comments

line 61: nonuniqueness is also going to be a problem if we employ newer frameworks like PINNs or dPL

line 110: it might be worthwhile to start with a reference to figure 1 and a down to earth explanation of the objective of your work, i.e. to calibrate model parameters across many sites, to capture the variation of parameters using neural networks, and to employ differentiable programming to speed up the identification process

line 118: please explain PFT again in this section

line 140: If you preserve eq. 4 and 5 in the paper, I think they should be presented in reverse order (f1 first, f2 second)

line 146-164: please include only methodological descriptions that are relevant for the results. of the julia implementation was not used, then it should not be described and discussed

line 183: you don't describe anywhere in your data how many PFTs you consider. it is therefore here also not clear how many dummy variables this model receives as input.

line 190-205: I think this information is not needed to understand the main message

eq. 10: why is psi_max replaced by psi_0? (missing explanation)

eq. 11: what is F_om?

line 232: the CLM4.5 data points should be documented in a dedicated data section. in general, i suggest they you separate the description of data and experiments

line 239: were all calculations performed only for the topsoil layer in all experiments?

Table 1: missing symbol explanations for means and standard deviations

line 383: please include time series for observations and model predictions

fig. 5: symbols in legend cannot be distinguished. are results shown for the test dataset?

line 426: i would add that you have identified parameter values that are optimized for the considered set of model equations and forcings. both of these have limitations. equations may be wrong, ERA5 is rather uncertain, and measurement principles can vary between stations. This is both a limitation and a strength of your framework. Parameter values will not be transferable to other inputs. On the other hand you can obtain optimized predictions for the given set of forcings.

The authors of the manuscript ‘A differentiable ecosystem modeling framework for large-scale inverse problems: demonstration with photosynthesis simulations’ describe the application of the ‘differentiable parameter learning’(dPL) framework to the photosynthesis module of FATES model. The framework, and concept, overcomes extrapolation limitations from site-by-site calibration approaches and allows leveraging information content in large-scale datasets towards a global parameterization of photosynthesis models. Neither the concept (Tsai et al., Nature Communications, https://www.nature.com/articles/s41467-021-26107-z, 2021; Bao et al., Authorea, https://www.authorea.com/doi/full/10.1002/essoar.10512186.3, 2022) nor the dPL framework (Tsai et al., 2021; Feng et al., 2022ab) are new. However, the framework is used in the FATES model for the first time and the results would be of interest for further model development, but also to the scientific community at large.

At this point, the experiment focuses on inverting two parameters, V_cmax25 and B, resulting in that the accuracy of the simulated net photosynthesis rate being slightly improved. The main concerns at this stage relate to apparently incorrect formulations of some key equations, to issues about the validation strategy, to the fact that the forcing data and the experiments are not described sufficiently, challenging the acceptance of the study, while hampering any reproducibility efforts. Please see below for details.

Major comments:

1. Two key equations are incorrect in the paper:

1) line 140: equation 5, C_i=C_a-A_n*P_atm*(1.4g_s+1.6g_b)/(g_s+g_b); 

2) line 505: equation A1, A_c=V_cmax*(C_i-Γ_*)/(C_i+K_c*(1+K_o/O_i)).

According to the user guide of the FATES model (https://fates-users-guide.readthedocs.io/projects/tech-doc/en/latest/fates_tech_note.html#fundamental-photosynthetic-physiology-theory), the equations should be:

1) C_i=C_a-A_n*P_atm*(1.4g_s+1.6g_b)/(g_s*g_b); 

2) A_c=V_cmax*(C_i-Γ_*)/(C_i+K_c*(1+O_i/K_o)).

Since the FATES model is reimplemented in Julia and PyTorch by the authors, the codes might be also wrong. If so, the unit of C_i will be incorrect, leading to errors in the inversion of V_cmax25 and B. The wrong computation of the effective Michaelis-Menten coefficient (=K_c*(1+O_i/K_o)) might only have a slight effect if the temperature is close to 25°C, but should be concerned if the temperature is too low or high (and I do see some points with low leaf temperature in the ‘Lin15’ database). Thus, I have doubts about the current results and relevant analysis.

1. As all the results are validated only once using the temporal holdout data or the random holdout data, the generalizability of the dPL (or NN_B+NN_v) is not clear. If the N-fold or leave-one-out cross-validation can be adopted, the statistical metrics can be more justifiable to reflect the model performance.
2. The forcing variables and parameters are not clearly differentiated in the paper. For example, is the leaf layer boundary conductance, g_b, a constant parameter across sites or a temporally changing variable? If it is a forcing variable for FATES, where is g_b from? is θ_ice a forcing variable or a parameter correlated with temperature and θ_liq? Is the C_a a constant value or variable? The model would be different if the spatial and temporal variability of all these factors are considered. If all these are parameters (i.e., scalars), what are the values?
3. Line 216-218: the reason for replacing saturated soil matric potential (Ψ_sat) with soil matric potential for closed stomata (Ψ_c) is not explained. Equation 10 shows that the Ψ_sat is replaced with soil matric potential for open stomata(Ψ_o), not Ψ_c. Furthermore, the Ψ_i was still calculated using Ψ_sat in Appendix A (equations A16-A18). I’m confused about which variable was used to calculate Ψ_i.
4. Line 218-220: is NN_B used to predict B_i or Ψ_i? B depends on only %clay and F_om according to equations A22-A23, while the authors add %sand, which is related to Ψ_sat and, therefore, Ψ_i. I didn’t find a direct relationship between B_i and %sand according to the original equations in the FATES model. If NN_B is used to predict Ψ_i, I think the equation can be Ψ_i=θ_liq*NN_B(%sand,%clay,PFT,F_om,T), where T represents the factors controlling θ_ice, e.g., temperature.
5. I think the neural networks (NN_B and NN_v) need constraints on V_cmax25 and Ψ_i. Although the authors declared that the predicted V_cmax25 without any constraints is within a rational range similar to the literature and measurement, the range of the predicted B is not discussed. If the predicted B_i is very large at some point, Ψ_i can be much higher than Ψ_o, leading to w_i being higher than 1 (i.e., exceeding the range defined in equation A15). Besides, the V_cmax25 is possibly to be inappropriate without any physical constraints at sites not considered in this study.
6. Line 235-236: ‘we tested retrieving both V_cmax25 and B, the latter of which varies spatially and temporally.’

If B varies temporally, it should be clarified how the training data is partitioned and how the ‘random holdout test’ is done. For example, is B changing per year or every N years? how many years/points per site are used to estimate B? Do the training points have to be in sequence or not?

1. Line 238-239: ‘For simplicity, the computations of B, Ψ_i, w_i, β_t were performed for the top soil layer only.’

In the synthetic experiment, only the top soil layer is considered. However, ‘B, Ψ_i, w_i‘ for the other layers are not clarified (=zero or default values in CLM?). Are the other soil layers considered in the real data experiment? If yes, how many ‘B’ was estimated (i.e., how many soil layers and how many temporally changing B_i)? If not, w_i can only represent the water availability at the top layer. The β_t is equal to w_i and the root distribution, r_i, at the top layer. What is r_i at the top layer (soil depth=0cm according to line 306)?

1. Line 239: ‘To retrieve B, we used NN_B but exclude the PFT term.’

I think it is not proper if the PFT is excluded from the training but included in the equation. If PFT is excluded, the term should be removed from equation 11. The sentence at line 222 ‘… along with categorical inputs (PFT), we used…’ should be rephrased.

1. Line 245: ‘The model passing the test of the synthetic case was then applied to a real dataset…’

The same NN was used for synthetic data and real data, but the NN information (layers, neurons activation functions) is not clear. As real data is much more complex, using a different NN structure from the synthetic test might have better performance.

1. Line 266-267: the loss function is very significant to evaluate the NN, but not explained in the paper. Without the equation of the loss function or the NN information, the dPL framework cannot be assessed by others, in other words, the experiment cannot be repeated. I think this doesn’t fulfil the requirement of Biogeosciences: ‘Is the description of experiments and calculations sufficiently complete and precise to allow their reproduction by fellow scientists (traceability of results)?’.
2. Line 268-272: the authors ‘hope to identify parameters that can generalize well in space’, so I think the readers would wonder if the parameters are estimated per site or per PFT. If parameters are estimated per site, how are they aggregated to parameters per PFT in figure 3a and 4a? If estimated per PFT, I’m afraid the spatial variability of the parameters is not fully captured by dPL.
3. Line 292-302: the sources of the soil moisture, stomatal conductance, meteorological forcings and the soil properties are mentioned, but the sources of C_a, g_b and P_atm are not clear.
4. The data source of ‘Lin15’ was not specified. I found a database at Lin et al., 2015, bud didn’t find the dates information on lines 296-300.
5. Line 304-305: the soil organic carbon content is collected, but the unit is not explained. Does the unit need to be transferred to get the soil organic matter fraction?
6. Line 410-411: the authors claim that the predicted V_cmax25 ‘were well correlated with’ literature values. However, the correlation coefficient or determination coefficient was never stated in the paper. Too few points are displayed in figure 6b, and the distribution pattern of only four PFT types (crop R, C3 grass, NET Boreal and BDS temperate) is similar to CLM.
7. Line 431-432: I cannot identify the C3 grass at the lower left corner of figure 5b. Maybe a violin plot per PFT can be helpful to show the difference between optimizing B or not for a specific plant type. The figures in the paper only show the net photosynthesis rate across all sites. However, the site-level comparison might be more meaningful to assess the four parameterization strategies: Vdef+B, Vdef+Bdef, V+Bdef, and V+B.
8. Line 445-450: I didn’t see any significant correlation between the estimated V_cmax25 and the PFT-mean from TRY database or other model default values. The authors should provide the scatter plots and the correlation coefficients between the V_cmax values to conclude that the dPL can get parameters correlated with literature values (line 490).

Minor comments:

1. Line 123: the right part looks very similar to the middle part in equation 3, but the subscript ‘W’ beside ‘argmin’ is not explained. As I understand, the ‘argmin’ in the right part is the same as the ‘argmin’ in the middle part.
2. Line 142: The short name for CO2 partial pressure at the leaf surface is ‘C_a’, but is ‘C_s’ in the appendix. Please use a uniform short name across the paper.
3. Line 187: equation 11 is cited at line 187 for the first time, but the full equation is placed at line 218. The equation should appear close to the first citation.
4. Line 193: does ‘i’ represent the soil layer number? I didn’t see the explanation around the equation.
5. Line 197: ‘across different soil different layers’ should be ‘across different soil layers’.
6. Line 203/equation 9: the second line should be T_i≤T_f-2 ‘or’ θ_liq≤0.
7. Line 205: the short name for the physical parameter at the second blue area should be θ but not θ_c.
8. Line 218/equation 11: B and F_om should have a subscript, i. B_i=NN_Bi(%sand,%clay,PFT,F_om,i).
9. Line 222: the ‘one-hot embedding’ was already stated at line 183. The definition should be explained where it is mentioned for the first time.
10. Line 228: the short name for ‘differentiable learning framework’ is defined but not used.
11. Line 310/Figure 2: the full names of the land cover types (e.g., BET tropical) are not explained before or around the figure.
12. Line 349: table 2 is mentioned for the first time, but the full table is placed after two pages.
13. Line 384: the CO2 should be CO_2(subscript).
14. Line 390/figure 5: I cannot understand the titles of the subplots. What is the meaning of ‘learning B’ and ‘learning V_cmax25’? The B is not optimized in figure 5a.
15. Line 514/equation A7: the Cs is not used.
16. Line 520-530: the three functions, Φ₁, Φ₂, and Φ₃, need to be clarified.
17. Appendix: the citations of equations are wrong (e.g, lines 503-504, 512, 520, 534…). The equations should be cited using A1-A23.