Comment on “Groundwater Affects the Geomorphic and Hydrologic Properties of Coevolved Landscapes” by Litwin et al.

The objective of this comment is to correct two sets of statements in Litwin et al. (2022, https://doi.org/10.1029/2021JF006239), which consider our research work (Bonetti et al., 2018, https://doi.org/10.1098/rspa.2017.0693; Bonetti et al., 2020, https://doi.org/10.1073/pnas.1911817117). We clarify here that (a) the specific contributing area is defined in the limit of an infinitesimal contour length instead of the product of a reference contour width (Bonetti et al., 2018, https://doi.org/10.1098/rspa.2017.0693), and (b) not all solutions obtained from the minimalist landscape evolution model of Bonetti et al. (2020, https://doi.org/10.1073/pnas.1911817117) are rescaled copies of each other. We take this opportunity to demonstrate that the boundary conditions impact the obtained solutions, which has not been considered in the dimensional analysis of Litwin et al. (2022, https://doi.org/10.1029/2021JF006239). We clarify this point by using dimensional analysis and numerical simulations for a square domain, where only one horizontal length scale (the side length l) enters the physical law.

wrote a very interesting paper on the problem of the co-evolution of groundwater hydrology and land-surface dynamics. The authors are to be commended for this very nice research contribution, which quantifies the feedback between the spatial patterns in the emerging topography and sub-surface properties.
This comment is intended only to rectify some statements in this paper that regard our work (Bonetti et al., 2018. The first one considers the definition of the specific contributing area a and its relation to the total contributing area A. The second one is related to the dimensional analysis and the scaling properties of the landscape evolution model (LEM) solutions. We feel that these corrections are necessary to avoid misinterpretations of our results as well as to clarify important aspects regarding the scaling properties of the obtained solutions.
The statements in question concern the LEM with governing equations Abstract The objective of this comment is to correct two sets of statements in Litwin et al. (2022, https://doi.org/10.1029/2021JF006239), which consider our research work (Bonetti et al., 2018(Bonetti et al., , https://doi. org/10.1098(Bonetti et al., /rspa.2017Bonetti et al., 2020Bonetti et al., , https://doi.org/10.1073. We clarify here that (a) the specific contributing area is defined in the limit of an infinitesimal contour length instead of the product of a reference contour width (Bonetti et al., 2018, https where z is the topographic elevation, t is time, D is a diffusion coefficient, K is an erosion coefficient, a is the specific contributing area, and U is the uplift rate. For constant K, the model-referred to as the NoHyd modelis a specific version of the minimalist LEM of Bonetti et al. (2020).
Regarding the definition of the specific contributing area, Litwin et al. (2022) write on page 5: "To make the conversion between A and a, we represent A as the product of a and a characteristic contour width v 0 , which is a chosen constant value." In reality, a is defined in the limit of an infinitesimal contour length (Bonetti et al., 2018;Gallant & Hutchinson, 2011), not as the product of a reference contour width. As a result, Equation 2 is only valid for the specific contributing area defined in this limit, as discussed in Bonetti et al. (2018).
With reference to similarity solutions and dimensional analysis, they write on page 4, "Additionally, the nondimensionalization generalizes our results and reconciles conflicting dimensional analyses provided by Theodoratos et al. (2018) and Bonetti et al. (2020)." Further, they mention on page 27, "We show that contrary to Bonetti et al. (2020) there is a single typology of the NoHyd model which can be rescaled to obtain all results the model may produce." These conclusions need to be revised because the role of the boundary conditions is not included in the dimensional analysis. While the dimensional analysis of the LEM is facilitated by the fact that some of the main variables and parameters are clearly listed in the governing equations, the obtained solutions depend on the initial and boundary conditions as much as they do on the equations (Bursten, 2021).
Considering the domain geometries for which only one horizontal length scale enters the physical law (more complicated geometries will add further governing quantities), the LEM solution can be formally written as For concreteness, we consider here a square domain of side length l with boundaries at fixed zero elevation values that is, where  is the "channelization index" derived in Bonetti et al. (2020).

of 4
Using the above-defined scales  ,  , and  , the non-dimensionalized form of the LEM (governing equations:

Equation 1 and Equation 2; boundary conditions: Equation 4 and Equation 5) can be written as
The absence of Π in the non-dimensionalized governing Equations 12 and 13 does not mean that the solutions for this model do not depend on it. Π l appears in the non-dimensionalized boundary conditions (Equation 14 and Equation 15) used to compute the steady-state solutions. While there is freedom in selecting the repeating variables, choosing one set over another for defining the relevant scales does not change the physics of the model at hand (Porporato, 2022). Selecting l, D, and U keeps Π l in the dimensionless governing equations , while Π l appears in the dimensionless boundary conditions when using K, D, and U for the non-dimensionalization.

of 4
The dependency of the non-dimensional solutions on the boundary conditions through Π l signifies that not all the solutions that the model can produce are rescaled copies of each other. This is true only for those solutions in which the non-dimensional group Π l remains invariant. A set of steady-state simulations for a square domain, represented in both dimensional and non-dimensional forms, demonstrates this point (Figure 1). In panels a, b, and c, the value of Π l is changed by increasing the erosion coefficient K for a domain of fixed side length l = 100 m. This is the case discussed in Bonetti et al. (2020), which shows different solutions ranging from the unchannelized case ( Figure 1a) to a case with multiple channels (Figure 1c). The non-dimensional solutions ( ) computed by employing suitable scales  and  are also dissimilar in these cases. On the contrary, panels d, b, and e have different values of K, but they preserve the same Π l because of the different length scale l, and therefore, are scaled versions of each other (similar non-dimensional forms labeled in the blue text). Theodoratos et al. (2018) consider this case only.
In reality, the reason why many solutions of LEMs look alike is that they have an interesting property of complete self-similarity that emerges at very large channelization index values (  = 3∕2 ) Hooshyar, Bonetti, et al., 2020;Hooshyar et al., 2021;Porporato, 2022). For these conditions, several spatially averaged quantities become invariant with respect to  , having reached a self-similar regime. This is similar to what happens to turbulent flows in the fully rough regime as they become independent of the Reynolds number or to regular polygons, which tends to circles as the number of sides tends to infinity (Barenblatt, 1996).

Data Availability Statement
The details of the numerical algorithm used for the simulations is described in Anand et al. (2020) and the well-commented Python code is available on GitHub https://github.com/ShashankAnand1996/LEM (Anand, 2022). Moore Science-to-Action Fund, and BP through the Carbon Mitigation Initiative (CMI) at Princeton University. S.K.A. and A.P. also acknowledge the support from the High Meadows Environmental Institute (HMEI). The simulations in this article were performed on computational resources provided by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.