The manuscript introduces the boundary integral method (BIM) to find equilibrium solutions, and their corresponding stability properties, in North-type energy balance models (EBMs). I fully support the use of EBMs to understand and characterize the main features of climate states under different conditions. The manuscript is well written and it is surely appropriate for NPG. I have some moderate suggestions to compare the results of this paper with previously obtained ones in a similar framework. 

Main comments

1. An interesting result is the "striped" pattern obtained in Figure 3. This pattern has been also reported in 1D EBMs including vegetation (e.g., Nevison et al., 1999; Adams et al., 2003; Alberti et al., 2015). I would suggest the authors to compare their results with those papers in terms of active feedbacks, especially those related to the ice-albedo feedback.
2. The authors consider a step-wise function for the albedo with a latitude dependence. However, they do not include in the feedback the extension of the continent (e.g., Wood et al., 2008; Rombouts and Ghil, 2015; Alberti et al., 2018). I would recommend to comment on this point and how the BIM can be generalized towards including additional contributions in the ice-albedo feedback.
3. I found the results on the larger number of equilibria very interesting. I would suggest to include more discussion on the bifurcation diagrams reported in Figure 2. What about inspecting how the extension of the stable and unstable regimes depend on ε? This would suggest and provide more information on the residence time in a specific state for the model that can be useful for a broader discussion on the role of ice-ocean-land distribution on the climate states. 

Minor comments

• Line 38: the values are slightly different from those usually reported in literature (Budyko, 1969; Rombouts and Ghil, 2015). Are these values corrected for current levels of greenhouse effect?
• Line 60: is in the definition of β a missing T_s or Eq. (3) is just written in dimensionless temperature (T → T/T_s)? Please clarify.
• Eq. (8): please clarify that this is valid in a linear approximation around an equilibrium solution.
• Eqs. (11)-(12): should these conditions valid for each selection of the ice-albedo feedback? Are there for boundary conditions at the poles, right?
• Line 144: why choosing l = π/4? Please also clarify that in the ice-albedo feedback there is no ice-ocean-land feedback. 
• Figure 3: I would recommend to increase the quality of the figure since it seems that θ labels are not fully resolved.

Suggested references

• Adams, B., et al. (2003) J. Theor. Biol., 223, 505.
• Adams, B., and Carr, J. (2003) Nonlinearity, 16, 1339.
• Alberti, T., et al. (2015) Phys. Rev. E, 92, 052717.
• Alberti, T., et al. (2017) ApJ, 844, 19.
• Alberti, T., et al. (2018) J. Phys. Comm., 2, 065018.
• Nevison, C., et al. (1999) Tellus B, 53, 288.
• Rombouts, J., and Ghil, M. (2015) NPG, 22, 275.
• Wood, A.J., et al. (2008) Rev. Geophys., 46, RG1001.

General

In the present manuscript the authors apply the boundary integral method (BIM) to find equilibrium  solutions of a North-type one dimensional energy balance model. The proposed method is described in detail and applied using some case studies concerning the effect of a continent.  The study demonstrates that BIM is a valuable method and may be applicable to conceptual models of similar type, which can still contribute significantly to our process understanding. The results of the showcase experiments indicate the existence a variety of equilibrium solutions depending on the location of the continent.

The manuscript is well written. The methodology is sound and sufficiently explained together with an appropriate application. Thus, in general, I can recommend publication. However, the authors may like to add some more discussion of the results (see Specific).

Specific

While I am happy with the methodology, there is (unfortunately) very little discussion of the (interesting) results. One may wonder why the authors perform (sort of) a sensitivity study with different continental set ups only to illustrate that the method is able to find equilibrium solutions. For this, one example would have been sufficient. In my view, a more thorough discussion of the results would add significant to the value of this paper.

Minor

1) I may have overlooked it, but the authors may state which root searching algorithm they have used.

2) It seems that starting from section 2 the variable T is non-dimensional (T/Ts)) which may be confusing given the dimensional T in Eq. (1). The authors may use different symbols for dimensional and non-dimensional variables, respectively.

3) All equations should be numbered.

This interesting paper investigates the steady states of a one-dimensional energy balance model with meridional heat diffusion using a so-called boundary integral method (BIM). The novelty lies in the application of the method to the problem at hand, which allows the treatment of arbitrary discontinuity, and which allows the authors to discover multistability in the system as well as a variety of unstable steady states with several ice edges and a continent.

