The flow of immiscible fluids in porous media is a common process in many natural and engineering settings, such as geological carbon sequestration [1, 2], enhanced oil recovery [3, 4], and groundwater contamination [5, 6]. Understanding and controlling the instability of the fluid–fluid interface are beneficial for increasing the efficiency of oil recovery and CO₂ sequestration [7–13]. Under quasi-static conditions, the geometry of the porous media is one of the key factors that control fluid invasion. The disorder effect has been sufficiently studied, and we now have a good understanding of the invasion process in disordered porous media [14–24]. Generally, decreasing disorder stabilizes the invasion front for drainage and imbibition conditions [17, 22, 24, 25]. In addition to the disordered porous media, fluid invasion in ordered porous media is also common in many practical applications, such as the removal of water in gradient gas diffusion layers for PEM fuel cells [26, 27] and flow in multi-layered soils [28]. The ordered porous medium signifies that the pore size remains constant in the transverse direction and changes linearly along the flow direction. In such a porous medium, it has been proven that the monotonous change of the pore size along the flow direction can suppress the viscous fingering and capillary fingering or trigger a single fingering [29–32]. The displacement patterns in the ordered porous media are remarkably different from those in the disordered porous media.

Wettability also significantly impacts the displacement pattern and is denoted by the contact angle θ of the invading phase [33–36]. To unravel the fundamentals of the invasion processes under various wetting conditions, Cieplak and Robbins [37, 38] proposed three instability events (burst, touch, and overlap) that describe the pore-scale meniscus motion. Generally, the burst and touch events solely depend on the local geometry of a single throat and are regarded as destabilizing the invasion front, while the overlap event involving the merging of the neighboring menisci suppresses branching of the fluid–fluid interface, thus stabilizing the invasion front [30, 39–45]. However, in an ordered porous medium with decreasing pore size, the burst event is reported to stabilize the invasion front and lead to a compact pattern [31]. This dramatic difference is due to the monotonical decrease of the pore size, which causes the burst event to be more likely to occur near the inlet. In addition, studies on the wettability effect provide a basic understanding that decreasing the contact angle of the invading phase in the range of 45°<θ < 180° will smooth the fluid–fluid interface [46–51]. However, for invasion in the narrow space with decreasing thickness between two parallel plates (a Hele–Shaw cell), inversely, decreasing contact angles will destabilize the invasion front [29]. Again, the impact of wettability on the invasion process in ordered porous media is significantly different from that in disordered porous media, which has not been well-understood.

In this study, we propose a phase diagram to directly predict displacement patterns with different pore size gradients and contact angles. For an ordered porous medium, the pore size remains constant in the same column and changes monotonously from the inlet to the outlet. This change in the pore size leads to the transition of the pore-filling events and eventually affects the displacement pattern. First, we calculate the transition of pore-filling events with different pore sizes and contact angles; second, we link the pore-filling events to local invasion morphologies and observe that, in the negative gradient porous media, the burst and overlap events lead to the compact front and tapered front, respectively. In the positive gradient porous media, however, the burst and touch events lead to the single finger or compact front, and the overlap event leads to the tapered front or widened fingering; third, the displacement patterns are obtained based on the combination of local invasion morphologies, and the theoretical model that predicts the phase diagram of displacement patterns is developed. The theoretical model is evaluated using pore network simulations in this work and experiments in the literature, proving that the theoretical model can reasonably predict displacement patterns in ordered porous media. Finally, we extend the predictive method to square-arranged porous media, indicating that this method can be extended to predict the displacement patterns in different lattice structures. This work improves the understanding of how pore size gradient and wettability impact the fluid invasion in ordered porous media and is of practical significance for multiphase flow control.

The pore-filling events proposed by Cieplak and Robbins [37, 38] can capture the meniscus motion in porous media and are the basis of the derivation of the theoretical model in this work. These modes are significantly affected by the geometry of porous media. For the ordered porous media considered in this work, the pore size remains constant in the same column and changes monotonously from the inlet to the outlet. This change in the pore size leads to the transition of the pore-filling events and eventually affects the displacement pattern.

Figure 1 shows the transition of pore-filling events and their critical radii of curvature on the d-θ plane; the derivation is shown in Supplementary Text S1. Here, we consider a simple three-post system with the same radius r ₀, which ranges from 0.1a to 0.4a, where a is the distance between the center of the neighboring posts. Therefore, the normalized pore size d = (a-2r ₀)/a ranges from 0.2 to 0.8. When the neighboring menisci form an angle with Φ = 180°, the burst event occurs at high contact angles, and the touch event and overlap event occur at low contact angles, where the touch event occupies the small pore size region and the overlap event occupies the large pore size region (Figure 1A). When the neighboring menisci form an angle with Φ = 120°, the touch event is completely covered by the overlap event since the overlap event always has a larger curvature radius. In this case, only the burst and overlap events occur (Figure 1B). To link pore-filling events to fluid invasion, we consider three specific contact angles θ = 50°, 90°, and 150° (dashed lines in Figure 1). For an ordered porous medium with a given pore size gradient along the flow direction, the pore size in each column can be obtained. If the contact angle is also determined, the type of the pore-filling events and their critical radius of curvature can be determined in each column, as shown in Figures 2A, 3A, 4A. Then, by analyzing the effect of pore-filling events on fluid invasion, the local invasion morphologies in each column can be determined, which will be discussed in the next subsection.

Since the pore size of ordered porous media remains constant in the same column, the same pore-filling events would occur in the local region and have similar impacts on local invasion morphologies. These impacts would change with the pore size from the inlet to the outlet. By analyzing the transition of pore-filling events and their impacts, the local invasion morphologies can be obtained.

We first link pore-filling events to local invasion morphologies at θ = 150°, where only the burst event occurs. Based on the contour plots in Figure 1, the variation of the normalized critical radius of curvature with the pore size is obtained (Figure 2A). Here, the critical radius of curvature of the burst event is calculated using Supplementary Equation S1, with r ₁ = r ₂ = a (1-d)/2. For a negative gradient porous medium, as shown in Figure 2B, the pore size decreases along the flow direction. Here, the decreasing radius of curvature causes the invading phase to fill all pores in a column before proceeding to the next column. Therefore, a compact front is obtained. For a positive gradient porous medium, as shown in Figure 2C, the pore size increases along the flow direction. In contrast, the critical radius of curvature of the meniscus increases from the inlet to the outlet, the invading phase preferentially filling the pore in the next column, leading to a single fingering.

We then link pore-filling events to local invasion morphologies at θ = 90°. Based on the contour plots in Figure 1, the pore-filling events can be divided into two zones. As Figure 3A shows, only the burst event occurs in zone Ⅰ (0.2 ≤ d ≤ 0.5), and both the burst event and overlap event occur in zone Ⅱ (0.5 < d ≤ 0.8). The critical radius of curvature of the burst event (black line in Figure 3A) and overlap event (red line in Figure 3A) can be determined using Supplementary Equations S1, S4, respectively. For a negative gradient porous medium, as Figure 3B shows, the invading phase invades from zone Ⅱ to zone Ⅰ. In the beginning, the fluid invasion is dominated by the burst event and overlap event, in which the overlap event always has a larger radius of curvature and will, therefore, occur first (Figure 3A). However, the overlap event needs the neighboring menisci to form an angle with Φ = 120°. This mechanism cannot solely cause the invading phase to proceed to the next column. When there is only one pore left in the given column, the burst event fills the last pore and leads to a compact front (Figure 3B1). Although both the burst event and overlap event can occur simultaneously, the local invasion morphology is determined by the burst event. After the invading phase reaches zone Ⅱ, only the burst event occurs. Again, the decreasing radius of curvature of the burst event leads to a compact front (Figure 3B2).

For a positive gradient porous medium, the invading phase invades from zone Ⅰ to zone Ⅱ. In the beginning, the increasing radius of curvature of the burst event leads to a single fingering (Figure 3C1). After the invading phase reaches zone Ⅱ, the occurrence of the overlap event widens the fingering, leading to a widened fingering (Figure 3C2). In contrast to negative gradient porous media, the impact of the burst event is covered by the overlap event in positive gradient porous media.

