Processing math: 33%
Research article

A regional macroeconomic approach to the link between immigration and wages in Spain: New insights from a spatial econometric perspective

  • Using an extended spatial wage equation that introduces novel elements in the weight matrix and the treatment of immigrants' human capital, we analyzed the impact of immigration on wages. Examining the Spanish case at the provincial level as a sort of laboratory over the period 2006–2021, the results indicated a small, negative effect that disappears when immigrants are highly skilled. Furthermore, the empirical evidence does not support an association between the unemployment rate and wages, confirming that Spain's well-known high degree of wage rigidity persisted in the aftermath of the financial and economic crisis of 2008. Our results are shown to be robust across a battery of tests that assess their reliability from different perspectives.

    Citation: Adolfo Maza, Maria Hierro. A regional macroeconomic approach to the link between immigration and wages in Spain: New insights from a spatial econometric perspective[J]. National Accounting Review, 2025, 7(2): 177-195. doi: 10.3934/NAR.2025008

    Related Papers:

    [1] Heesung Shin, Jiang Zeng . More bijections for Entringer and Arnold families. Electronic Research Archive, 2021, 29(2): 2167-2185. doi: 10.3934/era.2020111
    [2] Bin Han . Some multivariate polynomials for doubled permutations. Electronic Research Archive, 2021, 29(2): 1925-1944. doi: 10.3934/era.2020098
    [3] Shishuo Fu, Zhicong Lin, Yaling Wang . Refined Wilf-equivalences by Comtet statistics. Electronic Research Archive, 2021, 29(5): 2877-2913. doi: 10.3934/era.2021018
    [4] Jiafan Zhang . On the distribution of primitive roots and Lehmer numbers. Electronic Research Archive, 2023, 31(11): 6913-6927. doi: 10.3934/era.2023350
    [5] Dmitry Krachun, Zhi-Wei Sun . On sums of four pentagonal numbers with coefficients. Electronic Research Archive, 2020, 28(1): 559-566. doi: 10.3934/era.2020029
    [6] Massimo Grossi . On the number of critical points of solutions of semilinear elliptic equations. Electronic Research Archive, 2021, 29(6): 4215-4228. doi: 10.3934/era.2021080
    [7] Ji-Cai Liu . Proof of Sun's conjectural supercongruence involving Catalan numbers. Electronic Research Archive, 2020, 28(2): 1023-1030. doi: 10.3934/era.2020054
    [8] Hongjian Li, Kaili Yang, Pingzhi Yuan . The asymptotic behavior of the reciprocal sum of generalized Fibonacci numbers. Electronic Research Archive, 2025, 33(1): 409-432. doi: 10.3934/era.2025020
    [9] Hanpeng Gao, Yunlong Zhou, Yuanfeng Zhang . Sincere wide τ-tilting modules. Electronic Research Archive, 2025, 33(4): 2275-2284. doi: 10.3934/era.2025099
    [10] Taboka Prince Chalebgwa, Sidney A. Morris . Number theoretic subsets of the real line of full or null measure. Electronic Research Archive, 2025, 33(2): 1037-1044. doi: 10.3934/era.2025046
  • Using an extended spatial wage equation that introduces novel elements in the weight matrix and the treatment of immigrants' human capital, we analyzed the impact of immigration on wages. Examining the Spanish case at the provincial level as a sort of laboratory over the period 2006–2021, the results indicated a small, negative effect that disappears when immigrants are highly skilled. Furthermore, the empirical evidence does not support an association between the unemployment rate and wages, confirming that Spain's well-known high degree of wage rigidity persisted in the aftermath of the financial and economic crisis of 2008. Our results are shown to be robust across a battery of tests that assess their reliability from different perspectives.



    Poroelastic materials play an important role in various fields such as geophysics [1,2,3], biomechanics, chemical engineering [4,5], and materials science [6,7]. The physical process in porous media is usually described by a poroelasticity model. In homogeneous isotropic linear elastic porous media, the governing equations for the quasi-static poroelasticity model are given as

    2divσ(u)+αp=fin ΩT:=Ω×(0,T)Rd×(0,T), (1.1)
    (c0p+αdivu)t+divvf=ϕin ΩT, (1.2)

    where

    σ(u):=με(u)+λtr(ε(u))I,ε(u):=12(u+Tu), (1.3)
    vf:=Kμfp. (1.4)

    Here u and p represent the displacement vector of solid and the pressure of solvent, σ(u) is the (effective) stress tensor, vf is the volumetric solvent flux, and Eq (1.4) is so-called Darcy's law. I denotes the d×d identity matrix and ε(u) is the Green strain tensor. f is the body force and ϕ is the source term. The permeability tensor K=K(x) is assumed to be symmetric and uniformly positive definite in the sense that there exist positive constants K1 and K2 such that K1|ζ|2K(x)ζζK2|ζ|2 for a.e. xΩ and ζRd. The parameters μf,α and c0 are the solvent viscosity, Biot-Willis constant and constrained specific storage coefficient, respectively. λ and μ (Note that we use μ instead of 2μ in (1.3), as often seen in the literature, to simplify the notation) are Lamé constants, which can be computed from the Young's modulus E and Poisson ratio ν as

    λ=Eν(1+ν)(12ν),μ=E2(1+ν).

    Much effort has been carried out for the above poroelasticity model. Biot [1] proposed the three-dimensional consolidation theory, providing important theoretical support for the development of the poroelasticity model. Showalter [8] showed the well-posedness of a weak solution of a linear poroelasticity model. Phillips and Wheeler proposed and analyzed a continuous-in-time linear poroelasticity model [9] and a discrete-in-time linear poroelasticity model [10] with the continuous Galerkin mixed finite element method (MFEM) in space. After that, a discontinuous Galerkin MFEM for the model is also developed [11]. The nonconforming FEM was studied for the linear poroelasticity model [12,13], and the MINI and stabilized P1P1 elements were employed to solve the linear poroelasticity model [14]. Ju et al. [15] developed parameter-robust numerical algorithms for the Biot model and applied the algorithms to simulations of brain edema. Yi [16] studied two kinds of locking phenomena (Poisson locking and pressure oscillation) in linear poroelasticity and proposed a new family of MFEM to overcome them. To reveal the underlying deformation and diffusion processes of the original model, Feng et al. reformulated the linear poroelasticity model to a fluid-fluid coupled system and proposed a multiphysics MFEM to overcome the locking phenomena [17]. One can refer to more numerical methods for solving the poroelasticity model, such as the Galerkin least square methods [18,19], two-field finite element method [20,21], four-field finite element method[22,23], finite volume method [24], and Weak Galerkin finite element method [25,26].

    In practical problems, the model parameters have different ranges, such as the Poisson ratio ν=0.10.3 and the permeability K=10141016m2 in geophysics [6]. This requires us to construct numerical methods to be robust for the selection of model parameters. Following the idea of [17], we first reformulate the original problem into a four-field by introducing some auxiliary variables. After that, we discretize the four-field problem by the backward Euler method in time and design the coupled and decoupled algorithms by using MFEM with P2P1P1P1 element pairs in space. The inf-sup condition is proved to be uniformly independent of the physical parameters, and the optimal convergence rates both in space and time are also obtained for the proposed algorithms. It is pointed out that, for the decoupled algorithm [17], there exists the constraint τ=O(h2) on the mesh size h and time step τ while performing the error estimates. But there is no constraint on h and τ for the decoupled algorithm arising from our four-field problem. The reason is that the decoupled algorithm ensures that the variables p,ξ, and η are in the same time layer, thereby avoiding the use of inverse inequality. To improve computational efficiency and achieve high-accuracy in time, we further design a two-step backward differentiation formula (BDF2) Galerkin algorithm, which combines the MFEM with P2P1P1 element pairs for the space variables. The corresponding stability and error estimates of the BDF2 algorithm are performed, and the optimal convergence rate in time is also obtained in Theorem 7. Numerical tests are also provided to verify that our proposed algorithms are robust for selection of ν and K, overcome the Poisson locking in stress field and nonphysical oscillation in the pressure field, and have the optimal convergence rates. In addition, our algorithms are also used to simulate the brain edema. The numerical results on pressure are consistent with the experimental results [27]. The influence of the Poisson rate ν on the displacement of brain tissue and maximum pressure is investigated, which shows that the values of ν have little influence on the maximum pressure, and the maximum tissue deformation becomes smaller as ν gets bigger.

    The remainder of this paper is organized as follows: In Section 2, the original problem is equivalently transformed into a four-field problem, and the four-field problem is also proved to satisfy the inf-sup condition. In Section 3, the coupled and decoupled algorithms with the Euler method are proposed for the four-field problem, and the error estimates are also presented. In Section 4, the BDF2 algorithm for the poroelasticity problem is also designed, and the stability analysis and error estimates are established. In Section 5, numerical examples are provided to verify that the algorithms are robust with respect to the key physical parameters, reach the optimal convergence rates, and overcome the locking phenomenon. In Section 6, the algorithms are applied to simulate the brain edema. Finally, we draw a conclusion to summarize the main results of this paper.

    The standard function space notations are adopted; their precise definitions can be found in [28,29,30]. In particular, (,) and , denote respectively the standard L2(Ω) and L2(Ω) inner products, where ΩRd(d=2,3) is a bounded polygonal domain with the boundary Ω. For any Banach space B, we let B=[B]d, and use B to denote its dual space. We also introduce the function spaces

    L20(Ω):={qL2(Ω);(q,1)=0},X:=H1(Ω).

    Define the space of infinitesimal rigid motions by RM={r:=a+b×x;a,b,xRd}. As in [28,30,31], RM is the kernel of the operator ε, i.e., rRM if and only if ε(r)=0. Hence, we have

    ε(r)=0,divr=0,rRM. (2.1)

    Denote L2(Ω) and H1(Ω) respectively by the subspaces of L2(Ω) and H1(Ω), which are orthogonal to RM, and are given as,

    H1(Ω):={vH1(Ω);(v,r)=0rRM},L2(Ω):={gL2(Ω);g,r=0rRM}.

    It follows from [17,32] that there exists a constant ck>0 such that

    vL2(Ω)=infrRMv+rL2(Ω)rL2(Ω)ckε(v)L2(Ω),vH1(Ω). (2.2)

    To close the above system (1.1)–(1.4), here we set the following boundary and initial conditions:

    2σ(u)nαpn=f1on ΩT:=Ω×(0,T), (2.3)
    1μfKpn=ϕ1on ΩT, (2.4)
    u=u0,p=p0in Ω×{t=0}. (2.5)

    By introducing the variables q=divu,η=c0p+αq,ξ=αpλq, one has

    p=κ1ξ+κ2η,q=κ1ηκ3ξ, (2.6)

    with

    κ1=αα2+λc0,κ2=λα2+λc0,κ3=c0α2+λc0.

    Thus, the problem (1.1)–(1.4) can be equivalently written as

    2μdivε(u)+ξ=fin ΩT, (2.7)
    κ3ξ+divu=κ1ηin ΩT, (2.8)
    κ1ξ+κ2ηp=0in ΩT, (2.9)
    ηt1μfdiv(Kp)=ϕin ΩT. (2.10)

    Remark 1. We point out that the following results with respect to the pure Neumann condition can be extended to the case of the pure Dirichlet condition or mixed Dirichlet-Neumann condition. If the boundary condition (2.3) is replaced by a pure Dirichlet condition or the mixed Dirichlet-Neumann condition, there is no need to introduce the space H1(Ω) for the variable u.

    We here consider an equidistant partition of the time interval [0,T] with a time step τ=T/N and tn=nτ (n=1,2,,N), where NN. The focus of this section is on the proof of the inf-sup condition, which will hold uniformly independently of physical parameters. To do that, we first discretize the time-dependent problem (2.7)–(2.10) with (2.3)–(2.5) by the backward Euler method and BDF2 scheme in time to make it a static problem in each time step. Note that the discussion on inf-sup condition is similar for the backward Euler method and BDF2 scheme. Below, we mainly consider the backward Euler method. By the backward Euler method, the weak formulation of problem (2.7)–(2.10) with (2.3)–(2.5) reads: (v,φ,ω,ψ)H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω), (u,ξ,η,p)H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω) such that

    2μ(ε(u),ε(v))(ξ,divv)=(f,v)+f1,v, (2.11)
    κ3(ξ,φ)+(divu,φ)=κ1(η,φ), (2.12)
    (κ1ξ,ω)+(κ2η,ω)(p,ω)=0, (2.13)
    (η,ψ)+τμf(Kp,ψ)=(˜ϕ,ψ)+˜ϕ1,ψ, (2.14)

    where ˜ϕ(x,tn)=τϕ(x,tn)+η(x,tn1) and ˜ϕ1=τϕ1.

    Define

    A(u,ξ,η,p)=(μdivε00divκ3κ100κ1κ21001τμfdivK)(uξηp).

    Taking Λ=H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω) and ξ0 as the mean-value zero part of ξ, we define

    (u,ξ,η,p)Λ=μ12ε(u)L2(Ω)+κ123ξL2(Ω)+δ120ξ0L2(Ω)+κ122ηL2(Ω)+(τμf)12K12pL2(Ω), (2.15)

    where δ0 is a positive constant, uH1(Ω), ξL2(Ω), ηL2(Ω) and pH1(Ω).

    The following theorem shows that A is invertible and the inverse is a map from Λ to Λ.

    Theorem 1. There exists a constant β>0 for the problem (2.11)(2.14), independent of μ,κ2,κ3 and Kτ/μf, such that the following inf-sup condition holds:

    inf(u,ξ,η,p)Λsup(v,φ,ω,ψ)Λ(A(u,ξ,η,p),(v,φ,ω,ψ))(u,ξ,η,p)Λ(v,ϕ,ω,ψ)Λβ. (2.16)

    Proof. To prove (2.16), we need to prove the following inequalities: (0,0,0,0)(u,ξ,η,p)Λ, (v,φ,ω,ψ)Λ such that

    (v,φ,ω,ψ)Λβ1(u,ξ,η,p)Λ,(A(u,ξ,η,p),(v,φ,ω,ψ))β2(u,ξ,η,p)2Λ,

    where β1 and β2 are positive constants and independent of μ,κ1,κ2,κ3 and Kτ/μf.

    Suppose that (0,0,0,0)(u,ξ,η,p)Λ is given. It follows from the theory of the Stokes equation that there exists a constant β0>0, depending only on the domain Ω, and wH1(Ω) such that

    (divw,ξ)=ξ02L2(Ω),(ε(w),ε(w))β20ξ02L2(Ω). (2.17)

    Taking v=uδ0w, ϕ=ξ, ω=η, ψ=p in (A(u,ξ,η,p),(v,φ,ω,ψ)), we have

    (A(u,ξ,η,p),(v,φ,ω,ψ))=μ(ε(u),ε(u))δ0μ(ε(u),ε(w))(ξ,divu)+δ0(ξ,divw)+κ3(ξ,ξ)+(divu,ξ)κ1(η,ξ)+(κ1ξ+κ2η,η)(p,η)+(η,p)+τμfKp2L2(Ω)=με(u)2L2(Ω)δ0μ(ε(u),ε(w))+δ0(ξ,divw)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω). (2.18)

    Using Young's inequality and (2.17), we have

    (ε(u),ε(w))12θ0(ε(u),ε(u))+θ02(ε(w),ε(w))12θ0(ε(u),ε(u))+θ0β202(ξ0,ξ0),θ0>0. (2.19)

    Inserting (2.19) into (2.18), using (2.17) and choosing θ0=δ0=β20μ1, we obtain

    (A(u,ξ,η,p),(v,φ,ω,ψ))μ2ε(u)2L2(Ω)+δ02ξ02L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω). (2.20)

    Noting that

    (u,ξ,η,p)2Λ=(μ12ε(u)L2(Ω)+κ123ξL2(Ω)+κ122ηL2(Ω)+(τμf)12KpL2(Ω))24(με(u)2L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω))4(με(u)2L2(Ω)+δ0ξ02L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω)), (2.21)

    and using (2.21) in (2.20), we have

    (A(u,ξ,η,p),(v,φ,ω,ψ))18(u,ξ,η,p)2Λ. (2.22)

    It is easy to check that

    (v,φ,ω,ψ)Λ(1+μ12δ120β0)(u,ξ,η,p)Λ. (2.23)

    Combining (2.22) and (2.23), the (2.16) holds and β=1/(8(1+μ12δ120β0)). The proof is complete.

    Theorem 1 shows that the bilinear form (A(),) is coercive. Moreover, it is easy to check that (A(),) is bounded. According to the Lax-Milgram theorem, there exists a unique solution (u,ξ,η,p)Λ for the linear system (2.11)–(2.14). Thus, there holds the following theorem:

    Theorem 2. There exists a unique weak solution to problem (2.7)(2.10) with (2.3)(2.5). Likewise, there exists a unique weak solution to problem (1.1)(1.4) and (2.3)(2.5).

    Let Th be a quasi-uniform triangulation or rectangular partition of Ω with maximum mesh size h, and ˉΩ=KThˉK. In this part, we use the backward Euler method in time discretization. For a smooth function v defined on [0,T], we denote dtvn=(vnvn1)/τ.

    A number of stable mixed finite element spaces (Xh,Mh) have been presented in the literature [33]. A typical example is the following so-called Taylor-Hood element (cf. [33,34]):

    Xh={vhC0(¯Ω);vh|KP2(K)KTh},Mh={φhC0(¯Ω);φh|KP1(K)KTh}.

    Finite element approximation space Wh for p and η variables can be chosen independently; any piecewise polynomial space is acceptable provided that WhMh. A convenient choice is Wh=Mh.

    Define

    Vh:={vhXh;(vh,r)=0,rRM}. (3.1)

    If the problem (2.11)–(2.14) is discretized by the finite element spaces Vh,Mh, and Wh, then the discrete counterpart is to find (uh,ξh,ηh,ph)Vh×Mh×Wh×Wh such that

    2μ(ε(uh),ε(vh))(ξh,divvh)=(f,vh)+f1,vh,vhVh, (3.2)
    κ3(ξh,φh)+(divuh,φh)=κ1(ηh,φh),φhMh, (3.3)
    κ1(ξh,ωh)+κ2(ηh,ωh)(p,ωh)=0,ωhWh, (3.4)
    (ηh,ψh)+τμf(Kph,ψh)=(˜ϕ,ψh)+˜ϕ1,ψh,ψhWh. (3.5)

    Set Λh=Vh×Mh×Wh×Wh, and define the operator Ah:ΛhΛh arisen from (3.2)–(3.5). Then we have the following theorem:

    Theorem 3. For the problem (3.2)(3.5), there exists a constant β>0, independent of μ,κ2,κ3, and Kτ/μf, such that the following inf-sup condition holds:

    inf(uh,ξh,ηh,ph)Λhsup(vh,φh,ωh,ψh)Λh(Ah(uh,ξh,ηh,ph),(vh,φh,ωh,ψh))(uh,ξh,ηh,ph)Λh(vh,φh,ωh,ψh)Λhβ. (3.6)

    We omit the proof of Theorem 3, since it is completely analogous to the proof of Theorem 1. Moreover, Theorem 3 implies that Ah is invertible and the inverse is a map from Λh to Λh with operator norm independent of the parameters.

    Algorithm 1 A coupled algorithm
    Input: Evaluate u0hVh,p0hWh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      for n=0,1,2,, do
        solve for (un+1h,ξn+1h,ηn+1h,pn+1h)Vh×Mh×Wh×Wh such that
                μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηn+1h,φh),φhMh,κ1(ξn+1h,ωh)+κ2(ηn+1h,ωh)(pn+1h,ωh)=0,ωhWh,(dtηn+1h,ψh)+1μf(Kpn+1h,ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh,
      end for

    Algorithm 2 A decoupled algorithm
    Input: Evaluate u0hVh,p0hWh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      for n=0,1,2,, do
        (i) solve for (un+1h,ξn+1h)Vh×Mh such that
                μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηnh,φh),φhMh.
        (ii) Using (un+1h,ξn+1h) obtained in (i), solve for (ηn+1h,pn+1h)Wh×Wh by
                κ1(ξn+1h,ωh)+κ2(ηn+1h,ωh)(pn+1h,ωh)=0,ωhWh,(dtηn+1h,ψh)+1μf(Kpn+1h,ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh,
      end for

    Below, we assume f,f1, ϕ, and ϕ1 all are independent of t. Note that all the results can be easily extended to the case of time-dependent source functions.

    Now, we consider the fully discrete MFEM for the problem (2.7)–(2.10) with (2.3)–(2.5) based on the backward Euler method, which involves the coupled (i.e., Algorithm 1) and decoupled (i.e., Algorithm 2) algorithms. The coupled algorithm solves four unknowns simultaneously based on problem (2.7)–(2.10). The decoupled algorithm only solves two subproblems. Noting that problem (2.7) and (2.8) can be regarded as a generalized Stokes problem for a given η and problem (2.9) and (2.10) is a coupled problem associated with the diffusion problem, we can solve a generalized Stokes problem for u and ξ by using the solution of η at the previous time-step. After that, we solve the diffusion problem for η and p.

    To derive the error estimates of the MFEM, we define L2(Ω)-projection Qh:L2(Ω)Xkh by

    (Qhφ,ψh)=(φ,ψh),ψhXkh and φL2(Ω), (3.7)

    where Xkh:={ψhC0;ψh|KPk(K)KTh}, k is the degree of piecewise polynomial on K.

    Next, for any φH1(Ω), we define the elliptic projection Sh:H1(Ω)Xkh as

    (KShφ,φh)=(Kφ,φh),φhXkh, (3.8)
    (Shφ,1)=(φ,1), (3.9)

    and for any vH1(Ω), we define its elliptic projection Rh:H1(Ω)Vkh as

    (ε(Rhv),ε(wh))=(ε(v),ε(wh))whVkh, (3.10)

    where Vkh:={vhC0;vh|KPk(K),(vh,r)=0rRM}.

    It follows from [28] that for any 0sk, Qh,Sh and Rh satisfy

    QhφφL2(Ω)+h(Qhφφ)L2(Ω)Chs+1φHs+1(Ω),φHs+1(Ω), (3.11)
    ShφφL2(Ω)+h(Shφφ)L2(Ω)Chs+1φHs+1(Ω),φHs+1(Ω), (3.12)
    RhvvL2(Ω)+h(Rhvv)L2(Ω)Chs+1vHs+1(Ω),vHs+1(Ω). (3.13)

    Furthermore, we introduce the following error notations:

    Enu=u(tn)unh,Enξ=ξ(tn)ξnh,Enη=η(tn)ηnh,Enp=p(tn)pnh.

    Also, we denote

    Enu=u(tn)Rh(u(tn))+Rh(u(tn))unh:=Ynu+Znu,Eni=i(tn)Sh(i(tn))+Sh(i(tn))inh:=Yni+Zni,i=ξ,η,p,Eni=i(tn)Qh(i(tn))+Qh(i(tn))inh:=Fni+Gni,i=ξ,η,p.

    Then we present the error estimates for the coupled and decoupled algorithms by the following theorems.

    Theorem 4. Let {(unh,ξnh,ηnh,pnh)}n0 be solutions of the coupled algorithm (i.e., Algorithm 1). Then it holds

    max0nN[μ(u(tn)unh)L2(Ω)+κ2η(tn)ηnhL2(Ω)+κ3ξ(tn)ξnhL2(Ω)]+(τNn=01μfK12(p(tn)pnh)2L2(Ω))12C1τ+C2h2. (3.14)

    Theorem 5. Let {(unh,ξnh,ηnh,pnh)}n0 be solutions of the decoupled algorithm (i.e., Algorithm 2). Then for the case c0>0 and λ<, it holds

    max0nN[μ(u(tn)unh)L2(Ω)+κ2η(tn)ηnhL2(Ω)+κ3ξ(tn)ξnhL2(Ω)]+(τNn=01μfK12(p(tn)pnh)2L2(Ω))12ˉC1τ+ˉC2h2. (3.15)

    The proof of the aforementioned theorems follows a similar approach to that of [17, Theorem 3.7]; thus, further details are omitted here. Our primary focus is on the constraint condition τ=O(h2) in the error estimates for the algorithm with θ=0 [17]. This constraint arises because the choice of the test function ψh results in the variables p,ξ, and η being in different time layers, which in turn leads to the application of an inverse inequality in the estimation of the right-hand side term. To address this issue, a straightforward solution is to select an appropriate test function ψh that ensures the variables p,ξ, and η reside in the same time layer. To do that, we first examine the following error equations:

    2κ1(Gn+1ξ,ωh)+κ2(Gn+1η,ωh)(Gn+1p,ωh)=0,ωhWh, (3.16)
    (dtGn+1η,ψh)+1μf(KZn+1p,ψh)=(Rn+11,η,ψh),ψhWh, (3.17)

    which are derived from the fourth equation of Algorithm 2 and the continuous problem (2.10) the fourth equation of Algorithm 2 and the continuous problem (2.10). Based on Algorithm 2, we can define η1h by

    κ1(η1h,φh)=κ3(ξ0h,φh)+(divu0h,φh).

    Setting ωh=dtGnη,ψh=Znp in (3.16) and (3.17) separately after lowering the super index n+1 to n, applying the summation operator τln=0 to both sides of (3.17) and utilizing the fact that G1η=0, there holds

    κ22Glη2L2(Ω)+τln=0[κ1(dtGnη,Gn+1ξ)+κ2τ2dtGnη2L2(Ω)+1μfK12Znp2L2(Ω)]=τln=0[(dtGnη,YnpFnp)+(dtη(tn+1)ηt(tn+1),Znp)+κ1τ(dtGn+1ξ,dtGnη)]. (3.18)

    Note that the choice of ψh=Znp is critical. Unlike the setting in [17], where ψh=ˆZnp:=κ1Zn+1ξ+κ2Znη is employed, selecting ψh=Znp eliminates the need to apply the inverse inequality to bound a gradient term (as seen in [17, estimate (3.52)]). Instead, we only need to apply the Cauchy-Schwarz and Young inequalities to the last term on the right-hand side of (3.18) to derive

    κ1τ2ln=0(dtGn+1ξ,dtGnη)κ1τ2ln=0[12ϵdtGn+1ξ2L2(Ω)+ϵ2dtGnη)2L2(Ω)], (3.19)

    where ϵ(2κ2κ3,+) is a positive constant from Young inequality.

    Additionally, we point out that:

    (i) The aim of ϵ(2κ2κ3,+) is to ensure the signs of terms κ2κ1ϵ2τ2ln=0dtGnη2L2(Ω) and (κ34κ12ϵ)τ2ln=0dtGnη2L2(Ω) to be positive. In addition, for the case of λ+, there does not exist the term κ1τ2ln=0(dtGn+1ξ,dtGnη) in error estimates. For the case of c0=0, it has a similar discussion for the term κ1τ2ln=0(dtGn+1ξ,dtGnη).

    (ii) For the decoupled algorithm, it does not exist the constraint on the time step τ and mesh size h.

    We now consider the BDF2 scheme in time discretization for the problem (2.7)–(2.10) with (2.3)–(2.5). For a smooth function v defined on [0,T], we denote

    vn+12:=(vn+1+vn)/2andDtvn=3vn4vn1+vn22τ. (4.1)

    Based on (4.1), we introduce the BDF2 fully discrete scheme in Algorithm 3.

    Algorithm 3 BDF2 algorithm
    Input: Evaluate u0hVh,p0hMh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      Step 1: Solve for (u1h,ξ1h,η1h)Vh×Mh×Wh such that
            2μ(ε(u12h),ε(vh))(ξ12h,divvh)=(f,vh)+f1,vh,vhVh,(4.2)
            κ3(ξ12h,φh)+(divu12h,φh)=κ1(η12h,φh),φhMh,(4.3)
            (dtη1h,ψh)+1μf(K(κ1ξ12h+κ2η12h),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh.(4.4)
      Step 2:
      for n=1,2,, do
        (i) Solve for (un+1h,ξn+1h,ηn+1h)Vh×Mh×Wh such that
            2μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,(4.5)
            κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηn+1h,φh),φhMh,(4.6)
            (Dtηn+1h,ψh)+1μf(K(κ1ξn+1h+κ2ηn+1h),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh.(4.7)
        (ii) Update pn+1h and qn+1h by
            pn+1h=κ1ξn+1h+κ2ηn+1h,qn+1h=κ1ηn+1hκ3ξn+1h,(4.8)
      end for

    The next lemma establishes a discrete energy inequality for the BDF2 algorithm.

    Lemma 4.1. Let {(unh,ξnh,ηnh)}n0 be defined by the BDF2 algorithm. Then it holds

    μ8ε(ul+1h)2L2(Ω)+κ34ξl+1h2L2(Ω)+κ24ηl+1h2L2(Ω)+τ2μfln=1K12(κ1ξn+1h+κ2ηn+1h)2L2(Ω)C3(f2L2(Ω)+f12L2(Ω)+ε(u0h)2L2(Ω)+ξ0h2L2(Ω)+η0h2L2(Ω)+ϕ2L2(Ω)+ϕ12L2(Ω)),l1, (4.9)

    where C3=C3(ck,ca,cb,cp,μ,κ3,κ2,T,K,μf).

    Proof. Setting vh=dtu1h in (4.2), φh=dtξ1h in (4.3) and ψh=p12h in (4.4), we get

    2μ(ε(u12h),ε(dtu1h))(ξ12h,divdtu1h)=(f,dtu1h)+f1,dtu1h, (4.10)
    κ3(ξ12h,dtξ1h)+(divu12h,dtξ1h)=κ1(η12h,dtξ1h), (4.11)
    (dtη1h,p12h)+1μf(K(κ1ξ12h+κ2η12h),p12h)=(ϕ,p12h)+ϕ1,p12h. (4.12)

    Adding (4.10)–(4.12), there holds

    μ2τ(ε(u1h)2L2(Ω)ε(u0h)2L2(Ω))+κ32τ(ξ1h2L2(Ω)ξ0h2L2(Ω))+κ22τ(η1h2L2(Ω)η0h2L2(Ω))+1μfK12(κ1ξ12h+κ2η12h)2L2(Ω)=(f,dtu1h)+f1,dtu1h+(ϕ,p12h)+ϕ1,p12h+1τ((ξ0h,divu1h)(ξ1h,divu0h))+κ1τ((η0h,ξ1h)(η1h,ξ0h)). (4.13)

    Using the Cauchy-Schwarz inequality, triangle inequality, (2.2), (4.32), trace inequality, Poincaré inequality, and Young inequality, we can bound the right-side terms of (4.13) as follows:

    (f,dtu1h)1τ(5c2kμf2L2(Ω)+μ16ε(u1h)2L2(Ω)+μ4ε(u0h)2L2(Ω)), (4.14)
    f1,dtu1h1τ(5c2ac2kμf12L2(Ω)+μ16ε(u1h)2L2(Ω)+μ4ε(u0h)2L2(Ω)), (4.15)
    (ϕ,p12h)c2pμfK12ϕ2L2(Ω)+14μfK12p12h2L2(Ω), (4.16)
    ϕ1,p12hc2ac2pμfK12ϕ12L2(Ω)+14μfK12p12h2L2(Ω), (4.17)
    (ξ0h,divu1h)(ξ1h,divu0h)2c2bμξ0h2L2(Ω)+μ8ε(u1h)2L2(Ω)+κ38ξ1h2L2(Ω) (4.18)
    +2c2bκ3ε(u0h)2L2(Ω),κ1(η0h,ξ1h)κ1(η1h,ξ0h)2κ21κ3η0h2L2(Ω)+κ38ξ1h2L2(Ω)+κ24η1h2L2(Ω)+κ21κ2ξ0h2L2(Ω), (4.19)

    where ca is a positive constant related to Trace inequality.

    Substituting (4.14)–(4.19) into (4.13), we have

    μ4ε(u1h)2L2(Ω)+κ34ξ1h2L2(Ω)+κ24η1h2L2(Ω)+τ2μfK12(κ1ξ12h+κ2η12h)2L2(Ω)5c2kμf2L2(Ω)+5c2ac2kμf12L2(Ω)+c2pμfτK12ϕ2L2(Ω)+c2ac2pμfτK12ϕ12L2(Ω)+(μ+2c2bκ3)ε(u0h)2L2(Ω)+(κ32+2c2bμ+κ21κ2)ξ0h2L2(Ω)+(κ22+2κ21κ3)η0h2L2(Ω). (4.20)

    Next, setting vh=Dtun+1h,φh=ξn+1h and ψh=pn+1h in (4.5)–(4.7), we obtain

    2μ(ε(un+1h),ε(Dtun+1h))(ξn+1h,divDtun+1h)=(f,Dtun+1h)+f1,Dtun+1h, (4.21)
    κ3(Dtξn+1h,ξn+1h)+(divDtun+1h,ξn+1h)=κ1(Dtηn+1h,ξn+1h), (4.22)
    (Dtηn+1h,pn+1h)+1μf(K(κ1ξn+1h+κ2ηn+1h),pn+1h)=(ϕ,pn+1h)+ϕ1,pn+1h. (4.23)

    Adding (4.21)–(4.23) and using the fact of

    (Dtvn+1h,vn+1h)=14τ(vn+1h2L2(Ω)vnh2L2(Ω)+2vn+1hvnh2L2(Ω)2vnhvn1h2L2(Ω)+vn+1h2vnh+vn1h2L2(Ω)), (4.24)

    we have

    \begin{align} &\frac{\mu}{4\tau}\big(\|\varepsilon(\mathbf{u}_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(\mathbf{u}_{h}^{n})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(\mathbf{u}_{h}^{n+1})-\varepsilon(\mathbf{u}_{h}^{n})\|_{L^{2}(\Omega)}^{2}-\|2\varepsilon(\mathbf{u}_{h}^{n})-\varepsilon(\mathbf{u}_{h}^{n-1})\|_{L^{2}(\Omega)}^{2} \\ &+\|\varepsilon(\mathbf{u}_{h}^{n+1})-2\varepsilon(\mathbf{u}_{h}^{n})+\varepsilon(\mathbf{u}_{h}^{n-1})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4\tau}\big(\|\xi_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|\xi_{h}^{n}\|_{L^{2}(\Omega)}^{2}+\|2\xi_{h}^{n+1}-\xi_{h}^{n}\|_{L^{2}(\Omega)}^{2}\\ &-\|2\xi_{h}^{n}-\xi_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|\xi_{h}^{n+1}-2\xi_{h}^{n}+\xi_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{2}}{4\tau}\big(\|\eta_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|\eta_{h}^{n}\|_{L^{2}(\Omega)}^{2}+\|2\eta_{h}^{n+1}-\eta_{h}^{n}\|_{L^{2}(\Omega)}^{2}\\ &-\|2\eta_{h}^{n}-\eta_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|\eta_{h}^{n+1}-2\eta_{h}^{n}+\eta_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)+\frac{1}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}\\ & = (\mathbf{f}, D_{t}\mathbf{u}_{h}^{n+1})+\langle \mathbf{f}_1,D_{t}\mathbf{u}_{h}^{n+1}\rangle+(\phi, p_{h}^{n+1})+\langle \phi_1,p_{h}^{n+1}\rangle. \end{align} (4.25)

    Applying the summation operator \tau\sum_{n = 1}^{l} to both sides of (4.25), there holds

    \begin{align} &\frac{\mu}{4}\big(\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(\mathbf{u}_{h}^{l+1})-\varepsilon(\mathbf{u}_{h}^{l})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|\xi_{h}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2\xi_{h}^{l+1}-\xi_{h}^{l}\|_{L^{2}(\Omega)}^{2}\big) \\ &+\frac{\kappa_{2}}{4}\big(\|\eta_{h}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2\eta_{h}^{l+1}-\eta_{h}^{l}\|_{L^{2}(\Omega)}^{2}\big)+\sum\limits_{n = 1}^{l}\big(\frac{\mu}{2}\|\varepsilon(\mathbf{u}_{h}^{n+1})-2\varepsilon(\mathbf{u}_{h}^{n})+\varepsilon(\mathbf{u}_{h}^{n-1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\\ &\cdot\|\xi_{h}^{n+1}-2\xi_{h}^{n}+\xi_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|\eta_{h}^{n+1}-2\eta_{h}^{n}+\eta_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}\big)\\ & = \big(\mathbf{f}, \frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\big)+\langle \mathbf{f}_1,\frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\rangle+\tau\sum\limits_{n = 1}^{l}\left[(\phi, p_{h}^{n+1})+\langle \phi_1,p_{h}^{n+1}\rangle\right]\\ &+\frac{\mu}{4}\big(\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(\mathbf{u}_{h}^{1})-\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|\xi_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\|2\xi_{h}^{1}-\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}\big)\\ &+\frac{\kappa_{2}}{4}\big(\|\eta_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\|2\eta_{h}^{1}-\eta_{h}^{0}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.26)

    Similar to (4.13), we can bound the first four terms of the right side of (4.26) as follows:

    \begin{align} \big(\mathbf{f}, \frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\big)\leq&\frac{4c_{k}^{2}}{\mu}\|\mathbf{f}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{16}\|2\varepsilon(\mathbf{u}_{h}^{l+1})-\varepsilon(\mathbf{u}_{h}^{l})\|_{L^{2}(\Omega)}^{2} \end{align} (4.27)
    \begin{align} &+\frac{\mu}{16}\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{3\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2},\\ \langle \mathbf{f}_1,\frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\rangle\leq&\frac{4c_{a}^{2}c_{k}^{2}}{\mu}\|\mathbf{f}_{1}\|_{L^{2}(\partial\Omega)}^{2}+\frac{\mu}{16}\|2\varepsilon(\mathbf{u}_{h}^{l+1})-\varepsilon(\mathbf{u}_{h}^{l})\|_{L^{2}(\Omega)}^{2} \end{align} (4.28)
    \begin{align} &+\frac{\mu}{16}\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{3\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2},\\ \tau\sum\limits_{n = 1}^{l}(\phi, p_{h}^{n+1})\leq&\tau\sum\limits_{n = 1}^{l}\big( c_{p}^{2}\mu_{f}\|K^{-\frac{1}{2}}\phi\|_{L^{2}(\Omega)}^{2}+ \frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla p_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}\big), \end{align} (4.29)
    \begin{align} \tau\sum\limits_{n = 1}^{l}\langle \phi_1,p_{h}^{\frac{1}{2}}\rangle\leq& \tau\sum\limits_{n = 1}^{l}\big(c_{a}^{2}c_{p}^{2}\mu_{f}\|K^{-\frac{1}{2}}\phi_{1}\|_{L^{2}(\partial\Omega)}^{2}+\frac{1}{4\mu_{f}} \|K^{\frac{1}{2}}\nabla p_{h}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.30)

    Substituting (4.27)–(4.30) into (4.26) and using (4.20) for the resulting inequality, we can deduce that (4.9) holds. The proof is complete.

    In this subsection, we consider the optimal convergence analysis of the BDF2 algorithm. Before that, we have the following estimate for the BDF2 algorithm.

    Theorem 6. Let \{(\mathbf{u}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0} be defined by the BDF2 algorithm. Then there exists a positive constant C_{4} independent of mesh size h and time step \tau such that

    \begin{align} &\frac{\mu}{8}\|\varepsilon(Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2} \\ &+\frac{\tau}{2\mu_{f}}\sum\limits_{n = 1}^{l}\|K^{\frac{1}{2}} (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2}\leq C_{4}(h^{4}+\tau^{4}). \end{align} (4.31)

    Proof. First, it is easy to check that there exists a positive constant c_{b} such that

    \begin{align} \left\|\mathrm{div}\mathbf{w} \right\|_{L^{2}(\Omega)}\leq c_{b}\left\|\varepsilon (\mathbf{w}) \right\|_{L^{2}(\Omega)},\; \forall\; \mathbf{w}\in \mathbf{H}^{1}(\Omega). \end{align} (4.32)

    For n\geq 1 , at t = t_{n+1} , we have for (2.7)–(2.10) that

    \begin{align} {2} &\mu\bigl(\varepsilon(\mathbf{u}(t_{n+1})), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( \xi(t_{n+1}), {\rm div} \mathbf{v}_h \bigr) = (\mathbf{f}, \mathbf{v}_h)+\langle \mathbf{f}_1,\mathbf{v}_h\rangle, &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.33)
    \begin{align} &\kappa_3\bigl(\xi(t_{n+1}), \varphi_h \bigr) +\bigl({\rm div}\mathbf{u}(t_{n+1}), \varphi_h \bigr) = \kappa_1\bigl( \eta(t_{n+1}), \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.34)
    \begin{align} &\bigl(\eta_{t}(t_{n+1}), \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi(t_{n+1})+\kappa_{2}\eta(t_{n+1})),\nabla\psi_h \bigr) = (\phi, \psi_h)+\langle \phi_1,\psi_h\rangle, &&\quad \forall \psi_h\in W_h. \end{align} (4.35)

    Moreover, we have from (2.7)–(2.10) at t = t_{0} and t = t_{1} that

    \begin{align} {2} &\mu\bigl(\varepsilon(\mathbf{u}(t_{\frac{1}{2}})), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( \xi(t_{\frac{1}{2}}), {\rm div} \mathbf{v}_h \bigr) = (\mathbf{f}, \mathbf{v}_h)+\langle \mathbf{f}_1,\mathbf{v}_h\rangle, &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.36)
    \begin{align} &\kappa_3\bigl(\xi(t_{\frac{1}{2}}), \varphi_h \bigr) +\bigl({\rm div}\mathbf{u}(t_{\frac{1}{2}}), \varphi_h \bigr) = \kappa_1\bigl( \eta(t_{\frac{1}{2}}), \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.37)
    \begin{align} &\bigl(\eta_{t}(t_{\frac{1}{2}}), \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi(t_{\frac{1}{2}})+\kappa_{2}\eta(t_{\frac{1}{2}})),\nabla\psi_h \bigr) = (\phi, \psi_h)+\langle \phi_1,\psi_h\rangle, &&\quad \forall \psi_h\in W_h, \end{align} (4.38)

    where v(t_{\frac{1}{2}}) = (v(t_{0})+v(t_{1}))/2 .

    Subtracting (4.2)–(4.4) from (4.36)–(4.38), there holds

    \begin{align} {2} &\mu\bigl(\varepsilon(E_{\mathbf{u}}^{\frac{1}{2}}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( E_{\xi}^{\frac{1}{2}}, {\rm div} \mathbf{v}_h \bigr) = 0, &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.39)
    \begin{align} &\kappa_3\bigl(E_{\xi}^{\frac{1}{2}}, \varphi_h \bigr) +\bigl({\rm div}E_{\mathbf{u}}^{\frac{1}{2}}, \varphi_h \bigr) = \kappa_1\bigl( E_{\eta}^{\frac{1}{2}}, \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.40)
    \begin{align} &\bigl(d_{t}E_{\eta}^{1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}E_{\xi}^{\frac{1}{2}}+\kappa_{2}E_{\eta}^{\frac{1}{2}}),\nabla\psi_h \bigr) = (R_{1,\eta},\psi_{h})+(R_{2,\eta},\psi_{h}), &&\quad \forall \psi_h\in W_h, \end{align} (4.41)

    where

    \begin{align*} &R_{1,\eta}: = \frac{1}{2\tau}\int_{t_{0}}^{t_{1}}\Gamma^{2}(t)\eta_{ttt}dt,\qquad R_{2,\eta}: = \frac{1}{2}\int_{t_{0}}^{t_{1}}\Gamma(t)\eta_{ttt}dt,\\ &\Gamma(t): = \left\lbrace \begin{aligned} &t_{1}-t &&{\rm when}\; t\geq \frac{t_{0}+t_{1}}{2},\\ &t-t_{0}&&{\rm when}\; t < \frac{t_{0}+t_{1}}{2}. \end{aligned}\right. \end{align*}

    Using the definitions of the projection operators \mathcal{Q}_{h}, \mathcal{S}_{h}, \mathcal{R}_{h} for (4.39)–(4.41), we have

    \begin{align} {2} &\mu\bigl(\varepsilon(Z_{\mathbf{u}}^{\frac{1}{2}}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( G_{\xi}^{\frac{1}{2}}, {\rm div} \mathbf{v}_h \bigr) = \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} \mathbf{v}_h \bigr), &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.42)
    \begin{align} &\kappa_3\bigl(G_{\xi}^{\frac{1}{2}}, \varphi_h \bigr) +\bigl({\rm div}Z_{\mathbf{u}}^{\frac{1}{2}}, \varphi_h \bigr) = \kappa_1\bigl( G_{\eta}^{\frac{1}{2}}, \varphi_h \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.43)
    \begin{align} &\bigl(d_{t}G_{\eta}^{1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}}),\nabla\psi_h \bigr) = (R_{1,\eta},\psi_{h})+(R_{2,\eta},\psi_{h}), &&\quad \forall \psi_h\in W_h. \end{align} (4.44)

    Setting \mathbf{v}_{h} = d_{t}Z_{\mathbf{u}}^{1} in (4.42), \varphi_{h} = d_{t}G_{\xi}^{1} in (4.43) and \psi_{h} = Z_{p}^{\frac{1}{2}} in (4.44), we get

    \begin{align} {2} &\frac{\mu}{2\tau}\big(\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}\big)-\bigl( G_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr) = \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr), && \end{align} (4.45)
    \begin{align} &\frac{\kappa_{3}}{2\tau}\big(\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}-\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}\big) +\bigl({\rm div}Z_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr) = \kappa_1\bigl( G_{\eta}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr), && \end{align} (4.46)
    \begin{align} &\frac{\kappa_{2}}{2\tau}\big(\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}-\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}\big)+\kappa_{1}\bigl(d_{t}G_{\eta}^{\frac{1}{2}}, G_{\xi}^{\frac{1}{2}} \bigr) +\frac{1}{\mu_f} \|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2}\\ = &(R_{1,\eta},Z_{p}^{\frac{1}{2}})+(R_{2,\eta},Z_{p}^{\frac{1}{2}})+(d_{t}G_{\eta}^{1},Y_{p}^{\frac{1}{2}}-F_{p}^{\frac{1}{2}}). \end{align} (4.47)

    Adding (4.45)–(4.47), there holds

    \begin{align} &\frac{\mu}{2\tau}\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{2\tau}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{2\tau}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{1}{\mu_f} \|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2} \\ & = \frac{\mu}{2\tau}\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{2\tau}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{2\tau}\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}+ \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr)\\ &\quad+(R_{1,\eta},Z_{p}^{\frac{1}{2}})+(R_{2,\eta},Z_{p}^{\frac{1}{2}})+(d_{t}G_{\eta}^{1},Y_{p}^{\frac{1}{2}}-F_{p}^{\frac{1}{2}})+\frac{1}{\tau}\Big((G_{\xi}^{0},{\rm div}Z_{\mathbf{u}}^{1})-(G_{\xi}^{1},{\rm div}Z_{\mathbf{u}}^{0})\Big)\\ &\quad+\frac{\kappa_{1}}{\tau}\Big((G_{\eta}^{0},G_{\xi}^{1})-(G_{\eta}^{1},G_{\xi}^{0})\Big). \end{align} (4.48)

    Then, we bound the fourth to tenth terms on the right side of (4.48). Using the Cauchy-Schwarz inequality, triangle equality, and (4.32), Poincaré and Young inequalities, we have

    \begin{align} \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr) &\leq\frac{1}{2\tau}\big(\frac{2c_{b}^{2}}{\mu}\|F_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{2c_{b}^{2}}{\mu}\|F_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{2}\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2} \end{align} (4.49)
    \begin{align} &\quad+\frac{\mu}{2}\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}\big),\\ -\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr) &\leq\frac{1}{2\tau}\big(\frac{4c_{b}^{2}}{\kappa_{3}}\|\varepsilon(Y_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{4c_{b}^{2}}{\kappa_{3}}\|\varepsilon(Y_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.50)
    \begin{align} &\quad+\frac{\kappa_{3}}{4}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}\big),\\ (R_{1,\eta},Z_{p}^{\frac{1}{2}}) &\leq \frac{c_{p}^{2}\mu_{f}\tau^{3}}{640K_{1}}\|\eta_{ttt}\|_{L^{2}(t_{0},t_{1};L^{2}(\Omega))}^{2}+\frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla Z_{p}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.51)
    \begin{align} (R_{2,\eta},Z_{p}^{\frac{1}{2}}) &\leq \frac{c_{p}^{2}\mu_{f}\tau^{3}}{96K_1}\|\eta_{ttt}\|_{L^{2}(t_{0},t_{1};L^{2}(\Omega))}^{2}+\frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla Z_{p}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.52)
    \begin{align} (d_{t}G_{\eta}^{1},Y_{p}^{\frac{1}{2}}-F_{p}^{\frac{1}{2}}) &\leq\frac{1}{2\tau}\big(\frac{\kappa_{2}}{2}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{2}\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{4}{\kappa_{2}}\|Y_{p}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.53)
    \begin{align} &\quad+\frac{4}{\kappa_{2}}\|Y_{p}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{4}{\kappa_{2}}\|F_{p}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{4}{\kappa_{2}}\|F_{p}^{0}\|_{L^{2}(\Omega)}^{2}\big),\\ (G_{\xi}^{0},{\rm div}Z_{\mathbf{u}}^{1})-(G_{\xi}^{1},{\rm div}Z_{\mathbf{u}}^{0})&\leq\frac{2c_{b}^{2}}{\mu}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{8}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.54)
    \begin{align} &\quad+\frac{2c_{b}^{2}}{\kappa_{3}}\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2},\\ \kappa_{1}(G_{\eta}^{0},G_{\xi}^{1})-\kappa_{1}(G_{\eta}^{1},G_{\xi}^{0})&\leq\frac{2\kappa_{1}^{2}}{\kappa_{3}}\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{8}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{1}^{2}}{\kappa_{2}}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}. \end{align} (4.55)

    Substituting (4.49)–(4.55) into (4.48) and using (3.11)–(3.13), the facts of Z_{\mathbf{u}}^{0} = \mathbf{0}, \; G_{\xi}^{0} = 0 , and G_{\eta}^{0} = 0 , we find there exists a positive constant \widetilde{C} = \widetilde{C}\big(\|\mathbf{u}\|_{L^{\infty}(t_{0}, t_{1};H^{2}(\Omega))}, \|\xi\|_{L^{\infty}(t_{0}, t_{1};H^{2}(\Omega))}, \|\eta\|_{L^{\infty}(t_{0}, t_{1};H^{2}(\Omega))}, \\ \|\eta_{ttt}\|_{L^{2}(t_{0}, t_{1};L^{2}(\Omega))}, c_{b}, \mu, \kappa_{3}, c_{p}, \mu_{f}, K_{1}, \kappa_{2}, \kappa_{1}\big) such that

    \begin{align} &\mu\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\kappa_{3}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\kappa_{2}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{\mu_f} \|K\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2}\leq \widetilde{C}(h^{4}+\tau^{4}). \end{align} (4.56)

    Next, we show the error estimate at t = t_{n+1} . Subtracting (4.5)–(4.7) from (4.33)–(4.35), it holds

    \begin{align} {2} &\mu\bigl(\varepsilon(Z_{\mathbf{u}}^{n+1}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( G_{\xi}^{n+1}, {\rm div} \mathbf{v}_h \bigr) = \bigl( F_{\xi}^{n+1}, {\rm div} \mathbf{v}_h \bigr), &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.57)
    \begin{align} &\kappa_3\bigl(G_{\xi}^{n+1}, \varphi_h \bigr) +\bigl({\rm div}Z_{\mathbf{u}}^{n+1}, \varphi_h \bigr) = \kappa_1\bigl( G_{\eta}^{n+1}, \varphi_h \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{n+1}, \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.58)
    \begin{align} &\bigl(D_{t}G_{\eta}^{n+1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1}),\nabla\psi_h \bigr) = \bigl(R_{3,\eta}, \psi_h \bigr), &&\quad \forall \psi_h\in W_h, \end{align} (4.59)

    where

    \begin{align*} R_{3,\eta} = \frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_{n})^{2}\eta_{ttt}dt-\frac{1}{4\tau}\int_{t_{n-1}}^{t_{n+1}}(t-t_{n+1})^{2}\eta_{ttt}dt. \end{align*}

    Setting \mathbf{v}_{h} = D_{t}Z_{\mathbf{u}}^{n+1} in (4.57), \varphi_{h} = G_{\xi}^{n+1} in (4.58) after using the operator D_{t} to both sides of (4.58), \psi_{h} = Z_{p}^{n+1} in (4.59), using (4.24) and adding the result equation, we get

    \begin{align} &\frac{\mu}{4\tau}\big(\|\varepsilon(Z_{\mathbf{u}}^{n+1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(Z_{\mathbf{u}}^{n})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(Z_{\mathbf{u}}^{n+1})-\varepsilon(Z_{\mathbf{u}}^{n})\|_{L^{2}(\Omega)}^{2}-\|2\varepsilon(Z_{\mathbf{u}}^{n})-\varepsilon(Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2} \\ &+\|\varepsilon(Z_{\mathbf{u}}^{n+1})-2\varepsilon(Z_{\mathbf{u}}^{n})+\varepsilon(Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4\tau}\big(\|G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|G_{\xi}^{n}\|_{L^{2}(\Omega)}^{2}+\|2G_{\xi}^{n+1}-G_{\xi}^{n}\|_{L^{2}(\Omega)}^{2}\\ &-\|2G_{\xi}^{n}-G_{\xi}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|G_{\xi}^{n+1}-2G_{\xi}^{n}+G_{\xi}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{2}}{4\tau}\big(\|G_{\eta}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2}\\ &+\|2G_{\eta}^{n+1}-G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2}-\|2G_{\eta}^{n}-G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|G_{\eta}^{n+1}-2G_{\eta}^{n}+G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)\\ &+\frac{1}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2} = \bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr)-\bigl({\rm div}D_{t}Y_{\mathbf{u}}^{n+1}, G_{\xi}^{n+1}\bigr)+\bigl(R_{3,\eta}, Z_{p}^{n+1} \bigr)\\ &+\big(D_{t}G_{\eta}^{n+1},Y_{p}^{n+1}-F_{p}^{n+1}\big). \end{align} (4.60)

    Applying the summation operator \tau\sum_{n = 1}^{l} to both sides of (4.60), there holds

    \begin{align} &\frac{\mu}{4}\big(\|\varepsilon(Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(Z_{\mathbf{u}}^{l+1})-\varepsilon(Z_{\mathbf{u}}^{l})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|G_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\xi}^{l+1}-G_{\xi}^{l}\|_{L^{2}(\Omega)}^{2}\big) \\ &+\frac{\kappa_{2}}{4}\big(\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\eta}^{l+1}-G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2}\big)+\sum\limits_{n = 1}^{l}\big(\frac{\mu}{2}\|\varepsilon(Z_{\mathbf{u}}^{n+1})-2\varepsilon(Z_{\mathbf{u}}^{n})+\varepsilon(Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\\ &\cdot\|G_{\xi}^{n+1}-2G_{\xi}^{n}+G_{\xi}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|G_{\eta}^{n+1}-2G_{\eta}^{n}+G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2}\big)\\ & = \tau\sum\limits_{n = 1}^{l}\left[\bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr)-\bigl({\rm div}D_{t}Y_{\mathbf{u}}^{n+1}, G_{\xi}^{n+1}\bigr)+\bigl(R_{3,\eta}, Z_{p}^{n+1} \bigr)+\big(D_{t}G_{\eta}^{n+1},Y_{p}^{n+1}-F_{p}^{n+1}\big)\right]\\ &+\frac{\mu}{4}\big(\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(Z_{\mathbf{u}}^{1})-\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\xi}^{1}-G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}\big)\\ &+\frac{\kappa_{2}}{4}\big(\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\eta}^{1}-G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.61)

    Now, we begin to estimate the first four terms on the right side of (4.61). For the first term, it is easy to check that

    \begin{align} \tau\sum\limits_{n = 1}^{l}\bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr) = &\frac{1}{2}\Big(\big(F_{\xi}^{l+1},{\rm div}(3Z_{\mathbf{u}}^{l+1}-Z_{\mathbf{u}}^{l})\big)-\big(F_{\xi}^{l+1}-F_{\xi}^{l},{\rm div}Z_{\mathbf{u}}^{l}\big)-\big(3F_{\xi}^{2},{\rm div}Z_{\mathbf{u}}^{1}\big) \\ &+\big(F_{\xi}^{3},{\rm div}Z_{\mathbf{u}}^{1}\big)\Big)+\tau\sum\limits_{n = 3}^{l}\big(\frac{F_{\xi}^{n+1}-4F_{\xi}^{n}+3F_{\xi}^{n-1}}{2\tau},{\rm div}Z_{\mathbf{u}}^{n-1}\big). \end{align} (4.62)

    Applying the Cauchy-Schwarz inequality, (4.32), triangle and Young inequalities to (4.62), we have

    \begin{align} \tau\sum\limits_{n = 1}^{l}\bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr)\leq&\frac{1}{2}\Big(\frac{4c_{b}^{2}}{\mu}\|F_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|2\varepsilon (Z_{\mathbf{u}}^{l+1})-\varepsilon (Z_{\mathbf{u}}^{l})\|_{L^{2}(\Omega)}^{2} \\ &+\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{10c_{b}^{2}}{\mu}\|F_{\xi}^{l+1}-F_{\xi}^{l}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}\\ &+\frac{\mu}{8}\|2\varepsilon (Z_{\mathbf{u}}^{l+1})-\varepsilon (Z_{\mathbf{u}}^{l})\|_{L^{2}(\Omega)}^{2}+\frac{c_{b}^{2}}{\mu}\|F_{\xi}^{3}-3F_{\xi}^{2}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{4}\|\varepsilon (Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}\Big)\\ &+\tau\sum\limits_{n = 3}^{l}\Big(\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2}+\frac{2c_{b}^{2}}{\mu}\|\frac{F_{\xi}^{n+1}-4F_{\xi}^{n}+3F_{\xi}^{n-1}}{2\tau}\|_{L^{2}(\Omega)}^{2}\Big). \end{align} (4.63)

    For the second term, using the Canchy-Schwarz and Young inequalities, we obtain

    \begin{align} -\tau\sum\limits_{n = 1}^{l}\bigl({\rm div}D_{t}Y_{\mathbf{u}}^{n+1}, G_{\xi}^{n+1}\bigr)\leq\tau\sum\limits_{n = 1}^{l}\big(\frac{1}{\kappa_{3}}\|{\rm div}D_{t}Y_{\mathbf{u}}^{n+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.64)

    Applying the Canchy-Schwarz inequality, Poincaré inequality, Young inequality, and the fact of

    \begin{align*} \|R_{3,\eta}\|_{L^{2}(\Omega)}^{2}\leq\frac{\tau^{3}}{5}\|\eta_{ttt}\|_{L^{2}(t_{n},t_{n+1};L^{2}(\Omega))}^{2}+\frac{2\tau^{3}}{5}\|\eta_{ttt}\|_{L^{2}(t_{n-1},t_{n+1};L^{2}(\Omega))}^{2} \end{align*}

    to the third term, there holds

    \begin{align} \tau\sum\limits_{n = 1}^{l}\bigl(R_{3,\eta}, Z_{p}^{n+1} \bigr)&\leq\tau\sum\limits_{n = 1}^{l} \big(\frac{c_{p}^{2}\mu_{f}}{2}\|R_{3,\eta}\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\mu_{f}}\|\nabla Z_{p}^{n+1}\|_{L^{2}(\Omega)}^{2}\big)\\ &\leq\tau\sum\limits_{n = 1}^{l} \big(\frac{3\tau^{3}c_{p}^{2}\mu_{f}}{10K_{1}}\|\eta_{ttt}\|_{L^{2}(t_{n-1},t_{n+1};L^{2}(\Omega))}^{2}+\frac{1}{2\mu_{f}}\|K^{\frac{1}{2}}\nabla Z_{p}^{n+1}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.65)

    Similar to the first term, we can bound the fourth term as follows:

    \begin{align} \tau\sum\limits_{n = 1}^{l}\big(D_{t}G_{\eta}^{n+1},Y_{p}^{n+1}-F_{p}^{n+1}\big)\leq&\frac{1}{2}\Big(\frac{4}{\kappa_{2}}\|Y_{p}^{l+1}-F_{p}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|2G_{\eta}^{l+1}-G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2} \\ &+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{10}{\kappa_{2}}\|Y_{p}^{l+1}-Y_{p}^{l}-F_{p}^{l+1}+F_{p}^{l}\|_{L^{2}(\Omega)}^{2}\\ &+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|2G_{\eta}^{l+1}-G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}\\ &+\frac{2}{\kappa_{2}}\|Y_{p}^{3}-3Y_{p}^{2}-F_{p}^{3}+3F_{p}^{2}\|_{L^{2}(\Omega)}^{2}\Big)+\tau\sum\limits_{n = 3}^{l}\Big(\frac{\kappa_{2}}{8}\|G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}\\ &+\frac{2}{\kappa_{2}}\|\frac{Y_{p}^{n+1}-4Y_{p}^{n}+3Y_{p}^{n-1}}{2\tau}-\frac{F_{p}^{n+1}-4F_{p}^{n}+3F_{p}^{n-1}}{2\tau}\|_{L^{2}(\Omega)}^{2}\Big). \end{align} (4.66)

    Substituting (4.63)–(4.66) into (4.61), using the facts of Z_{\mathbf{u}}^{0} = \mathbf{0}, \; G_{\xi}^{0} = 0, \; G_{\eta}^{0} = 0 , (4.56), (3.11)–(3.13), there exists a positive constant \widehat{C} = \widehat{C}(c_{b}, \mu, \kappa_{3}, \mu_{f}, K_{1}, c_{p}, \kappa_{2}) such that

    \begin{align} &\frac{\mu}{8}\|\varepsilon(Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2} \\ &\quad+\frac{\tau}{2\mu_{f}}\sum\limits_{n = 1}^{l}\|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2}\\ &\leq \tau\sum\limits_{n = 1}^{l}\big(\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{n+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{n+1}\|_{L^{2}(\Omega)}^{2}\big)+ \widehat{C}\big(h^{4}\|\xi\|_{L^{\infty}(0,T;H^{2}(\Omega))}\\ &\quad+h^{4}\|\xi_{t}\|_{L^{2}(0,T;H^{2}(\Omega))}+h^{4}\|\nabla\mathbf{u}_{t}\|_{L^{2}(0,T;H^{2}(\Omega))}+\tau^{4}\|\eta_{ttt}\|_{L^{2}(0,T;L^{2}(\Omega))}+h^{4}\|\eta\|_{L^{\infty}(0,T;H^{2}(\Omega))}\\ &\quad+h^{4}\|\eta_{t}\|_{L^{2}(0,T;H^{2}(\Omega))}+\widetilde{C}(h^{4}+\tau^{4})\big). \end{align} (4.67)

    Applying the Gronwall inequality to (4.67), one has that (4.31) holds. The proof is complete.

    Furthermore, we state the main result of this part by a theorem.

    Theorem 7. There exists a positive constant C_{5} independent of h and \tau such that the solution of the BDF2 algorithm satisfies the following error estimates:

    \begin{align} &\max\limits_{0\leq n\leq N}\left[\|\nabla(\mathbf{u}(t_{n})-\mathbf{u}_{h}^{n})\|_{L^{2}(\Omega)}+\sqrt\kappa_{2}\|\eta(t_{n})-\eta_{h}^{n}\|_{L^{2}(\Omega )}\right. \end{align} (4.68)
    \begin{align} &\left.+\sqrt{\kappa_{3}}\|\xi(t_{n})-\xi_{h}^{n}\|_{L^{2}(\Omega )} \right]\leq C_{5}(h^{2}+\tau^{2}),\\ &\max\limits_{0\leq n\leq N}\|p(t_{n})-p_{h}^{n}\|_{L^{2}(\Omega )}\leq C_{5}(h^{2}+\tau^{2}). \end{align} (4.69)

    Proof. Using (3.11)–(3.13), the relation p = \kappa_{1}\xi+\kappa_{2}\eta , Theorem 6, (4.56), and the triangle inequality on

    \begin{align*} &\mathbf{u}(t_{n})-\mathbf{u}_{h}^{n} = Y_{\mathbf{u}}^{n}+Z_{\mathbf{u}}^{n},\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \xi(t_{n})-\xi_{h}^{n} = Y_{\xi}^{n}+Z_{\xi}^{n} = F_{\xi}^{n}+G_{\xi}^{n},\\ &\eta(t_{n})-\eta_{h}^{n} = Y_{\eta}^{n}+Z_{\eta}^{n} = F_{\eta}^{n}+G_{\eta}^{n},\; \; \; \; \; \; \; p(t_{n})-p_{h}^{n} = Y_{p}^{n}+Z_{p}^{n} = F_{p}^{n}+G_{p}^{n}, \end{align*}

    one can see that (4.68) and (4.69) hold. The proof is complete.

    Here the \mathbf{P}_{2}-P_{1}-P_{1}-P_{1} element pairs are employed for the coupled algorithm and decoupled algorithms, and \mathbf{P}_{2}-P_{1}-P_{1} element pairs are employed for the BDF2 algorithm. To verify the robustness of our algorithms with respect to parameters, we compute the spatial errors and convergence rates of the variables \mathbf{u} and p for different values of Poisson rates and permeability tensor separately. In the meantime, we show the effectiveness of our algorithms in overcoming the Poisson locking as \nu\rightarrow 0.5 . After that, to illustrate that there is no nonphysical oscillation in the pressure field, we use the decoupled algorithm to test a Barry-Mercer problem. Finally, we use the following examples to compute the convergence rates for the above proposed algorithms in subsections 5.1, 5.2, and 5.4.

    Example 1. Let \Omega = (0, 1)^2 , \Gamma_1 = \{(1, x_2); 0\leq x_2\leq 1\} , \Gamma_2 = \{(x_1, 0); 0\leq x_1\leq 1\} , \Gamma_3 = \{(0, x_2); 0\leq x_2\leq 1\} , \Gamma_4 = \{(x_1, 1); 0\leq x_1\leq 1\} and T = 1 . The boundary and initial conditions are

    \begin{alignat*} {2} p = t\sin(\pi x_{1})\cos(\pi x_{2}) &\qquad\mbox{on }\partial\Omega_T,\\ u_1 = e^{-t}\sin(2\pi x_{2})(-1+\cos(2\pi x_{1}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}) &\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 1,3,\\ u_2 = e^{-t}\sin(2\pi x_{1})(1-\cos(2\pi x_{2}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}) &\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 2,4,\\ \sigma({\mathbf{u}})\mathbf{n}-\alpha p \mathbf{n} = \mathbf{f}_1 &\qquad \mbox{on } \partial\Omega_T,\\ \mathbf{u}(\mathbf{x},0) = \mathbf{0}, \quad p(\mathbf{x},0) = 0.& \end{alignat*}

    The exact solutions of (1.1)–(1.4) and (2.3)–(2.5) are p = t\sin(\pi x_{1})\cos(\pi x_{2}) and

    \begin{alignat*} {2} u_1 & = e^{-t}\sin(2\pi x_{2})(-1+\cos(2\pi x_{1}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}),\\ u_2 & = e^{-t}\sin(2\pi x_{1})(1-\cos(2\pi x_{2}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}). \end{alignat*}

    We first test the robustness of the coupled and decoupled algorithms on the Poisson ratio \nu by fixing \alpha = 1 , \mu_{f} = 1 , E = 2.8\text{e}3 , K = 1 e-7 \mathbf{I} , and c_{0} = 0.2 for \nu = 0.49 and \nu = 0.4999999 , respectively.

    Tables 1 and 2 show the spatial errors and convergence rates of the coupled algorithm for \nu = 0.49 and \nu = 0.4999999 , respectively. Numerical results show that the H^{1} convergence rates of \mathbf{u} and the L^{2} convergence rates of p are all approaching 2 . The L^{2} convergence rates of \mathbf{u} are around 3 . The H^{1} convergence rates of p are around 1 . It can be concluded that the convergence rates of \mathbf{u} and p are optimal, which is consistent with Theorem 4.

    Table 1.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm \nu = 0.49 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5755e-2 4.2637e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7262e-2 2.31 2.0452e-1 1.06
    h=1/16, \tau=1/256 6.5936e-4 3.10 1.1349e-2 1.96 3.8013e-3 2.18 9.9924e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 9.8977e-4 1.94 4.9462e-2 1.01

     | Show Table
    DownLoad: CSV
    Table 2.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm for \nu = 0.4999999 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7651e-2 1.6055e-1 8.5700e-2 4.2574e-1
    h=1/8, \tau=1/64 5.6669e-3 3.07 4.4078e-2 1.86 1.7255e-2 2.31 2.0459e-1 1.06
    h=1/16, \tau=1/256 6.5948e-4 3.10 1.1349e-2 1.96 3.7745e-3 2.19 9.9988e-2 1.03
    h=1/32, \tau=1/1024 7.9662e-5 3.05 2.8614e-3 1.99 8.7440e-4 2.11 4.9501e-2 1.01

     | Show Table
    DownLoad: CSV

    For the decoupled algorithm, we also take \nu = 0.49 and \nu = 0.4999999 separately to observe the convergence rates of \mathbf{u} and p . Tables 3 and 4 show that the convergence rates of all variables are optimal.

    Table 3.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for \nu = 0.49 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5755e-2 4.2637e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7262e-2 2.31 2.0452e-1 1.06
    h=1/16, \tau=1/256 6.5935e-4 3.10 1.1349e-2 1.96 3.8013e-3 2.18 9.9924e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 9.8977e-4 1.94 4.9462e-2 1.01

     | Show Table
    DownLoad: CSV
    Table 4.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for \nu = 0.4999999 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7651e-2 1.6055e-1 8.5700e-2 4.2574e-1
    h=1/8, \tau=1/64 5.6669e-3 3.07 4.4078e-2 1.86 1.7255e-2 2.31 2.0459e-1 1.06
    h=1/16, \tau=1/256 6.5948e-4 3.10 1.1349e-2 1.96 3.7745e-3 2.19 9.9988e-2 1.03
    h=1/32, \tau=1/1024 7.9662e-5 3.05 2.8614e-3 1.99 8.7440e-4 2.11 4.9501e-2 1.01

     | Show Table
    DownLoad: CSV

    Thus, Tables 14 imply the coupled and decoupled algorithms are robust on the Poisson rate \nu .

    Table 5 reports the CPU time for both the coupled and decoupled algorithms. The results indicate that the decoupled algorithm significantly outperforms the coupled one in terms of runtime, demonstrating its superior efficiency.

    Table 5.  The CPU time of coupled and decoupled algorithms for \nu = 0.49 .
    h, \tau the CPU time of coupled algorithm the CPU time of decoupled algorithm
    h=1/4, \tau=1/16 0.19 s 0.13 s
    h=1/8, \tau=1/64 2.61 s 1.63 s
    h=1/16, \tau=1/256 51.82 s 30.32 s
    h=1/32, \tau=1/1024 915.23 s 503.59 s

     | Show Table
    DownLoad: CSV

    We now test the robustness of the algorithms with respect to the permeability tensor K by fixing \alpha = 1 , \mu_{f} = 1 , E = 2.8\text{e}3 , \nu = 0.49 , and c_{0} = 0.2 for K = 1\mathbf{I} and K = 1 e- 4\mathbf{I} , respectively.

    For the cases of K = 1\mathbf{I} and K = 1 e- 4\mathbf{I} , Tables 6 and 7 exhibit the spatial errors and convergence rates of \mathbf{u} and p for the coupled algorithm. The convergence rates of \mathbf{u} with respect to H^{1} norm and L^{2} norm are all around 3 and 2 separately. The convergence rates of p with respect to H^{1} norm and L^{2} norm are all around 2 and 1, respectively. This implies that the coupled algorithm is robust for K .

    Table 6.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm for K = 1\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 1.3174e-1 3.7941e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 3.5460e-2 1.89 1.9461e-1 0.96
    h=1/16, \tau=1/256 6.5936e-4 3.10 1.1349e-2 1.96 9.0399e-3 1.97 9.7954e-2 0.99
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 2.2712e-3 1.99 4.9060e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 7.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm for K = 1 e- 4\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5753e-2 4.2449e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7287e-2 2.31 2.0376e-1 1.06
    h=1/16, \tau=1/256 6.5936e-4 3.10 1.1349e-2 1.96 3.8309e-3 2.17 9.9511e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 1.0018e-3 1.93 4.9279e-2 1.01

     | Show Table
    DownLoad: CSV

    Again, for the cases of K = 1\mathbf{I} and K = 1 e- 4\mathbf{I} , we compute the convergence rates of \mathbf{u} and p by the decoupled algorithm. Tables 8 and 9 show that the convergence rates of \mathbf{u} and p are optimal, which is consistent with the results in Theorem 5 and shows that the decoupled algorithm is robust for K .

    Table 8.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for K = 1\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 1.3174e-1 3.7941e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 3.5460e-2 1.89 1.9461e-1 0.96
    h=1/16, \tau=1/256 6.5935e-4 3.10 1.1349e-2 1.96 9.0399e-3 1.97 9.7954e-2 0.99
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 2.2712e-3 1.99 4.9060e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 9.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for K = 1 e- 4\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5753e-2 4.2449e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7287e-2 2.31 2.0376e-1 1.06
    h=1/16, \tau=1/256 6.5935e-4 3.10 1.1349e-2 1.96 3.8309e-3 2.17 9.9511e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 1.0018e-3 1.94 4.9279e-2 1.01

     | Show Table
    DownLoad: CSV

    We now verify that there is no nonphysical oscillation in the pressure field for the proposed algorithms. For brevity, we only consider the decoupled algorithm. The benchmark problem for the poroelasticity model is given as follows:

    Example 2. (Barry-Mercer problem) Again we take \Omega = (0, 1)^2 and T = 1 . The boundary segments \Gamma_{j}(j = 1, 2, 3, 4) are defined in Example 1. The source functions \mathbf{f} = 0 and \phi = 0 , and the boundary and initial conditions are

    \begin{alignat*} {2} p & = 0 &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 1,3,4,\\ p& = p_{1} &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 2,\\ u_1 & = 0 &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 2,4,\\ u_2 & = 0 &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 1,3,\\ \sigma(\mathbf{u})\mathbf{n}-\alpha p \mathbf{n} & = \mathbf{f}_1: = (0,\alpha p)^{T} &&\qquad \mbox{on } \partial\Omega_T,\\ \mathbf{u}(\mathbf{x},0) = \mathbf{0}, \quad p(\mathbf{x},0) & = 0 &&\qquad\mbox{in } \Omega, \end{alignat*}

    where

    \begin{align*} &p_{1} = \left\lbrace \begin{aligned} &\sin t &&{\rm when}\; ( x_{1},t)\in[0.2,0.8)\times(0,T),\\ &0&&{\rm others}. \end{aligned}\right. \end{align*}

    Set \alpha = 1 , \mu_{f} = 1 , E = 2.8 e 3 , \nu = 0.4 , K = 1 e- 8\mathbf{I} , and c_{0} = 0 . For the case of c_{0} = 0 , low permeability K , and small time step \tau , it will lead to \mathrm{div}\mathbf{u}^{1}\approx 0 for the discrete form of Eq (1.2). If using \mathbf{P}_2-P_1 element pairs to solve the original problem, there exists the nonphysical oscillation in the pressure field at a very short time see Figure 1. But Figure 2 shows that there is no nonphysical oscillation in the pressure field for the decoupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs. In addition, the arrows near the boundary match very well with the inner arrows in Figure 3. This implies that the decoupled algorithm is locking-free and stable.

    Figure 1.  The numerical solution p_{h}^{n+1} by using the \mathbf{P}_2-P_1 element pairs for the original problem at t = 0.00001 .
    Figure 2.  The numerical solution p_{h}^{n+1} by using the decoupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs at t = 0.00001 .
    Figure 3.  Arrow plot of numerical solution \mathbf{u}_h by using the decoupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs at t = 0.00001 .

    We here test the convergence rates for the BDF2 algorithm by fixing \alpha = 1 , \mu_{f} = 1 , E = 2.8 e 3 , K = 1 e- 7\mathbf{I} and c_{0} = 0.2 for \nu = 0.49 and \nu = 0.4999999 , respectively.

    Tables 10 and 11 show the spatial errors and convergence rates of the BDF2 algorithm for \nu = 0.49 and \nu = 0.4999999 , respectively. The numerical results show that the H^{1} convergence rates of \mathbf{u} and the L^{2} convergence rates of p are all approaching 2 . The L^{2} convergence rates of \mathbf{u} are around 3 . The H^{1} convergence rates of p are around 1 . This implies that the BDF2 algorithm has the optimal convergence rates for the different values of \nu .

    Table 10.  Spatial errors and convergence rates of \mathbf{u} and p of the BDF2 algorithm for \nu = 0.49 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/4 4.8482e-2 1.6013e-1 5.7104e-2 4.0384e-1
    h=1/8, \tau=1/8 5.7075e-3 3.09 4.4041e-2 1.86 1.3284e-2 2.10 1.9872e-1 1.02
    h=1/16, \tau=1/16 6.6098e-4 3.11 1.1349e-2 1.96 3.2745e-3 2.02 9.8508e-2 1.01
    h=1/32, \tau=1/32 7.9759e-5 3.05 2.8615e-3 1.99 9.3353e-4 1.81 4.9128e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 11.  Spatial errors and convergence rates of \mathbf{u} and p of the BDF2 algorithm for \nu = 0.4999999 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/4 4.8493e-2 1.6012e-1 5.7040e-2 4.0369e-1
    h=1/8, \tau=1/8 5.7096e-3 3.09 4.4040e-2 1.86 1.3273e-2 2.10 1.9874e-1 1.02
    h=1/16, \tau=1/16 6.6116e-4 3.11 1.1349e-2 1.96 3.2409e-3 2.03 9.8517e-2 1.01
    h=1/32, \tau=1/32 7.9771e-5 3.05 2.8615e-3 1.99 8.0517e-4 2.01 4.9130e-2 1.00

     | Show Table
    DownLoad: CSV

    In this section, the coupled and decoupled algorithms with \mathbf{P}_2-P_1-P_1-P_1 element pairs and BDF2 algorithm with \mathbf{P}_2-P_1-P_1 element pairs are used to simulate brain swelling caused by brain injury. Since the coupled and decoupled algorithms present almost the same result, we only report the simulation result based on the coupled algorithm. For more details about \Omega, \; \Gamma = \Gamma_{1}+\Gamma_{2} , source functions, boundary conditions, and the values of parameters, one can see [15]. Here we list the required information, which is obtained from [15]. The \Omega is a 2D cross-section of a 3D model of the brain (see Figure 4), and the values of parameters are shown in Table 12. The injured area is marked with a yellow line in the right-hand figure of Figure 4.

    Figure 4.  A slice of the magnetic resonance imaging for a human brain (left) and the finite-element mesh (right).
    Table 12.  Values of parameters.
    Parameters Values Parameters Values
    c_{0} 4.5e-7 {\rm Pa}^{-1} K 1.4e-9 {\rm mm}^{2}
    c_{b} 3e-5 mm/min/Pa \alpha 1
    p_{\text{ SAS}} 1070 Pa \nu 0.35
    \mu_{f} 1.48e-5 {\rm Pa\cdot min} E 9010 Pa

     | Show Table
    DownLoad: CSV

    Moreover, the source functions are \mathbf{f} = \mathbf{0}, \; \phi = 9 e- 3\; \text{ml}/{\rm min} , and the boundary conditions are

    \begin{alignat*} {2} \mathbf{u} & = \mathbf{0} &&\qquad\mbox{on }\Gamma_1,\\ (K\nabla p)\cdot\mathbf{n}& = c_{b}(p_{\text{ SAS}}-p) &&\qquad\mbox{on }\Gamma_1,\\ p & = 1100\; {\rm Pa} &&\qquad\mbox{on }\Gamma_{2},\\ (\sigma-\alpha p\mathbf{I})\mathbf{n} & = -p\mathbf{n} &&\qquad \mbox{on } \Gamma_{2}. \end{alignat*}

    Figures 5 and 6 show the simulation results of the coupled and BDF2 algorithms for an injured brain with respect to \nu = 0.35, \; 0.49 . From Figure 5, one can see that the maximum pressures of the coupled algorithm are 3025.17 Pa and 3025.69 Pa separately, which are consistent with that in [27]. And Figure 6 illustrates the similar facts for the BDF2 algorithm by the values of maximum pressures ( 3028.13 Pa and 3025.87 Pa). In addition, the brain tissue can deform due to the influence of the total stress. In our simulation, the values of maximum tissue deformation are 0.99342 mm and 0.048082 mm for the coupled algorithm and 1.0018 mm and 0.048174 mm for the BDF2 algorithm, respectively. It is easy to find that the change of \nu has little influence on the maximum pressure, and the maximum tissue deformation becomes smaller as \nu gets bigger. When the Poisson ratio is approaching 0.5 , the brain tissue is nearly incompressible such that the tissue deformation is very small. Moreover, from Table 13, the execution time of the CPU of the coupled algorithm is longer than the BDF2 algorithm; thus the latter has higher computational efficiency.

    Figure 5.  Simulation results of the coupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs.
    Figure 6.  Simulation results of the BDF2 algorithm with \mathbf{P}_2-P_1-P_1 element pairs.
    Table 13.  Execution time of CPU.
    \nu Coupled algorithm BDF2 algorithm
    0.35 1547.65s 1239.13s
    0.496 126.97s 87.21s

     | Show Table
    DownLoad: CSV

    In this paper, we transform the poroelasticity model into a four-field problem by introducing the new variables and discuss the inf-sup condition about the four-field problem. Besides, we design the coupled and decoupled algorithms for the four-field problem and perform the error estimates. Optimal convergence is obtained in space and time. It can be observed that there is no constraint on mesh size and time step in error estimates. Then we introduce the BDF2 algorithm and perform stability analysis and error estimates. Furthermore, some numerical examples are presented to verify the theoretical results. More importantly, we successfully apply the algorithms to simulate the brain edema and observe the influence on the maximum tissue deformation and maximum pressure.

    Peizhen Wang: Conceptualization, investigation, writing-original draft; Wenlong He: Validation; Yajuan Di: Writing-review, editing. All authors have read and approved the final version of the manuscript for publication.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work was supported by the National Natural Science Foundation of China (Grant No.11901197) and the Young Teacher Foundation of Henan Province (Grant No.2021GGJS080). The numerical simulations have been done on the supercomputing system in the Supercomputing Center of Wuhan University.

    The authors declare that they have no competing interests.



    [1] Amuedo-Dorantes C, de la Rica S (2013) The immigration surplus and the substitutability of immigrant and native labor: Evidence from Spain. Empir Econ 44: 945–958. https://doi.org/10.1007/s00181-011-0534-4 doi: 10.1007/s00181-011-0534-4
    [2] Atienza-Montero P, Romo-Calixto CE (2021) El efecto de la población inmigrante sobre el empleo y salarios de los trabajadores nativos. Una investigación empírica para España, 2009–2018. Estud Econ Apl 39: 133–3197. https://doi.org/10.25115/eea.v39i2.3706 doi: 10.25115/eea.v39i2.3706
    [3] Baltagi BH, Başkaya YS (2022) Spatial wage curves for formal and informal workers in Turkey. J Spat Econometrics 3. https://doi.org/10.1007/s43071-022-00021-y doi: 10.1007/s43071-022-00021-y
    [4] Blanco MR, Golik M (2024) Why moving there? Spanish SIEs: Factors and motivations involved in the choice of the host destination. J Glob Mobil 12: 520–544. https://doi.org/10.1108/JGM-08-2023-0059 doi: 10.1108/JGM-08-2023-0059
    [5] Blau FD, Kahn LM (2015) Immigration and the distribution of incomes, In: Handbook of the Economics of International Migration, 1: 793–843. https://doi.org/10.1016/B978-0-444-53768-3.00015-1
    [6] Borjas GJ (2003) The labor demand curve is downward sloping: Reexamining the impact of immigration on the labor market. Q J Econ 118: 1335–1374. https://doi.org/10.1162/003355303322552810 doi: 10.1162/003355303322552810
    [7] Borjas GJ (2016) The wage impact of the Marielitos: Additional evidence. Working Paper 21850, National Bureau of Economic Research.
    [8] Card D, Dustmann C, Preston I (2012) Immigration, wages, and compositional amenities. J Eur Econ Assoc 10: 78–119. https://doi.org/10.1111/j.1542-4774.2011.01051.x doi: 10.1111/j.1542-4774.2011.01051.x
    [9] Card D, Peri G (2016) Immigration economics by George J. Borjas: A review essay. J Econ Lit 54: 1333–1349. https://doi.org/10.1257/jel.20151248 doi: 10.1257/jel.20151248
    [10] Carrasco R, Jimeno JF, Ortega AC (2008) The effect of immigration on the labour market performance of native-born workers: Some evidence for Spain. J Popul Econ 21: 627–648. https://doi.org/10.1007/s00148-006-0112-9 doi: 10.1007/s00148-006-0112-9
    [11] Cavalleri MC, Luu N, Causa O (2021) Migration, housing and regional disparities: A gravity model of inter-regional migration with an application to selected OECD countries. OECD Economics Department Working Papers 1691, OECD Publishing. https://doi.org/10.1787/421bf4aa-en
    [12] Cheng W, Lee LF (2017) Testing endogeneity of spatial and social networks. Reg Sci Urban Econ 64: 81–97. https://doi.org/10.1016/j.regsciurbeco.2017.03.005 doi: 10.1016/j.regsciurbeco.2017.03.005
    [13] Cuéllar-Martín J, Martín-Román ÁL, Moral A (2019) An empirical analysis of natural and cyclical unemployment at the provincial level in Spain. Appl Spat Anal Policy 12: 647–696. https://doi.org/10.1007/s12061-018-9262-x doi: 10.1007/s12061-018-9262-x
    [14] Daouli J, Demoussis M, Giannakopoulos N, et al. (2017) The wage curve before and during the Greek economic crisis. Empir Econ 52: 59–77. https://doi.org/10.1007/s00181-016-1073-9 doi: 10.1007/s00181-016-1073-9
    [15] Devicienti F, Maida A, Pacelli L (2008) The resurrection of the Italian wage curve. Econ Lett 98: 335–341. https://doi.org/10.1016/j.econlet.2007.05.013 doi: 10.1016/j.econlet.2007.05.013
    [16] Dolado J, Felgueroso F, Jimeno J (2021) Past, present and future of the Spanish labour market: When the pandemic meets the megatrend. Appl Econ Anal 29: 21–41. https://doi.org/10.1108/AEA-11-2020-0154 doi: 10.1108/AEA-11-2020-0154
    [17] Domingo A, Bayona-i-Carrasco J (2024) Second Latin American migratory boom in Spain: From recovery to COVID-19. Migr Stud 12: 93–113. https://doi.org/10.1093/migration/mnad039 doi: 10.1093/migration/mnad039
    [18] Dustmann C, Frattini T, Preston IP (2013) The effect of immigration along the distribution of wages. Rev Econ Stud 80: 145–173. https://doi.org/10.1093/restud/rds019 doi: 10.1093/restud/rds019
    [19] Dustmann C, Preston I (2007) Racial and economic factors in attitudes to immigration. B.E. J Econ Anal Policy 7. https://doi.org/10.2202/1935-1682.1655 doi: 10.2202/1935-1682.1655
    [20] Dustmann C, Schönberg U, Stuhler J (2016) The impact of immigration: Why do studies reach such different results? J Econ Perspect 30: 31–56. https://doi.org/10.1257/jep.30.4.31 doi: 10.1257/jep.30.4.31
    [21] Elhorst JP, Blien U, Wolf K (2007) New evidence on the wage curve: A spatial panel approach. Int Reg Sci Rev 30: 173–191. https://doi.org/10.1177/0160017606298426 doi: 10.1177/0160017606298426
    [22] Espinosa AM, Díaz-Emparanza I (2021) The long-term relationship between international labour migration and unemployment in Spain. J Int Migr Integr 22: 145–166. https://doi.org/10.1007/s12134-019-00716-6 doi: 10.1007/s12134-019-00716-6
    [23] Facchini G, Lodigiani E (2014) Attracting skilled immigrants: An overview of recent policy developments in advanced countries. Natl Inst Econ Rev 229: R3–R21. https://doi.org/10.1177/002795011422900102 doi: 10.1177/002795011422900102
    [24] Felgueroso F, Hidalgo-Pérez M, Jiménez-Martín S (2016) The puzzling fall of the wage skill premium in Spain. Manch Sch 84: 390–435. https://doi.org/10.1111/manc.12116 doi: 10.1111/manc.12116
    [25] Fellini I (2018) Immigrants' labour market outcomes in Italy and Spain: Has the Southern European model disrupted during the crisis? Migr Stud 6: 53–78. https://doi.org/10.1093/migration/mnx029 doi: 10.1093/migration/mnx029
    [26] Fingleton B, Palombi P (2013) The wage curve reconsidered: Is it truly an 'Empirical law of economics'?. Région et Développement 38: 49–92.
    [27] Fundación CYD (2024) Informe CYD 2024. Available from: https://www.fundacioncyd.org/publicaciones-cyd/informe-cyd-2024/.
    [28] Ghosh D, Dickey H (2024) The wage impact of immigration into the UK after the Great Recession. J Int Migr Integr 25: 1943–1961. https://doi.org/10.1007/s12134-024-01152-x doi: 10.1007/s12134-024-01152-x
    [29] Golgher AB, Voss PR (2016) How to interpret the coefficients of spatial models: Spillovers, direct and indirect effects. Spat Demogr 4: 175–205. https://doi.org/10.1007/s40980-015-0016-y doi: 10.1007/s40980-015-0016-y
    [30] González L, Ortega F (2011) How do very open economies adjust to large immigration flows? Evidence from Spanish regions. Labour Econ 18: 57–70. https://doi.org/10.1016/j.labeco.2010.06.001 doi: 10.1016/j.labeco.2010.06.001
    [31] Gutiérrez-Portilla M, Villaverde J, Maza A, et al. (2018) A spatial approach to the impact of immigration on wages: Evidence from Spain. Reg Stud 54: 505–514. https://doi.org/10.1080/00343404.2018.1424994 doi: 10.1080/00343404.2018.1424994
    [32] Kyriazi A, Mendes MS, Rone J, et al. (2023) The politics of emigration in Europe: A research agenda. J Common Mark Stud 61: 563–575. https://doi.org/10.1111/jcms.13392 doi: 10.1111/jcms.13392
    [33] LeSage JP, Pace RK (2009) Introduction to spatial econometrics, Boca Raton: CRC Press. https://doi.org/10.1201/9781420064254
    [34] Lin KH, Weiss I (2019) Immigration and the wage distribution in the United States. Demography 56: 2229–2252. https://doi.org/10.1007/s13524-019-00828-9 doi: 10.1007/s13524-019-00828-9
    [35] Longhi S (2010) Job competition and the wage curve. Reg Stud 46: 611–620. https://doi.org/10.1080/00343404.2010.521145 doi: 10.1080/00343404.2010.521145
    [36] Longhi S, Nijkamp P, Poot J (2005) A meta-analytic assessment of the effect of immigration on wages. J Econ Surv 19: 451–477. https://doi.org/10.1111/j.0950-0804.2005.00255.x doi: 10.1111/j.0950-0804.2005.00255.x
    [37] Longhi S, Nijkamp P, Poot J (2006) Spatial heterogeneity and the wage curve revisited. J Reg Sci 46: 707–731. https://doi.org/10.1111/j.1467-9787.2006.00474.x doi: 10.1111/j.1467-9787.2006.00474.x
    [38] Longhi S, Nijkamp P, Poot J (2010) Joint impacts of immigration on wages and employment: Review and meta-analysis. J Geogr Syst 12: 355–387. https://doi.org/10.1007/s10109-010-0111-y doi: 10.1007/s10109-010-0111-y
    [39] Martín-Román ÁL, Cuéllar-Martín J, Moral A (2020) Labor supply and the business cycle: The "Bandwagon Worker Effect". Pap Reg Sci 99: 1607–1642. https://doi.org/10.1111/pirs.12542 doi: 10.1111/pirs.12542
    [40] Maza A, Gutiérrez-Portilla M, Hierro M, et al. (2019) Internal migration in Spain: Dealing with multilateral resistance and nonlinearities. Int Migr 57: 75–93. https://doi.org/10.1111/imig.12472 doi: 10.1111/imig.12472
    [41] Maza A, Moral-Arce I (2006) An analysis of wage flexibility: Evidence from the Spanish regions. Ann Reg Sci 40: 621–637. https://doi.org/10.1007/s00168-005-0033-7 doi: 10.1007/s00168-005-0033-7
    [42] Maza A, Villaverde J (2009) Provincial wages in Spain: Convergence and flexibility. Urban Stud 46: 1969–1993. https://doi.org/10.1177/0042098009106018 doi: 10.1177/0042098009106018
    [43] Nedoncelle C, Marchal L, Aubry A, et al. (2024) Does immigration affect native wages? A meta-analysis. Available from: https://hdl.handle.net/10419/281775.
    [44] Nickell S, Saleheen J (2015) The impact of immigration on occupational wages: Evidence from Britain. Staff working paper, Bank of England, 2015. http://dx.doi.org/10.2139/ssrn.2706493
    [45] Nieto S, Ramos R (2017) Overeducation, skills and wage penalty: Evidence for Spain using PIAAC data. Soc Indic Res 134: 219–236. https://doi.org/10.1007/s11205-016-1423-1 doi: 10.1007/s11205-016-1423-1
    [46] Norris P, Inglehart R (2019) Cultural Backlash: Trump, Brexit, and Authoritarian Populism, Cambridge: Cambridge University Press. https://doi.org/10.1017/9781108595841.002
    [47] Porras-Arena S, Martín-Román AL (2023) The heterogeneity of Okun's law: A metaregression analysis. Econ Model 128: 106490. https://doi.org/10.1016/j.econmod.2023.106490 doi: 10.1016/j.econmod.2023.106490
    [48] Porras-Arena S, Martín-Román ÁL, Dueñas-Fernández D, et al. (2024) Okun's law: The effects of the Covid-19 pandemic and the temporary layoffs procedures (ERTEs) on Spanish regions. Investig Reg 59: 105–125. https://doi.org/10.38191/iirr-jorr.24.013 doi: 10.38191/iirr-jorr.24.013
    [49] Ramos R, Nicodemo C, Sanromá E (2015) A spatial panel wage curve for Spain. Lett Spat Resour Sci 8: 125–139. https://doi.org/10.1007/s12076-014-0118-y doi: 10.1007/s12076-014-0118-y
    [50] Rincke J (2010) A commuting-based refinement of the contiguity matrix for spatial models, and an application to local police expenditures. Reg Sci Urban Econ 40: 324–330. https://doi.org/10.1016/j.regsciurbeco.2010.04.002 doi: 10.1016/j.regsciurbeco.2010.04.002
    [51] Roos C, Nagel M, Kieschnick H, et al. (2025) The emigration conundrum: EU countries of origin of migrants between integration and demarcation. J Common Mark Stud 63: 160–178. https://doi.org/10.1111/jcms.13637 doi: 10.1111/jcms.13637
    [52] Sanchis-Llopis JA, Cutanda A (2020) The Spanish cyclicality of the user cost of labour. Econ Bull 40: 1893–1899.
  • This article has been cited by:

    1. Shishuo Fu, Jiaxi Lu, Yuanzhe Ding, A skeleton model to enumerate standard puzzle sequences, 2021, 30, 2688-1594, 179, 10.3934/era.2022010
    2. Kanasottu Anil Naik, Rayappa David Amar Raj, Chepuri Venkateswara Rao, Thanikanti Sudhakar Babu, Generalized cryptographic image processing approaches using integer-series transformation for solar power optimization under partial shading, 2022, 272, 01968904, 116376, 10.1016/j.enconman.2022.116376
    3. Sen-Peng Eu, Tung-Shan Fu, Springer Numbers and Arnold Families Revisited, 2024, 10, 2199-6792, 125, 10.1007/s40598-023-00230-9
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(691) PDF downloads(68) Cited by(0)

Figures and Tables

Figures(1)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog