{'fraction_non_alphanumeric': 0.04280296500603344, 'fraction_numerical': 0.011308395104292364, 'mean_word_length': 4.003967911670836, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Capillary condensation in one-dimensional irregular confinement

A lattice-gas model with heterogeneity is developed for the description of fluid condensation in finite sized one-dimensional pores of arbitrary shape. An exact solution of the model is presented for zero-temperature that reproduces the experimentally observed dependence of the amount of fluid adsorbed in the pore on external pressure. Finite-temperature Metropolis dynamics simulations support analytical findings in the limit of low temperatures. The proposed framework is viewed as a fundamental building block of the theory of capillary condensation necessary for reliable structural analysis of complex porous media from adsorption-desorption data.


I. INTRODUCTION
Physical systems which consist of networks of pores, such as Vycor 1 , Silica aerogels 2 , porous rocks 3 , soil 4 and others, have a wide spectrum of applications, ranging from molecular filters and catalysts 5 to fuel storage 6 . Capillary condensation is an important and peculiar physical phenomenon occurring in many such systems. Depending upon the pore structure and the materials involved, the adsorbed fluid density can exhibit both hysteresis and avalanches (abrupt changes in density) as the pressure varies. In recent years, a lot of experimental and theoretical work has been undertaken in order to understand this dependence 2,7,8 . However, the understanding is not fully complete 8 and, for example, routinely used classical theories fail to explain the appearance and shape of hysteresis in sorption curves of capillary condensation even in simple one-dimensional (1D) isolated pores [9][10][11] . Indeed, according to classical theories, pores closed at one end should not have hysteresis. In contrast, the adsorption of N 2 into MCM-41 mesoporous silica consisting of pores closed at one end at the temperature T = 78K reveals a hysteresis loop of the so-called H2-type 12 , i.e. involving a smooth increase in density for adsorption and sharp drop in density for desorption. An accurate theory for capillary condensation in 1D pores is ultimately necessary to ensure that, for instance, this phenomenon can be used as a reliable technique to probe the structure of generic porous media consisting of a network of interlinked 1D pores 2,13 .
A number of techniques have been developed to study capillary condensation theoretically, including microscopic molecular dynamics 14 , density functional theory 15,16 and lattice-gas mean-field theory 1,17,18 . Such studies have been conducted on a variety of different porous media. For one-dimensional pores, mean-field theory 19 and multi-scale molecular dynamics studies 20,21 have previously been used to test the hypothesis that hysteresis is caused by heterogeneity 11 , which includes variations in pore diameter, chemical heterogeneity in the pore walls and roughness of internal surfaces 22-25 . Additionally, these numerical studies have revealed the occurrence of avalanches in the amount of adsorbed fluid during condensation and evaporation. Such avalanches bare similarity to avalanches in magnetisation found in the random-field Ising model (RFIM) 1,17,26,27 . Employing this similarity, we map RFIM to the lattice-gas model and demonstrate that: (i) a heterogeneous lattice-gas model is a minimally sufficient model to reproduce experimental observations of variations of fluid density with pressure in finite-sized pores; (ii) this model can be solved exactly analytically at zero temperature (T = 0) by a novel technique, with the solution being fully supported by numerical simulations; (iii) such an analytical solution leads to simple physical explanations and interpretations of experimental results for condensation in 1D pores that remain qualitatively valid at sufficiently low finite temperatures.
More specifically, our findings are as follows. (i) The effects of the closed and open ends of the pores, considered important in the classical theories, are significantly reduced for large disorder in strengths of the interactions of the fluid with the pore walls. (ii) A positive-skewed distribution in strengths of such interactions can lead to sorption curves of the same form as those found experimentally 10,11 . (iii) It is predicted that the mechanism for adsorption depends crucially on the length of pores and their geometry. For short pores, adsorption isotherms may depend on whether pores are open or close at the ends. In contrast, the dependence of adsorption on the characteristics at the pore ends is lost for long enough pores. (iv) Sorption is shown to have several different mechanisms, leading to different forms of sorption curves. However, sorption in long pores, or pores characterised by a large degree of heterogeneity, shows universal features typical of a disorder-controlled regime. (v) In cylindrical pores consisting of two parts of different diameter, condensation and evaporation in one part can induce condensation and evaporation in the rest of the pore for low disorder, but for high disorder, the two parts of the pore behave independently.
The structure of this paper is as follows. The model for condensation and pore geometries studied are intro-duced in Sec. II. Sec. III presents the exact solution of the proposed model at zero temperature. The results for adsorption/desorption are presented in sections IV and V for T = 0 and T > 0, respectively. Finally, the conclusions are given in Sec. VI.

II. MODEL AND PORE GEOMETRIES
A. Lattice-gas model with disorder in matrix-fluid interactions The proposed model is based on a standard lattice-gas model of capillary condensation 1,2,9,17 . In this model, the 3D space is split into cells, and each cell can be filled either by matrix (the solid substrate), liquid or vapour. If the cell i is occupied by matrix (expressed by setting parameter η i = 0), then it cannot become occupied by fluid. It is assumed that the variables η i are quenched for the whole system and thus the matrix state cannot change during the condensation. If cell i is not occupied by the matrix (η i = 1), then it can be occupied by either liquid (τ i = 1) or vapour (τ i = 0). The variables τ i can vary with the change in chemical potential, µ. The Hamiltonian which describes the lattice gas model is given by 1,17,18 , where the summations run over all the cells in the system in the first term and over all nearest-neighbour pairs ij in the other terms. The fluid-fluid interaction parameter, w ff > 0, is assumed to be the same for all pairs of cells. The matrix-fluid interaction strength between the matrix at cell j and fluid at the neighbouring cell i is described by the parameter w mf ij . The values of w mf ij are considered to be independently distributed quenched random variables with the probability density function, ρ i (w mf ij ), which can be cell dependent. The random distribution of this parameter has been studied previously in the context of chemical heterogeneity of the pore walls 28 , but below it is assumed to characterise all types of heterogeneity. For concreteness, we focus on two forms of disorder in w mf ij , representing heterogeneity on different scales ranging from local fluctuations at a single point on the pore wall to a variable diameter of the pore. More specifically, we consider (i) a normal distribution, and (ii) an exponential distribution with correlations ensuring that all w mf ij are the same for the same cell i (i.e. w mf ij is independent of j), and distributed according to, where x = (w mf ij − w mf i ) + ∆ i and Θ(x) is the Heaviside step function. In Eqs. (2) and (3), w mf i and ∆ 2 i refer to the mean and variance, respectively, of the corresponding distributions. Gaussian uncorrelated heterogeneity describes local fluctuations of the matrix-fluid interaction. In contrast, correlations in the exponential heterogeneity are introduced to describe heterogeneity in pore diameter.
B. Mapping to the random-field Ising model The lattice gas model can be mapped to the RFIM 1,17 which gives an intuitively simpler representation and is technically more convenient for the calculations presented below. In the RFIM representation, the Hamiltonian described by Eq. (1) is given by, where the variable s i = (2τ i − 1)η i (s i = ±1) represents a spin state, H = µ/2 refers to an external magnetic field, and J = w ff /4 describes the spin-spin interaction. The fields h i at cell i are given by with the sum running over all nearest neighbours of i.
Since w mf ij are randomly distributed, the fields h i defined by (5) are also random variables, with distribution ρ h (h i ) which depends on the pore geometry chosen. More explicitly, ρ h (h i ) depends on the numbers, n m i = j/i (1 − η j ) and n u i = j/i η j , of neighbouring cells which are occupied and unoccupied by matrix, respectively. When w mf ij is normally distributed, the random field at cell i is distributed according to ρ hi (h i ) = N ( h i i , ∆ hi 2 ) with h i i = (n m i /2) w mf i + (n u i /4)w ff and ∆ 2 hi = (n m i /2)∆ 2 i . When w mf ij follows a correlated exponential distribution, the random field at cell i is dis- , with the mean h i i being the same as for the normal distribution and standard deviation ∆ hi = n m i ∆ i /2. In both cases, the values of h i i and ∆ hi 2 can differ between the cells.

C. Pore geometries
Three types of 1D pore geometries consisting of N cells embedded in a simple cubic lattice are analysed below. They are linear pores of (i) type I with both ends closed by matrix (see Fig. 1(a)), (ii) type II with both ends open and bounded by vapour (see Fig. 1(b)), and (iii) type III with one open end and one closed end. In the first two cases, we consider w mf to be identically and independently distributed random variables with w mf i = w mf (a) Type I (b) Type II (c) Type III and ∆ i = ∆. In the last case, we consider two distinct sections of lengths N 1 and N − N 1 characterised by different values of mean matrix-fluid interaction w mf 1 and w mf 2 , but the same variance ∆, respectively (see Fig. 1(c)). Pores of type III are intended to model inkbottle and funnel pore geometries, with weaker matrixfluid interaction representing a larger diameter.

D. Dynamics
We study adsorption and desorption isotherms obtained by sweeping the chemical potential at a rate r from µ = −∞ to ∞ and back again. The state of the system is described in terms of the mean volume of the absorbed liquid, V = N −1 i τ i and the variance of this quantity, When the system is driven in this way, it evolves through a rugged free energy landscape corresponding to the Hamiltonian H which, due to the presence of random fields, consists of an exponentially large number of local minima for given µ and V 1,29,30 . The evolution of the system (i.e. changes of the occupation numbers, {τ i }) is caused either by changes in the energy landscape due to variations of µ, or by thermal effects. Below the critical temperature, T c , where condensation occurs in a discontinuous manner 2 , random fields introduce glass-like behaviour to the system and thermally activated transitions are unlikely on the time-scale of real experiments 1 . That is why the results of mean-field lattice-gas models 1,18 , which ignore thermal fluctuations, and experiments 2 are in good agreement. Temperature is known to affect hysteresis 2 but this has been properly accounted for by mean-field theories in terms of entropic contributions that represent a quenched modification to the energy landscape 1,18 . Motivated by these results, we first obtain the exact analytical solution to the proposed model at zero temperature and then study the effects of non-zero temperature by means of Monte-Carlo numerical simulations. These simulations include thermally activated events but their effect is expected to play a secondary role on hysteresis at low temperatures and are not analysed in detail.
At T = 0, we employ a single-spin flip metastable dynamics that has been widely used within the context of the zero-temperature RFIM (zt-RFIM) 26 . According to this dynamics, for adsorption (desorption) each cell is initially empty (occupied), s i = −1 (s i = +1), and can become occupied by liquid (gas), s i = +1 (s i = −1), once its local field, is positive, f i > 0 (negative, f i < 0). Here, n f i = j/i η j τ j is the number of fluid cells neighbouring cell i. According to the above rule, a configuration of phases is stable if all the fluid cells satisfy the condition s i f i > 0. The system is driven quasistatically under the assumption that the rate of relaxation of spins is much larger than the rate r of variation of µ 26,27 . In practise, this is achieved by sweeping µ until at least one spin becomes unstable. At this point, an avalanche starts and µ is kept constant until a new stable configuration is reached. After that, activity can only resume if µ is varied and a new avalanche is induced by this variation.
Adsorption/desorption processes at T > 0 are simulated numerically using Metropolis dynamics implemented as follows 31 . For given µ, a cell is chosen at random and a change of its state, τ i (or s i ), is proposed. The change is accepted with probability p = min{1, exp(−β∆H)}, where β = (k B T ) −1 and ∆H = f i s i is the change in energy of the proposed change of state. This process is referred to as a single Monte Carlo step, with the unit of time in the simulation being Monte-Carlo Steps per Spin (MCSS). As the simulation progresses the value of µ is incremented at a fixed rate r measured in MCSS −1 so that the increment at each step is δµ = r/N . This dynamics reduces to the T = 0 dynamics in the limit of r → 0 and β → ∞.

III. EXACT SOLUTION OF THE MODEL AT
In this section, we obtain exact analytical solutions for sorption curves for a 1D pore. In order to do that, we develop a novel method based on the use of the generating function formalism 32,33 . Within this formalism, the volume of fluid in the finite pore is a random value and can be expressed by a generating function, where P (n) is the probability that n = N i=1 τ i cells in the finite pore are occupied by fluid. The mean volume of fluid is given by, with volume of a cell set to V 0 = 1, and the variance is, where ∂ x G(1) and ∂ xx G(1) refer to the first and second derivatives of G(x) with respect to x evaluated at x = 1.
We derive a form for these expression by establishing a recursion relation for G(x) in the following way. We assume particular sequences for occupancy of cells that are convenient from a mathematical viewpoint. This is possible owing to the abelian property of the zt-RFIM which ensures that the final state of the system is independent of the order in which the cells are occupied by fluid (the order spins flip) as long as the system ends up in a stable or metastable state 26,34 . The sequence we chose assumes that the relaxation of the system into a metastable state takes place in a series of N time-steps, t, 1 ≤ t ≤ N . Initially, during time-step t = 1, all the cells with i ≥ 2 are artificially prevented from being occupied, while the first cell i = 1 is allowed to change its state. Cell i = 1 can either change state from unoccupied to occupied if the local field given by Eq. (6) is positive, f 1 > 0, or remain unoccupied if the local field is negative. This process with two possible outcomes is called relaxation of cell 1. Next, during time-step t = 2, we allow cell i = 2 to relax, while cells i ≥ 3 are still held in the unoccupied state, and cell i = 2 can become occupied if f 2 > 0. If cell i = 2 does become occupied, then the local field f 1 at cell i = 1 will increase. This can cause cell i = 1 to become occupied if it was not occupied already, i.e. an avalanche can pass from cell i = 2 to cell i = 1 during time-step t = 2. Similarly, we allow the next cell in the pore, i.e. i = 3, to relax and if it becomes occupied an avalanche can pass back along the pore towards cell i = 1 if those cells with i < 3 were not occupied. This method is recursively applied N times until all the cells in the system are relaxed.
Let us consider the cell i = N at the end of the pore, which is the last cell to be allowed to relax in the above procedure. At the start of time-step t = N , the neighbouring pore cell i = N − 1 can be occupied or unoccupied, i.e. s N −1 = ±1. If the neighbouring cell i = N − 1 is occupied (s N −1 = +1), and the random field at cell i = N is h N > −J − H, then the local field f N > 0 and cell i = N will become occupied. This occurs with cell-dependent probability p ′ N,1 where, with m being the number of occupied neighbours of cell i = N , i.e. m = 1 in this case. If, however, the neighbouring cell i = N − 1 is unoccupied (s N −1 = −1), then the local random field must be above a higher threshold for cell N to become occupied, h N > J −H, which occurs with probability, p ′ N,0 , given by Eq. (10) with m = 0. In this case, the field at cell i = N − 1 will increase, and an avalanche can propagate back along the pore.
The probability that there are n occupied cells in the lattice at the end of the relaxation process can therefore be written as, in terms of the probabilities, P A (i, n ′ ), P B (i, n ′ ) and P C (i, n ′ ), which can be recursively determined. The quantity, is the probability that at the end of time is the probability that cell i is unoccupied at end of timestep t = i, s i (t = i) = −1, and at the end of the relaxation process (end of time-step t = N ) there are n i = n ′ occupied cells in the range 1 ≤ j ≤ i, given that cell i + 1 becomes occupied, s i+1 (t = N ) = +1, during some timestep t ′ , i < t ′ ≤ N (which causes avalanches to pass back along the pore towards cell 1, changing the occupation number n i ). The third quantity, is the probability that at the end of time-step t = i there are n i (t = i) = n ′ occupied cells in the range 1 ≤ j ≤ i and that cell i is unoccupied at this time-step, s i (t = i) = −1.
Using Eqs. (7) and (11), the generating function G(x) for the total number of occupied cells in the pore can be written in terms of the generating functions A i (x), B i (x) and C i (x) for the corresponding probabilities defined by Eqs. (12)- (14) as, and C i (x) are defined according to the same relation with A replaced by B and C, respectively.
Expressions for the generating functions A i (x), B i (x) and C i (x) for i > 1 can be found recursively in the following way. If cell i − 1 is occupied at the start of timestep t = i and the local field at cell i is positive, i.e. f i = h i + H > 0, then cell i will become occupied during time-step t = i. This occurs with probability p i,1 , where If cell i − 1 is unoccupied then the random field has to be at a higher threshold, h i > −H + 2J, in order for cell i to become occupied during time-step i. The random field will be above this higher threshold with probability p i,0 , given by Eq. (16). If cell i does become occupied during time-step t = i and s i−1 = −1 at this time, then an avalanche might pass back along the pore towards cell i = 1. This gives the expressions for P A (i, n ′ ) and P C (i, n ′ ) defined in Eqs. (12) and (14), If cell i does not become occupied at time step t = i, but cell i + 1 becomes occupied at some later time-step t ′ (i < t ′ ≤ N ) then cell i can also become occupied during time-step t ′ , i.e. an avalanche can propagate back down the pore. The probability of cell i becoming occupied in this way at time step t ′ depends upon whether cell i − 1 is occupied or not. In fact, if cell i − 1 is occupied and cell i + 1 becomes occupied on time step t ′ then cell i will also become occupied at time step t ′ only if the random field at cell i is in the range −H < h i < −H + 2J, which occurs with probability p i,2 − p i,1 . On the other hand, if cell i − 1 is unoccupied then the random field must be in the range −H − 2J < h i < −H in order for cell i to become occupied at time step t ′ when cell i + 1 becomes occupied, which occurs with probability p i,1 − p i,0 . In the case that cell i − 1 is unoccupied, an avalanche of spin flips can propagate back down the pore from cell i towards cell 1 during time step t ′ . This gives the expression for P B (i, n ′ ) defined in Eq. (13), Eqs. (17) and (18) lead to the following recursive relations for the generating functions, The boundary values of P A (1, n ′ ) and P C (1, n ′ ) can be found using the following relations, where p ′ 1,0 is the probability that cell 1 has a positive field (and thus becomes occupied) at the first time step (when all other cells are unoccupied). The value of P B (1, n ′ ) is given by the relation, is the probability that cell 1 has a negative local field during time-step t = 1, but the field becomes positive when cell 2 becomes occupied, and 1−p ′ 1,1 is the probability that cell 1 still has a negative local field after cell 2 becomes occupied. Eqs. (20) and (21) result in the following expression for the boundary generating functions, The generating function given by Eqs. (15), (19) and (22) can be written as a matrix equation, Eq. (23) is the main analytical result of our analysis allowing exact evaluation of ∂ x G(1) and ∂ xx G(1). The mean and variance of the volume of fluid in the pore can be found for both adsorption and desorption regimes using Eqs. (8) and (9) along with the derivatives of G(x). Technically, the derivatives of G(x) can be calculated by numerical iteration, i.e. the derivatives of The results presented in this section correspond to the numerical solution for G(x) given by Eq. (23). In order to test the validity of the exact solution, we performed Monte-Carlo simulations of condensation within the framework of the lattice-gas model. During each simulation, the value of µ was fixed, and hysteresis plots were obtained by running separate simulations for each value of µ. The fluid occupation number of the matrix-free cells was changed following single spin-flip zero-temperature Metropolis dynamics 26,35 . The mean and variance of V were calculated by averaging over 10 4 realisations of the disorder in random-fields for a fixed matrix structure and values of the parameters. Figs. 2, 3, 7 and 8 show that the analytical calculations (lines) agree with the numerical simulations (symbols).
The shape of the hysteresis loop and the distribution of ρ for a given value of µ depends on the pore type and are governed by several parameters. These parameters are the degree of disorder in matrix-fluid interactions, ∆, the system size, N , the relative strengths of matrix-fluid, w mf , and fluid-fluid interactions, w ff , and the value of the chemical potential, µ. The values of the parameters influence the dynamics of the avalanche, i.e. they control the location of the avalanche nucleation points, being either at end points or internal points of the pore, and the number of points at which the avalanches are pinned. As a result, there are several distinct regimes for sorption in linear pores of all types.
We start out analysis with a detailed description of the different regimes in pores of types I and II. First, the different possible forms for the sorption curves in the case of a normal distribution in w mf ij are described. Second, we analyse the new features which an exponential distribution brings to the sorption curves. Finally, a description of the sorption curves in pores of type III is presented.

A. Nucleation and Pinning
When the increasing value of µ reaches a certain value one of the cells in the system can be filled with liquid, i.e. an avalanche can be nucleated. This avalanche can propagate either through the whole system or stop at some cell, which becomes a pinning point. When the avalanche is pinned, a further increase in µ is required for de-pinning, i.e. for adsorption to continue. Therefore, there exist two typical values of µ, corresponding to nucleation and de-pinning. Let us first estimate the value of µ at which nucleation occurs in pores of both types I and II.
Nucleation of an avalanche occurs at a particular cell i if it is an energetically favourable process, i.e. when the local field given by Eq. (6) at that cell is positive. This happens if the value of chemical potential becomes greater than µ i , a local nucleation potential, given by, where n m i and n f i are defined in Secs. II B and II D, respectively. In pores of type I, there are two kinds of cell, end cells, (i = 1 and i = N ) and interior cells (1 < i < N ). For end cells, number of surrounding matrix cells (n m 1 = n m N = 5) is greater than for inner cells (n m i = 4) favouring nucleation at the end cells of type-I pores. In contrast, all cells in pores of type-II are equivalent for nucleation events because n m i = 4 for 1 ≤ i ≤ N . The values of µ i are random as a consequence of the disorder in matrix-fluid interaction strengths w mf ij . First, we analyse the case of a normal distribution of w mf ij given by Eq. (2). For a quenched configuration of matrixfluid interaction there will be a cell with a minimal value of µ i , and the first nucleation event will occur at that cell. In a pore of type II, this minimal value has a mean, µ o min , which can be estimated using the relation Prob[µ i ≤ µ o min ] ≃ 1/N , for the mean minimum of N independently and identically distributed random values of µ i (see e.g. Ref. 36). For a normal distribution of w mf (see Eq. (2)), the estimate is given by the following equation, where n m with n m 1 = n m N = 5 (see the black lines in Fig. 4(a) and (b) for ∆ ∆ * ). For inner cells, the expression for the mean value of the minimum, µ inner min , is given by Eq. (26) with N replaced by N − 2, i.e. for large values of N , If µ end min < µ inner min which is true for ∆ < , nucleation starts typically at the end points, and the mean of the minimal value of µ i in type-I pores, µ c min , coincides with µ end min . Otherwise, if ∆ > ∆ * , nucleation starts at any cell in the pore with approximately equal probability and µ c min ≃ µ inner min < µ end min (see the coinciding dashed and solid lines Fig. 4(a) and (b) for ∆ ∆ * ).
Once an avalanche is nucleated, it starts propagating. This propagation can be stopped by unfavourable variations in the matrix-fluid interaction strength, i.e. it can be pinned at a certain cell. Pinning will occur at an interior cell i if it is not energetically favourable for cell i to become occupied, even when the avalanche causes one of the cells neighbouring cell i to become occupied by liquid, i.e. pinning occurs at cell i when the random value µ i (4, 1) > µ. Typically, pinning will occur when µ < µ pin max where µ pin max is the mean maximal value of µ i for 1 < i < N , which can be found using arguments similar to those for Eq. (26) as, with n m i = 4 and n f i = 1. The value of µ pin max does not depend on the type of the pore and increases with disorder (see the dot-dashed lines in Fig. 4(a) and (b)).

B. Regimes for Avalanches
First, we analyse the different adsorption regimes of adsorption for pores of type I. The relative values of µ pin max , µ c min and µ inner min (with µ c min ≤ µ inner min ) depend on w ff , ∆ and N , and four different regimes of adsorption exist. The first regime is defined by the following sequence of characteristic chemical potentials, An example of adsorption in this regime is shown by the vertical dashed line marked with (1) in Fig. 4(a), which first crosses the solid line, then the dot-dashed  (26)) occurs at ∆ < ∆ * , and for a certain range of ∆, e.g. at ∆ = 0.05 in Fig. 4(a). The adsorption in this regime (see solid line in Fig. 2(a) upper panel) begins for µ ≃ µ c min < µ inner min , implying that nucleation typically occurs at the end of the pore. The variance takes a relatively small value as compared to its maximum possible value of Var[V ] max = 0.25 corresponding to a bi-modal distribution of V with equally probable values V = 0 and V = 1 37 (see solid line in Fig. 2(a) lower panel), meaning that the avalanches nucleated at the ends of the pore progress gradually along the pore between several pinning points with no large sudden jumps in density. This is in contrast to cases when the varience takes its maximum value, which correspond to a single large avalanche (see regimes (3) and (4) below).
In the regime of large disorder in w mf ij (regime (2)), pinning dominates over nucleation, and adsorption takes place in the form of many small avalanches nucleated mainly in the inner cells. An example of adsorption in regime (2) is shown in Fig. 4(a) by a vertical dashed line, with the following sequence of crossing points, µ c min < µ inner min < µ pin max .
In this example, adsorption can be nucleated first at the closed ends, when µ ≃ µ c min . The avalanche that occurs as a result of such a nucleation is small because of the large number of pinning points which exist for µ < µ pin max . At a higher µ µ inner , many more small avalanches are nucleated in the middle of the pore, in contrast to regime (1). The resulting filling process is gradual, occurring mainly over the range of µ, µ inner µ µ pin max (see solid line in upper panel of Fig. 2 then the dominant effect for adsorption is pinning, meaning that adsorption is in regime (2). Such a regime thus occurs for large enough ∆ or N for any value of w ff , i.e. when the value of ∆ is to the right of the crossing point between the dashed and dot-dashed lines in Fig. 4(a) or (b). The third and fourth regimes (see Fig. 4(b)), can be achieved by increasing the value of w ff , corresponding to a downward shift of the dot-dashed line so that the intersection point of the dot-dashed and dashed line occurs at ∆ > ∆ * . In regimes (3) and (4), pinning is not important, i.e. µ pin max is the smallest relevant value of µ, and adsorption occurs in a single avalanche. This is demonstrated by the large peak in Var[V ] (≃ 0.25) for both regimes, see right-hand solid peak-shaped curves in the lower panels of Fig. 3(a) and (b) corresponding to regimes (3) and (4), respectively. The boundary between these two regimes occurs at ∆ = ∆ * . For regime (3), ∆ < ∆ * meaning that, In this case, adsorption occurs in a single avalanche nucleated at one of the end-cells at µ ≃ µ c min (see right-hand solid curve in Fig. 3(a) upper panel). For regime (4), ∆ > ∆ * and, so that adsorption occurs in a single avalanche nucleated at one of the inner cells at µ ≃ µ inner min (see right-hand solid curve in Fig. 3(b) upper panel).
The length, N , of the pore influences the configuration of the boundaries shown in Fig. 4 by changing the slope of the dashed and dot-dashed lines. The solid line is independent of N for ∆ < ∆ * . When N increases, the magnitude of both slopes increases proportionally to √ ln N . As such, for very large N , only regime (2) can be observed, with nucleation first occurring at very small values of µ and adsorption taking place as a series of small avalanches until µ is very large.
The four regimes for adsorption in type-I pores described above can be accessed at constant values of ∆ and N by varying the values of w mf and w ff . A diagram showing the boundaries between the regimes (1)-(4) in the parameter space of w mf and w ff is shown in Fig. 5(a). As seen from this diagram, all of these boundaries meet at a point A w mf , w ff located at, For values of w ff lower than its value at point A the adsorption is in regime (2), i.e. for such values of w ff there will be many small avalanches nucleated at inner cells. For larger values of w ff the regime depends on the value of w mf . Regimes (1) and (3), corresponding to avalanches nucleated only at the end-cells of the pore, appear at higher values of both w ff and w mf than the point A (the region bounded by the dashed line in Fig. 5(a)). Regimes (3) and (4), corresponding to avalanches that are not affected by pinning and occur in a single jump, appear at values of w ff which are higher than both the value at point A, and the value on a line w ff = w mf + ∆ √ 8 erfc −1 (2/N ) + 5/π which passes through A with gradient 1 (the shaded region in Fig. 5(a)). To summarise, regime (1) . It can be shown that only two regimes exist, (2) and (4). The boundary between these two regimes in the parameter space ( w mf , w ff ) is a line of constant w ff passing through the same point A as in Fig. 5(a). This means that pores of both types are in regime (2) for the same range of parameters. Conversely, when type-I pores are in regimes (1), (3) or (4), type-II pores are in regime (4).
The shapes of the adsorption curves in different regimes reveal the similarities and differences between the pores of types I and II. For the same set of parameters, adsorption can be in regime (1) for a type-I pore while it is in regime (4) for a type-II pore. In this case, adsorption in the type-I pore is nucleated at the endcells and progresses slowly between many pinning points along the pore while adsorption in the type-II pore is nucleated in an inner cell and happens in a single avalanche because there are no pinning points. The difference between the two pore types is due to the additional inter-  34) and (35), respectively. Adsorption in pores of type II exhibits two regimes, (2) and (4). Regime (2) corresponds to the region (2) in panel (a) and regime (4) spans over the regions (1), (3) and (4) in (a). Similarly, desorption in pores of type II occurs only in two regimes, (1) and (2). Regime (2) corresponds to the region (2) in panel (b) and regime (1) spans over the regions (1), (3) and (4)  actions between the fluid and the end-cap matrix cell of the type-I pores, which are absent in pores of type II and ensure that the value of µ at which adsorption starts in the type-I pores is lower than in the type-II pores ( µ c min ≃ −5.1 µ o min ≃ −4.2, compare solid and dashed lines in Fig. 2(a) upper panel). At the low values of µ when adsorption starts in the pore of type-I, there are many pinning points, which prevent the occurrence of a single avalanche. This is supported by the fact that the maximum value of the variance, Var[V ], is smaller than 0.25 for the type-I pore, but is approximately equal to 0.25 for a type-II pore (see the lower panel of Fig. 2(a)).
Adsorption in regime (2) is governed by disorder for both types of pore and their adsorption curves practically coincide (see Fig. 2(b) both panels). The relatively small values of variance (see Fig. 2(b) lower panels) confirm the presence of many pinning points leading to many small avalanches.
When the adsorption is in regime (3) for the pores of type I, it is in regime (4) for the pores of type II. This means that both pores exhibit a large avalanche, but it is nucleated at the end-cells in a type-I pore and at an inner cell in a type-II pore. Adsorption in the type-I pore therefore occurs at a lower value of µ than in a pore of type II (see Fig. 3(a) upper panel). The variance reaches a peak of 0.25 for both (see Fig. 3(a) lower panel), confirming that adsorption occurs in a single avalanche. At lower values of w mf , adsorption is in regime (4) for both types of pore, so that it is nucleated in an inner cell in both cases and occurs in a single avalanche. This causes the adsorption curves to practically coincide for both pore types (see Fig. 3(b)).
For desorption, a similar analysis can be undertaken, and the results of this analysis are presented by the set of lines to the left of the adsorption curves in Figs. 2 and 3. Four regions, similar to those for adsorption, exist for desorption in pores of type I (see Fig. 5(b)). These are separated by several boundaries, which cross at the point D( w mf , w ff ) located at, For values of w ff lower than that at point D (regime (2)), desorption exhibits many small avalanches nucleated in the inner cells, as in regime (2) for adsorption. For w ff higher than that at the point D, there are several regimes. Regimes (3) and (4) correspond to desorption taking place in a single avalanche, which occurs when both w ff and w mf are greater than their values at point D (shaded region in Fig. 5(b)). Regimes (1) and (3) correspond to desorption which is nucleated at the end of the cell, and occur for values of w ff greater than that at D and greater than a line w ff = w mf + ∆ √ 8erfc −1 (2/N ) − 5/π which passes through point D with gradient 1 (the region bounded by the dashed line in Fig. 5(b)). Note that the boundaries between the regimes are, in general, different from those for adsorption. Larger values of w mf encourage (discourage) nucleation of adsorption (desorption) at the end-cells, and also encourage (discourage) pinning to occur by making the adsorption (desorption) process start at a lower value of µ.
For pores of type-II, the desorption can be either in regime (1) or (2) only. These two regimes are separated in parameter space ( w mf , w ff ) by the line at constant w ff passing through points A and D. As such, when the parameters are chosen in such a way that desorption in type-II pores is in regime (2), a type-I pore with the same parameters will exhibit desorption in regime (2) also. In this case, the desorption curves and variances coincide for the two pore types (see Fig. 2(b)). For parameters such that desorption in a pore of type II is in regime (1), a type-I pore can be in either types (1), (3) or (4) and the desorption curves, in general, do not coincide.