Finally, we link pore-filling events to local invasion morphologies at θ = 50°. In this case, the invasion can also be divided into two zones. As Figure 4A shows, the touch event and overlap event (Φ = 120°) occur in zone Ⅰ (0.2 ≤ d ≤ 0.36), and the overlap events with Φ = 120° and 180° occur in zone Ⅱ (0.36 < d ≤ 0.8). Here, the critical radius of curvature of the touch event (black solid line in Figure 4A) and overlap event (red line and black-dashed line in Figure 4A) can be determined using Supplementary Equations S2, S4, respectively. In negative gradient porous media, the invading phase invades from zone Ⅱ to zone Ⅰ. In the beginning, the invasion is dominated by the overlap events with Φ = 120° and 180°. For multiphase flow in the porous media, there are many menisci at the fluid–fluid interface. When the adjacent menisci in some places form an angle with Φ = 180°, the adjacent meniscus in other places may form an angle with Φ = 120°. Therefore, the Φ = 120° and Φ = 180° cases can simultaneously exist at the fluid–fluid interface. Since the overlap events (Φ = 120°) need the neighboring menisci to form an angle with Φ = 120°, it leads to a tapered front of decreasing width (Supplementary Figure S3). On the other hand, the radius of curvature of the overlap event with Φ = 180° does not change with the pore size (Figure 4A), which will not affect the tapered front determined by the overlap event with Φ = 120° (Figure 4B1). When the invading phase reaches zone Ⅰ, the fluid invasion is dominated by the touch event and overlap event. Since the overlap event has a larger radius of curvature than the touch event (Figure 4A), the pores are filled by the overlap event, leading to a tapered front (Figure 4B2). As a result, the overlap event with Φ = 120° can suppress the impacts of both the touch event and overlap event with Φ = 180° in negative gradient porous media.

In positive gradient porous media, the invading phase invades from zone Ⅰ to zone Ⅱ. In zone Ⅰ, the fluid invasion is dominated by the touch event and overlap event. Similar to the invasion process in Figure 3B1, the pores are first filled by the overlap event, and then the touch event fills the last pore, leading to a compact front (Figure 4C1). Therefore, the impact of the overlap event with Φ = 120° is covered by the touch event. After the invading phase reaches zone Ⅱ, again, the overlap event with Φ = 180° does not affect the local invasion morphology and the overlap event with Φ = 120° leads to a tapered front (Figure 4C2).

In summary, the local invasion morphologies in ordered porous media follow the following two rules:

1) The local invasion morphologies are determined by the pore-filling events. When the corresponding pore-filling events can occur at a given contact angle, first, the burst event leads to a compact front in negative gradient porous media (Figures 2B, 3C) and single fingering in positive gradient porous media (Figures 2D, 3C1). This dramatic change in the local invasion morphology is due to the opposite variation of the radius of curvature along the flow direction. Second, the touch event is suppressed by the overlap event (Φ = 120°) in negative gradient porous media. However, in positive gradient porous media, the touch event leads to a compact front (Figure 4C1). Finally, the overlap event (Φ = 120°) has similar impacts in both negative and positive gradient porous media. It leads to a tapered front when adjacent to a compact front and a widened fingering when adjacent to a single fingering channel.

2) When more than one pore-filling events occur simultaneously, the local invasion morphology is determined by a dominant pore-filling event. In negative gradient porous media, the local invasion morphology is first determined by the burst event (Figure 3B1) and then by the overlap event (Φ = 120°) (Figure 4B2). In contrast, in positive gradient porous media, it is first determined by the touch event (Figure 4C1), then by the overlap event (Φ = 120°) (Figure 3C2), and finally by the burst event (Figures 2C, 3C1).

Based on the aforementioned two conclusions and the transition of pore-filling events on the d-θ plane (Figure 1), the transition of the local invasion morphology on the d-θ plane is obtained, as shown in Figure 5. In negative gradient porous media, the burst event leads to a compact front (zone Ⅰ), and the overlap event leads to a tapered front (zone Ⅱ) (Figure 5A). In positive gradient porous media, however, the burst event leads to a single fingering (zone Ⅴ), while the touch event leads to a compact front (zone Ⅲ) (Figure 5B). For the overlap event, it leads to a tapered front or a widened fingering, depending on which local invasion morphology it is adjacent to (zone Ⅳ). Moreover, the displacement patterns in ordered porous media can be determined based on the combination of the local invasion morphologies, which will be discussed in the next section.

Section 3 describes how the pore-filling events affect the fluid invasion in the ordered porous media, providing a basic understanding of the relationship between the pore-filling events and local invasion morphologies. In this section, we link local invasion morphologies to displacement patterns and establish a theoretical model by combining the pore-filling events and local invasion morphologies to predict the displacement patterns in the ordered porous media.

