Wortmann et al. present the Earth Science Box Modeling Toolkit (ESBMTK), which is a Python library designed for building and analyzing box models in Earth science. It uses a modular, object-oriented approach to study topics like the long-term carbon cycle and the impact of atmospheric CO2 changes on seawater chemistry. ESBMTK separates model geometry from the underlying numerical implementation, and thus allows users to focus on the conceptual challenges, rather than mathematical theory. Such a tool is very useful for teaching and research requiring fast conceptual models. In addition to predefined setup, the user can customize rather easily his/her own model such as the number and volume of boxes/reservoirs, the flux between them, the isotope species, ... This tool will be very useful for the climate community, and I therefore recommend publication after minor revisions.

General remarks

• Several species are already defined in species_definition.py such as stable water isotope ²H and ¹⁸O, ¹³C, ...  I would recommend to add the 14C because it is quite usual to use ocean box modeling to model Δ¹⁴C, especially in the framework of IntCal. See for example Bard et al. (1997, https://doi.org/10.1016/S0012-821X(97)00082-4).
• In my opinion, it would be beneficial to give more details on the modeling of stable water isotopes ²H and ¹⁸O, and not only on the carbon-related species. The transport between the ocean box models are quite simple in that case (one to one), but the authors could mention how to set up a simulation when considering fractionation effect between atmosphere and ocean boxes e.g, the evaporation from the ocean to the atmosphere.
• The authors discuss the ESBMTK results with Boudreau et al. (2010a) model. They show the good ability of ESBMTK to replicate other models, which is important to encourage potential users to switch to ESBMTK. One other important aspect is to show how realistic model results are. Is there any way for the authors to compare the ESBMTK results with observations for a typical simulation (or with the setup of Boudreau et al. (2010a))?

Minor and technical comments

• Line 49: a point is missing.
• Legend of Figure 1: remove ". caption".
• Lines 62-63: In the following, we only describe the pertinent implementation details.
• Table 2: I would show the information for the all the fluxes, and in the order F1, F2, F3, ..., F9.
• Legend of Table 3: all parameters are after Boudreau et al. (2010a).

This manuscript introduces the Earth Science Box Modeling Toolkit (ESBMTK), a Python library designed to streamline the build and operation of box models. The toolkit incorporates many commonly used processes such as air-sea exchange, marine carbonate chemistry, and isotope calculations. Its feasibility and robustness are validated through replicating the work by Boudreau et al., 2010.

Given the significant utility of box models in the Earth Science, ESBMTK has considerable potential for applications in both teaching and research settings. Below are several suggestions that, hopefully, could help further enhance the quality and impact of this project:

1. The toolkit’s ability to separate model geometry from its underlying numerical implementation is highly convenient, particularly for users without extensive programming experience. However, this abstraction also hides the core functions, making debugging more challenging. One potential enhancement could be enabling the model to output the codes for governing equation it generates during runtime. Such a feature would not only improve the model transparency but also serve as a valuable resource for educational purposes, allowing users to better understand the mechanics of the model.
2. Several processes and their parameterizations are described in Section 2. However, I found it somewhat challenging to clearly link these processes to the toolkit framework as depicted in Figure 1. It would be helpful if the authors could integrate these processes explicitly within the framework illustration or provide some example codes.
3. I believe perturbation experiments are crucial for exploring system behavior under various conditions. I was happy to learn from the documentation that the model accepts external forcings. I suggest highlighting this feature more prominently and giving brief introduction on its implementation in the main text
4. Several things might need further clarifications.
   • Line 99, the authors state that “the above parameters at the beginning of each run and assumes that they are constant over the integration interval.” However, it is unclear what “the above parameters” specifically refer to. My interpretation is that this might relate to thermodynamic and kinetic constants. However, the carbonate system is typically time-dependent and calculated in each integration.
   • Again, from the same paragraph, the authors state that modeling temperatures remain constant throughout the runtime. However, this seems inconsistent with examples where transient and long-term carbon cycles are simulated. Since temperature can affect thermodynamic constants in the carbon system, this assumption may not be realistic. If my understanding is incorrect, additional explanation would be appreciated to resolve this confusion.
   • Line 105, It is noted that [H+] is initially calculated using thepyco2sys library, and subsequently, the iterative approach of Follows et al. (2006) is applied. I’m confused by the rationale behind using two different methods for this calculation.
5. Carbon isotopes are crucial tracers in carbon cycle modeling. In Example 4 (https://github.com/uliw/ESBMTK-Examples/blob/main/Examples_from_the_manual/po4_4.png), the steady-state surface box δ¹³C reaches approximately 8‰, which is not realistic. Since the Boudreau et al. (2010) model does not incorporate isotope modeling, I recommend that the authors conduct additional validation work for this component.