C. Exponential disorder in matrix-fluid interaction strength
The above analysis has been done for the case of a normal distribution in w mf ij . The overall picture is qualitatively the same for a correlated exponential distribution of w mf ij given by Eq. (3). However, the presence of a well-defined lower bound for w mf ij associated with the sharp cut-off in the distribution at w mf min = w mf − ∆ leads to two important differences between the two types of disorder. Indeed, the sorption curves for correlated exponential disorder display a number of cusp singularities (discontinuities in the derivative of V with respect to µ) that contrast with the smooth curves for normal disorder in w mf ij (see Fig. 7). The second difference is that, for large system size, the limiting behaviour of the sorption is not the same for both types of heterogeneity, i.e. the hysteresis curves for the exponential disorder are asymmetric and of the H2-type 12 , in contrast to the parallel-sided H1-type seen for a normal distribution of w mf ij .
Similarly to the case of a normal distribution, certain values of µ can be found using Eq. (25) at which nucleation and de-pinning typically occur for both adsorption and desorption. For adsorption, the values of µ c min , µ o min and µ pin max can be found in a similar way as for a normal distribution and the resulting dependence of these characteristic values of µ on ∆ is shown in Fig. 6(a).
For the case of desorption (see Fig. 6(b)), nucleation will typically occur at the end-cells for µ ≃ µ end max , where n m 1 = n m N = 5 for type-I pores (see solid line in Fig. 6(b)) and n m 1 = n m N = 4 for type-II pores (see dashed line in Fig. 6(b)). For the inner cells, desorption will typically first nucleate at µ ≃ µ inner max (see dotted line in Fig. 6 with the additional factor 2 before w ff corresponding to two neighbouring cells being occupied by fluid as opposed to 1 at the ends and n m i = 4 for 1 < i < N . For long pores (N → ∞), the gradient of µ inner max tends to a limiting value of n m i = 4. The values of µ end max can be, depending on the values of w ff , N and ∆, either greater than or less than µ inner max for pores of both types I and II. This is in contrast to adsorption, when µ inner max < µ end max always for pores of type II.
Pinning during desorption can occur only for µ approximately greater than the mean minimal value of µ i (4, 1), µ µ pin min = min i=2,...,N −1 µ i (4, 1) , i.e. when there is some inner cell which remains occupied when one of its neighbours is unoccupied, thus impeding the propagation of the avalanche. The value of µ pin min can be calculated and it is equal to, with n m i = 4 (see dot-dashed line in Fig. 6(b)). There are several regimes for desorption in the case of exponential disorder depending on relative values of the characteristic chemical potentials, µ c max , µ inner max , µ inner max and µ pin min (see Fig. 6(b)). If the point where  in Fig. 6(b). In fact, similarly to the case of normally distributed disorder in w mf ij , the parameter space, w mf , w ff for fixed values of ∆ and N , can be split into four regions corresponding to four regimes for either adsorption or desorption. The boundaries between the different regimes have the same configuration as in Fig. 5 for normal distribution, but are characterised by different lo-cations of the points A w mf , w ff and D w mf , w ff , Below, we analyse two representative sets of parameters. For the first set of parameters, µ o min > µ c min > µ pin max , and the type-I and II pores are in regimes (3) and (4), respectively, for adsorption, and in regimes (3) and (1), respectively, for desorption. The sorption curves for this choice of parameters are presented in Fig. 7(a). It is confirmed that the adsorption in both pore types takes place in a single avalanche. Typically this avalanche takes place at a smaller value of µ in type-I pores (µ ≃ µ c min ) than in type-II pores (µ ≃ µ o min ). The desorption also agrees with the prediction, in that it takes place over a range of µ for the type-II pore around µ ≃ µ o max and in a single large avalanche for the type-I pore at a lower value µ ≃ µ c max . For the second set of parameters, for which µ o > µ c ≃ µ inner > µ pin , both pore types are in the disorder-controlled regime (2) for adsorption and desorption. For this set of parameters, the adsorption curves are practically identical (see Fig. 7(b)). The desorption curves, however, differ slightly. There are several cusps which appear in these sorption curves as a consequence of the lower bound for matrix-fluid interaction which are w mf ij ≥ w mf min = w mf − ∆. The existence of the lower bound to w mf ij leads to appearance of upper bounds to the random variables µ end max , µ inner max and µ pin max such that, which can be derived from Eq. (25) by substituting w mf ij = w mf min = w mf − ∆. Limits are given for desorption at the ends of pores of type I by µ c max = M end with n m 1 = 5 and for desorption at the ends of pores of type II by  Fig. 7(b) by arrows. The value of M pin corresponds to a cusp in the adsorption curves for pores of both types, at which V rapidly approaches unity and Var[V ] displays a sharp decay to zero. This cusp occurs because the number of pinning points rapidly approaches zero as the increasing value µ → M pin . For µ > M pin , there are no pinning points and pore will be fully occupied with probability 1 if any nucleation event has occurred. The mean occupied volume is then close to 1 in this region.
For desorption, nucleation can only occur if µ ≤ max M end , M inner , and as a consequence, the occupied volume remains equal to 1 on decreasing µ until the above condition is satisfied. When µ passes from above to below M end , the probability of nucleation at the ends of the pore begins to increase, giving rise to cusps in the solid and dashed desorption curves at µ = µ c max and µ = µ o max , respectively (see Fig. 7(b)). At these cusps the value of V begins to decrease. The cusps are only significant if the number of pinning points is small enough for the desorption at the end of the cell to cause fluid to desorb in a large part of the pore, meaning that the rate of decrease of V is large. A similar effect occurs when µ passes from above to below M inner , except that the probability of nucleation increases more rapidly with reducing µ, and the resulting cusp in the solid and dashed desorption curves is much sharper, because there are a greater number of inner cells than end cells. Typically, the cusp at M inner can only be seen when desorption is in regimes (2) or (4), when nucleation occurs in the inner cells. Outside of these regimes, desorption fully occurs at higher values of µ than M inner , and the cusp is insignificant. For a type-II pore there can be two cusps visible in regime (2), corresponding to nucleation at the end cells and at inner cells. The desorption curve of a type-I pore can also show two cusps in regime (2) if M inner < µ c max , i.e. if w ff ≤ w mf − ∆. However, if this condition is not satisfied, only one cusp will be visible for a type-I pore, because most desorption will occur as a result of nucleation in the inner part of the pore, i.e. there is no significant cusp at µ c max .
The sorption curves for large exponential correlated disorder in w mf are of the, so-called, H2-type 12 (see upper panel of Fig. 7(b)). Such a shape of sorption curves contrasts with that for Gaussian disorder which are of the parallel-sided H1-type (see upper panels of Figs. 2 and 3). It should be noted that the asymmetry in the distribution of h i = j/i w mf ij for the case of exponential disorder is a consequence of a correlated asymmetric distribution of w mf ij . Indeed, if the values of w mf ij are uncorrelated for a given value of i then the central limit theorem ensures that their sum, j/i w mf ij , is approximately distributed according to a normal distribution, which is symmetric. Therefore correlations in w mf ij play a significant role in achieving a skewed distribution of h i and H2-type hysteresis, the effect being maximal when they are fully correlated, i.e. when all w mf ij are equal to each other for a given i. This implies that H2-type hysteresis in heterogeneous pores might arise due to variations in pore diameter (represented by correlated disorder in w mf ij ) rather than due to individual defects, in agreement with previous numerical studies 19 . On the other hand, symmetric (normal) disorder in local fields (i.e. uncorrelated disorder in w mf ij ) can only cause a parallel sided H1-type hysteresis ( Fig. 2 and 3), and may represent the effect of uncorrelated structural defects in the pore surface on small length scales.

D. Pores of type III
In this section, we present results on the effects of interaction between the two parts of a pore, each of length N 1 = N 2 = 50, characterised by different mean matrixfluid interaction strength w mf 1 and w mf 2 (as shown in Fig. 1(c)). The difference in the interaction strength of the fluid with the matrix can represent variable pore diameter for a pore of either a funnel or an ink-bottle shape. A major feature of sorption in such pores is that, if the difference between w mf 1 and w mf 2 is large enough, it occurs in the narrow part of the pore at a lower value of µ than it does in the wider part. This gives rise to two steps in the sorption curves (see upper panels of Fig. 8), in agreement with experimental observations 11,38 for both pore shapes. In order to understand the form of the sorption curves in such pores, it is helpful to compare it with the mean of the sorption curves for two independent pores of type II of length N 1 = N 2 = N = 50, one with w mf = w mf 1 and the other with w mf = w mf 2 (see dotted lines in upper panel of Fig. 8). The sorption curves (solid and dashed lines) differ from the dotted lines representing the behaviour of two independent type-II pores. These differences are due to the interaction between fluid in the two parts of the pore and of fluid with the closed end-cap of the type-III pore.
We analyse these differences, first, for a pore with nor-mal disorder in w mf ij (see dashed lines in Fig. 8(a)) and a funnel shape, characterised by w mf 1 < w mf 2 , i.e. the matrix-fluid interaction is weaker in the part of the pore with the larger diameter. Adsorption in each part of the pore is enhanced in comparison with a type-II pore (the dashed line corresponding to adsorption is shifted to the left by around 0.25 with respect to the dotted line for the two steps). The shift of the lower step (at µ ≃ −5.5) is due to the increased matrix-fluid interaction between the narrow section and the closed end and the shift of the upper step (at µ ≃ −4.75) is due to the fluid-fluid interaction between the wide section and the narrow section. Conversely, desorption in the funnel pore occurs later than in two independent open-ended pores (the dashed line corresponding to desorption is shifted to the left by around 0.1 with respect to the dotted line). This is because the funnel pore has only one open end at which desorption can nucleate and thus it occurs later than in a pore open at both ends.
Second, we compare the sorption curves for an inkbottle pore with the sorption curves of two independent pores (cf. solid and dotted lines in Fig. 8(a), respectively). Both the adsorption and desorption curves for the narrow part of the ink-bottle pore closely match those for an open ended pore of the same diameter (the dotted and solid lines practically coincide in the lower step until the shoulder develops). Qualitatively, for this part of the adsorption and desorption curves, the wider part of the pore is not occupied by the fluid, meaning that the narrow part effectively has two open ends. The upper step in the sorption curves, after the shoulder, represent adsorption and desorption in the wide part of the ink-bottle pore and differ from an open ended pore (cf. the upper steps in the solid and dotted lines in Fig. 8(a)). This is because sorption in the wide part occurs when the narrow part is fully occupied with fluid, meaning that the wide part of the ink-bottle pore behaves like a pore closed at both ends. As such, adsorption occurs for smaller values of µ (the upper step in the solid adsorption curve is shifted to the left by around 0.25 compared to the dotted line), since adsorption is nucleated either by the closed end or by the fluid which has already condensed in the narrow part. Similarly, desorption occurs at a lower value of µ, since it is not nucleated at either end (the upper step in the solid desorption curve is also shifted by around 0.25 to the left compared to the dotted line). It can be noted that the behaviour on desorption observed here for ink-bottle pore and funnel pore topologies is similar to that experimentally observed in Ref. 38. The adsorption curves observed here differ from the experimental results of Ref. 38, in that condensation in the wider part of an ink-bottle pore is observed at higher values of µ than in the wider part of a funnel pore. This might be related to the fact that, in our model, condensation is allowed to occur in the inner part of the ink-bottle pore regardless of the fact that gas has no way to flow into the wider part of the pore (the pore is blocked by fluid in the narrow part) 39 .
Results for exponential disorder in w mf ij are shown in Fig 8(b). The dotted line corresponds to a pore with two open ends and the dashed line corresponds to a funnel pore, which effectively has a single open end. The solid line shows the isotherms for the ink-bottle pore where, effectively, the wide part has two closed ends and the narrow part has two open ends. As such, the main differences between the different sorption curves are similar to the differences observed between open-and closedended pores in Fig. 7(b), i.e. the adsorption curves are broadly similar for both funnel and ink-bottle pores. The main differences between the isotherms for different pore geometries is observed in two intervals of µ: −5.6 < µ < −4.6 and −6.5 µ −6. In the first interval, desorption can be nucleated at the open end of the wide part of the funnel pore, but nucleation cannot occur in the wide part of the ink-bottle pore because this part effectively has both ends closed. This is why the dashed line is below the solid line in this region. In the second interval, the wide parts of both the funnel and ink-bottle pores are empty. This means that the narrow part of the funnel-shaped pore effectively has one open end while the narrow part of the ink-bottle pore effectively has two open ends. Desorption is therefore more likely to be nucleated in the narrow part of the ink-bottle shaped pore than in the narrow part of the funnel pore. Correspondingly, the dashed line is above the solid line in this region.
In the disorder-controlled regime (see Fig. 8(c)), the sorption curves for both funnel and ink-bottle pores are identical to those found by adding together the sorption curves for two independent open-ended pores of different diameters obtained by cutting the funnel and ink-bottle pores in half, i.e. the effects of the ends of the pores are small for disorder on this scale. This behaviour appears to correspond to the experimental observations in Ref. 11.
The effect of finite temperature on sorption processes is analysed here for both a funnel and ink-bottle pore. In particular, the conditions corresponding to Figs. 8(a) and 8(c) were repeated for several values of temperature and with a constant rate, r, of change of µ. The results are shown in Fig. 9.
For finite sufficiently low temperature, β = 50, and Gaussian disorder of strength ∆ = 0.1, it can be seen by comparing Figs. 8(a) and 9(a) that the mean volume of condensed fluid closely matches that at zero-temperature for both ink-bottle and funnel pores. The main difference from the zero-temperature and small r limit is that the meniscus propagates slowly along the pore rather than the whole pore becoming occupied simultaneously, reflected in a reduction in the maximum value of Var [V ]. For instance, the maximum value of the variance for an ink-bottle shaped pore is Var[V ] ≃ 0.04 (see peak of solid line in Fig. 9(a) middle panel) in comparison with Var[V ] ≃ 0.06 at zero-temperature (see solid line in Fig. 8(a), lower panel). For larger temperatures, β = 5, hysteresis loops become narrower (cf. dot-dashed curve with solid curve in Fig. 9(a)).
For exponential disorder of strength ∆ = 0.25, the main effect of the finite temperature is in the smoothing of the cusps in the hysteresis loops (see Fig. 9(b) upper panel). As the temperature is increased, first to β = 15 (double-dot dashed curves) and then to β = 5 (dotdashed curve), the area of the hysteresis loop reduces gradually. This observation with increasing temperature agrees with the behaviour observed for 1D pores 38 and is similar to the decreasing width of sorption hysteresis loops in 3D porous media 2 . Although the behaviour in 1D and 3D systems is similar, the proposed model and its mapping to the RFIM suggests that the origin of hysteresis is not identical in both cases. Indeed, the disordertemperature phase diagram of a 3D lattice gas (or the 3D RFIM) consists of a phase at low temperature and disorder with ferromagnetic order and a paramagnetic phase at high temperature and disorder 40 . In the ferromagnetic phase, the free energy consists of two global minima such that gas and liquid (or states with positive and negative magnetisation in the RFIM) could in principle coexist in the thermodynamic limit 1 . Hysteresis is associated in this case with both the global minima of the free energy and the existence of many local minima where the system can remain trapped for very long times. The mapping of the proposed model to the 1D-RFIM indicates that ferromagnetic order does not exist as a stable macroscopic phase at any finite temperature and/or disorder 40 . Instead, the free energy exhibits a single global minimum. Therefore, hysteresis in linear pores at non-zero disorder and temperature is expected to be only associated with the rugged character of the energy landscape which consists of many local minima.

VI. CONCLUSIONS
To conclude, a heterogeneous lattice-gas model has been proposed to describe fluid condensation in 1D pores of different shapes and rough surfaces. Heterogeneity, missed in classical theories, is the key and sufficient feature of the model which allows it to reproduce the main experimental findings. We demonstrate that a simple coarse-grained representation of pores consisting of 1D chains of cells is a minimal model sufficient to account for the effects of heterogeneity. Within a single cell, liquid interacts with the elements of the surface. These interactions can be either identical or variable (random) within the cell, and also can vary for different cells. Accounting for such different types of heterogeneity in the model, results in the different shapes of hysteresis loop found experimentally 12 , including both the H2-type 11 and the H1-type 41 hysteresis loops.
In addition, the model is able to reproduce the shape of the sorption curves for some more complex experi- mentally studied systems, including ink-bottle and funnel shaped pores. The physical phenomena inherent for this model include nucleation of adsorption and desorption (cavitation) and propagation of a meniscus through the pore, which are known 38 to be two main effects observed in such systems. In this respect, our model could be easily extended to account for the fluid blocking effect that prevents flow of gas to the wider part of the ink-bottle structure 38,39 .
Besides reproducing experimental observations, our model also suggests interesting predictions that motivate new experiments. For instance, we have demonstrated that the sorption mechanisms (i.e. whether it starts at the ends of pores or in the interior part) might depend on the length and diameter of the pores and the degree of heterogeneity. However, for large heterogenity, or long pores, the sorption tends to a limiting, and apparently universal, disorder-controlled regime.
In physical systems exhibiting a more complex topology, e.g. a 3D maze-like network of 1D channels, an analytical solution may be nontrivial, however numerical simulations within our model can be performed straightforwardly for porous media of arbitrary topology 42 .

VII. ACKNOWLEDGEMENTS
TPH acknowledges the support of EPSRC.