To investigate the displacement pattern in the ordered porous media, we consider ordered porous media with a constant pore size at the inlet, with d _inlet = 0.2 for negative gradient porous media (Figure 6A) and d _inlet = 0.8 for positive gradient porous media (Figure 6B). The pore size changes linearly from the inlet to the outlet, with a gradient λ = (d _outlet-d _inlet)a/L, where d _outlet is the pore size in the outlet, a is the distance between the center of neighboring posts, and L is the length of the porous media.

Then, the displacement patterns in such ordered porous media are obtained by enumerating all combinations of local invasion morphologies. For negative gradient porous media, there are four combinations of local invasion morphologies, as shown in Figure 7A. Among all combinations, if the compact front occurs in the end part of the porous media (Figures 7A1,A2), the compact pattern (CP) is obtained (Figure 7B1), and if the tapered front occurs in the end part (Figures 7A3,A4), then the taper shape pattern (TS) is obtained (Figure 7B2). Here, the taper shape pattern (TS) represents the morphology, in which the width of the invaded region decreases along the flow direction. Since the compact front is caused by the burst event (Figure 5A), the transition from the compact pattern (CP) to the taper shape pattern (TS) in negative gradient porous media corresponds to the critical condition for the burst event to occur in the last column: θ C P _ T S ( λ ) :   R b , N − 1 = R N − 1 / 2 , (1) where R _b,N-1 is the critical radius of curvature of the burst event in the N-1th column, which can be determined by substituting r 1 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 1 ) / 4 and r 2 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 2 ) / 4 into Supplementary Equation S1, where r ₁ and r ₂ are the radii of posts 1 and 2 in Supplementary Figure S1A, respectively. R _N-1/2 is the critical radius of curvature for menisci in the tip of the invasion front, and N is the maximum number of columns. The derivation of the critical radius of curvature for menisci in the tip of the invasion front is presented in SI Text S2. As shown in Figure 8A, two lines (θ _{CP_TS1}(λ) and θ _{CP_TS2}(λ)) are obtained using Eq. 1 because different pore-filling events occur at the tip of the invasion front.

For positive gradient porous media, we have five combinations of local invasion morphologies, as shown in Figure 7C. First, if only the compact front occurs (Figure 7C1), the compact pattern (CP) is obtained (Figure 7D1), and second, if the tapered front occurs in the end part of the porous media (Figure 7c2, c3), the taper shape pattern (TS) is obtained (Figure 7D2). Since the compact front is caused by the touch event (Figure 5B), the transition from the compact pattern (CP) to the taper shape pattern (TS) in positive gradient porous media corresponds to the critical condition for the touch event to occur in the last column, which leads to θ C P _ T S ′ ( λ ) :   R t , N − 1 = R N − 1 / 2 (2) where R _t,N-1 is the critical radius of curvature of the touch event in the N-1th column, which can be determined by substituting r 1 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 1 ) / 4 and r 2 = r 3 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 2 ) / 4 into Supplementary Equation S2. Similarly, two lines ( θ C P _ T S 1 ′ ( λ ) and θ C P _ T S 2 ′ ( λ ) , Figure 8B) are obtained due to different pore-filling events occurring at the tip of the invasion front.

Third, when the single fingering starts to occur in the first column (Figure 7C4), displacement patterns translate from the taper shape pattern (TS, Figure 7D2) to the kite shape pattern (KS, Figure 7d3). This transition corresponds to the critical condition that the overlap event occurs in the first column: θ T S _ K S ( λ ) :   R o , 1 / 2,1 = R 3 / 2 , (3) where R _o,1/2,1 is the critical radius of curvature of the overlap event between the 1/2th column and 1st column, which is determined by substituting r 1 = r 2 = a ( 1 − d i n l e t ) / 2 , r 4 = a ( 1 − d i n l e t ) / 2 − 3 a λ / 4 , and Φ = 120° into Supplementary Equation S4 for θ ≤ 90° and into Supplementary Equation S5 for θ > 90°.

Finally, when only the single fingering occurs (Figure 7C5), the single-fingering pattern (SF) is obtained (Figure 7D4). The transition from the kite shape pattern (KS) to the single-fingering pattern (SF) can be captured under the critical condition for the overlap event to occur in the last column: θ K S _ S F ( λ ) :   R o , N − 1 , N − 2 = R N − 1 / 2 , (4) where R _o,N-1,N-2 is the critical radius of curvature of the overlap event between the N-1th column and N-2th column, which is determined by substituting r 1 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 3 ) / 4 , r 2 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 2 ) / 4 , r 4 = a ( 1 − d i n l e t ) / 2 − 3 a λ ( N − 1 ) / 4 , and Φ = 120° into Supplementary Equation S4 for θ ≤ 90° and into Supplementary Equation S5 for θ > 90°.