The paper is well-written and the results are presented clearly. I recommend publication after a few minor comments are addressed.

1. At the beginning of Section 2, it may be helpful for readers unfamiliar with the method to explicitly introduce what constitutes the boundary integral method. Either in general terms, or in the context of the presented equations.

2. For extra claricty, in Section 2.1, perhaps explicitly mention how the solutions in this work are constructed. I.e., first make an "ansatz" for the albedo function and then solve the previously obtained relations in the relevant latitudinal domains.

3. In the Discussion, it could be helpful to put the problem of multistability into a broader context. It seems like BIM is a good tool to determine steady states analytically, and especially the possibility to find a variety of unstable states is worth mentioning. But I am wondering how the method fairs in the presence of a high number of co-existing states? It is not so clear whether Fig. 3 constitutes a "complete" bifurcation diagram, or just a continuation of several (out of potentially many more) steady states that have been found via a set of boundary conditions. Is it possible that the model shows a high degree of multistability, akin to a "fragmented" or "hierarchical" stability landscape, as has been found in other climate models (Bastiaansen et al 2022, Lohmann et al 2024)? In this case, can the BIM with its reliance on prescribing an "ansatz" for the albedo function help, or does one need to resort to numerical simulation in order to find the "global" stability landscape?

References:

Bastiaansen et al., "Fragmented tipping in a spatially heterogeneous world", Environ. Res. Lett. 17 045006 (2022)

Lohmann et al., "Multistability and intermediate tipping of the Atlantic Ocean circulation", Sci. Adv. 10, eadi4253 (2024)

This article is concerned with zonal energy balance models (EBMs) for Earth’s climate system.  EBMs were first considered by Budyko and Sellers in the 1960s and studied extensively by North and coworkers in the 1970s. Needless to say, much analytical and computational work has been done since then to extend the original results to more complicated configurations.  The purported novelty of the article is the application of boundary integral methods. 

The configuration is kept relatively simple.  The Earth is modeled as a solid sphere, symmetric with respect to the equatorial plane, with one or more zonal regions and latitudinal energy transfer (diffusion). The authors focus on steady state solutions of Eq. (1).  The nondimensional form is given in Eq. (3).  The equations are linear, but they allow for ice-albedo feedback.

The authors construct a Green’s function for the perturbed diffusion operator in terms of Legendre polynomials.  The Green’s function is then used to produce solutions on the relevant subintervals (polar ice caps, land mass, and equatorial regions) and patched together to find a global representation.  The main results are given in Figure 2 (bifurcation diagrams) for three slightly different configurations (distinguished by the value of the parameter \epsilon). It shows that the solution of an EBM can be highly sensitive to the parameter values. Also, the bifurcation diagrams show some unexpected unstable branches.

The article is self-contained and would be a nice contribution to the literature on EBMs.  The problem is that EBMs have been around for so long that a full-blown introduction with a slew of references seems like overkill.  (I missed references to some articles by Ester Widiashi et al., who focused on the ice lines.)  So, I recommend that the authors cut most of the introduction and give us their essential findings, leaving the technical details to the reader. After all, the problem is linear, and most readers are (or ought to be) familiar with integration by parts and linear stability analysis.

In summary: I am somewhat hesitant to recommend publication in its present form. A short summary of the new elements and findings would be acceptable.

This paper proposes the analytical method, the boundary integral method (BIM), to solve a type of simplified climate model known as a North-type energy balance model (EBM). The authors demonstrate the effectiveness of the BIM by solving the EBM with various configurations, including a zonally symmetric continent with an altered ice-albedo feedback mechanism. They show that the BIM can handle complex scenarios with multiple ice edges and ice belts, making it a powerful tool for understanding the dynamics of these simplified climate models.

The argument is surely relevant and the methodology is well-explained, so that I recommend the publication after minor corrections. In the following, I provide some suggestions for further improving the manuscript:

1. It would be useful to provide some examples of the application and efficacy of such methodology in different research fields (e.g. how it works in the cited paleoclimatology field).
2. Readers would appreciate knowing all the limitations of the North-type EBMs as well as the advantages.
3. In my opinion, the authors should broaden the conclusions by summarizing the key points and further highlighting the importance of North-type EBMs.