## 2 Description of Model and Simulation Parameters
It can be seen as a kind of advection equation of the probability distribution function of fictitious particles, \(f\) (\(t\), \(x\), \(c\)), with an external body force term, that can be written for each \(i\) of \(Q_{m}\) directions (lattice linkages) in velocity space as follows:  
\[\frac{\partial f_{i}}{\partial t}+c_{i}\,\cdot\
abla_{x}f_{i}+\frac{F}{\rho} \,\cdot\,\
abla_{c}f_{i}=\Omega_{i} \tag{1}\]  
where \(\mathbf{x}\) [\(x_{x}\), \(x_{y}\), \(x_{z}\)] is a spatial point coordinate vector, \(\mathbf{c_{i}}\) [\(c_{x}\), \(c_{y}\), \(c_{z}\)] is the lattice velocity vector: \(c_{i}=\frac{\delta_{x}}{\delta_{i}}e_{i}\), where \(\mathbf{e}_{i}\) is the unit vector along the \(i\) lattice direction, \(\frac{\delta_{x}}{\delta_{i}}\) is the lattice speed = 1 as \(\delta_{x}\) and \(\delta_{t}\), that are respectively the lattice spacing and the lattice time step, and both equal 1, \(t\) represents time, \(\mathbf{F}\) [\(F_{x}\), \(F_{y}\), \(F_{z}\)] is an external body force vector, \(\Omega\) is the collision operator--a source term, which represents the local redistribution of \(f\) due to collisions, and \(\rho\) is a LB mass density:  
\[\rho=\sum_{i=0}^{i=m}f_{i} \tag{2}\]  
\(\mathbf{u}\) [\(u_{x}\),\(u_{y}\),\(u_{z}\)] for incompressible fluid can be calculated from:  
\[u=\sum_{i=0}^{i=m}f_{i}c_{i}+\frac{1}{2}F \tag{3}\]  
where \(\
abla_{x}:=\left[\frac{\partial}{\partial x_{x}},\frac{\partial}{\partial x _{y}},\frac{\partial}{\partial x_{z}}\right]\) and \(\
abla_{c}:=\left[\frac{\partial}{\partial c_{x}},\frac{\partial}{\partial c _{y}},\frac{\partial}{\partial c_{z}}\right]\). Figure 2: The numerical grid: (**a**) View from the top; (**b**) View from the side; (**c**) The cross-section through the reservoir. Grid nodes colored by node type: 1—fluid, 36—north pressure, 41–46—“velocity 0”, 61—full-way bounce-back. (A, B) Cross-section planes. A-A cross-section through A plane. (C, D) Lines placed in axes. The channel length is on the X-axis, the channel width is on the Y-axis, and the channel height is on the Z-axis. The collision operator, \(\Omega_{i}\), was approximated with the most popular Bhatnagar-Gross-Krook model [24]:  
\[\Omega_{i}=-\frac{1}{\tau}\Big{[}f_{i}-f_{i}^{eq}\Big{]} \tag{4}\]  
and the external body force term in Equation (1) with the He-Shan-Doolen scheme (HSD) [25];  
\[\frac{F}{\rho}\ \cdot\ \
abla_{c}f_{i}=\frac{\frac{F}{\rho}\cdot(c_{i}-u)}{c_{s}^ {2}}f_{i}^{eq} \tag{5}\]  
where \(\tau\) is the LB relaxation time related to the lattice fluid viscosity: \(\upsilon=1/3(\tau-1/2)\), and \(f_{i}^{eq}\) is the equilibrium distribution function along the \(i\) lattice direction given by the following formula for incompressible fluid [26]:  
\[f_{i}^{eq}=\omega_{i}\Bigg{[}\rho+\frac{c_{i}\cdot u}{c_{s}^{2}}+\frac{(c_{i} \cdot u)^{2}}{2c_{s}^{4}}-\frac{u^{2}}{2c_{s}^{2}}\Bigg{]} \tag{6}\]  
where \(c_{s}\) is the lattice constant: \(c_{s}=\frac{1}{\sqrt{3}}\), \(u\) is the macroscopic fluid velocity vector, and \(\omega_{i}\) is a weighting scalar for the \(i\) lattice direction. In MF3D, all solution parameters are collected in case_params.cfg and thread_params.cfg (Supplementary Files S2 and S3). The flow simulations were performed for spontaneous gravitationally induced fluid flow in CC. The force density vector was calculated for the microsystem sloped at an angle of 30\({}^{\circ}\): \(F=\rho\ \times g=1000\ \times\ 9.81=9810\ (\mathrm{N/m^{3}})\) (\(F_{x}=4905\), \(F_{y}=0\), \(F_{z}=-8495.7\)), and its constant value was fixed for the whole region. A free solution parameter--the relaxation time, \(\tau\)--was set to 1.0. The solution termination criteria were set for scaled velocity residue to 1.0 \(\times\) 10\({}^{-7}\) and 0.1% for the mass flow balance error. The MF3D configuration files necessary to run the simulation are available in additional resources of the article (Supplementary File S4).