Based on the theoretical model (Eqs 1–4), the phase diagram of displacement patterns is obtained. We evaluate this phase diagram in the next section.

To evaluate the theoretical model (Eqs 1–4), we use our previous pore network model [44], which can reasonably simulate the fluid invasion processes in the porous media [30, 44]. The pore network model adapts the geometries as shown in Figure 6. For the negative gradient, the pore size remains constant at the inlet, with d _inlet = 0.8, and decreases from the inlet to the outlet. For the positive gradient, the pore size also remains constant at the inlet, with d _inlet = 0.2, and increases along the flow direction. The simulations were performed at Ca = 10^–15, which can be regarded as the quasi-static condition, i.e., each time the pore adjacent to the meniscus with the largest critical radius of curvature gets filled.

We perform pore network simulations with 45°≤θ ≤ 120° and with -1.5 × 10^–2<λ < -5 × 10^–4 for negative gradient porous media or 5 × 10^–4<λ < 1.5 × 10^–2 for positive gradient porous media. In total, λ×θ = 14 × 12 = 168 cases are simulated. The simulated displacement patterns are shown in Figure 9 and Supplementary Figures S4, S5. Visually, for negative gradient porous media, TS occurs at low contact angles (θ ≤ 60°) and translates to CP as θ increases (Figure 9; Supplementary Figure S4), indicating that increasing θ destabilizes the fluid invasion. The critical θ that separates TS and CP increases with λ. For positive gradient porous media, CP and TS occur at low contact angles (θ < 70°), while KS and SF occur at high contact angles (θ ≥ 70°). Contrary to negative gradient porous media, increasing θ stabilizes the fluid invasion in positive gradient porous media.

We also evaluate the theoretical model (Eqs 1–4) with invading fluid saturation S _br and normalized fluid–fluid interface L _inter/L _max at the time of breakthrough, as shown in Figure 10. The dots represent the simulated displacement patterns in Figure 9 and Supplementary Figures S4, S5, with green circles for the compact pattern (CP), purple squares for the taper shape pattern (TS), yellow diamonds for the kite shape pattern (KS), and blue triangles for the single fingering pattern (SF). We also present the experimental results in the literature [30], which were performed at θ = 61.5° ± 5° in the positive gradient porous media (stars in Figures 10B,D). The displacement patterns are the compact pattern (CP) at λ = 5 × 10^–4 (gray star) and the taper shape pattern (TS) at λ = 5 × 10^–3 and λ = 1 × 10^–2 (white stars). From Figure 10, we observe that the compact pattern (CP) has the highest S _br (S _br>0.95) and the lowest L _inter/L _max (L _inter/L _max<0.25). In contrast, the single-fingering pattern (SF) has the lowest S _br (S _br≤0.05) and the highest L _inter/L _max (L _inter/L _max = 1). It is worth noting that our theoretical model is initially derived from the three same size post system (Figure 1). However, in an ordered porous medium, the radii of these three posts are different. This difference may lead to a decrease in the critical radius of curvature of the touch event along the flow direction in the negative gradient porous media at certain contact angles and gradients (θ = 65°, λ = 1.5 × 10^–2 and 1.2 × 10^–2 in Figures 10A,C), and eventually, lead to the compact pattern (CP). Except for these two cases, the theoretical model (Eqs 1–4) well-predicts the regimes of CP, TS, KS, and SF in the λ-θ plane.

