{"id":"c3d5f554-c6cd-4260-bb4f-6b1c04a454be","arxiv_id":"2504.17729","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"A lowest-order four-field virtual element method with strongly symmetric stress approximates Biot poroelasticity at first order on general 3D polyhedral meshes.","lead":"This paper presents a numerical method for poroelasticity, the coupled deformation and fluid flow in saturated porous materials such as soil and rock. The method handles complicated 3D meshes and computes stress and fluid velocity directly, with first-order accuracy and stability in hard material limits.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The robustness proof in Theorem 3.1 Step 3 assumes 'without loss of generality' that s0 ≤ C3/κ; this fails for fixed positive s0 as λ→∞, so the claimed s0/λ-independence is not established.","rationale":"The reader identified the weakest point accurately: the proof of Theorem 3.1 Step 3 uses an unjustified bound s0 ≤ C3κ^{-1} exactly where robustness with respect to s0 and λ is claimed. This is not a cosmetic gap. The storage term s0 ph(t) must be bounded pointwise to control ∥∇·wh∥, and the energy available from (29) contains only ∥s0^{1/2} ph∥ and κ^{-1/2} terms. The assumption is what makes the pointwise bound possible, and it is not a consequence of the physical definitions. The failure mode is concrete: fixed positive s0 with λ→∞, a physically admissible regime and precisely the near-incompressible limit. The paper's own numerical tests avoid this corner by taking s0 = 0 in the high-λ test and λ = 1 in the positive-s0 test. I therefore agree with the reader's conditional verdict: the method is plausible and the proofs are mostly detailed, but the advertised robustness claim needs either a corrected proof or a revised statement. No change to the reader's verdict is warranted; the conditional acceptance remains the appropriate assessment.","tokens_in":24365,"tokens_out":5036,"duration_ms":52380,"concrete_test":"Rerun the 3D convergence experiment of Test 2 (λ = 10^6, nearly incompressible) with s0 = 0.002 and also s0 = 1 instead of s0 = 0, on the Cube mesh, keeping the same manufactured solution and time step. Compare Eu, Eσ,Π, Ep, and Ew,Π at λ = 1, 10^6, and 10^8. If the errors, or the discrete stability constants and condition numbers, deteriorate with λ for fixed s0 > 0, the claimed independence from s0 and λ is false. If they do not deteriorate, the central claim may still be true, but the proof of Step 3 still needs an independent derivation of (34) that does not invoke s0 ≤ C3κ^{-1}; such a derivation would settle whether the 'without loss of generality' statement can be justified.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The load-bearing point is in Theorem 3.1, Step 3. After integrating (22d) in time and testing with qh = div wh(t), the proof must control the pointwise term ∥s0 ph(t)∥. It then states: 'Owing to the definitions of the storativity coefficient s0 and bulk modulus κ, we can assume without loss of generality that s0 ≤ C3 κ^{-1} ≤ C3 μ^{-1}.' This inequality is doing real work: it converts s0^2∥ph∥² into (C3/μ)∥s0^{1/2} ph∥², the only way the energy (29) can absorb the storage term. No such relation follows from the model; s0 and κ are independent physical parameters. For fixed s0 > 0 and λ→∞, κ = (2μ + 3λ)/3 → ∞, so the inequality fails exactly in the near-incompressible limit that the paper claims to cover. The subsequent bound (34) and the choice ε4 = 1/(4C4) inherit this assumption. Theorem 3.4 says it proceeds 'as in the second and third steps of the proof of the stability estimate in Theorem 3.1,' so the convergence estimate is likewise not proven independent of s0 and λ as stated. The numerics do not close the gap: Test 2 uses s0 = 0, making the assumption trivial, and Test 1 uses λ = 1 with s0 = 0.002; no run combines fixed positive s0 with large λ.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a lowest-order four-field virtual element method for the quasi-static Biot poroelasticity problem. The unknowns are symmetric stress, displacement, Darcy velocity, and pore pressure; stress symmetry is imposed in the discrete VEM space, avoiding an extra Lagrange multiplier, and the method is formulated on polyhedral meshes with piecewise-constant coefficients. After describing the discrete spaces, computable bilinear forms, and the semi-discrete formulation, the authors prove a stability estimate (Theorem 3.1) and an a priori error estimate (Theorem 3.4) with constants stated to be independent of mesh size, storage coefficient s0, and Lamé parameter λ. The scheme is then combined with backward Euler and tested in 3D on cube, tetrahedral, CVT, and random Voronoi meshes, including a nearly incompressible case with s0 = 0, together with a footing benchmark. The central claim is first-order convergence robust to limiting material parameters.","tokens_in":24714,"tokens_out":14781,"duration_ms":146809,"significance":"If the parameter-robustness claim is established, the paper is a useful contribution to mixed VEM methods for poroelasticity: it removes a Lagrange multiplier for stress symmetry, supports general polyhedral meshes, and directly addresses the nearly incompressible and zero-storage limits relevant in geomechanics. The paper's strengths are the explicit construction of the discrete spaces and degrees of freedom, the detailed step-by-step stability proof, the treatment of non-homogeneous boundary conditions, and numerical experiments on four mesh families plus a standard benchmark. The main limitation is a single unsupported 'without loss of generality' assumption in Theorem 3.1, Step 3, which controls the simultaneous behavior of s0 and λ; until that point is repaired, the headline independence claim is not proven.","major_comments":[{"comment":"The sentence 'Owing to the definitions of the storativity coefficient s0 and bulk modulus κ, we can assume without loss of generality that s0 ≤ C3 κ^{-1} ≤ C3 μ^{-1}' is not a consequence of the model. Since κ = (2μ + 3λ)/3, for any fixed s0 > 0 the inequality fails as λ → ∞. This inequality is doing essential work: it converts the pointwise term ∥s0 p_h(t)∥ into a multiple of the energy term ∥s0^{1/2} p_h(t)∥ in (29), leading to the bound (34). Without it, the proof does not establish the claimed independence of the constant in (25) from s0 and λ. The numerical tests do not cover the missing regime: Test 2 uses s0 = 0 and Test 1 uses λ = 1 with s0 = 0.002. Please either replace the 'without loss of generality' claim by a proof that covers all s0 and λ, or explicitly restrict the robustness statement, for example to s0 ≤ Cκ^{-1} or to sequential limits taken in a specified order.","section":"Theorem 3.1, Step 3 (bound for ∇·w̄_h, Eq. (34))"},{"comment":"The proof of the error estimate says it proceeds 'as in the second and third steps of the proof of the stability estimate in Theorem 3.1'; it therefore inherits the unsupported assumption s0 ≤ C3 κ^{-1} from Theorem 3.1, Step 3. Since Theorem 3.4 claims a constant C independent of h, λ, and s0, the convergence result is not fully proven for fixed positive s0 as λ → ∞. The same repair as for Theorem 3.1 is required.","section":"Theorem 3.4 (proof, paragraph after Eq. (44))"}],"minor_comments":[{"comment":"The symbol tf denotes both the final time and, in the definition of Th(f), a tangent vector on a face; please use a different symbol for the tangent vector.","section":"Section 3.2, Eq. (7)"},{"comment":"The boundary-condition notation is inconsistent: after defining ∂Ω = ∂wΩ ∪ ∂pΩ, the flux condition in (4) is written on ∂qΩ instead of ∂wΩ.","section":"Section 2, boundary conditions (4)"},{"comment":"The trace of the fourth-order tensor A used in ξ1,E = 1/2 tr(A|E) is not defined; please state explicitly that it is the trace of A as a map on the space of symmetric tensors.","section":"Section 3.3, Eq. (16)"},{"comment":"The estimated convergence rates are printed as several numbers per curve, which is hard to read; a table listing mesh level, h, and the computed rates for each error would improve readability.","section":"Figures 3 and 4"},{"comment":"The fully discrete scheme is not accompanied by a time-discretization error estimate; the paper states only that first-order convergence is observed numerically. A sentence clarifying that the a priori analysis is for the semi-discrete problem would be helpful.","section":"Section 4.1, Eq. (45)"}],"recommendation":"major_revision","confidential_remarks":"The paper relies substantially on [12] (submitted, with author overlap) for inf-sup and interpolation estimates, and on [21] for the Hellinger–Reissner VEM spaces. I did not treat this as a technical flaw, but the editor may wish to confirm the status of [12] before acceptance. The decisive issue is the s0–κ assumption in Theorem 3.1; if the authors can supply a clean argument covering the case of fixed positive s0 as λ → ∞, or honestly restrict the claimed robustness regime, the paper would meet the bar."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The paper is the first VEM for Biot in a fully mixed four-field form with the stress symmetry enforced strongly in the discrete space, and it deserves the attention of anyone working on poromechanics with polytopal meshes. The lowest-order choice, the computability arguments, and the three-dimensional tests on hexahedral, tetrahedral, and Voronoi meshes are all handled carefully. The convergence rates are around one, matching the stated order, and the footing benchmark is a sensible sanity check. I believe the central construction is sound and the writing is honest.\n\nWhere the paper is soft is in the robustness proof. Theorem 3.1, Step 3, binds ∥s0 p_h(t)∥ by assuming, 'without loss of generality,' that s0 ≤ C3 κ^{-1}. That is not a generality: s0 and κ are independent material parameters, and for fixed s0 > 0 the inequality fails when λ → ∞. The step is load-bearing because it is the only way the storage term is absorbed into the energy (29), and Theorem 3.4 explicitly inherits it. So as written the paper does not establish independence from s0 and λ, despite the abstract's emphasis on robustness. The numerics do not cover the missing corner: Test 2 sets s0 = 0 and Test 1 takes λ = 1. This is a genuine gap, but it is a gap in a proof, not in the formulation; the method may well be robust, and the paper should be revised to either supply a corrected argument or state the theorem under a condition like s0 κ bounded.\n\nTwo smaller points. The proofs rely on the submitted companion [12] for key interpolation and stability results; since [12] has author overlap, the referee should ask for the full statements or an appendix. And the convergence plots show some superconvergence on coarse meshes for displacement and pressure; the explanation given (preasymptotic regime) is plausible but brief.\n\nI would send this to peer review. The idea and the numerical study are solid enough that a serious referee time is warranted, and the robustness claim can be fixed or clarified. If I worked on VEM for poroelasticity, I would cite it for the four-field construction, though with a note about the proof gap.","headline":"A competent four-field VEM paper for Biot whose main robustness theorem overreaches; the construction is genuinely new and the numerics back the convergence, but the claimed s0/λ-independence is not proven.","tokens_in":25222,"tokens_out":3857,"would_cite":true,"duration_ms":36813,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65M12","65M60","74F10","76S05"],"pacs":[],"model":"deepseek-v4-flash","headline":"A lowest-order four-field virtual element method for Biot poroelasticity is stable and first-order convergent with constants independent of the mesh size, the storage coefficient, and the Lamé parameter λ.","keywords":["Virtual Element Method","Biot poroelasticity","mixed formulation","polyhedral meshes","stress symmetry","nearly incompressible","robust error analysis","poromechanics"],"falsifier":"Run the manufactured-solution test with a fixed positive storage coefficient, say $s_0 = 10^{-3}$, and a very large Lamé parameter, say $\\lambda = 10^8$, refining both $h$ and $\\Delta t$; if the error constants grow or the convergence rate falls below first order, the claimed robustness with respect to $s_0$ and $\\lambda$ is not borne out. Separately, one can check the hypothesis directly: with $s_0>0$ fixed and $\\lambda \\to \\infty$, $\\kappa^{-1} \\to 0$, so the inequality $s_0 \\le C_3 \\kappa^{-1}$ cannot hold for any finite $C_3$.","tokens_in":1609,"feed_emoji":"💧","tokens_out":4484,"duration_ms":60815,"temperature":0.7,"pith_summary":"The paper proposes and analyzes a lowest-order four-field virtual element discretization of Biot's poroelasticity equations, enforcing stress symmetry directly inside the discrete space rather than through a Lagrange multiplier. It claims the scheme is stable and first-order convergent on general polyhedral meshes, with error constants independent of mesh size, the storage coefficient $s_0$, and $\\lambda$. If correct, the method yields simultaneous $O(h)$ accuracy for stress, displacement, velocity, and pressure, including nearly incompressible and zero-storage limits. Three-dimensional tests on hexahedral, tetrahedral, and Voronoi meshes support the predicted rates.","feed_headline":"Biot poroelasticity solved at first order on polyhedral meshes","feed_subtitle":"A four-field virtual element scheme keeps stress and flow errors stable as materials approach incompressible limits.","key_machinery":"The load-bearing construction is the three-dimensional symmetric Hellinger–Reissner virtual element stress space, where the normal traction on each face lies in the six-dimensional space $T_h(f)$ and the divergence is a rigid body motion. This space is coupled to a piecewise-rigid displacement space, a lowest-order virtual Raviart–Thomas velocity space, and a piecewise-constant pressure space. Local bilinear forms are made computable through $L^2$ projections onto constants plus face-based stabilization, and the discrete stress projection, the coupling terms, and the mixed terms are all recovered from the degrees of freedom; the stability analysis then runs on a time-integrated energy identity with weighted norms that expose the dependence on $s_0$ and $\\lambda$.","core_discovery":"The central claim is stated in Theorem 3.4: for the semi-discrete problem, the errors in displacement, stress, velocity, and pressure satisfy a bound of the form $O(h)$ with a constant independent of $h$, $\\lambda$, and $s_0$. The lowest-order fully discrete method, obtained with backward Euler in time, is then confirmed numerically to converge at first order in space. The distinctive feature is a four-field formulation in which the stress tensor lives in a symmetric $H(\\mathrm{div})$-conforming virtual element space, so symmetry is strongly imposed and no Lagrange multiplier is needed; the divergence of the discrete stress lies in the space of rigid body motions, which makes the mixed terms computable from the degrees of freedom.","pith_inferences":["The proof of robustness in Theorem 3.1, Step 3, assumes $s_0 \\le C_3 \\kappa^{-1}$ 'without loss of generality'; for fixed positive $s_0$ and $\\lambda \\to \\infty$, this assumption fails, so the claimed independence from $s_0$ and $\\lambda$ would be fully established only under a small-storage or large-$\\kappa$ regime, and a numerical test with fixed $s_0>0$ and very large $\\lambda$ would probe this","Because the stress space is symmetric and divergence-conforming, the same degree-of-freedom layout could be extended to hybrid-dimensional models where fractures are lower-dimensional flow domains, a direction the paper names as future work; the strongly symmetric stress space would then also serve contact and frictional interface conditions.","The structure of the coupling block suggests that the same four-field VEM could be paired with multigrid or preconditioned saddle-point solvers inherited from the two mixed subproblems, although no preconditioner is analyzed here."],"forward_implications":["On general polyhedral meshes, the method gives first-order accuracy for stress, displacement, velocity, and pressure without an extra Lagrange multiplier for stress symmetry.","The error and stability constants do not degrade when the storage coefficient tends to zero or the material approaches incompressibility, so the scheme is intended to cover the limiting cases that standard formulations struggle with.","Local mass and momentum conservation follow from the mixed form of both the flow and the mechanical subproblems, which is useful in subsurface and fracture-scale applications.","The fully discrete system retains a symmetric saddle-point structure with a symmetric coupling block, which is directly relevant for designing block preconditioners for practical poroelasticity simulations."],"supporting_citations":[{"why":"Supplies the three-dimensional Hellinger–Reissner virtual element stress space, its degrees of freedom, and the discrete inf-sup condition used for the stress-displacement coupling.","marker":"[21]"},{"why":"Provides the stability and interpolation estimates of Hellinger–Reissner virtual element spaces used in the error analysis and in the stabilization bounds.","marker":"[12]"},{"why":"Supplies the basic principles of mixed virtual element methods, including the lowest-order virtual Raviart–Thomas velocity space and its projections.","marker":"[14]"},{"why":"Gives the robust error analysis of coupled mixed methods for Biot's consolidation model that the present work extends to a four-field strongly symmetric setting.","marker":"[28]"},{"why":"Establishes the differential-algebraic-equation framework for existence and uniqueness of the semidiscrete solution and presents an earlier four-field mixed formulation for Biot's model.","marker":"[40]"}],"fun_headline_variants":["Four-field VEM solves Biot at first order","First-order Biot solver: four-field VEM, no multiplier","Four-field VEM: first-order Biot on polyhedral meshes","Biot poroelasticity: four-field VEM, first order","Stress-symmetric four-field VEM gives first-order Biot"],"cache_read_input_tokens":27264,"weakest_assumption_plain":"The stability proof assumes, in Step 3 of Theorem 3.1, that $s_0 \\le C_3 \\kappa^{-1}$ holds 'without loss of generality'; this is not a consequence of the model definitions and fails for fixed positive $s_0$ as $\\kappa \\to \\infty$, so the claimed independence of the stability constant from $s_0$ and $\\lambda$ rests on that assumption.","fun_headline_variants_meta":{"raw":{"variants":["Four-field VEM solves Biot at first order","First-order Biot solver: four-field VEM, no multiplier","Four-field VEM: first-order Biot on polyhedral meshes","Biot poroelasticity: four-field VEM, first order","Stress-symmetric four-field VEM gives first-order Biot"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001411,"raw_usage":{"total_tokens":5656,"prompt_tokens":854,"completion_tokens":4802,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":470,"completion_tokens_details":{"reasoning_tokens":4713}},"tokens_in":470,"tokens_out":4802,"duration_ms":32511,"temperature":1.0,"reasoning_tokens":4713,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T10:32:58.144436+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the manufactured-solution test with a fixed positive storage coefficient, say $s_0 = 10^{-3}$, and a very large Lamé parameter, say $\\lambda = 10^8$, refining both $h$ and $\\Delta t$; if the error constants grow or the convergence rate falls below first order, the claimed robustness with respect to $s_0$ and $\\lambda$ is not borne out. Separately, one can check the hypothesis directly: with $s_0>0$ fixed and $\\lambda \\to \\infty$, $\\kappa^{-1} \\to 0$, so the inequality $s_0 \\le C_3 \\kappa^{-1}$ cannot hold for any finite $C_3$.","supporting_citations":[{"cited_title":"Dassi, C","cited_arxiv_id":null,"evidence_quote":"Supplies the three-dimensional Hellinger–Reissner virtual element stress space, its degrees of freedom, and the discrete inf-sup condition used for the stress-displacement coupling."},{"cited_title":"Botti, L","cited_arxiv_id":null,"evidence_quote":"Provides the stability and interpolation estimates of Hellinger–Reissner virtual element spaces used in the error analysis and in the stabilization bounds."},{"cited_title":"Brezzi, R","cited_arxiv_id":null,"evidence_quote":"Supplies the basic principles of mixed virtual element methods, including the lowest-order virtual Raviart–Thomas velocity space and its projections."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the robust error analysis of coupled mixed methods for Biot's consolidation model that the present work extends to a four-field strongly symmetric setting."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Establishes the differential-algebraic-equation framework for existence and uniqueness of the semidiscrete solution and presents an earlier four-field mixed formulation for Biot's model."}],"review_version":1}