Electronic Journal of Differential Equations, Vol. 2020 (2020), No. 13, pp. 1–25. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ORLICZ ESTIMATES FOR GENERAL PARABOLIC OBSTACLE PROBLEMS WITH p(t, x)-GROWTH IN REIFENBERG DOMAINS HONG TIAN, SHENZHOU ZHENG Abstract. This article shows a global gradient estimate in the framework of Orlicz spaces for general parabolic obstacle problems with p(t, x)-Laplacian in a bounded rough domain. It is assumed that the variable exponent p(t, x) sat- isfies a strong log-Hölder continuity, that the associated nonlinearity is mea- surable in the time variable and have small BMO semi-norms in the space variables, and that the boundary of the domain has Reifenberg flatness. 1. Introduction We devote this article to obtaining a nonlinear Calderón-Zygmund type esti- mate in the framework of Orlicz spaces for general parabolic obstacle problems of nonstandard growths with weaker regularity assumptions imposed on given datum. First, let us review recent studies on the related topic. The Calderón-Zygmund estimate for elliptic p-Laplacian in the scalar setting N = 1 had been first obtained by Iwaniec [17], while the vectorial setting N > 1 was treated by DiBenedetto and Manfredi [14]. An extension to general elliptic equations with VMO lead- ing coefficients was achieved by Kinnunen and Zhou [22]. Recently, a nonlinear Calderón-Zygmund estimate for parabolic obstacle problems involving possibly de- generate operators of p-growth was obtained by Bögelein, Duzaar and Mingione [6]. Byun and Cho [8] also established a local Calderón-Zygmund estimate for par- abolic variational inequalities of general type degenerate and singular operators in divergence form, and they proved that for any q ∈ (1,∞) it holds |ψt|p ′ , |Dψ|p, |F |p ∈ Lqloc(ΩT ) =⇒ |Du|p ∈ Lqloc(ΩT ). A local regularity version in Lorentz spaces for the gradients of weak solution to parabolic obstacle problems has been also achieved by Baroni [3]. Later, Byun and Cho in [9] showed a global regularity in Orlicz spaces for the gradients of weak solution to parabolic variational inequalities of p-Laplacian type under weak assumptions that the nonlinearities are merely measurable in the time-variable and have small BMO semi-norms in the spatial variables, while the underlying domain is a Reifenberg flatness. Tian and Zheng [27] also derived a global weighted Lorentz estimate to nonlinear parabolic equations with partial regular nonlinearity in a nonsmooth domain. On the other hand, we would like to mention that Zhang and 2010 Mathematics Subject Classification. 35B65, 35K86, 46E30. Key words and phrases. Parabolic obstacle problems; discontinuous nonlinearities; p(t, x)-growth; Orlicz spaces; Reifenberg flat domains. c©2020 Texas State University. Submitted August 21, 2019. Published January 27, 2020. 1 2 H. TIAN, S. ZHENG EJDE-2020/13 Zheng [28] got Lorentz estimates for asymptotically regular elliptic equations in quasiconvex domains, and Liang and Zheng [19] established the gradient estimate in Orlicz spaces for elliptic obstacle problems with partially BMO nonlinearities. Very recently, Liang, Zheng and Feng [20] showed a global Calderón-Zygmund type estimate in Lorentz spaces for a variable power of the gradients of weak solution pair (u, P ) to generalized steady Stokes system in a bounded Reifenberg domain. Nonlinear elliptic and parabolic problem under consideration with a variable growth naturally originates from some mathematical modeling of fluid dynamics, such as certain models for non-Newtonian fluids and electrorheological fluids. In- deed, there are also various phenomena involved some energy functionals, for exam- ple, elastic mechanics, porous media problems, and thermistor problems (cf. [26]). Therefore, it is a rather interesting topic in the fields of analysis and PDEs to get nonlinear Calderón-Zygmund theory for general elliptic and parabolic equations with variable growths. In recent decades, a lot of attention has been paid to a sys- tematic study on the Calderón-Zygmund theory for nonlinear elliptic and parabolic problems with nonstandard growths. For instance, some regularities regarding gen- eral elliptic equations of p(x)-growth have been treated by Acerbi and Mingione [1]. Naturally, there also have been many interesting theoretic developments involving more general obstacle problems since this kind of problems of variable growths al- ways appeared in various phenomena of physical applications. It was observed by Bögelein and Duzaar in [5] that it holds a higher integrability for the gradients of weak solutions to possibly degenerate parabolic systems with nonstandard growth. Later Baroni and Bögelein in [4] showed nonlinear Calderón-Zygmund estimate for evolutionary p(t, x)-Laplacian system in requiring the variable exponent p(t, x) being a logarithmic Hölder continuity and the coefficients a(t, x) satisfying VMO condition in the spatial variables. Erhardt [16] considered an interior Lq-estimate of |Du|p(t,x) for general parabolic variational inequality in the weak form as 〈φt, φ− u〉ΩT + ∫ ΩT a(t, x)|Du|p(t,x)−2Du ·D(φ− u) dx dt + 1 2 ‖φ(a, ·)− ua‖2L2(Ω) ≥ ∫ ΩT |f|p(t,x)−2f ·D(φ− u) dx dt, (1.1) and he showed that |Du|p(t,x) belongs to a local integrability with the same index as an assigned obstacle |Dψ|p(t,x), |ψt|γ ′ 1 as well as |f|p(t,x), which implies that |ψt|γ ′ 1 , |Dψ|p(t,x), |f|p(t,x) ∈ Lqloc(ΩT ) =⇒ |Du|p(t,x) ∈ Lqloc(ΩT ) for any q ∈ (1,∞). On the other hand, Li [24] handled a higher integrability for the derivatives of very weak solutions to parabolic systems of p(t, x)-Laplacian type with the inhomogeneity being different growths, respectively. Furthermore, Bui and Duong [7] derived global weighted estimate in Lorentz spaces for nonlinear parabolic equations of p(t, x)-growth in a Reifenberg flat domain with the non- linearities a(t, x; ξ) being small BMO in the spatial variables, while the variable growth p(t, x) satisfying a strong log-Hölder continuity. Byun and Ok [10] reached a global Ls(t,x)-integrability with s(t, x) > p(t, x) for the gradients of weak solution to general parabolic equations of p(t, x)-growth in Reifenberg flat domains by im- posing the same weak regular assumptions as shown in [7] on a(t, x; ξ), p(t, x) and EJDE-2020/13 ORLICZ ESTIMATES 3 the boundary of the underlying domains. Li, Zhang and Zheng [18] established a local Orlicz estimate for nondivergence linear elliptic equations with partially BMO coefficients, and Chlebicka in [12] provided the Lorentz and Morrey estimates for the gradients of solution to general nonlinear elliptic equations with the datum of Orlicz growths. Byun and Park [11] considered global weighted Orlicz estimate to nonlinear parabolic equation with measurable nonlinearity in a bounded nonsmooth domain while the right-hand side is of finite signed Radon measure. This article is inspired by these above-mentioned recent progresses. The aim of this article is to show a global Calderón-Zygmund type estimate in Orlicz spaces for nonlinear parabolic obstacle problems of nonstandard growth with weaker regularity assumptions on the given datum, which means an implication that |ψt|γ ′ 1 , |Dψ|p(t,x), |f|p(t,x) ∈ Lφ(ΩT ) =⇒ |Du|p(t,x) ∈ Lφ(ΩT ) (1.2) for Young’s function φ ∈ ∆2 ∩ ∇2 defined below. As we know, the Orlicz space is a generalization of Lebesgue spaces. Jia, Li and Wang [21] recently obtained a global Orlicz estimate to linear elliptic equations of divergence form with small BMO coefficients in Reifenberg flat domains. Byun and Cho [9] obtained Orlicz estimates to parabolic obstacles problems of p-Laplacian type for 2d d+2 < p < ∞, they proved that |ψt|p ′ , |Dψ|p, |F |p ∈ Lφ(ΩT ) =⇒ |Du|p ∈ Lφ(ΩT ) for φ ∈ ∆2 ∩ ∇2 while the nonlinearity is small BMO in spatial variables and the domain is Reifenberg flatness. A key ingredient under consideration is the power p(t, x) being a variable function with respect to the independent variables (t, x). In this way, the Hardy-Littlewood maximal operators technique does not work well for parabolic equations of p(t, x)- growth since the usual scaling arguments used for p = 2 do not work smoothly. The main difficulty for parabolic setting comes from the nonhomogeneous scaling behavior for variational inequalities so that any solution multiplied by a constant is in general no longer a solution of original problem. We here employ the technique of the so-called intrinsic parabolic cylinder first introduced by DiBenedetto and Friedman [13], which applies the time-space scaling dependent on a local behavior of the solution itself to re-balance the nonhomogeneous scaling for parabolic problem of p-Laplacian. Another point is that we adapt the so-called large-M-inequality principle from Acerbi and Mingione’s work [2] to our situation with non-trivial modifications and significant improvements. In order to get a suitable power decay for the following upper level we set{ (t, x) ∈ ΩT : |Du|p(t,x) > κ } with the scaling parameter κ > 0 sufficiently large, we make use of the so-called stop-time argument and the modified Vitali type covering with a countable covering by the intrinsic parabolic cylinder {Qκri(τi, yi)} ∞ i=1 satisfying − ∫ Qκri (τi,yi) |Du|p(t,x) dx dt ≈ κ, which will be discussed in Section 3. The rest of this article is organized as follows. In the next section we present the weaker regular assumptions on the datum, and state our main result. Section 3 is to give necessary preliminary lemmas, in which shows various comparison estimates to 4 H. TIAN, S. ZHENG EJDE-2020/13 the reference problems. Finally, we devote Section 4 to the proof of main Theorem 2.5. 2. Minimal assumptions on the datum and main result Let Ω be a bounded domain in Rd for d ≥ 2 with its rough boundary ∂Ω specified later. For a fixed a ∈ R and 0 < T < ∞, let ΩT = (a, a + T ) × Ω denote the parabolic cylinder in R×Rd, and the typical parabolic boundary ∂ΩT =( (a, a + T ) × ∂Ω ) ∪ ( {t = a} × Ω ) be the typical parabolic boundary of ΩT . We suppose that the main nonlinearity a(t, x; ξ) = ( a1(t, x; ξ), a2(t, x; ξ), . . . , ad(t, x; ξ) ) : ΩT × Rd −→ Rd is a Carathéodory vectorial-valued function with the following basic structural con- ditions: for a.a. (t, x) ∈ ΩT and all ξ, η ∈ Rd, there exist constants 0 < λ ≤ 1 ≤ Λ and 0 ≤ µ ≤ 1 such that λ ( µ2 + |ξ|2 ) p(t,x)−2 2 |η|2 ≤ Dξa(t, x; ξ)η · η, |a(t, x; ξ)|+ ( µ2 + |ξ|2 )1/2 |Dξa(t, x; ξ)| ≤ Λ ( µ2 + |ξ|2 ) p(t,x)−1 2 . (2.1) Let the given obstacle function ψ : ΩT → R satisfy ψ ∈ C0([a, a+ T ];L2(Ω)) ∩W p(t,x)(ΩT ), ψt ∈ Lγ ′ 1(a, a+ T ;W−1,γ′1(Ω)), ψ ≤ 0 a.e. on (a, a+ T )× ∂Ω, ψ(a, ·) ≤ 0 a.e. on Ω; (2.2) and let an initial value ua be such that ua = u(a, ·) ∈ L2(Ω) and ua ≥ ψ(a, ·) a.e. on Ω. We introduce an admissible set defined by A(ΩT ) = { φ ∈ C0([a, a+ T ];L2(Ω)) ∩W p(t,x) 0 (ΩT ) : φ ≥ ψ a.e. on ΩT } . (2.3) Note that a minimizing the energy functional with certain constraint in A(ΩT ) immediately leads to the following form: for u = u(t, x) ∈ A(ΩT ) it holds in the weak form of the parabolic variational inequality 〈φt, φ− u〉ΩT + ∫ ΩT a(t, x,Du) ·D(φ− u) dx dt+ 1 2 ‖φ(a, ·)− ua‖2L2(Ω) ≥ ∫ ΩT |f|p(t,x)−2f ·D(φ− u) dx dt (2.4) for all test functions φ ∈ A′(ΩT ) = { φ ∈ A(ΩT ) : φt ∈ (W p(t,x)(ΩT ))′ } , where f ∈ Lp(t,x)(ΩT ) is a given inhomogeneous term. For convenience, throughout this paper we assume that R ≤ 1 is an arbitrary given positive number, while δ ∈ (0, 1/8) is to be determined later. Let us now endow the variable exponent p(t, x) : ΩT → R with the regularity of the so- called strong log-Hölder continuity. We write dp(z1, z2) the parabolic distance by dp(z1, z2) := max { |x− y|, √ |τ − t| } for any z1 = (t, x), z2 = (τ, y) ∈ Rd+1. We say that p(t, x) is locally strong log-Hölder continuous, denoted by p(t, x) ∈ SLH(ΩT ), if for some constant ρ̄ > 0 such that for all z1 = (t, x), z2 = (τ, y) ∈ ΩT with EJDE-2020/13 ORLICZ ESTIMATES 5 0 < dp(z1, z2) < ρ̄, one has that there exists a nondecreasing continuous function ω(·) : [0,∞) −→ [0, 1] satisfying ω(0) = 0 such that |p(z1)− p(z2)| ≤ ω(dp(z1, z2)) (2.5) with lim sup ρ→0 ω(ρ) log (1 ρ ) = 0. It is easy to check that if p(t, x) is a strong logarithmic Hölder continuity, then for any given δ ∈ (0, 1/8) there exists a small R > 0 such that sup 0<ρ 2d d+2 is unavoidable even in the constant exponent setting p(z) ≡ p, for more details see [10, Section 2]. With the assumptions (2.1) (2.2) (2.6) and (2.7) in hand, the existence of such weak solution is ensured by the result from Erhardt [16], which leads to that there exists a unique weak solution u ∈ A(ΩT ) to the parabolic variational inequality (2.4) with the estimate sup t∈[a,a+T ] ∫ Ω |u(t, x)|2 dx+ ∫ ΩT |Du|p(t,x) dx dt ≤ C (∫ ΩT ( |ψt|γ ′ 1 + |Dψ|p(t,x) + |f|p(t,x) + 1 ) dx dt ) , (2.8) where C is a positive constant depending only on d, γ1, γ2, λ,Λ and ‖ua‖L2(Ω), see also [16, Theorem 7.1]. We now recall that the space Lp(t,x)(ΩT ) is defined to be the set of these mea- surable functions g(t, x) : ΩT → Rk for k ∈ N, which satisfies |g|p(t,x) ∈ L1(ΩT ), i.e. Lp(t,x)(ΩT ) := { g(t, x) : ΩT → Rk is measurable in ΩT : ∫ ΩT |g|p(t,x) dx dt < +∞ } , which is a Banach space equipped with the Luxemburg norm ‖g‖Lp(t,x)(ΩT ) := inf { λ > 0 : ∫ ΩT ∣∣ g λ ∣∣p(t,x) dx dt ≤ 1 } . (2.9) The Sobolev spaces W p(t,x)(ΩT ) is defined by W p(t,x)(ΩT ) := { g ∈ Lp(t,x)(ΩT ) : Dg ∈ Lp(t,x)(ΩT ) } endowed with the norm ‖g‖Wp(t,x)(ΩT ) := ‖g‖Lp(t,x)(ΩT ) + ‖Dg‖Lp(t,x)(ΩT ). (2.10) It would be worthwhile to mention that for g ∈ W p(t,x) 0 (ΩT ) it indicates that g(t, x) = 0 in the sense of trace on the boundary of Ω. For 1 < p(t, x) < ∞, we 6 H. TIAN, S. ZHENG EJDE-2020/13 also denote the dual space of W p(t,x) 0 (ΩT ) by (W p(t,x)(ΩT ))′, which means that for g ∈ (W p(t,x)(ΩT ))′ there exist functions gi ∈ Lp ′(t,x)(ΩT ) with p′(t, x) = p(t,x) p(t,x)−1 for i = 0, 1, . . . , d such that the dual parting 〈g, w〉ΩT = ∫ ΩT ( g0w + d∑ i=1 giDiw ) dx dt for all w ∈W p(t,x) 0 (ΩT ). In particular, if p(t, x) = γ1 it yields that W γ1(ΩT ) = Lγ1 ( a, a+ T ;W 1,γ1(Ω) ) . Consequently, the dual space of W γ1 0 (ΩT ) is given by( W γ1 0 (ΩT ) )′ = ( Lγ1(a, a+ T ;W 1,γ1 0 (Ω)) )′ = Lγ ′ 1(a, a+ T ;W−1,γ′1(Ω)), where 1 γ1 + 1 γ′1 = 1. Now we impose some regularity assumptions on the nonlinearities a(t, x; ξ) and on the boundary ∂Ω of domain. For this, let ρ, θ > 0, Bρ(y) = {x ∈ Rd : |x−y| < ρ}, and the local parabolic cylinders Q(θ,ρ)(τ, y) = (τ − θ, τ + θ)×Bρ(y) with any (τ, y) ∈ R×Rd. For the abbreviations, Bρ = Bρ(0) , Q(θ,ρ) = Q(θ,ρ)(0, 0) and Qρ = Q(ρ2,ρ), we measure the oscillation of a(t, x; ξ)/(µ2 + |ξ|2) p(t,x)−1 2 in the x-variables over the ball Bρ(y) by Θ[a;Bρ(y)](t, x) := sup ξ∈Rd ∣∣∣ a(t, x; ξ) (µ2 + |ξ|2) p(t,x)−1 2 − ( a(t, ·, ξ) (µ2 + |ξ|2) p(t,·)−1 2 ) Bρ(y) ∣∣∣, where ( a(t, ·, ξ) (µ2 + |ξ|2) p(t,·)−1 2 ) Bρ(y) := 1 |Bρ(y)| ∫ Bρ(y) a(t, x; ξ) (µ2 + |ξ|2) p(t,x)−1 2 dx represents an integral average of a(t, x; ξ)/(µ2 + |ξ|2) p(t,x)−1 2 in the x-variables over Bρ(y) for any fixed ξ ∈ Rd and t ∈ R. Assumption 2.1. Let δ ∈ (0, 1/8) to be specified later. We say that (a,ΩT ) is a (δ,R)-vanishing in the spatial variables, if for every point (τ, y) ∈ ΩT there exists a constant 0 < R ≤ 1 such that for any ρ ∈ (0, R) the following relation holds: (i) If dist(y, ∂Ω) = min x∈∂Ω dist(y, x) > √ 2ρ, then there exists a coordinate system depending on (τ, y) and ρ, whose variables are still denoted by (t, x) such that in this new coordinate system (τ, y) is the origin, and for every θ ∈ (0, ρ2) one has − ∫ Q(θ,ρ)(τ,y) ∣∣Θ[a;Bρ(y)](t, x) ∣∣2 dx dt ≤ δ2; (ii) while dist(y, ∂Ω) = min x∈∂Ω dist(y, x) = dist(y, ȳ) ≤ √ 2ρ EJDE-2020/13 ORLICZ ESTIMATES 7 for some ȳ ∈ ∂Ω, there exists a new coordinate system depending on (τ, y) and ρ, whose variables are denoted by (t, x), such that in this new coordinate system (τ, ȳ) is the origin, and for any θ ∈ ( 0, (3ρ)2 ) it holds B3ρ(ȳ) ∩ {x1 > 3δρ} ⊂ B3ρ(ȳ) ∩ Ω ⊂ B3ρ(ȳ) ∩ {x1 > −3δρ} (2.11) and − ∫ Q(θ,3ρ)(τ,ȳ) ∣∣Θ[a;B3ρ(ȳ)](t, x) ∣∣2 dx dt ≤ δ2. Remark 2.2. Roughly speaking, the nonlinearity a(t, x; ξ)/(µ2 + |ξ|2) p(t,x)−1 2 is assumed to be a small BMO semi-norm in the x-variables, but there is no regular requirement in the t-variable, uniformly in ξ ∈ Rd; while the domain Ω is assumed to be the (δ,R)-Reifenberg flatness as a necessary geometric condition if (2.11) holds, which leads to the following measure density conditions: sup 0 0 with φ(2ν) ≤ K̄φ(ν). (2.14) On the other hand, the Young function φ is said to be global ∇2-condition, denoted by φ ∈ ∇2, if there exists a constant ā > 1 such that for every ν > 0 one has φ(ν) ≤ φ(āν) 2ā . (2.15) Remark 2.3. Actually, φ ∈ ∆2 implies that for any β1 > 1 there exists α1 = log2K̄ such that φ(β1ν) ≤ K̄βα1 1 φ(ν), which describes the growth for φ(ν) near ν = ∞. Meanwhile, the condition φ ∈ ∇2 means that for any 0 < β2 < 1, there exists α2 = logā2 + 1 such that φ(β2ν) ≤ 2āβα2 2 φ(ν), and it describes the growth for φ(ν) near ν = 0. The simplest example for φ(ν) satisfying the ∆2 ∩∇2 condition is the power function φ(ν) = νp with p > 1. Moreover, we also remark that for p > 1, φ(ν) = |ν|p(1 + |log|ν||) ∈ ∆2 ∩∇2. Definition 2.4. Let D be an open subset in Rd+1 and φ be a Young function. The Orlicz class Kφ(D) is called to be the set of all measurable functions g : D → R satisfying ∫ D φ (|g|) dx dt <∞. 8 H. TIAN, S. ZHENG EJDE-2020/13 Orlicz space Lφ(D) is just a linear hull of Kφ(D). It consists of all measurable functions f such that η̂f ∈ Kφ(D) for some η̂ > 0. Moreover, the norm ‖ · ‖Lφ(D) is denoted by ‖g‖Lφ(D) = inf { λ > 0 : ∫ D φ ( |g| λ ) dx dt ≤ 1 } . If D is bounded, then Lα1(D) ⊂ Lφ(D) ⊂ Lα2(D) ⊂ L1(D) with the constants α1 and α2 as Remark 2.3, for more details see [25]. We are now in a position to state the main result of this paper. Theorem 2.5. Let the Young function φ ∈ ∆2 ∩∇2, and p(t, x) ∈ SLH(ΩT ) with its range in [γ1, γ2] shown as (2.7). Assume that u ∈ A(ΩT ) is a weak solution of the variational inequality (2.4) with the given datum |ψt|γ ′ 1 , |Dψ|p(t,x), |f|p(t,x) ∈ Lφ(ΩT ). Then, there exists a small positive constant δ = δ(d, λ,Λ, γ1, γ2, ∂Ω) such that if (a,ΩT ) satisfies (δ,R)-vanishing as Assumption 2.1, then we have |Du|p(t,x) ∈ Lφ(ΩT ) with the estimate∫ ΩT φ ( |Du|p(t,x) ) dx dt ≤ C [ φ (( − ∫ ΩT Ψ(t, x) dx dt )m) + ∫ ΩT φ(Ψ(t, x)) dx dt ] , where C = C(d, γ1, γ2, λ,Λ, α1, α2, δ, R, T, |Ω|, ‖ua‖L2(Ω)), Ψ(t, x) = |ψt|γ ′ 1 + |Dψ|p(t,x) + |f|p(t,x) + 1, and m ≥ 1 with m = sup (τ,y)∈ΩT m(τ, y), (2.16) m(τ, y) = { p(τ,y) 2 if p(τ, y) ≥ 2, 2p(τ,y) p(τ,y)(d+2)−2d if 2d d+2 < p(τ, y) < 2. (2.17) 3. Comparison estimates to the reference problems We start this section with introducing some related notations and basic facts which will be useful in the paper. Throughout the paper, we always use Ci and ci for i = 1, 2, . . . , to denote positive constants that only depend on d, λ,Λ, γ1, γ2, . . . , but whose values may differ from line to line. For any fixed point z = (τ, y) ∈ Rd+1 with τ ∈ R and y ∈ Rd, we denote the spatial open ball Bρ(y) ⊂ Rd with center y and the radius ρ > 0. For any κ > 1 we write the intrinsic parabolic cylinder by Qκρ(z) = Qκρ(τ, y) = ( τ − κ 2−p(z) p(z) ρ2, τ + κ 2−p(z) p(z) ρ2 ) Bρ(y). (3.1) We also set Ωρ = Ω ∩Qρ, Kκ ρ (z) = Qκρ(z) ∩ ΩT , ∂Qκρ(z) = ( τ − κ 2−p(z) p(z) ρ2, τ + κ 2−p(z) p(z) ρ2 ) ∂Bρ(y), ∂Kκ ρ (z) = ( Qκρ(z) ∩ ( (a, a+ T )× ∂Ω )) ∪ ( ∂Qκρ(z) ∩ ΩT ) EJDE-2020/13 ORLICZ ESTIMATES 9 for a ∈ R and T > 0. For the sake of convenience, while z = (τ, y) = (0, 0) we simply write Qκρ = Qκρ(0), Kκ ρ = Kκ ρ (0) and ∂Kκ ρ = ∂Kκ ρ (0). In the following we write B+ ρ (y) = Bρ(y) ∩ {x1 > 0} , Qκ+ ρ (z) = ( τ − κ 2−p(z) p(z) ρ2, τ + κ 2−p(z) p(z) ρ2 ) B+ ρ (y), Tκρ (z) = ( τ − κ 2−p(z) p(z) ρ2, τ + κ 2−p(z) p(z) ρ2 )( Bρ(y) ∩ {x1 = 0} ) . Also, we briefly denote B+ ρ = B+ ρ (0), Qκ+ ρ = Qκ+ ρ (0) and Tκρ = Tκρ (0). We use the following localizing technique, which is first used by Bögelein and Duzaar in [5]. As we know, an interior estimate for parabolic obstacle problems with nonstandard growth had been obtained by Erhardt in [16]. Owing to the measure density (2.12), this readily allows an obvious extension to the Reifenberg flat domain. More precisely, we state the following boundary estimate by setting Kκ ρ (z) = Qκρ(z) ∩ ΩT for a fixed z = (τ, y) ∈ (a, a+ T )× ∂Ω. Lemma 3.1. Suppose that p(t, x) ∈ SLH(ΩT ) with its range in [γ1, γ2] shown as (2.7), and M := ∫ ΩT |Du|p(t,x) dx dt+ ∫ ΩT Ψ(t, x) dx dt (3.2) with Ψ(t, x) := |ψt|γ ′ 1 + |Dψ|p(t,x) + |f|p(t,x) + 1. (3.3) For any fixed δ ∈ (0, 1/8), M > 1 and α := min { 1, γ1 d+2 4 − d 2 } ∈ (0, 1], let ρ1 = Γ− 2 α with Γ := 4 (( 2 1− δ )dMM 2δωd + 1 )1/2 ≥ 4, (3.4) where ωd denotes the measure of the unit ball of Rd. If Ω is a (δ,R)-Reifenberg flat domain; moreover, for any fixed κ > 1, z = (τ, y) ∈ (a, a + T ) × ∂Ω and for any 0 < ρ < ρ1 we have κ ≤M ( − ∫ Kκ ρ (z) |Du|p(t,x) dx dt+ 1 δ − ∫ Kκ ρ (z) Ψ(t, x) dx dt ) , (3.5) then there exists ca := exp ( γ2 ( δ + δ(d+2) α )) > 1 such that p2 − p1 ≤ ω(Γ ρα), κ 2 p(z) ≤ Γ2 ρ−(d+2), κp2−p1 ≤ ca, (3.6) where p1 = p(z1) = inf Kκ ρ (z) p(t, x), p2 = p(z2) = sup Kκ ρ (z) p(t, x). (3.7) Proof. For a fixed point z = (τ, y) ∈ (a, a+T )×∂Ω, it suffices to prove our estimate in the setting ( τ−κ 2−p(z) p(z) ρ2, τ+κ 2−p(z) p(z) ρ2 ) ⊂ (a, a+T ). Otherwise, if Qκρ(z) touches the bottom or the top of ΩT , i.e. ( τ − κ 2−p(z) p(z) ρ2, τ + κ 2−p(z) p(z) ρ2 ) 6⊂ (a, a + T ), then we may consider an extended variational inequality (2.4) in (a− T, a+ 2T )×Ω in terms of an argument from [10, Remark 2.6], which results in that( τ − κ 2−p(z) p(z) ρ2, τ + κ 2−p(z) p(z) ρ2 ) ⊂ (a− T, a+ 2T ). Consequently, it yields the same process as follows. 10 H. TIAN, S. ZHENG EJDE-2020/13 Now, by the measure density (2.12) we know that 1 |Kκ ρ (z)| = 1 |Qκρ(z)| |Qκρ(z)| |Kκ ρ (z)| = 1 2ωdρd+2κ 2−p(z) p(z) |Bρ(y)| |Bρ(y) ∩ Ω| ≤ 1 2ωd ρd+2κ 2−p(z) p(z) ( 2 1− δ )d . Hence, from (3.5) it follows that κ ≤ M |Kκ ρ (z)| (∫ Kκ ρ (z) |Du|p(t,x) dx dt+ 1 δ ∫ Kκ ρ (z) Ψ(t, x) dx dt ) ≤ M 2ωdρd+2κ 2−p(z) p(z) ( 2 1− δ )d(∫ ΩT |Du|p(t,x) dx dt+ 1 δ ∫ ΩT Ψ(t, x) dx dt ) , which implies that κ 2 p(z) ≤ M 2ωdρd+2 ( 2 1− δ )d(∫ ΩT |Du|p(t,x) dx dt+ 1 δ ∫ ΩT Ψ(t, x) dx dt ) ≤ MM 2ωdρd+2 ( 2 1− δ )d , (3.8) where we have used (3.2) in the last inequality. Recalling the definitions of p1 and p2, by (2.5) it yields p2 − p1 ≤ |p2 − p1| ≤ ω(dp(z1, z2)) ≤ ω ( 2ρ+ √ 2κ 2−p(z) p(z) ρ ) . So, if 2 ≤ p(z) ≤ γ2 <∞, then p2 − p1 ≤ ω(4ρ); (3.9) if 2d d+2 < γ1 ≤ p(z) < 2 then by (3.8) and (3.4) we obtain p2 − p1 ≤ ω(Γργ1 d+2 4 − d 2 ). (3.10) Combining (3.9) and (3.10), we obtain the first estimate of (3.6). Further, putting (3.8) and (3.4) together, we also obtain the second estimate of (3.6). Finally, recalling p(t, x) ∈ SLH(ΩT ) and 0 < ρ < ρ1 = Γ− 2 α , we have Γp2−p1 ≤ exp(δ), ρ−(p2−p1) ≤ exp (2δ α ) , which implies κp2−p1 ≤ ( Γ ρ− d+2 2 )(p2−p1)γ2 ≤ exp ( γ2 ( δ + δ(d+ 2) α )) = ca. (3.11) This concludes the proof. � Let us recall the modified Vitali type covering lemma with a covering of intrinsic parabolic cylinders, see [10, Lemma 3.5]. EJDE-2020/13 ORLICZ ESTIMATES 11 Lemma 3.2. For κ > 1, we set that F = {Qκρi(zi)}i∈J is a family of intrinsic par- abolic cylinders with zi = (τi, yi) ∈ Rd+1 and ρi > 0, which satisfy that ∪i∈JQκρi(zi) is bounded in Rd+1 and κp + i −p − i ≤ ca for all i ∈ J , where ca > 1 is the same as Lemma 3.1. Let p+ i = sup Qκρi (zi) p(t, x) and p−i = inf Qκρi (zi) p(t, x), then there exists a countable sub-collection G ⊂ F of disjoint parabolic cylinders such that ∪Qκρi (zi)∈FQ κ ρi(zi) ⊂ ∪Qκρi (zi)∈GχQ κ ρi(zi), where χ ≥ { 5, ( 8c 4 γ2 1 a + 1 )1/2} , and χQi denotes the χ-time enlarged cylinder Qi. To obtain the interior and boundary comparison estimates with the reference problems on the intrinsic parabolic cylinders, respectively, we suppose that u ∈ A(ΩT ) is a weak solution of (2.4) under the regularity assumptions that p(t, x) ∈ SLH(ΩT ) with its range (2.7), and (a,R×Ω) is (δ,R)-vanishing with the specified δ ∈ (0, 1/8) and R ∈ (0, 1). It is clearly checked that the condition (2.1) easily leads to the following monotonicity( a(t, x; ξ)− a(t, x, η) ) (ξ − η) ≥ C1 ( |µ|2 + |ξ|2 + |η|2 ) p(t,x)−2 2 |ξ − η|2 if 2d d+ 2 < p(t, x) < 2,( a(t, x; ξ)− a(t, x, η) ) (ξ − η) ≥ C2|ξ − η|p(t,x) if p(t, x) ≥ 2 (3.12) for all ξ, η ∈ Rd and a.a. (t, x) ∈ ΩT , where C1 and C2 are positive constants depending only on d, γ1, γ2, λ and Λ, see [16, Section 2] or [7, Formula (10)]. Setting W (ΩT ) := { g ∈W p(t,x)(ΩT ) : gt ∈ ( W p(t,x)(ΩT ) )′} . We recall the following comparison principle, which is useful to construct a com- parison that it almost everywhere satisfies an obstacle constrain ψ ≤ k, see [16, Lemma 3.15]. Lemma 3.3. Let ΩT be an open subset of Rd+1. Assume that p(t, x) ∈ SLH(ΩT ) satisfying (2.7), and ψ, k ∈ W (ΩT ) satisfy the following relations with a(t, x; ξ) such that (3.12) holds, ψt − div(a(t, x,Dψ)) ≤ kt − div(a(t, x,Dk)) in ΩT , ψ ≤ k on ∂ΩT . (3.13) Then ψ ≤ k a.e. on ΩT . We set a fixed point z = (τ, y) ∈ ΩT , κ > 1 and a sufficiently small r > 0 specified later. Without loss of generality, we assume that y = 0, i.e., z = (τ, 0). We only consider the boundary case of B+ 6r ⊂ Ω6r := B6r ∩Ω ⊂ {x1 > −12rδ} and( τ −κ 2−pz pz (6r)2, τ +κ 2−pz pz (6r)2 ) ⊂ (a, a+T ) with pz = p(z) since the interior case 12 H. TIAN, S. ZHENG EJDE-2020/13 is simpler for Qκ6r(z) = Kκ 6r(z) ⊂ ΩT . By an argument of normalization we can assume that for suitable r > 0 such that − ∫ Kκ 6r(z) |Du|p(t,x) dx dt+ 1 δ − ∫ Kκ 6r(z) Ψ(t, x) dx dt ≤ c∗κ (3.14) for some c∗ > 1, where Ψ(t, x) is as (3.3). Let k ∈W (Kκ 6r(z)) be any weak solution of the following local initial-boundary problem kt − div(a(t, x,Dk)) = ψt − div(a(t, x,Dψ)) in Kκ 6r(z), k = u on ∂Kκ 6r(z). (3.15) Then, by Lemma 3.3 we immediately conclude the following, cf. [16, Lemma 8.2]. Lemma 3.4. Under the normalization assumption of (3.14), for any ε1 ∈ (0, 1) there exists a small constant δ = δ(d, λ,Λ, γ1, γ2, ε1) > 0 such that − ∫ Kκ 4r(z) |Du−Dk|p(t,x) dx dt ≤ ε1κ and − ∫ Kκ 4r(z) |Dk|p(t,x) dx dt ≤ c1κ (3.16) for some c1 = c1(d, λ,Λ, γ1, γ2, ∂Ω) > 1. Let w ∈W (Kκ 4r(z)) be the weak solution of wt − div(a(t, x,Dw)) = 0 in Kκ 4r(z), w = k on ∂Kκ 4r(z). (3.17) Lemma 3.5. Under the normalization assumption of (3.14), for any ε2 ∈ (0, 1) there exists a small δ = δ(d, λ,Λ, γ1, γ2, ε2) > 0 such that − ∫ Kκ 4r(z) |Dk −Dw|p(t,x) dx dt ≤ ε2κ and − ∫ Kκ 4r(z) |Dw|p(t,x) dx dt ≤ c2κ (3.18) for some c2 = c2(d, λ,Λ, γ1, γ2, ∂Ω) > 1, see [10, Lemma 4.1]. Now let us recall a self-improving integrability of Dw to (3.18). For 0 < ρ = 6r < ρ1, p1 and p2 shown as in (3.7), we assume that p2 − p1 ≤ ω(Γ (6r)α), κ 2 p(z) ≤ Γ2(6r)−(d+2), κp2−p1 ≤ ca (3.19) for some α ∈ (0, 1), Γ ≥ 4 and ca > 1 defined by Lemma 3.1. By Lemma 3.5 it holds − ∫ Kκ 4r(z) |Dw|p(t,x) dx dt ≤ c2κ with c2 > 1. Then, thanks to [4, Corollary 5.2] we conclude that there exist ε0 > 0 and ρ2 > 0 such that for 0 < 4r < ρ2 it holds − ∫ Kκ 2r(z) |Dw|p(t,x)(1+ε0) dx dt ≤ cκ1+ε0 , (3.20) where c is a positive constant depending only on d, λ,Λ, µ, γ1, γ2, δ, R, ω(·). As in [10, 7], let p2 − p1 ≤ min { λ 4Λ , 1, ε0(γ1 − 1) 4 } (3.21) and the vector-valued function b(t, x; ξ) : Kκ 2r(z)× Rd → Rd is introduced by b(t, x; ξ) = a(t, x; ξ) ( µ2 + |ξ|2 ) pz−p(t,x) 2 . EJDE-2020/13 ORLICZ ESTIMATES 13 By using (2.1) and (3.21), we obtain (λ/2) ( µ2 + |ξ|2 ) pz−2 2 |η|2 ≤ Dξb(t, x; ξ) η · η, |b(t, x; ξ)|+ ( µ2 + |ξ|2 )1/2 |Dξb(t, x; ξ)| ≤ 3Λ ( µ2 + |ξ|2 ) pz−1 2 (3.22) for a.a. (t, x) ∈ Kκ 2r(z) and all ξ, η ∈ Rd, see [10, Eq. (4.18)] or [7, Lemma 3.6]. For the interior case, we define b(t, ξ) : ( τ −κ 2−pz pz (2r)2, τ +κ 2−pz pz (2r)2 ) ×Rd → Rd by b(t, ξ) = − ∫ B2r(y) b(t, x; ξ) dx. Then, by 2.1 it yields − ∫ Qκ2r(z) sup ξ∈R |b(t, ξ)− b(t, x, ξ)| (µ2 + |ξ|2) pz−1 2 dx dt = − ∫ Qκ2r(z) Θ[a;B2r(y)](t, x) dx dt ≤ δ. For the boundary case, we define b̃(t, ξ) : ( τ−κ 2−pz pz (2r)2, τ+κ 2−pz pz (2r)2 ) ×Rd → Rd by b̃(t, ξ) := { b(t, ξ) = − ∫ B+ 2r(y) b(t, x; ξ) dx (t, x) ∈ B+ 2r(y), b(t, ξ) (t, x) ∈ Ω2r(y)\B+ 2r(y). Again by Assumption 2.1 we see that − ∫ Qκ+ 2r (z) sup ξ∈R |b(t, ξ)− b(t, x, ξ)| (µ2 + |ξ|2) pz−1 2 dx dt = − ∫ Qκ+ 2r (z) Θ[a;B+ 2r(y)](t, x) dx dt ≤ 4δ. Moreover, for both cases we see that b(t, ξ) satisfies (3.22) with b(t, x; ξ) replaced by b(t, ξ). With b̃(t, ξ) in hand, we further recall the following two comparisons with the so-called limiting problems. Let h ∈W 1,pz (Kκ 2r(z)) be a weak solution of ht − div(b̃(t,Dh)) = 0 in Kκ 2r(z), h = w on ∂Kκ 2r(z). (3.23) Lemma 3.6. Let 0 < r ≤ min {ρ1 6 , ρ2 4 , (4e)−1Γ−( d+3 α +2), (Γ−1R) 1 α } , (3.24) where ρ1, ρ2 are the radi appearing in (3.7) and (3.20), respectively. For any given ε3 ∈ (0, 1) there exists a small constant δ = δ(d, λ,Λ, γ1, γ2, ∂Ω, ε3) > 0 such that − ∫ Kκ 2r(z) |Dw −Dh|pz dx dt ≤ ε3κ and − ∫ Kκ 2r(z) |Dh|pz dx dt ≤ c3κ (3.25) for some c3 = c3(d, λ,Λ, γ1, γ2, ∂Ω) > 0, see [10, Lemma 4.2]. Lemma 3.7. For each ε4 ∈ (0, 1), there exists a small constant δ > 0, δ = δ(d, λ,Λ, γ1, γ2, ε4), such that for the weak solution v ∈W 1,pz (Qκ+ 2r (z)) of vt − div(b(t,Dv)) = 0 in Qκ+ 2r (z), v = 0 on Tκ2r(z), (3.26) it holds − ∫ Qκ+ 2r (z) |Dv|pz dx dt ≤ c3κ and − ∫ Kκ r (z) |Dh−Dv̄|pz dx dt ≤ ε4κ, 14 H. TIAN, S. ZHENG EJDE-2020/13 where c3 is defined in Lemma 3.6. Here, we extend v from Qκ+ 2r (z) to Kκ 2r(z) by zero-extension denoted it by v̄, see [10, Lemma 4.3]. We also recall the L∞-estimate for the gradients of weak solution to the limiting problem of general p-Laplacian type with the nonlinearity independent of the spatial variable. Indeed, DiBenedetto showed an interior gradient bound for parabolic systems, see [15, Theorems 5.1 and 5.2], and Lieberman [23] extended it up to the boundary case for parabolic equations. Lemma 3.8. (i) (interior case) For a fixed κ > 1 and r > 0, we suppose that v ∈W 1,pz (Qκ2r(z)) is any weak solution of vt − div(b(t,Dv)) = 0 in Qκ2r(z) ⊂ ΩT with − ∫ Qκ2r(z) |Dv|pz dx dt ≤ c∗κ for some c∗ > 1. Then ‖Dv‖pzL∞(Qκr (z)) ≤ Cκ, (3.27) where C = C(d, λ,Λ, γ1, γ2, c∗) > 0. (ii) (boundary case) Let κ > 1 and r > 0, we suppose that v ∈W 1,pz (Qκ+ 2r (z)) is a weak solution of vt − div(b(t,Dv)) = 0 in Qκ+ 2r (z), v = 0 on Tκ2r(z) (3.28) with − ∫ Qκ+ 2r (z) |Dv|pz dx dt ≤ c∗κ for some c∗ > 1, then ‖Dv‖pz L∞(Qκ+ r (z)) ≤ Cκ, (3.29) where C = C(d, λ,Λ, γ1, γ2, c∗, ∂Ω) > 0. We finish this section by recalling the following two lemmas. Lemma 3.9. Let φ ∈ Φ be a Young function with φ ∈ ∆2 ∩ ∇2 and g ∈ Lφ(ΩT ). Then ∫ ΩT φ(|g|) dx dt = ∫ ∞ 0 ∣∣{(t, x) ∈ ΩT : |g| > k} ∣∣ dφ(k). Lemma 3.10. Let φ ∈ Φ be a Young function as shown in Lemma 3.9. Then, for any â, b̂ > 0 one has I = ∫ ∞ 0 1 κ (∫ {(t,x)∈ΩT :|g|>âκ} |g| dx dt ) dφ(b̂κ) ≤ C ∫ ΩT φ(|g|) dx dt, where C = C(â, b̂, φ), see [8, Lemma 3.4]. EJDE-2020/13 ORLICZ ESTIMATES 15 4. Proof of Theorem 2.5 Let us assume that p(t, x) ∈ SLH(ΩT ) with its range [γ1, γ2] shown as (2.7), (a,ΩT ) is (δ,R)-vanishing for R ∈ (0, 1) with a small δ ∈ (0, 1/8) such that the validity of Lemmas 3.4–3.7. Let the given datum |ψt|γ ′ 1 , |Dψ|p(t,x), |f|p(t,x) ∈ Lφ(ΩT ) for Young’s function φ ∈ ∆2 ∩ ∇2, and u ∈ A(ΩT ) be the weak solution of variational inequality (2.4) with the constants M , α, Γ, ca as in Lemma 3.1, m := sup(τ,y)∈ΩT m(τ, y) as (2.16), and R0 > 0 chosen as 0 < 2R0 ≤ min {ρ1 6 , ρ2 4 , (4e)−1Γ−( d+3 α +2), (Γ−1R) 1 α } , (4.1) ω(4R0) ≤ min { λ 4Λ , 1, ε0(γ1 − 1) 4 } , (4.2) where ρ1, ρ2 are shown in Lemma 3.1 and (3.20), ε0 > 0 as in (3.20). For any κ > 0, we set κ0 = ( − ∫ ΩT |Du|p(t,x) dx dt+ 1 δ − ∫ ΩT Ψ(t, x) dx dt )m , the upper-level set E(κ) = { (t, x) ∈ ΩT : |Du|p(t,x) > κ } , and for fixed (τ, y) ∈ ΩT and ρ > 0, J ( Kκ ρ (τ, y) ) = − ∫ Kκ ρ (τ,y) |Du|p(t,x) dx dt+ 1 δ − ∫ Kκ ρ (τ,y) Ψ(t, x) dx dt. Without loss of generality, we take a suitable positive constant K such that |ΩT | < |QKR0 |, where R0 > 0 is defined by (4.1) and (4.2). Step 1. We prove the modified Vitali covering for the major portion of E(κ) by a family of countably many disjoint cylinders. To this end,we have the following. Lemma 4.1. For κ ≥ κ1 := (( 2 1−δ )d (48χK)d+2 )m κ0, there exists a family of disjoint cylinders {Kκ ρi(τi, yi)}i≥1 with (τi, yi) ∈ E(κ) and 0 < ρi < min { κ pi−2 2pi , 1 } R0 48χ such that E(κ) ⊂ ( ∪i≥1 χK κ ρi(τi, yi) ) ∪ a negligible set, where the constant χ is shown as in Lemma 3.2, pi = p(τi, yi), and for each i ≥ 1 it holds J ( Kκ ρi(τi, yi) ) = κ, J ( Kκ ρ (τi, yi) ) < κ for all ρ ∈ ( ρi, min{κ pi−2 2pi , 1}R0 2 ] . Proof. For every fixed point z0 = (τ0, y0) ∈ E(κ), we consider the radius ρ with min{κ p0−2 2p0 , 1}R0 48χ ≤ ρ ≤ min{κ p0−2 2p0 , 1}R0 2 , (4.3) 16 H. TIAN, S. ZHENG EJDE-2020/13 where χ is as in Lemma 3.2 and p0 = p(z0). It is clear that for any z0 ∈ E(κ) it holds J (ρ) < κ. Indeed, it follows from the measure density conditions (2.12) that J(Kκ ρ ) = 1 |Kκ ρ | (∫ Kκ ρ |Du|p(t,x) dx dt+ 1 δ ∫ Kκ ρ Ψ(t, x) dx dt ) ≤ |Qκρ | |Qκρ ∩ ΩT | |ΩT | |Qκρ | ( − ∫ ΩT |Du|p(t,x) dx dt+ 1 δ − ∫ ΩT Ψ(t, x) dx dt ) < ( 2 1− δ )d |QKR0 | |Qκρ | κ 1 m 0 = ( 2 1− δ )d(KR0 ρ )d+2 κ p0−2 p0 κ 1 m 0 . We now divide it into the cases 2 ≤ p0 < γ2 and γ1 ≤ p0 < 2. If 2 ≤ p0 ≤ γ2, we ob- tain that m = sup(τ,y)∈ΩT m(τ, y) ≥ m(τ0, y0) = p0 2 by (2.16) and min{κ p0−2 2p0 , 1} = 1. Therefore, J(Kκ ρ ) < ( 2 1− δ )d (48χK)d+2κ p0−2 p0 κ 1 m 0 ≤ κ p0−2 p0 κ 2 p0 = κ; If γ1 ≤ p0 < 2, one gets that m = sup(τ,y)∈ΩT m(τ, y) ≥ m(τ0, y0) = 2p0 p0(d+2)−2d by (2.16) and min{κ p0−2 2p0 , 1} = κ p0−2 2p0 . This implies that J(Kκ ρ ) < ( 2 1− δ )d (48χK)d+2κ 2−p0 2p0 (d+2)κ p0−2 p0 κ 1 m 0 ≤ κ (2−p0)d 2p0 κ p0(d+2)−2d 2p0 = κ. In summary, J ( Kκ ρ ) < κ for all ρ ∈ [ min{κ p0−2 2p0 , 1}R0/(48χ),min{κ p0−2 2p0 , 1}R0/2 ] . (4.4) On the other hand, by the Lebesgue differentiation theorem we infer that lim ρ→0 J ( Kκ ρ ) ≥ |Du(z0)|p0 > κ. Consequently, one can select a maximal radius ρ0 ∈ ( 0,min{κ p0−2 2p0 , 1}R0/(48χ) ] by the intermediate value theorem such that J ( Kκ ρ0 ) = κ and J ( Kκ ρ ) < κ for all ρ ∈ ( ρ0,min{κ p0−2 2p0 , 1}R0/2 ] . Now, let us take {Kκ ρz (z) : z = (τ, y) ∈ E(κ)} as a covering of E(κ), and note that κ (48χ)d+2 ≤ − ∫ 48χKκ ρz (z) |Du|p(t,x) dx dt+ 1 δ − ∫ 48χKκ ρz (z) Ψ(t, x) dx dt ≤ κ. (4.5) Therefore, by taking M1 = (48χ)d+2 > 1 as in Lemma 3.1, we have κp + z −p − z ≤ ca for all z ∈ E(κ), where ca > 1 is as in Lemma 3.1, p+ z = supKκ ρz (z) p(t, x) and p−z = infKκ ρz (z) p(t, x). Finally, by employing the Vitali’s covering lemma 3.2 we can find a family of disjoint cylinders {Kκ ρi(τi, yi)}i≥1 with (τi, yi) ∈ E(κ) and ρi ∈ ( 0,min{κ pi−2 2pi , 1}R0/(48χ) ] , which reached the desired result. � Step 2. We are now in a position to show a suitable decay estimate to each of the above-mentioned covering {Kκ ρi(τi, yi)}i≥1. EJDE-2020/13 ORLICZ ESTIMATES 17 Lemma 4.2. Under the same hypotheses as in Lemma 4.1, we have |Kκ ρi(zi)| ≤ 2 κ (∫ Kκ ρi (zi)∩{z∈ΩT :|Du|p(z)>κ 4 } |Du|p(z) dx dt + 1 δ ∫ Kκ ρi (zi)∩{z∈ΩT :Ψ(t,x)> δκ 4 } Ψ(t, x) dx dt ) . Proof. By Lemma 4.1 we have − ∫ Kκ ρi (τi,yi) |Du|p(t,x) dx dt+ 1 δ − ∫ Kκ ρi (τi,yi) Ψ(t, x) dx dt = κ, which implies that κ|Kκ ρi(τi, yi)| = ∫ Kκ ρi (τi,yi) |Du|p(t,x) dx dt+ 1 δ ∫ Kκ ρi (τi,yi) Ψ(t, x) dx dt. Now we split the two integrals above in two parts, κ|Kκ ρi(τi, yi)| ≤ ∫ Kκ ρi (τi,yi)∩{(t,x)∈ΩT :|Du|p(t,x)>κ 4 } |Du|p(t,x) dx dt+ κ 4 |Kκ ρi(τi, yi)| + 1 δ ∫ Kκ ρi (τi,yi)∩{(t,x)∈ΩT :Ψ(t,x)> δκ 4 } Ψ(t, x) dx dt+ κ 4 |Kκ ρi(τi, yi)|, which yields the desired result. � Based on Lemma 4.1, we constructed a family of disjoint cylinders {Kκ ρi(zi)}i≥1 for zi = (τi, yi) ∈ E(κ) with 0 < ρi < min { κ pi−2 2pi ,1 } R0 48χ . We denote K0 zi = Kκ ρzi (zi), K1 zi = χKκ ρi(zi), K2 zi = 2χKκ ρi(zi), K3 zi = 4χKκ ρi(zi), K4 zi = 6χKκ ρi(zi), Q4 zi = 6χQκρi(zi); and consider the following estimates by parting the settings of Q4 zi ⊂ ΩT and Q4 zi 6⊂ ΩT . Case 1. For the interior case Q4 zi ⊂ ΩT , let k, w and h be the unique solution to the initial-boundary value problems (3.15), (3.17) and (3.23) with Q4 zi , Q 3 zi and Q2 zi instead of Kκ 6r, K κ 4r and Kκ 2r, respectively. With the same argument as the estimate (4.5), it holds κ (6χ)d+2 ≤ − ∫ Qjzi |Du|p(t,x) dx dt+ 1 δ − ∫ Qjzi Ψ(t, x) dx dt ≤ κ (4.6) for j ∈ {0, 1, 2, 3, 4}. Note that 0 < 6χρi ≤ min{κ pi−2 2pi , 1}R0 8 ≤ R0, (4.7) we take M2 = (6χ)d+2 > 1 in Lemma 3.1, and obtain p+ i − p − i ≤ ω(Γ(6χρi) α), κ 2 pi ≤ Γ2(6χρi) −(d+2) and κp + i −p − i ≤ ca (4.8) for zi = (τi, yi) ∈ E(κ) and 0 < ρi < min{κ pi−2 2pi ,1}R0 48χ , where pi = p(zi) = p(τi, yi), p − i = infQ4 zi p(t, x), p+ i = supQ4 zi p(t, x) with Γ ≥ 4, α ∈ (0, 1) and ca > 1 18 H. TIAN, S. ZHENG EJDE-2020/13 being the same as in Lemma 3.1. Let us now replace (3.14) by (4.6), and make use of the argument of Proposition 3.10 in [7] by Bui and Duong for general parabolic equation of p(t, x)-growth. Then we only make very slightly modifications to our problem, and immediately conclude Corollary 4.3. Assume that u ∈ A(ΩT ) is a weak solution of (2.4), and hi ∈ W 1,pi(Q2 zi) is a weak solution of (3.23). Then, for any ε > 0 there is a small δ = δ(d, λ,Λ, γ1, γ2, ε) > 0 with (4.6) such that − ∫ Q1 zi |Du−Dhi|p(t,x) dx dt ≤ εκ and ‖Dhi‖piL∞(Q1 zi ) ≤ N1κ, (4.9) where N1 = N1(d, λ,Λ, γ1, γ2, c∗) > 1 is the positive constant independent of the index i. We omit its proof, which is very similar to that of following boundary setting, but a simple process. Case 2. For the boundary case Q4 zi 6⊂ ΩT , we suppose that dist{yi, ∂Ω} ≤ 6χρi, and take y′i ∈ ∂Ω with |yi − y′i| ≤ 8χρi. Since 0 < 48χρi < min{κ pi−2 2pi , 1}R0 ≤ R0, by Assumption 2.1 there exists a new spatial coordinate system, still denoting x = {x1, . . . , xd}-coordinate, with the origin at y′i = 0′ such that B48χρi(0 ′) ∩ {x1 > 48χρiδ} ⊂ B48χρi(0 ′) ∩ Ω ⊂ B48χρi(0 ′) ∩ {x1 > −48χρiδ}. Since 0 < δ < 1/8, it leads to B40χρi(48χρiδe1) ⊂ B48χρi(0 ′) for e1 = {1, 0, . . . , 0}. We then translate the spatial coordinate system to the x1-direction by 48χρiδ, and denote the new origin by 48χρiδe1 = 0, so that it yields B+ 40χρi (0) ⊂ B40χρi(0) ∩ Ω ⊂ B40χρi(0) ∩ {x1 > −96χρiδ}. By considering this transformation is composed of only the translation and the rotation, it leads to that the basic structure of the problems (2.4) and the main assumptions are invariant. Here, we will continuously use the original symbols and notations in this new coordinate system. Since |yi| ≤ |yi − y′i| + |y′i| ≤ 6χρi + 48χρiδ ≤ 12χρi, for zi = (τi, 0) we obtain K1 zi ⊂ K 5 zi ⊂ K 6 zi ⊂ K 7 zi ⊂ K 8 zi , (4.10) where K5 zi = 15χKκ ρi(zi), K6 zi = 28χKκ ρi(zi), K7 zi = 38χKκ ρi(zi), K8 zi = 48χKκ ρi(zi). Similarly, we now let that k, w, h and v is a unique solution to the initial-boundary value problems (3.15), (3.17), (3.23) and (3.26) with K8 zi , K 7 zi , K 6 zi and Q6+ zi instead of Kκ 6r, K κ 4r, K κ 2r and Qκ+ 2r , respectively, where Q6+ zi = 28χQκ+ ρi (zi). With the same argument as for (4.5), it holds κ (48χ)d+2 ≤ − ∫ Kj zi |Du|p(t,x) dx dt+ 1 δ − ∫ Kj zi Ψ(t, x) dx dt ≤ κ (4.11) for each j = 0, 1, . . . , 8. By taking M = (48χ)d+2 > 1 in Lemma 3.1, we obtain p+ i − p − i ≤ ω(Γ (48χρi) α), κ 2 pi ≤ Γ2 (48χρi) −(d+2), κp + i −p − i ≤ ca (4.12) EJDE-2020/13 ORLICZ ESTIMATES 19 for every zi ∈ E(κ), 0 < ρi < min{κ pi−2 2pi ,1}R0 48χ , pi = p(zi), p − i = infK8 zi p(t, x), p+ i = supK8 zi p(t, x); where Γ ≥ 4, α ∈ (0, 1) and ca > 1 are as in Lemma 3.1. Similar to the argument in [7, Corollary 3.9], we have the following result. Corollary 4.4. Assume that u ∈ A(ΩT ) is the weak solution of (2.4). For any ε > 0, there exist small constant δ = δ(d, λ,Λ, γ1, γ2, ε) > 0 with (4.11) and a weak solution vi ∈W 1,pi(Q6+ zi ) of (3.26) such that − ∫ K1 zi |Du−Dv̄i|p(t,x) dx dt ≤ εκ and ‖Dv̄i‖piL∞(K1 zi ) ≤ N2κ, (4.13) where N2 = N2(d, λ,Λ, γ1, γ2, c∗, δ, R) > 1 is the constant independent of i. Here, we extend vi from Q6+ zi to Q6 zi by zero extension, and denote it by v̄i. Proof. Let us begin with the fact that |Dw −Dv̄i|pi ≤ 2γ2−1 (|Dw −Dh|pi + |Dh−Dv̄i|pi) , then, by Lemmas 3.6 and 3.7 we obtain − ∫ K5 zi |Dw −Dv̄i|pi dx dt ≤ C1 ( − ∫ K5 zi |Dw −Dh|pi dx dt+− ∫ K5 zi |Dh−Dv̄i|pi dx dt ) ≤ C1 (ε3 + ε4) κ with ε3 and ε4 being the same as in Lemmas 3.6 and 3.7, respectively, where C1 is a positive constant depending only on d, λ,Λ, γ1 and γ2. This together with Hölder’s inequality implies − ∫ K5 zi |Dw −Dv̄i|p(t,x) dx dt = − ∫ K5 zi |Dw −Dv̄i| pi 2 |Dw −Dv̄i|p(t,x)− pi2 dx dt ≤ ( − ∫ K5 zi |Dw −Dv̄i|pi dx dt )1/2( − ∫ K5 zi |Dw −Dv̄i|2p(t,x)−pi dx dt )1/2 ≤ C2 ( (ε3 + ε4)κ )1/2( − ∫ K5 zi |Dw|2p(t,x)−pi dx dt +− ∫ K5 zi |Dv̄i|2p(t,x)−pi dx dt )1/2 . (4.14) Next, by using (4.2) and (4.12) we have 2p(t, x)− pi ≤ p(t, x) ( 1 + p+ i − p − i ) ≤ p(t, x) ( 1 + ω ( Γ (48χρi) α )) ≤ p(t, x) (1 + ε0) in K5 zi where ε0 is the same as in the inequality (3.20). Therefore, from (3.20) it follows that − ∫ K5 zi |Dw|2p(t,x)−pi dx dt ≤ − ∫ K5 zi |Dw|p(t,x)(1+ω(Γ(48χρi) α)) dx dt+ 1 ≤ κ1+ω(Γ (48χρi) α) + 1. (4.15) 20 H. TIAN, S. ZHENG EJDE-2020/13 By (4.12), κω(Γ(48χρi) α) ≤ ca for ca > 1, which leads to − ∫ K5 zi |Dw|2p(t,x)−pi dx dt ≤ C3κ. (4.16) Using Lemma 3.8-(ii) for the weak solution vi ∈W 1,pi(Q6+ zi ) mentioned in Corollary 4.4, and formula (4.12) for κ > 1 we obtain − ∫ K5 zi |Dv̄i|2p(t,x)−pi dx dt = 1 |K5 zi | ∫ K5 zi |Dv̄i|2p(t,x)−pi dx dt = 1 |K5 zi | (∫ Q5+ zi |Dvi|2p(t,x)−pi dx dt+ |K5 zi\Q 5+ zi | ) ≤ C4 |K5 zi | ( sup K5 zi κ 2p(t,x)−pi pi |Q5+ zi |+ |K 5 zi\Q 5+ zi | ) ≤ C5 |K5 zi | ( κκ 2(p + i −p− i ) γ1 |Q5+ zi |+ |K 5 zi\Q 5+ zi | ) ≤ C6κ. (4.17) Now, combining (4.16) and (4.17) with (4.14) yields − ∫ K5 zi |Dw −Dv̄i|p(t,x) dx dt ≤ C7 ( ε3 + ε4 )1/2 κ. (4.18) Then, it follows from Lemmas 3.4 and 3.5 that − ∫ K5 zi |Du−Dv̄i|p(t,x) dx dt ≤ 3γ2−1 ( − ∫ K5 zi |Du−Dk|p(t,x) dx dt+− ∫ K5 zi |Dk −Dw|p(t,x) dx dt +− ∫ K5 zi |Dw −Dv̄i|p(t,x) dx dt ) ≤ C ( ε1 + ε2 + (ε3 + ε4)1/2 ) κ with εi > for i = 1, 2, 3, 4 being the same as in Lemmas 3.4–3.7, respectively, which means that − ∫ K1 zi |Du−Dv̄i|p(t,x) dx dt ≤ C10 ( ε1 + ε2 + (ε3 + ε4)1/2 ) κ. Therefore, by taking εi > 0, i = 1, 2, 3, 4 sufficiently small we ensure the validity of first inequality of (4.13). The second inequality of (4.13) is proved by following from Lemma 3.8 (ii). � Step 3. We now prove a decay estimate of the upper-level set for a variable power of the gradients of weak solution to (2.4). To this end, let N = max{N1, N2} > 1 with N1, N2 as shown in Corollary 4.3 and Corollary 4.4, and let A = N γ2 γ1 c 1 γ1 a 2γ2 > 1. Lemma 4.5. Let κ ≥ κ1 be as in Lemma 4.1. For any ε ∈ (0, 1), there exists a small constant δ = δ(d, λ,Λ, γ1, γ2) > 0 with δ ∈ (0, 1/8) such that if (a,ΩT ) EJDE-2020/13 ORLICZ ESTIMATES 21 satisfies (δ,R)-vanishing for some 0 < R < 1, p(t, x) ∈ SLH with the range [γ1, γ2], and u ∈ A(ΩT ) is the weak solution of (2.4). Then we have the estimate |E(Aκ)| ≤ Cε { 2 κ (∫ {(t,x)∈ΩT :|Du|p(t,x)>κ 4 } |Du|p(t,x) dx dt + 1 δ ∫ {(t,x)∈ΩT :Ψ(t,x)> δκ 4 } Ψ(t, x) dx dt )} , where C = C(d, λ,Λ, γ1, γ2, R, ca). Proof. Considering the fact E(Aκ) ⊂ E(κ) for A > 1, by Lemma 4.1 we obtain that the family {K1 zi}i≥1 can cover almost all E(Aκ), which implies that |E(Aκ)| = ∣∣{(t, x) ∈ ΩT : |Du|p(t,x) > Aκ }∣∣ ≤ ∞∑ i=1 |{(t, x) ∈ K1 zi : |Du|p(t,x) > Aκ}| = ∑ i:interior case |{(t, x) ∈ K1 zi : |Du|p(t,x) > Aκ}| + ∑ i:boundary case |{(t, x) ∈ K1 zi : |Du|p(t,x) > Aκ}|. (4.19) For the interior estimate, by Corollary 4.3 it yields sup K1 zi |Dhi|p(t,x) ≤ sup K1 zi N p(t,x) pzi 1 κ p(t,x) pzi ≤ N γ2 γ1 1 κ p + 1 −p − 1 γ1 κ ≤ N γ2 γ1 1 c 1 γ1 a κ. Recalling 2d d+2 < γ1 ≤ p(t, x) ≤ γ2 <∞ in ΩT leads to |Du|p(t,x) ≤ 2γ2−1 ( |Du−Dhi|p(t,x) + |Dhi|p(t,x) ) for all (t, x) ∈ K1 zi ⊂ ΩT . Therefore, |{(t, x) ∈ K1 zi : |Du|p(t,x) > Aκ}| ≤ |{(t, x) ∈ K1 zi : ( |Du−Dhi|p(t,x) + |Dhi|p(t,x) ) > 21−γ2Aκ}| ≤ |{(t, x) ∈ K1 zi : |Du−Dhi|p(t,x) > N γ2 γ1 1 c 1 γ1 a κ}| ≤ 1 N γ2 γ1 1 c 1 γ1 a κ ∫ K1 zi |Du−Dhi|p(t,x) dx dt ≤ |K1 zi |ε N γ2 γ1 1 c 1 γ1 a , which implies |{(t, x) ∈ K1 zi : |Du|p(t,x) > Aκ}| ≤ C1ε|Kκ ρi(zi)|. (4.20) For the boundary case, we carry out the same procedure as the estimate of (4.20), and use Corollary 4.4 to discover that |{(t, x) ∈ K1 zi : |Du|p(t,x) > Aκ}| ≤ C2ε|Kκ ρi(zi)|. (4.21) 22 H. TIAN, S. ZHENG EJDE-2020/13 Putting (4.20) and (4.21) into (4.19) yields |E(Aκ)| ≤ C3ε ∞∑ i=1 |Kκ ρi(zi)|. Thanks to Lemma 4.2, |E(Aκ)| ≤ C3ε ∞∑ i=1 { 2 κ (∫ Kκ ρi (zi)∩{z∈ΩT :|Du|p(t,x)>κ 4 } |Du|p(t,x) dx dt + 1 δ ∫ Kκ ρi (zi)∩{z∈ΩT :Ψ(t,x)> δκ 4 } Ψ(t, x) dx dt )} . Note that {Kκ ρi(zi)} are non-overlapping in ΩT , then the required result follows. Step 4. The step is devoted to Orlicz estimate for the derivatives of weak solution. Using Lemma 3.9, we have∫ ΩT φ(|Du|p(t,x)) dx dt = ∫ ∞ 0 |{(t, x) ∈ ΩT : |Du|p(t,x) > Aκ}| dφ(Aκ) = ∫ κ1 0 |{(t, x) ∈ ΩT : |Du|p(t,x) > Aκ}| dφ(Aκ) + ∫ ∞ κ1 |{(t, x) ∈ ΩT : |Du|p(t,x) > Aκ}| dφ(Aκ) := J1 + J2. (4.22) Frist we estimate of J1. Recalling the above definitions of κ0 and κ1, and using that φ ∈ ∆2 ∩∇2, we have J1 ≤ |ΩT |φ(Aκ1) = |ΩT |φ ( N γ2 γ1 c 1 γ1 a 2γ2 (( 2 1− δ )d (48χK)d+2 )m κ0 ) ≤ C1|ΩT |φ (κ0) = C1|ΩT |φ (( − ∫ ΩT |Du|p(t,x) dx dt+ 1 δ − ∫ ΩT Ψ(t, x) dx dt )m) . Then by (2.8) it follows that J1 ≤ C2φ (( − ∫ ΩT Ψ(t, x) dx dt )m) , where C2 = C2(d, γ1, γ2, λ,Λ, α1, α2, δ, R, T, |Ω|), and Ψ(t, x) is defined by (3.3). Now we estimate J2. From Lemmas 4.5 and 3.10 we observe that J2 ≤ C3ε ∫ ∞ 0 ( 2 κ (∫ {z∈ΩT :|Du|p(z)>κ 4 } |Du|p(z) dx dt + 1 δ ∫ {z∈ΩT :Ψ(z)> δκ 4 } Ψ(z) dx dt )) dφ(Aκ) ≤ C4ε (∫ ΩT φ(|Du|p(t,x)) dx dt+ ∫ ΩT φ(Ψ(t, x)) dx dt ) . Inserting the estimates for J1 and J2 into (4.22),∫ ΩT φ(|Du|p(t,x)) dx dt EJDE-2020/13 ORLICZ ESTIMATES 23 ≤ C4ε ∫ ΩT φ(|Du|p(t,x)) dx dt+ C2φ (( − ∫ ΩT Ψ(t, x) dx dt )m) + C5 ∫ ΩT φ(Ψ(t, x)) dx dt. If ∫ ΩT φ(|Du|p(t,x)) dx dt <∞, then we can select ε > 0 sufficiently small such that 0 < C4ε < 1/2, then∫ ΩT φ(|Du|p(t,x)) dx dt ≤ C6 ( φ (( − ∫ ΩT Ψ(t, x) dx dt )m) + ∫ ΩT φ(Ψ(t, x)) dx dt ) , (4.23) where C6 = C6(d, γ1, γ2, λ,Λ, α1, α2, δ, R, T, |Ω|). Otherwise, the integral on the left may be +∞, and wee need to refine the estimate for |Du|. Let us consider the truncated gradients |Du(t, x)|p(t,x) λ̄ := min { |Du(t, x)|p(t,x), λ̄ } for (t, x) ∈ ΩT and λ̄ ∈ [κ1,∞), and set Eλ̄(κ) = {(t, x) ∈ ΩT : |Du|p(t,x) λ̄ > κ}. By Lemma 4.5 we obtain |Eλ̄(Aκ)| ≤ Cε { 2 κ (∫ {(t,x)∈ΩT :|Du|p(t,x) λ̄ >κ 4 } |Du|p(t,x) λ̄ dx dt + 1 δ ∫ {(t,x)∈ΩT :Ψ(t,x)> δκ 4 } Ψ(t, x) dx dt )} . In the case λ̄ ≤ Aκ, we have |Eλ̄(Aκ)| = 0 so that (4.23) holds trivially. Otherwise, while λ̄ > Aκ, working exactly as in the previous lines, we obtain the inequality∫ ΩT φ(|Du|p(t,x) λ̄ ) dx dt ≤ C6 { φ (( − ∫ ΩT Ψ(t, x) dx dt )m) + ∫ ΩT φ(Ψ(t, x)) dx dt } , (4.24) instead of (4.23). Taking λ̄→∞ and using the lower semi-continuity of Orlicz norm with respect to almost everywhere convergence, we obtain (4.23). This completes the proof. � As a direct consequence of Theorem 2.5, by taking φ(ν) = νq for q ∈ (1,∞) we conclude the classical Calderón-Zygmund theory for parabolic obstacle problems with p(t, x)-growth. Corollary 4.6. Let q ∈ (1,∞). Assume that p(t, x) ∈ SLH with the range [γ1, γ2], and (a,R × Ω) satisfies (δ,R)-vanishing. If u ∈ A(ΩT ) is a weak solution of the variational inequality (2.4) with |ψt|γ ′ 1 , |Dψ|p(t,x), |f|p(t,x) ∈ Lq(ΩT ), then we have |Du|p(t,x) ∈ Lq(ΩT ) with the estimate∫ ΩT |Du|p(t,x)q dx dt ≤ C (∫ ΩT ( Ψ(t, x) )q dx dt+ 1 )m , 24 H. TIAN, S. ZHENG EJDE-2020/13 where C = C(d, γ1, γ2, λ,Λ, δ, R, T, |Ω|, ‖ua‖L2(Ω)), Ψ(t, x) and m ≥ 1 are the same as in Theorem 2.5. Acknowledgements. H. Tian was supported by the National Natural Science Foundation of China Youth Fund, grant no. 11901429. References [1] Acerbi, E.; Mingione, G.; Gradient estimates for the p(x)-Laplacian system, J. Reine Angew. Math., 584 (2002), 117–148. [2] Acerbi, E.; Mingione, G.; Gradient estimates for a class of parabolic systems, Duke Math. J., 136 (2007), 285–320. [3] Baroni, P.; Lorentz estimates for obstacle parabolic problems, Nonlinear Anal., 96 (2014), 167–188. [4] Baroni, P.; Bögelein, V.; Calderón-Zygmund estimates for parabolic p(x, t)-Laplacian Sys- tems, Rev. Mat. Iberoam., 30 (2014), 1355–1386. [5] Bögelein, V.; Duzaar, F.; Higher integrability for parabolic systems with nonstandard growth and degenerate diffusions, Publicacions Matemátiques, 55 (2011), 201–250. [6] Bögelein, V.; Duzaar, F.; Mingione, G.; Degenerate problems with irregular obstacles, J. Reine Angew. Math., 650 (2011), 107–160. [7] Bui, T. A.; Duong, X. T.; Weighted Lorentz estimates for parabolic equations with non- standard growth on rough domains, Calc. Var., 56 (177) (2017), https: //doi.org/ 10.1007 /s00526-017-1273-y. [8] Byun, S. S.; Cho, Y.; Nonlinear gradient estimates for parabolic problems with irregular obstacles, Nonlinear Anal., 94 (2014), 32–44. [9] Byun, S. S.; Cho, Y.; Nonlinear gradient estimates for parabolic obstacle problems in non- smooth domains, Manuscripta Math., 146 (2015), 539–558. [10] Byun, S. S.; Ok, J.; Nonlinear parabolic equations with variable exponent growth in non- smooth domains, SIAM J. Appl. Math., 48 (2016), 3148–3190. [11] Byun, S. S.; Park, J. T.; Global weighted Orlicz estimates for parabolic measure data prob- lems: Application to estimates in variable exponent spaces, J. Math. Anal. Appl., 467 (2018), 1194–1207. [12] Chlebicka, I.; Gradient estimates for problems with Orlicz growth, Nonlinear Anal., (2018), https://doi.org/10.1016/j.na.2018.10.008. [13] DiBenedetto, E.; Friedman, A.; Regularity of solutions of nonlinear degenerate parabolic systems, J. Reine Angew. Math., 349 (1984), 83–128. [14] DiBenedetto, E.; Manfredi, J.; On the higher integrability of the gradient of weak solutions of certain degenerate elliptic systems, Amer. J. Math., 115 (1993), 1107–1134. [15] DiBenedetto, E.; Degenerate parabolic equations, Universitext Springer, New York, 1993. [16] Erhardt, A.; Existence and gradient estimates in parabolic obstacle problems with nonstan- dard growth, Dissertationsschrift, Universitat Erlangen, 2013. [17] Iwaniec, T.; Projections onto gradient fields and Lp-estimates for degenerate elliptic equa- tions, Studia Math., 75 (1983), 293–312. [18] Li, H. Z.; Zhang, J. J.; Zheng, S. Z.; Orlicz estimates for nondivergence linear elliptic equa- tions with partially BMO coefficients, Complex Var. Elliptic Equ., 63 (6) (2018), 871–885. [19] Liang, S.; Zheng, S. Z.; Gradient estimate in Orlicz spaces for elliptic obstacle problems with partially BMO nonlinearities, Electron. J. Differential Equations, 2018 (58) (2018), 1–15. [20] Liang, S.; Zheng, S. Z.; Feng, Z.; Variable Lorentz estimate for generalized Stokes systems in non-smooth domains, Electron. J. Differential Equations, 2019 (109) (2019), 1–29. [21] Jia, H. L.; Li, D. S.; Wang, L. H.; Global regularity for divergence form elliptic equations in Orlicz spaces on quasiconvex domains, Nonlinear Anal., 74 (2011), 1336–1344. [22] Kinnunen, J.; Zhou, S. L.; A local estimate for nonlinear equations with discontinuous coef- ficients, Comm. Partial Differential Equations, 24 (1999), 2043–2068. [23] Lieberman, G. M., Boundary regularity for solutions of degenerate parabolic equations, Non- linear Anal., 14 (1990), 501–524. [24] Li, Q. F.; Very weak solutions of subquadratic parabolic systems with non-standard p(t, x)- growth, Nonlinear Anal., 156 (2017), 17–41. [25] Rao, M. M.; Ren, D. Z.; Applications of Orlicz Spaces, New York, Marcel Dekker Inc, 2002. EJDE-2020/13 ORLICZ ESTIMATES 25 [26] Růžička, M., Electrorheological Fluids: Modeling and Mathematical Theory, Springer Lecture Notes in Math, Vol. 1748, Springer-Verlag, Berlin, Heidelberg, New York, 2000. [27] Tian, H.; Zheng, S. Z.; Global weighted Lorentz estimates to nonlinear parabolic equations over nonsmooth domains, J. Math. Anal. Appl., 456 (2017), 1238–1260. [28] Zhang, J. J.; Zheng, S. Z.; Lorentz estimates for asymptotically regular elliptic equations in quasiconvex domains, Electron. J. Differential Equations, 2016 (142) (2016), 1–13. Hong Tian College of Science, Tianjin University of Technology, Tianjin 300384, China. Department of Mathematics, Beijing Jiaotong University, Beijing 100044, China Email address: 14118404@bjtu.edu.cn Shenzhou Zheng (corresponding author) Department of Mathematics, Beijing Jiaotong University, Beijing 100044, China Email address: shzhzheng@bjtu.edu.cn 1. Introduction 2. Minimal assumptions on the datum and main result 3. Comparison estimates to the reference problems 4. Proof of Theorem ?? Acknowledgements References