Based on the contour plots in Figure 10, we further present the details of the wettability impact on S _br and L _inter/L _max at high gradient (|λ| = 1.25 × 10^–2, Figures 11A,B) and low gradient (|λ| = 5 × 10^–4, Figures 11C,D). The boundaries separating different displacement patterns predicted by the theoretical model (Eqs 1–4) are also presented with vertical dashed lines. Red, yellow, green, and blue represent the predicted displacement patterns of CP, TS, KS, and SF, respectively. For negative gradient porous media, increasing θ increases S _br and decreases L _inter/L _max at both gradients (Figures 11A,C), thus stabilizing the invasion front and causing the displacement patterns to translate from TS to CP at the critical contact angle, θ _{CP_TS}. For positive gradient porous media, on the contrary, a destabilizing effect occurs as θ increases (Figures 11B,D). At a high gradient (λ = 1.25 × 10^–2), increasing θ decreases S _br and increases L _inter/L _max, causing the displacement patterns to translate from TS to KS at θ _{TS_KS} and then to SF at θ _{KS_SF} (Figure 11B). At low gradient (λ = 1.25 × 10^–2), S _br and L _inter/L _max first, respectively, slightly increases and decreases with θ in the range of θ < 60° because the displacement patterns translate from TS to CP. Then, increasing θ destabilizes the invasion front, with displacement patterns translating from CP to TS at θ C P _ T S , ′ , and to KS at θ _{TS_KS}, and to SF at θ _{KS_SF}.

We also present the details of the gradient impact on S _br and L _inter/L _max at θ = 70° for negative gradient porous media (Figure 11E) and at θ = 55° and θ = 90° (Figures 11C,D) for positive gradient porous media. Again, the boundaries of λ _{CP_TS} and λ _{KS_SF} can capture the transitions from CP to TS (Figures 11E,F) and from SF to KS (Figure 11G). Here, λ _{CP_TS} are obtained by θ _{CP_TS}(λ) = 70° and θ C P _ T S ′ ( λ ) = 55°, and λ _{KS_SF} is obtained by θ _{KS_SF} (λ) = 90°. Generally, the theoretical model (Eqs 1–4) can capture the boundaries of displacement patterns.

We establish a link between pore-filling events and displacement patterns in ordered porous media and propose a predictive method to directly predict displacement patterns based on the pore structure. This method is originally derived based on the pore-filling event in the triangular lattice porous media, with a pore size gradient along the flow direction. For a different lattice structure, the pore-filling events will change. Whether the predictive method can be extended to other lattice structures should be discussed. In this section, we take the square-arranged porous media as an example to discuss how to extend this method in porous media with different lattice structures.

To begin with, we investigate the transition of pore-filling events in square-arranged porous media. Compared with the triangular arrangement, the square arrangement is a loose packing and has significant porosity for the same pore size. Since the touch event is more likely to occur in the geometry with dense packing, it is less likely to occur in the square-arranged porous media. We calculate the critical condition of the touch event in square-arranged porous media, as shown in SI Supplementary Text S3. The touch event can only occur when d < d _c, where the maximum value of d _c is d _cmax = 0.105. For the range considered in this work (45°<θ < 180° and 0.2 < d < 0.8), d < d _c cannot be satisfied, indicating that the touch event will not occur. On the other hand, as discussed in Section 3 in the main text, the overlap event with Φ = 180° will not affect the local invasion morphologies, so we only consider the burst event and the overlap event with Φ = 90° in square-arranged porous media.

Similarly, we consider a simple four-post system with the same radius. Figure 12 presents the transition of pore-filling events and their critical radius of curvature in square-arranged porous media. The boundary between the burst event and the overlap event is determined by substituting Φ = 90° into Supplementary Equations S8, S9. The critical radius of curvature of the burst event is determined by substituting r ₁ = r ₂ = a (1-d)/2 into Supplementary Equation S1. The critical radius of curvature of the overlap event is determined by substituting r ₁ = r ₂ = r ₃ = a (1-d)/2 and Φ = 90° into Supplementary Equation S4 for θ ≤ 90° and into Supplementary Equation S5 for θ > 90°. By comparing Figure 12 and Figure 1B, we find that the transition of pore-filling events is similar in both pore structures, with only a slight difference in the boundary curves and magnitude of the radius of curvature. In both pore structures, the critical radius of curvature of the burst event increases with the pore size, leading to the burst event being more likely to occur in the region with a larger pore size. The burst event leads to a compact front for the negative gradient porous media and a single fingering for the positive gradient porous media. On the other hand, the overlap event with Φ = 90° in the square-arranged porous media involves the merging of the neighboring menisci and its effect is related to the filling state of the invading phase. Similar to the overlap event with Φ = 120° in the triangle-arranged porous media, the overlap event with Φ = 90° leads to a tapered front when it is adjacent to the compact front and a widened fingering when it is adjacent to the single fingering. Therefore, the relation between the pore-filling events and the local invasion morphologies will not change with lattice structures. In a general porous media with a pore size gradient, as long as the pore structure is determined, the local invasion morphologies can be determined by the calculation of the pore-filling event.

