{"id":"860210d0-0f31-413a-8e88-e3eb1c5bb88c","arxiv_id":"2508.11154","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"A Gauss-Legendre quadrature discretization of the geothermal g-function integral equations, with regularization for stability, computes g-functions 20 to 200 times faster than the stacked finite line source model.","lead":"This paper introduces a faster numerical method for computing g-functions, the standard curves used to design geothermal borehole heating and cooling systems. The method uses Gauss-Legendre quadrature and regularization to cut computation times by 20 to 200 times versus the current stacked finite line source approach while keeping similar or better accuracy.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Accuracy comparison is anchored to the proposed method's own regularized high-order solution, so the claimed error advantage over SFLS is not yet established against the physical g-function.","rationale":"The paper's efficiency contribution is credible: replacing segment-to-segment integral response factors with analytical point-to-point factors is a genuine algorithmic improvement, and the measured speedups of 20-160x are consistent with the complexity reduction. The convergence trends shown in Section 3.1 are also internally coherent. The weakest part of the central claim is the accuracy comparison, because every error number is computed relative to the proposed method's own high-order regularized solution. That is a convergence self-check, not validation against the physical g-function. The reader's weakest_assumption identifies the fixed regularization parameter, which is part of the problem, but the more load-bearing issue is that the benchmark itself is the regularized proposed solution. Even a perfect lambda sensitivity study would not fix the circularity unless an independent reference is used. The proposed and SFLS solutions converge to nearby but not identical values, so the reported 'accuracy improvement' is sensitive to which solution is treated as truth. An independent FEM or SBM comparison is therefore the single check that would settle whether the claimed accuracy advantage is real. This concern does not require changing the reader's conditional verdict; it sharpens the condition that must be met before the paper's headline accuracy claim can be accepted.","tokens_in":19126,"tokens_out":6531,"duration_ms":73790,"concrete_test":"Run an independent finite-element (COMSOL) or SBM reference for the 2-borehole and 3x3 borehole cases with the same geometry and thermal parameters as Table 1 and the same uniform borehole wall temperature boundary condition, computing g(t) at the same ten time points. Then recompute the relative errors in Tables 12-14 for the proposed method and for the SFLS nonuniform-opt scheme against this external reference. If the proposed method is not consistently closer to the FEM/SBM g-function at matched discretization order, the central accuracy claim fails.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central accuracy claim is supported only by errors measured against a self-referential benchmark. Section 3.2 defines g_exact(tk) as the proposed method's own solution at Nheat_quad=1024, Ntemp_quad=50, and Eq. (42) uses this same solution as the reference for all error values in Tables 10 and 12-14. This reference is not an independent value of the UBWT g-function, and it is itself produced with the fixed regularization lambda=1e-8 from Eq. (41). Because Eq. (37) alters the Fredholm first-kind equation by adding an extra lambda*f term, the high-order solution used as 'exact' can carry a regularization bias whose magnitude is never bounded or studied. If that bias is comparable to the reported errors (roughly 1e-4), the claimed 'nearly 10 times improvement in accuracy' at Nheat_quad=512 could shrink or reverse when all methods are compared against an unregularized or external reference. The convergence tables do not resolve this: at N=1024 the proposed method gives g(20a)=7.2889 while the optimal nonuniform SFLS gives 7.2904 (Tables 6 and 9), so the choice of reference matters at exactly the scale of the accuracy claim.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper reformulates the computation of g-functions for vertical geothermal boreholes under the uniform borehole wall temperature (UBWT) boundary condition as a system of spatio-temporal integral equations. It shows that the stacked finite line source (SFLS) model is a stepwise-constant discretization of these integral equations, and proposes an alternative Nyström-type discretization using Gauss-Legendre quadrature in space and point-to-point analytical response factors. The authors further identify the governing equations as Fredholm integral equations of the first kind, explain the observed divergence at high discretization orders as ill-conditioning, and add a diagonal regularization term lambda*I to stabilize the linear systems. Numerical experiments on 2, 2x2, and 3x3 borehole fields are used to claim that the proposed method achieves comparable or better accuracy than optimized SFLS discretizations while being 20-200 times faster.","tokens_in":19482,"tokens_out":5016,"duration_ms":53525,"significance":"If the accuracy and efficiency claims are substantiated, the paper makes a useful contribution: it provides a principled quadrature-based discretization that avoids empirical segment selection, replaces segment-to-segment integrals with cheaper point-to-point response factors, and gives a plausible explanation for the previously observed divergence of g-functions at high discretization orders. The derivation from the Green's function to the discretized linear systems is coherent, and the complexity reduction is plausible from the stated O(N_b^2 N_heat_quad^2 N_temp_quad) structure. However, the central accuracy claims currently rest on a self-referential benchmark, and the regularization parameter is fixed without a sensitivity study. These issues must be addressed before the reported error reductions can be considered established.","major_comments":[{"comment":"The 'exact' g-function used as the reference in the error metric is the proposed method's own regularized solution at Nheat_quad=1024 and Ntemp_quad=50. This reference is not independent of the method being evaluated, and it inherits any bias introduced by the fixed regularization parameter lambda=1e-8. The convergence tables do not resolve the ambiguity: at t=20a and N=1024, the proposed method gives g=7.2889 (Table 6) while the optimal nonuniform SFLS method gives g=7.2904 (Table 9), a relative difference of about 2e-4, which is the same order as several of the errors reported in Table 12. Without comparison to an external reference (for example, a high-resolution finite-element or super-position-borehole-model solution, or an independently converged unregularized solution), the claim of 'nearly 10 times improvement in accuracy' is not established.","section":"Section 3.2, Eq. (42), Tables 10 and 12-14"},{"comment":"The regularization parameter lambda is fixed at 1e-8 for every numerical result, with no sensitivity study and no error bound relating the solution of the regularized system to the solution of the original Fredholm first-kind equation. Since Eq. (37) modifies the physical boundary condition to T_reg = T - lambda*Q, and since the condition number at N=1024 remains as large as 1.48e10 (Table 6), one cannot exclude a regularization bias at the 1e-4 level. The paper should report g-function values for a range of lambda (for example, 1e-12 to 1e-4) at each discretization order and demonstrate that the chosen lambda leaves the results unchanged within the claimed tolerance.","section":"Section 2.4, Eqs. (37) and (41)"},{"comment":"The scalar error defined in Eq. (42) measures the relative difference between sums of g-function values over all ten time points, rather than a pointwise or norm-based error. Errors with different signs at different time points can cancel in this aggregate, so Tables 10 and 12-14 may understate local inaccuracies. Reporting per-time-point relative errors or an L2/max norm over the time points would make the accuracy comparison more robust and would let the reader verify the claimed accuracy improvement at the scale of individual response times.","section":"Section 3.2, Eq. (42)"}],"minor_comments":[{"comment":"The first equality in Eq. (5) states that J(d,0,t2) equals the integral from t1 to t2, but the definition and the erfc expression correspond to an integral from 0 to t2; the integration limits should be corrected.","section":"Eq. (5)"},{"comment":"The segment bounds use z_{(i-1)j}, which is not a well-defined quantity for the i-th borehole's quadrature points; this appears to be a typo for z_{i(j-1)} or a similar index, and should be clarified so the segmentation scheme is unambiguous.","section":"Eqs. (32)-(33)"},{"comment":"The abstract and the introduction claim speedups of 20-200 times, but the reported numerical results show 20-160 times in Table 12 and 20-150 times in Table 14; the magnitudes should be reconciled or the claim should be stated as 'up to 20-160 times' based on the actual results.","section":"Abstract and Section 3.2"},{"comment":"There are several typos, including 'thay' in Section 3.1.1, 'seperate' in Section 2.3, 'segementation' in Section 3, and 'discretimzation' in Section 3.2, which should be corrected in a revised version.","section":"Throughout"},{"comment":"The column header in Table 7 reads 'Nheat_quad' for the SFLS model, but the row variable is the number of segments Ns; using a consistent header would avoid confusion.","section":"Table 7"},{"comment":"The paper does not include a data or code availability statement; since the timing comparisons depend on implementation details, a link to the implementation or a detailed description of the numerical libraries and hardware would strengthen reproducibility.","section":"Reproducibility"}],"recommendation":"major_revision","confidential_remarks":"The core derivation and efficiency argument are plausible, but the accuracy benchmark is self-referential and the regularization parameter is fixed without sensitivity analysis. I would be willing to accept after the authors compare against an independent reference and report a lambda-sensitivity study; both are feasible within the scope of the manuscript. The claimed novelty of being 'the first time' the convergence issue is discussed in an integral-equation framework should also be softened unless a more thorough literature check supports it."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here's the thing you should know: this paper gives a clean integral-equation reformulation of the UBWT g-function problem, and the 20-200x speedup it reports is credible. The accuracy claim, though, is anchored to a self-referential benchmark, so I wouldn't take the '10x more accurate' headline at face value yet.\n\nWhat's new: they write the g-function as Fredholm first-kind integral equations in space and time, then apply Nyström/Gauss-Legendre quadrature. SFLS falls out as a stepwise special case. The point-to-point response factors replace the expensive segment-to-segment integrals of pygfunction, which explains the speedup. The ill-conditioning diagnosis is also genuinely useful: condition numbers explode at high discretization order, and a small diagonal regularization (lambda=1e-8) stabilizes it. This explains why some high-order SFLS runs diverge, which is real and worth knowing.\n\nSoft spots, in proportion: the accuracy comparison is not yet established. Section 3.2 defines g_exact as the proposed method's own solution at Nheat=1024, Ntemp=50 with lambda=1e-8, and Eq. 42 measures everything against that. If you instead anchor errors to the SFLS-opt 1024 value (7.2904 vs their 7.2889), the proposed method's advantage at N=512 basically vanishes, and it can look worse. That's because the reference is on their trajectory, not on an independent physical answer. There is no sensitivity study for lambda, and no bound on regularization bias. Also, the divergence story is oversold: uniform SFLS at N=1024 has condition number 1116 and converges fine, so divergence is not generic to SFLS, it's specific to certain nonuniform schemes and their own high-order method. No code or data are released, so the numbers can't be reproduced.\n\nOverall: the method is worth engaging. The efficiency improvement is likely real and important for design optimization. The accuracy claim needs a proper external benchmark (FEM or legacy SBM) and a lambda sensitivity analysis. I'd send it to peer review, not desk reject, and would ask for those two additions before acceptance.","headline":"Genuinely faster g-function computation via Gauss-Legendre quadrature, but the accuracy improvement is measured against the method's own regularized high-order solution and needs independent validation.","tokens_in":19890,"tokens_out":3045,"would_cite":false,"duration_ms":31153,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65R20","45B05","80A19"],"pacs":[],"model":"deepseek-v4-flash","headline":"Geothermal borehole g-functions can be computed 20–200× faster by solving the underlying Fredholm integral equation with Gauss-Legendre quadrature.","keywords":["geothermal energy","borehole heat exchanger","g-function","discretization method","Gauss-Legendre quadrature","Fredholm integral equation","regularization","numerical simulation"],"falsifier":"Take the two-borehole test case at $N_{\\rm heat}^{\\rm quad}=512$ and compute the regularized g-function for $\\lambda$ ranging from $10^{-12}$ to $10^{-4}$; if the value moves by more than the quoted relative error of roughly $10^{-5}$ as $\\lambda$ varies, the accuracy claim depends on an unexamined parameter choice. Cross-check the regularized value against an independent, finely meshed finite-element simulation with the same uniform borehole wall temperature boundary condition.","tokens_in":18970,"feed_emoji":"♨️","tokens_out":14218,"duration_ms":127611,"temperature":0.7,"pith_summary":"This paper argues that the dimensionless thermal response function of a geothermal borehole field, the g-function used to design ground-source heat pumps and seasonal thermal storage, should be computed by solving the spatio-temporal integral equations that enforce a uniform borehole wall temperature rather than by the usual stacked finite line source (SFLS) approximation. The SFLS model is shown to be one special discretization of those equations, in which the heat extraction rate is taken constant within each depth segment. The paper replaces that stepwise approximation with Gauss-Legendre quadrature, so depth integrals become weighted sums at prescribed points and the costly segment-to-segment response integrals become analytic point-to-point factors; it also shows the discretized system is a Fredholm equation of the first kind that can become ill-conditioned enough at high order for the g-function to diverge, and that a small diagonal regularization restores convergence. If the numerical tests are representative, the method computes g-functions with comparable or better accuracy at 20 to 200 times the speed of optimized SFLS, reducing the tested 3×3 bore-field case from about 31 minutes to about 14 seconds.","feed_headline":"Gauss-Legendre rule cuts borehole g-function time 20–200x","feed_subtitle":"Replaces slow segment-to-segment integrals with analytic point-to-point factors and fixes high-order divergence.","key_machinery":"The load-bearing object is the Nyström discretization (a quadrature-based way of solving integral equations) applied to the governing equations: integrals over depth are replaced by quadrature sums, so the heat-extraction profile enters only through values at Gauss-Legendre nodes $z^{\\rm heat}_{ij}$ with weights $w^{\\rm heat}_{ij}$, and the wall-temperature constraint is evaluated at temperature quadrature points $z^{\\rm temp}_{mnl}$ with weights $w^{\\rm temp}_{mnl}$. The coefficient matrix is built from analytic point-to-point response factors $J(d,t_1,t_2)=\\frac{1}{4\\pi k d}[\\operatorname{erfc}(d/2\\sqrt{\\alpha t_2})-\\operatorname{erfc}(d/2\\sqrt{\\alpha t_1})]$, which replace the double-integral segment-to-segment factors of SFLS. The second load-bearing piece is the diagnosis that the governing equations form a Fredholm integral equation of the first kind, meaning the unknown heat extraction rate appears only inside the integral and the inverse problem is ill-posed; the paper's remedy is to add $\\lambda I$ to the coefficient block (Eq. 41), turning the unstable first-kind system into a stable regularized one.","core_discovery":"The central claim is that the g-function under the uniform borehole wall temperature boundary condition is governed by a closed system of integral equations (Eqs. 11–14) coupling the depth- and time-dependent heat extraction rate $Q_i(z,t_k)$ to the uniform wall temperature $T(t_k)$, with the total extraction rate normalized to 1. In this formulation the SFLS model is exactly the discretization that approximates $Q_i$ by a step function constant on each segment, and its segment-to-segment response factors are double integrals of the kernel $J(d,t_1,t_2)$. The proposed method instead applies a Nyström discretization: the depth integrals of heat extraction are replaced by weighted sums at Gauss-Legendre points, the average wall temperature is evaluated with a second Gauss-Legendre rule over temperature points, and the resulting linear system couples only point-to-point response factors $J(r(\\cdot),t_k-t_p,t_k-t_{p-1})$, which are analytic. Because the infinite-dimensional problem is a Fredholm integral equation of the first kind, the discretized matrix is ill-conditioned and its condition number grows rapidly with the number of unknowns, reaching about $10^{18}$ at 1024 points in the tested case; adding $\\lambda I$ with $\\lambda=10^{-8}$ to the coefficient matrix stabilizes the system. With regularization, the method converges at $N_{\\rm heat}^{\\rm quad}=512$ to a relative error of $2.8\\times 10^{-5}$, about an order of magnitude below the best nonuniform SFLS scheme tested, while taking 20–200 times less computation.","pith_inferences":["Editorial inference: a sweep of $\\lambda$ from $10^{-12}$ to $10^{-4}$ at $N_{\\rm heat}^{\\rm quad}=512$ would reveal whether the fixed $\\lambda=10^{-8}$ is load-bearing; if the solution is insensitive across that range, the accuracy claim is robust, and if not, the durable contribution is the ill-conditioning diagnosis rather than the specific speed numbers.","Editorial inference: the Fredholm-first-kind formulation invites applying standard inverse-problem parameter selection rules, such as the discrepancy principle or L-curve, to choose $\\lambda$ per time step and per bore field, which could make the regularization adaptive rather than fixed.","Editorial inference: the same Nyström-plus-regularization view may extend to other borehole boundary conditions, inclined boreholes, or boreholes with internal thermal resistance, and whether the first-kind ill-posedness survives those generalizations is a testable question.","Editorial inference: if the reported speed and accuracy hold in practice, on-the-fly g-function generation in optimization loops becomes feasible, replacing precomputed libraries as the standard way to size ground-source heat pump fields."],"forward_implications":["If the claim is right, fine-grained g-functions for large bore fields can be produced in seconds instead of minutes, making simulation-based design optimization and uncertainty quantification practical.","The divergence at high discretization order is explained by matrix conditioning rather than by a failure of the physical model, so users of fine SFLS discretizations now know when to mistrust the output and what to fix.","Gauss-Legendre point placement removes the need to hand-tune or empirically optimize segment lengths, so the accuracy claims transfer across borehole configurations without re-tuning.","At the highest tested resolution the quadrature method's relative error is about $2.8\\times 10^{-5}$, roughly ten times smaller than the optimized nonuniform SFLS value of $2.6\\times 10^{-4}$, so the accuracy advantage coexists with the speed advantage.","The linear system is built by pairing $N_b N_{\\rm heat}^{\\rm quad}$ heat points with $N_b N_{\\rm heat}^{\\rm quad} N_{\\rm temp}^{\\rm quad}$ temperature points, a structure whose cost grows quadratically with the number of unknowns but is still cheaper than SFLS because each factor is analytic."],"supporting_citations":[{"why":"Defines the g-function and the superposition borehole model whose uniform wall temperature boundary condition this paper adopts.","marker":"[4]"},{"why":"Introduces the stacked finite line source model under uniform wall temperature, the baseline identified as a stepwise special case and compared against.","marker":"[21]"},{"why":"Documents the quadratic cost of segment-pair response factors and the need for finer discretization as borehole count grows.","marker":"[30]"},{"why":"Provides the optimized nonuniform discretization schemes used as the accuracy and speed comparison targets.","marker":"[32]"},{"why":"Supplies the point-heat-source Green's function kernel used in the integral equations.","marker":"[33]"},{"why":"Reduces segment-to-segment response factors to a single integral, the object the proposed point-to-point factors replace.","marker":"[34]"},{"why":"Supplies the Nyström method for solving integral equations by quadrature, the basis of the proposed discretization.","marker":"[35]"},{"why":"States the ill-posed nature of Fredholm integral equations of the first kind, motivating the divergence diagnosis.","marker":"[36]"},{"why":"Supplies the regularization technique of adding a small term to convert the first-kind equation to a stable second-kind form.","marker":"[37]"}],"fun_headline_variants":["Gauss-Legendre rule speeds borehole g-functions 20–200x","Integral reformulation makes borehole g-functions 200x faster","Point-to-point factors cut borehole g-function time 20–200x","Regularized Nyström method handles borehole g-function divergence"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"All numerical results rest on a single fixed regularization parameter $\\lambda=10^{-8}$ being simultaneously small enough to leave the physical g-function unchanged and large enough to suppress the ill-conditioning, but the paper gives no sensitivity study or error bound connecting the regularized solution to the true solution of the Fredholm equation.","fun_headline_variants_meta":{"raw":{"variants":["Gauss-Legendre rule speeds borehole g-functions 20–200x","Integral reformulation makes borehole g-functions 200x faster","Point-to-point factors cut borehole g-function time 20–200x","Regularized Nyström method handles borehole g-function divergence"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000375,"raw_usage":{"total_tokens":2089,"prompt_tokens":1124,"completion_tokens":965,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":740,"completion_tokens_details":{"reasoning_tokens":883}},"tokens_in":740,"tokens_out":965,"duration_ms":9198,"temperature":1.0,"reasoning_tokens":883,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T17:27:21.946895+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take the two-borehole test case at $N_{\\rm heat}^{\\rm quad}=512$ and compute the regularized g-function for $\\lambda$ ranging from $10^{-12}$ to $10^{-4}$; if the value moves by more than the quoted relative error of roughly $10^{-5}$ as $\\lambda$ varies, the accuracy claim depends on an unexamined parameter choice. Cross-check the regularized value against an independent, finely meshed finite-element simulation with the same uniform borehole wall temperature boundary condition.","supporting_citations":[{"cited_title":"Eskilson, Thermal analysis of heat extraction boreholes (1987)","cited_arxiv_id":null,"evidence_quote":"Defines the g-function and the superposition borehole model whose uniform wall temperature boundary condition this paper adopts."},{"cited_title":"Cimmino, M","cited_arxiv_id":null,"evidence_quote":"Introduces the stacked finite line source model under uniform wall temperature, the baseline identified as a stepwise special case and compared against."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Documents the quadratic cost of segment-pair response factors and the need for finer discretization as borehole count grows."},{"cited_title":"Cimmino, J","cited_arxiv_id":null,"evidence_quote":"Provides the optimized nonuniform discretization schemes used as the accuracy and speed comparison targets."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the point-heat-source Green's function kernel used in the integral equations."},{"cited_title":"Claesson, S","cited_arxiv_id":null,"evidence_quote":"Reduces segment-to-segment response factors to a single integral, the object the proposed point-to-point factors replace."},{"cited_title":"Hackbusch, Integral equations: theory and numerical treatment, Vol","cited_arxiv_id":null,"evidence_quote":"Supplies the Nyström method for solving integral equations by quadrature, the basis of the proposed discretization."},{"cited_title":"Wazwaz, Linear and nonlinear integral equations, Vol","cited_arxiv_id":null,"evidence_quote":"States the ill-posed nature of Fredholm integral equations of the first kind, motivating the divergence diagnosis."},{"cited_title":"Wazwaz, The regularization method for fredholm integral equations of the first kind, Computers & Mathematics with Applications 61 (10) (2011) 2981–2986","cited_arxiv_id":null,"evidence_quote":"Supplies the regularization technique of adding a small term to convert the first-kind equation to a stable second-kind form."}],"review_version":2}