{"id":"b38eb2fa-f5e2-4691-a8e3-64ecfdb02c8b","arxiv_id":"2509.26284","paper_version":3,"verdict":"ACCEPT","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"For the six-equation single-velocity two-phase model, different discretizations of non-conservative products converge to different shock solutions, and the Bassi-Rebay treatment fits the path-conservative framework.","lead":"Numerical schemes for the six-equation two-phase flow model are compared, showing that different treatments of non-conservative terms converge to different shock solutions at extreme pressure ratios. The paper also shows that a Bassi-Rebay treatment of those terms fits the path-conservative framework, supporting a robust HLLC-based strategy.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Missing mesh-convergence measurement undermines the claim that the inter-scheme spread is a converged property of the model; the fine-mesh comparisons are visual only.","rationale":"The reader's weakest_assumption identified the same issue, so agreement is 'agree'. The concern is load-bearing because the paper's main novelty claim—distinct schemes converge to different shock solutions—contains the word converge. The only support offered is fine-mesh profiles and a qualitative statement of mesh convergence. The proposed refinement sequence directly tests whether the fine mesh is in the convergent regime and whether inter-scheme differences persist or shrink, which is the crux. I recommend CONDITIONAL rather than REJECT because the theoretical analysis of non-unique jump conditions is coherent and the numerical evidence is suggestive; the missing piece is a quantitative convergence check. If the check passes, ACCEPT is appropriate; if not, the central claim should be weakened to an observation about finite-resolution behavior.","tokens_in":27390,"tokens_out":7377,"duration_ms":64082,"concrete_test":"Reproduce the epoxy-spinel case (Tables 1–3) at p_L/p_R = 10^6 with the three schemes used in Fig. 12 (HLLC wave-propagation, HLLC+BR-2023, Rusanov+BR-2023) on N = 8192, 16384, 32768, 65536, 131072. For each run, record the intermediate-state volume fraction α1 (e.g., the plateau value behind the shock) and compute (i) the intra-scheme differences |α1(N) − α1(2N)| and (ii) the maximum inter-scheme spread at each N. If (i) → 0 while (ii) stays bounded away from 0 as N doubles, the central claim is confirmed; if (ii) also → 0, the reported spread is an under-resolution artifact and Section 7.3's conclusion must be revised.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central claim is that at extreme pressure ratios different discretizations of the 6-equation hyperbolic operator converge to different shock solutions after instantaneous mechanical relaxation (Sections 1, 4.3, 7.3). The word 'converge' is load-bearing, and the paper explicitly asserts that mesh convergence is established: 'we employ the fine mesh, so as to achieve mesh convergence' (Section 4) and 'we have achieved mesh convergence' (Section 8). Yet no mesh-refinement sequence, error norm, or convergence rate is reported; all fine-mesh comparisons are made by visual inspection of profiles (e.g., Figures 6, 8, 11, 12). For a non-conservative system, this is not a cosmetic gap: first-order shock computations have O(h) numerical layers, and different numerical-viscosity models can have different convergence paths toward possibly different weak limits. Without measuring that each scheme's own solution is stationary as h → 0, the spread in Fig. 12 could be an under-resolution artifact rather than a property of the non-unique jump conditions. The conclusion 'different numerical schemes may tend towards different numerical solutions' is exactly what requires this measurement.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper studies a six-equation single-velocity two-phase flow model with instantaneous mechanical relaxation, used as a numerical path to the five-equation Kapila model. The authors compare several first-order finite-volume discretizations of the homogeneous hyperbolic operator, including an HLLC wave-propagation scheme, Rusanov/HLLC fluxes combined with Bassi–Rebay or Crouzet-type treatments of non-conservative products. A central contribution is an algebraic equivalence result in Section 3.2 showing that two HLLC energy closures used in the literature are different discrete realizations of the same approximate jump condition, and a formal path-conservative reformulation of the BR approach in Section 5. Numerical tests cover homogeneous problems and cases with instantaneous mechanical relaxation. The paper's main claim is that, while most test cases give nearly identical results after relaxation, at extreme pressure ratios different schemes converge to different solutions, reflecting the absence of uniquely defined jump conditions for the model.","tokens_in":27641,"tokens_out":2346,"duration_ms":22050,"significance":"If the central claim is supported, the paper makes an important and cautionary contribution: shock solutions of the six-equation model, and hence some solutions of the Kapila model, are discretization-dependent at extreme pressure ratios even when global conservation is enforced. The algebraic analysis in Section 3.2 is a genuine clarification of the literature, and the path-conservative interpretation of BR-type discretizations in Section 5 is a useful conceptual bridge. The numerical study covers a wide range of benchmarks, including cavitation, strong shocks, and a 2D problem, and uses open-source code and no fitted parameters, which are methodological strengths. However, the paper's principal conclusion is explicitly tied to a claim of mesh convergence that is asserted but never quantitatively demonstrated, and this gap is load-bearing because the model is non-conservative and different schemes may converge to different weak limits.","major_comments":[{"comment":"The paper repeatedly asserts mesh convergence ('we employ the fine mesh, so as to achieve mesh convergence'; 'we have achieved mesh convergence'), but no mesh-refinement sequence, error norm, or convergence rate is reported anywhere. All fine-mesh comparisons are visual (e.g., Figs. 6, 8, 11, 12). This matters for the central claim of Section 1 and the conclusions: for a non-conservative system, first-order schemes have O(h) numerical layers, and different numerical viscosity models can have different convergence paths toward possibly different weak limits. Without measuring that each scheme's own solution is stationary as h→0, the spread in Fig. 12 could be an under-resolution artifact rather than a property of the non-unique jump conditions. I request a quantitative convergence study — e.g., L1 or L∞ errors of key variables (volume fraction, velocity, phasic pressures) for at least N =","section":"Section 4 (introductory paragraph), Section 8"},{"comment":"The phrasing 'different numerical schemes converge to different solutions' is stronger than what the evidence shows. In Section 4.3 the fine-mesh results still show visible differences, but the text only says different schemes 'can converge to different solutions.' The conclusion (Section 8) correctly softens to 'may tend towards different numerical solutions.' The central claim should be stated consistently in its falsifiable form, and the proposed convergence study should be used to decide whether 'converge' is warranted. Without that measurement, the unique-solution claim is not established.","section":"Section 1 and Section 4.3"},{"comment":"The derivation linking the two HLLC closures is algebraically sound and a strength of the paper. However, the claim just before Eq. (30) that assuming (28) 'owing to (30) is now shown to also yield (24)' should be made explicit: the step from the sum of the two bracketed terms being zero and the second bracket being zero to the first bracket being zero is valid only if the prefactor ρk(u−S) is nonzero. For a stationary contact or a transonic configuration this degenerates. The authors should state this nondegeneracy condition or note that the closure is intended for genuinely supersonic shock waves.","section":"Section 3.2, Eqs. (23)–(30)"}],"minor_comments":[{"comment":"Typo: 'well-know' should be 'well-known'.","section":"Section 1"},{"comment":"The notation for the jump identity is confusing: the formula JϕψK = {{ϕ}}JψK + {{ψ}}JϕK is correct, but the surrounding text uses J·K with inconsistent subscripts (e.g., 'J ˜mkuK { {u} }' should be 'J ˜mkuK {{u}}' — check the placement of the average braces). Minor typographical issue.","section":"Section 3.2, display before Eq. (30)"},{"comment":"The comparison with the analytical Euler solution (Fig. 7) is interesting but should be qualified: the 6-equation model with finite-rate (here zero) relaxation is not identical to the Euler equations, even for quasi-pure phases, as the phasic energy equations and non-conservative terms change the shock structure. The text notes this informally, but a sentence making the caveat explicit would avoid a misleading reading.","section":"Section 4.3"},{"comment":"The reformulation of BR as a path-conservative scheme is a nice contribution, but the claim that this is 'the first rigorous analysis' linking these approaches would benefit from a more precise literature comparison: many works use double integration by parts for non-conservative terms (the paper cites some), and the novelty is the trapezoidal-rule path interpretation. A clearer statement of what exactly is new would help.","section":"Section 5"}],"recommendation":"major_revision","confidential_remarks":"The paper is technically solid in its derivations and the numerical comparison is extensive, but the central claim—that at extreme pressure ratios different discretizations converge to different solutions—hinges on an asserted but unmeasured mesh convergence. The missing quantitative convergence study is fixable within the manuscript's scope and would substantially strengthen the paper. I therefore recommend major revision rather than rejection."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Bottom line: the paper is genuinely useful, but its central claim is stronger than its evidence. The algebraic work is real; the mesh-convergence claim is asserted, not measured.\n\nWhat is new: Section 3.2 shows the two HLLC closures used in [60] and [51] are two discretizations of the same approximate jump condition (24). That is clean and worth having. Section 5 makes the case that the Bassi-Rebay treatment of non-conservative products is a path-conservative scheme with segment paths and a trapezoidal rule; that explains why the BR variants behave consistently and why the [19] approach does not fit the framework. The systematic comparison across Rusanov, HLLC, wave-propagation, and multiple non-conservative-product treatments, both for the homogeneous model and with instantaneous relaxation, is exactly the kind of test practitioners need. The authors are also honest about the model's lack of unique jump conditions and about the fact that the observed spread only shows up at extreme pressure ratios.\n\nThe soft spot is the one the stress-test flags. The paper says it uses the fine mesh 'so as to achieve mesh convergence' and later 'we have achieved mesh convergence,' but no mesh-refinement sequence, error norm, or convergence rate is reported. The comparisons are visual. For a non-conservative system this is not cosmetic: first-order shock computations have O(h) numerical layers, and different numerical-viscosity models can be at different stages of convergence at 65,536 cells. The observation that different discretizations do not settle on the same solution is worth publishing, but without measuring each scheme's own convergence as h goes to zero, the spread in the epoxy-spinel figure could be partly an under-resolution artifact. The authors' comment that 65,536 cells is 'unfeasible for practical applications' does not fix this; the claim is about the numerical limit, not about practice.\n\nA smaller practical gap: no pinned code, input files, or output data are provided. The samurai framework is open, but exact setups would be hard to reproduce.\n\nWho this is for: people computing compressible two-phase flows with relaxation methods, and anyone designing Riemann solvers for non-conservative systems. It deserves a serious referee, but the referee should push for measured convergence and archived test cases.","headline":"Useful comparison and a clean algebraic link between HLLC closures, but 'mesh convergence' is asserted, not measured, and that is the load-bearing point.","tokens_in":28145,"tokens_out":3053,"would_cite":true,"duration_ms":25177,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65M08","76T10","35L65"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper establishes that shock solutions of a six-equation two-phase flow model become discretization-dependent at extreme pressure ratios, and that a wave-propagation HLLC scheme with carefully chosen approximate jump conditions is a ro","keywords":["six-equation two-phase model","non-conservative products","HLLC Riemann solver","wave-propagation scheme","path-conservative schemes","instantaneous mechanical relaxation","Kapila model","shock jump conditions"],"falsifier":"A grid-convergence study on the epoxy-spinel shock (for instance, 131,072 cells with measured L1 errors against a reference solution) that shows whether the spread between HLLC wave-propagation, HLLC plus double-integration-by-parts, and Rusanov plus double-integration-by-parts shrinks to zero. If the profiles collapse, the observed spread is under-resolution; if they persist while errors plateau, the non-uniqueness is a property of the model. A complementary test is to compute the intermediate volume fraction using a different path in a path-conservative scheme and check whether the value cha","tokens_in":27286,"feed_emoji":"💥","tokens_out":6274,"duration_ms":48955,"temperature":0.7,"pith_summary":"The paper tackles a practical question: when a compressible two-phase flow is simulated with a six-equation model that is then relaxed to mechanical equilibrium, do different numerical methods give the same answer? Because the model contains non-conservative products, it has no complete set of jump conditions, so the answer is not guaranteed. The authors compare several first-order schemes on very fine meshes and find that for most configurations the methods agree, but under extreme pressure ratios they converge to different shock profiles even though total energy is conserved. They trace this to different partitions of phasic energies across the shock layer and argue that an HLLC wave-propagation scheme, together with a double-integration-by-parts treatment of non-conservative terms, forms a robust framework. The practical message is that for moderate conditions the five-equation model's solutions are insensitive to the discretization, while at extreme conditions the model itself does not single out one shock profile.","feed_headline":"At extreme pressures, two-phase shock solutions depend on the solver","feed_subtitle":"The six-equation model lacks unique jump conditions; solver choices pick different weak solutions at fine resolution.","key_machinery":"The key object is the HLLC approximate Riemann solver for the six-equation system, whose intermediate states are closed by approximate jump conditions for the phasic total energies. The paper shows that two distinct closures found in the literature are two discretizations of the same condition, and that the closure adopted here gives a consistent representation of contact discontinuities, including the non-conservative term u·Σ. The fluctuation form of the wave-propagation method then treats non-conservative terms only at contacts, which is what makes it robust. A second mechanism is the reformulation of the double-integration-by-parts treatment as a path-conservative scheme, which places it","core_discovery":"The central claim is that the six-equation single-velocity two-phase model, and therefore the five-equation model obtained from it by instantaneous pressure relaxation, does not determine a unique shock solution. At extreme pressure ratios, different discretizations of the non-conservative phasic-energy terms converge to different weak solutions, corresponding to different partitions of phasic energies across the shock, even though mixture momentum and total energy remain conserved. The paper also proves that the double-integration-by-parts treatment of non-conservative products, commonly used for diffusion operators, is precisely a path-conservative scheme with a linear path, which explains","pith_inferences":["The findings imply that a 'converged' shock profile for the five-equation model is actually one point in a family of weak solutions; code certification for high-pressure-ratio applications should include a sensitivity study over the discretization of the non-conservative product, not just mesh refinement.","A natural extension would be to parameterize the path in a path-conservative scheme and compute the resulting shock-layer partition of phasic energies, giving a quantitative map of the non-uniqueness without needing very fine meshes.","For finite-rate relaxation, where the six-equation model is a standalone model, the choice of interfacial pressure and the spatial discretization interact, so combined robustness tests are needed before trusting predictions near material interfaces.","Comparing converged six-equation solutions with those of the full seven-equation model under identical extreme conditions would reveal how much of the spread is an artifact of the reduced model's missing jump conditions."],"forward_implications":["For most practical two-phase flow configurations, all tested schemes agree after instantaneous mechanical relaxation, so the non-uniqueness has limited practical impact at moderate pressure ratios.","At extreme pressure ratios, the six-equation model without relaxation is not predictive at the level of shock profiles: different solvers select different weak solutions even on fine meshes.","HLLC-based schemes handle the water-air shock tube at standard Courant numbers, while Rusanov-based central-upwind schemes require smaller time steps, making HLLC wave-propagation the most robust of the tested strategies.","Because the double-integration-by-parts treatment is a path-conservative scheme, the known limitations of path-conservative methods, such as possible convergence to non-entropic weak solutions, carry over to that treatment.","For out-of-equilibrium flows where pressure relaxation is finite or absent, the six-equation model's shock non-uniqueness may become problematic, which argues for moving to better-posed seven-equation models."],"fun_headline_variants":["Shock solutions not unique in two-phase flow model","Solver choice decides shock outcome in two-phase model","Six-equation model fails to pin down shock solutions","Non-conservative terms break uniqueness of two-phase shocks"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The conclusions about solver-dependent convergence rest on the assumption that the fine-mesh solutions have genuinely converged: the paper asserts mesh convergence from visual inspection of profiles on meshes up to 65,536 cells and reports no convergence rates or error norms.","fun_headline_variants_meta":{"raw":{"variants":["Shock solutions not unique in two-phase flow model","Solver choice decides shock outcome in two-phase model","Six-equation model fails to pin down shock solutions","Non-conservative terms break uniqueness of two-phase shocks"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000409,"raw_usage":{"total_tokens":1994,"prompt_tokens":813,"completion_tokens":1181,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":557,"completion_tokens_details":{"reasoning_tokens":1117}},"tokens_in":557,"tokens_out":1181,"duration_ms":7127,"temperature":1.0,"reasoning_tokens":1117,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-04T13:33:29.467917+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"A grid-convergence study on the epoxy-spinel shock (for instance, 131,072 cells with measured L1 errors against a reference solution) that shows whether the spread between HLLC wave-propagation, HLLC plus double-integration-by-parts, and Rusanov plus double-integration-by-parts shrinks to zero. If the profiles collapse, the observed spread is under-resolution; if they persist while errors plateau, the non-uniqueness is a property of the model. A complementary test is to compute the intermediate volume fraction using a different path in a path-conservative scheme and check whether the value cha","supporting_citations":[],"review_version":1}