Since the square arrangement is a loose packing, the touch event will not occur in the range considered in this work (45°<θ < 180° and 0.2 < d < 0.8). Therefore, the compact front will not occur in the square-arranged porous media with a positive gradient. Except for this difference, the effects of the burst event and overlap event on the local invasion morphologies are similar to those in the triangle-arranged porous media, that is, the burst event leads to a compact front in the negative gradient porous media and a single fingering in the positive gradient porous media, and the overlap event leads to a tapered front or a widened fingering. The displacement patterns can be predicted by combining the local invasion morphologies, as shown in Supplementary Figure S7. Compared to the triangle-arranged porous media (Figures 7B,D), the displacement patterns are also observed in the square-arranged porous media, except for the compact pattern (CP) under the positive gradient condition.

The transition of the displacement patterns can be determined using the critical condition of pore-filling events. In the negative gradient porous media, the burst event leads to a compact front and the overlap event leads to a tapered front. The transition from the compact pattern (CP) to the taper shape pattern (TS) corresponds to the critical condition for the burst event to occur in the last column: θ s , C P _ T S ( λ ) :   R s , b , N − 3 / 2 = R s , N − 1 / 2 , (5) where R _s,b,N-3/2 is the critical radius of curvature of the burst event in the N-3/2th column, which can be determined by substituting r 1 = a ( 1 − d i n l e t ) / 2 − 2 a λ ( N − 2 ) / 4 and r 2 = a ( 1 − d i n l e t ) / 2 − 2 a λ ( N − 3 ) / 4 into Supplementary Equation S1. R _s,N-1/2 is the critical radius of curvature for menisci on the tip of the invasion front in square-arranged porous media. Here, the derivation of the critical radius of curvature of menisci on the tip of the invasion front in square-arranged porous media is presented in Supplementary Text S4. The subscript s represents the square lattice.

In the positive gradient porous media, the burst event leads to a single fingering, and the overlap event leads to a tapered front or a widened fingering. The displacement patterns translate from the taper shape pattern (TS) to the kite shape pattern (KS) when the burst event just occurs in the first column: θ s , T S _ K S ( λ ) :   R s , o , 3 / 2 = R s , b , 3 / 2 , (6) where R _s,o,3/2 is the critical radius of curvature of the overlap event in the 3/2th column, which is determined by substituting r 1 = r 4 = a ( 1 − d i n l e t ) / 2 , r 2 = a ( 1 − d i n l e t ) / 2 − 2 a λ / 4 , and Φ = 90° into Supplementary Equation S4 for θ ≤ 90° and into Supplementary Equation S5 for θ > 90°; and R _s,b,3/2 is the critical radius of curvature of the burst event in the 3/2th column, which can be determined by substituting r 1 = a ( 1 − d i n l e t ) / 2 and r 2 = a ( 1 − d i n l e t ) / 2 − 2 a λ / 4 into Supplementary Equation S1.

Finally, the transition from the kite shape pattern (KS) to the single-fingering pattern (SF) can be captured under the critical condition that the overlap event just occurs in the last column: θ s , K S _ S F ( λ ) :   R s , o , N − 1 / 2 = R s , b , N − 1 / 2 , (7) where R _s,o,N-1,N-2 is the critical radius of curvature of the overlap event in the N-1/2th column, which is determined by substituting r 1 = r 4 = a ( 1 − d i n l e t ) / 2 − 2 a λ ( N − 1 ) / 4 , r 2 = a ( 1 − d i n l e t ) / 2 − 2 a λ ( N − 2 ) / 4 , and Φ = 90° into Supplementary Equation S4 for θ ≤ 90° and into Supplementary Equation S5 for θ > 90°; and R _s,b,N-1/2 is the critical radius of curvature of the burst event in the N-1/2th column, which can be determined by substituting r 1 = a ( 1 − d i n l e t ) / 2 − 2 a λ ( N − 1 ) / 4 and r 2 = a ( 1 − d i n l e t ) / 2 − − 2 a λ ( N − 2 ) / 4 into Supplementary Equation S1.

We investigate the displacement patterns in square-arranged porous media with a constant pore size at the inlet, and the pore size changes linearly along the flow direction. We consider d _inlet = 0.55 for negative gradient porous media and d _inlet = 0.1 for positive gradient porous media, which are consistent with the geometry in the experiments in Lu et al. [31]. They have performed many experiments at different gradients, and the pore structure (d _inlet) can only be determined in the two sets of experiments with the largest gradient. Here, we select these two experiments for evaluation [31]. The length of the porous media is L = 25 mm, and the maximum number of columns is N = 30. The phase diagram of the displacement patterns can be determined using Eqs 5–7 and is shown in Figure 13. The phase diagram in the square-arranged porous media (Figure 13) is similar to that in the triangle-arranged porous media (Figure 8), except that the compact pattern (CP) does not occur under the positive gradient condition and the boundary curves of displacement patterns are slightly different because the lattice structure only affects the occurrence of pore-filling events and the magnitude of the radius of curvature, but not the relation between the pore-filling events and the displacement patterns.

We evaluate this phase diagram using the experiments in [31]. The contact angle is determined to be θ = 118.5 ± 5°, according to the fluids and the material of the porous media (SI Supplementary Text S5) [52–55]. The compact pattern (CP) is observed at λ = -2.04 × 10^–2 (solid circle in Figure 13A), and the single-fingering pattern (SF) is observed at λ = 2.04 × 10^–2 (open circle in Figure 13B) [31], in agreement with the theoretical model. It is worth noting that the experiments are performed with a controlled amount of disorder. The difference between the maximum and minimum pore size in a column is d _max–d _min = 5.5 × 10^–3, much smaller than the gradient (|λ| = 2.04 × 10^–2), indicating that the predictive method is applicable when the effect of pore size gradient can suppress the effect of disorder. We also evaluate the predictive method using numerical simulations in the literature [32], proving that our predictive method can reasonably predict the displacement patterns when the effect of viscous force is suppressed by the pore size gradient (Supplementary Text S6).

The theoretical model also shows that decreasing contact angle in square-arranged porous media leads to the transition from the compact pattern (CP) to the taper shape (TS) in the negative gradient porous media and from the single-fingering pattern (SF) to the kite shape pattern (KS) to the taper shape (TS) in the positive gradient (Figure 13). Since the existing experiments and numerical simulations are performed under drainage conditions [31, 32], the transition of displacement patterns with wettability has not been observed in these studies and needs to be further evaluated.

In general, the predictive method can be extended to different lattice structures. The lattice structures only affect the occurrence of pore-filling events and the magnitude of the curvature of radius, but not the relation between the pore-filling events and displacement patterns. For a generally arranged porous medium with a pore size gradient along the flow direction, the pore-filling events can be determined through the pore structure, and then the displacement patterns can be predicted as long as the pore size gradient can suppress the effect of disorder. This method directly predicts displacement patterns through pore structures and significantly improves the efficiency of predicting displacement patterns, without performing experiments and simulations.

We propose a method to directly predict the displacement patterns in ordered porous media with different pore size gradients and contact angles. The pore size in the ordered porous media remains constant in the same column and changes monotonously along the flow direction, leading to the transition of the pore-filling events and eventually affecting the displacement patterns. We first link the pore-filling events to the local invasion morphologies, providing a basic understanding of how pore-filling events impact the local invasion process in the ordered porous media. In contrast to the well-accepted picture that the burst and touch events always play a role in destabilizing the fluid invasion in disordered porous media, we observe that they can lead to a compact front in negative and positive gradient porous media, respectively. Then, the theoretical models that predict displacement patterns are developed based on the combination of local invasion morphologies. We find that decreasing the contact angle in the range of 45°<θ < 120°, which is regarded as smoothing the invasion front, destabilizes the fluid invasion in the negative gradient porous media. Finally, the theoretical models are evaluated using pore network simulations and experiments in the literature. The phase diagram of displacement patterns predicted by the theoretical models is in good agreement with simulations and experiments in the literature.

This work is an extension of our previous work [30] and aims to present the boundaries of the different displacement patterns. In summary, there are three major differences between these two studies.

1) Our previous work can only predict the displacement efficiency and not displacement patterns. In this work, we investigate the pattern transition in ordered porous media and propose a predictive model to directly predict displacement patterns.

This work improves the understanding of how pore size gradient and contact angles impact the invasion process and is of practical significance for controlling fluid invasion in ordered porous media, such as multi-layered soils.