Corentin Jouault
email: jouaultc0@gmail.com]

Gi-Soo Nam

André Nel

The first Geinitziidae (Polyneoptera: Reculida) from the Upper Triassic Amisan Formation of Korea

Keywords: Shurabia, Reculida, new species, fossil record, taxonomy Shurabia, Reculida, nouvelle espèce, registre fossile, taxinomie

Shurabia taewani sp. nov., is first representative of the reculid family Geinitziidae described from the Upper Triassic of South Korea. The preservation of numerous wing venation characters on this new fossil allows for its placement in the genus Geinitzia and a deep comparison with other genera of Geinitziidae. This discovery suggests that the insect paleofauna of the Amisan Formation is rich and needs to be further investigated. Shurabia taewani sp. nov. differs from the other species currently included in the genus, inter alia, because its vein ScP ends near wing mid-length, its fork of M is at the same level as the fork of R, its RP has five branches, its MA four branches, its MP two branches, and its vein m-cua 'M5' originates from the stem of M and ends into CuA.

Résumé

Shurabia taewani sp. nov., est le premier représentant de la famille de Geinitziidae (Reculida) décrit du Trias supérieur de Corée du Sud. La préservation de nombreux caractères de nervation alaire sur ce nouveau fossile permet de le ranger dans le genre Shurabia par une comparaison complète avec les autres genres de Geinitziidae. Cette découverte suggère que la paléo-faune de la Formation d'Amisan est riche et doit faire l'objet d'une étude plus approfondie. Shurabia taewani sp. nov. diffère des autres espèces actuellement incluses dans le genre, inter alia, par sa nervure ScP se terminant près du milieu de l'aile, la fourche de M au même niveau que celle de R, la RP avec cinq branches, la MA avec quatre branches, sa MP avec deux branches et la vein m-cua 'M5' provenant de la tige de M et se terminant dans CuA.

Introduction

The 'Grylloblattodea' is a group, likely polyphyletic or paraphyletic, of insects composed of the extant apterous order Grylloblattida Walker, 1914, and numerous fossil families and genera of uncertain affinities. To date, there is still no evidence (synapomorphy) supporting the relationships within this group or the monophyly of the 'Grylloblattodea'. Nevertheless, [START_REF] Storozhenko | Sistematika, filogeniya i evolyutsiya grilloblattidovykh nasekomykh (Insecta: Grylloblattida) [Systematics, phylogeny and evolution of the grylloblattids (Insecta: Grylloblattida)[END_REF] proposed a first revision and phylogenetic analysis of the group providing crucial state of art and information for the following decade. Unfortunately, his phylogenetic analysis was not polarized (i.e., no outgroup) and needs to be updated. Recently, Aristov (2015a) revised the group but additional work is to be conducted to fully disentangle the relationships within this lineage. Aristov (2015a) notably transferred some former 'Grylloblattodea' families, such as the speciose Geinitziidae, into the order Reculida.

Despite the recent progress in phylogenetic analyses, the systematic of many clades within the 'Grylloblattodea' and the Reculida remains blurry and the current delineations of the families are poorly justified. This issue stems from the fossil material -often poorly preserved and partial (i.e., wings only) -that often consists of a single isolated wing that does not allow for the circumscription of the variability of the venation (a similar problem is known for other Polyneoptera see [START_REF] Jouault | A new stonefly species (Plecoptera: Perlodidae) from Eocene Baltic amber and questions on the wing venation potential for species diagnostic of fossil Plecoptera[END_REF]. It is to be hoped that future descriptions will bring additional information on the morphology of this group, and its evolution to refine the limit of its constitutive families.

The fossil representatives of the 'Grylloblattodea' including the Reculida span from late Carboniferous to mid-Cretaceous with numerous species known from the Permian and Triassic epochs [START_REF] Cawood | The first 'Grylloblattida' of the family Liomopteridae from the middle Permian in the Onder Karoo, South Africa (Insecta: Polyneoptera)[END_REF]. Similarly, they are worldwide known -even in Gondwanan deposits [START_REF] Lara | First records of Geinitziidae (Insecta: Grylloblattida) from the Upper Triassic of Argentina (Mendoza)[END_REF], which are often poorly documentedsuggesting that they could be used to investigate the impact of at least three extinctions events: the end-Permian mass extinction, the Triassic-Jurassic extinction, and the Toarcian anoxic event [START_REF] Zhang | Widespread grylloblattid insects after the End Permian mass extinction[END_REF]. However, before proposing hypotheses on the potential impact of these extinctions on the 'Grylloblattodea' and the Reculida in particular, it is necessary to clarify their systematic. Embracing this vision, we describe the first species and representative of Geinitziidae from the Korean Peninsula and provide a detailed justification for its placement. This species belongs to Shurabia Martynov, 1937 one of the most speciose genera of the family, known between the Middle Permian and the Late Jurassic, with a maximum of species in the Jurassic (Table . 1).

Material and methods

A single forewing with preserved venation was discovered in the middle shale unit of the Amisan Formation exposed in the Myeongcheon section (36°20'21"N, 126°37'34"E; see [START_REF] Park | A new titanopteran Magnatitan jongheoni n. gen. n. sp. from southwestern Korean Peninsula[END_REF] fig. 1). The Amisan Formation is known to 1,000 m thick and consists of the lower sandstone units, the lower shale unit, the middle sandstone unit, the middle shale unit, and the upper sandstone unit [START_REF] Yang | The lower Mesozoic strata[END_REF]. To date, the Amisan Formation is thought to be of Triassic age but the floral analyses of [START_REF] Kimura | Geological age of the Daedong flora in the Korean Peninsula and its phytogeographical significance in Asia[END_REF] have to be challenged because of the life span extension (to the early Jurassic or Cretaceous) of numerous genera initially assumed to be restricted to the Late Triassic [START_REF] Kim | Lobatannularia nampoensis Kawasaki from the Upper Triassic Baegunsa Formation, Nampo Group, Korea[END_REF][START_REF] Kim | Three new Equisetites species found from the Upper Triassic Amisan Formation, Nampo Group, Korea[END_REF][START_REF] Kim | Species diversity and leaf form of ginkgoaleans from the Mesozoic and Cenozoic strata in Korea[END_REF][START_REF] Kim | Weltrichia sp. from the Late Triassic Amisan Formation of Nampo Group, Korea[END_REF][START_REF] Kim | A new species of Leptostrobus from the Upper Triassic Amisan Formation of the Nampo Group in Korea[END_REF][START_REF] Lee | A revision of Mesozoic Equisetales Annuriopsis bunkeiensis Kimura et Kim from the Amisan Formation of Nampo Group, Korea[END_REF][START_REF] Kim | Organ fossils of Neocalamites carrerei from the Amisan Formation of the Nampo Group, Korea[END_REF]. The analysis of clam shrimp fauna supports a Triassic age [START_REF] Kim | Margaritifera cf. isfarensis (Chernishev) from the Amisan Formation, Nampo Group, Korea[END_REF], but the only bivalve genus Margaritifera Schumacher, 1816 known from the Formation is possibly indicative of a post-Triassic age [START_REF] Kim | Margaritifera cf. isfarensis (Chernishev) from the Amisan Formation, Nampo Group, Korea[END_REF], although a recent phylogenetic analysis of the Margaritiferidae supports a Triassic origin for this family [START_REF] Araujo | The origin and phylogeny of Margaritiferidae (Bivalvia, Unionoida): a synthesis of molecular and fossil data[END_REF]. Radiometric analyses established a Jurassic age for the Nampo Group, encompassing the Amisan Formation, questioning the Triassic age derived from the paleontological evidences [START_REF] Koh | Tectonic implication of the Mungyeong-Jeongseon tectonic line, the Yeongweol Nappe and the Bansong Group in the Ogcheon belt[END_REF][START_REF] Jeon | Early Archean to Middle Jurassic evolution of the Korean Peninsula and its correlation with Chinese cratons: SHRIMP U-Pb zircon age constraints[END_REF]. Finally, the composition of insect fauna, with Titanoptera and Triadophlebiomorpha, indicates a Triassic age for the Formation [START_REF] Park | A new titanopteran Magnatitan jongheoni n. gen. n. sp. from southwestern Korean Peninsula[END_REF][START_REF] Jouault | Koreaphlebia gen. nov. (Odonatoptera: Triadophlebiomorpha): new evidence of a Triassic age for the Amisan Formation in Korea[END_REF]. Herein, we consider the age of the Amisan Formation to be Late Triassic until clear-cut analyses clarify the age of the formation.

The samples were photographed using a Canon EOS 6D camera with an attached Canon EF 100 mm f/2.8 USM macro lens. Images were cropped and enhanced using Adobe

Systematic paleontology

Class Insecta Linnaeus, 1758 Clade Polyneoptera Martynov, 1938Order Reculida Handlirsch, 1906Family Geinitziidae Handlirsch, 1906Genus Shurabia Martynov, 1937 Type species: Shurabia ovata [START_REF] Martynov | Liassic insects from Shurab and Kisyl-Kiya, Part I, Various orders e x c e p t B l a t t o d e a a n d C o l e o p t e r a . A k a d e m i y a N a u k S S S R[END_REF] Other species included in the genus: Shurabia aleda [START_REF] Wappler | New orthopteroid insects from the Upper Triassic Molteno Formation, Karoo Basin, southern Africa[END_REF], Shurabia angustata [START_REF] Martynov | Liassic insects from Shurab and Kisyl-Kiya, Part I, Various orders e x c e p t B l a t t o d e a a n d C o l e o p t e r a . A k a d e m i y a N a u k S S S R[END_REF], Shurabia anomala Rasnitsyn, 1982, Shurabia australis Rasnitsyn, 1982, Shurabia bashkuevi Aristov, 2011, Shurabia creta Aristov, 2020, Shurabia elegans [START_REF] Fujiyama | Mesozoic insect faunas of east Asia, part I. Introduction and Upper Triassic faunas[END_REF], Shurabia ferganensis [START_REF] Rasnitsyn | Triassic and Jurassic Insects of the genus Shurabia (Grylloblattida, Geinitziidae)[END_REF], Shurabia grandis [START_REF] Huang | New "Grylloblattida" related to the genus Prosepididontus Handlirsch, 1920 in the Middle Jurassic of China (Insecta: Geinitziidae)[END_REF], Shurabia hissarica [START_REF] Aristov | New and little-known grylloblattids of the family Geinitziidae (Insecta: Grylloblattida) from the Triassic and Jurassic of Europe, Asia, and South Africa[END_REF], Shurabia inferior Aristov, 2015, Shurabia izyumica Aristov, 2022, Shurabia kapokkraalensis (Wappler, 2001), Shurabia lukashevichae [START_REF] Aristov | New and little known Grylloblattida (Insecta) from intertrappean deposits of the Tunguska Basin of Siberia[END_REF], Shurabia magna Rasnitsyn, 1982, Shurabia minuta Rasnitsyn, 1982, Shurabia minutissima Aristov, 2022, Shurabia ovata Martynov, 1937, Shurabia parvula Rasnitsyn, 1982, Shurabia permiana Aristov, 2013, Shurabia postiretis (Huang et al., 1991), Shurabia serrata [START_REF] Aristov | New and little-known grylloblattids of the family Geinitziidae (Insecta: Grylloblattida) from the Triassic and Jurassic of Europe, Asia, and South Africa[END_REF], Shurabia sogutensis Rasnitsyn, 1982, Shurabia tanga Aristov, 2018, and Shurabia taewani sp. nov.

Shurabia taewani sp. nov. (Figs. 1-2) Zoobank xxxx

Etymology

The specific epithet honors Mr. Kim Taewan, a mentor of (G-SN) and renowned fossil collector in Korea; and is to be treated as a noun in a genitive case.

Material

Holotype GNUE112003 (a nearly complete wing, with basal and distal part slightly damaged), housed in the collection of the Gongju National University of Education, Gongju, Republic of Korea.

Locality, unit and age

Amisan Formation, Upper Triassic, Myeongcheon Section, Seongju-myeon, Boryeong City, Chungcheongnam-do, Republic of Korea.

Diagnosis

Space between each vein relatively broad; costal space broad; ScP ending near wing midlength, slightly s-shaped; fork of M and fork of R into RA and RP at same level; RP with five branches, not pectinate; MA with four branches (three directed anteriorly), distinct from RP along its entire length; first fork of MA slightly anteriad fork of MP; fork of the anterior branch of MA close to wing apex; MP with two branches; m-cua ('M5') originating from stem of M and ending into CuA; fork of CuA1 located near wing mid-length, space between posterior branch of CuA1 and anterior branch of CuA2 broad; CuA1 with two branches.

Description

Macropterous forewing, at least 10 mm long (likely ca. 12-13 mm if fully preserved) and ca.

4.5 mm wide, sides convex; crossveins present between main veins (number impossible to certify because of preservation); ScP short and slightly sigmoidal, with at least 10 branches ending into anterior wing margin; costal space more than twice as broad as subcostal space; R with two distinct branches (i.e., no fusion); RA branches ending into anterior wing margin, first fork of RA located well before first fork of RP, second fork of RA slightly anteriad first fork of RP; space between RA and RP slightly wider near first fork of RA; RP field nearly twice as broad as RA field; RP branches ending into anterior wing margin, first fork of RP closer to MA that RA; M divided into two distinct branches (i.e., no fusion); MA encompassing wing apex; MA fork located slightly anteriad fork of MP, second fork of MA distad first fork of MP; MP field reduced, nearly half the size of MA field; space between M (or MA) and CuA broad, broader than space between R (RP) and M (MA); m-cua ('M5') vertical; fork of CuA located distad m-cua; CuA1 fork located closer to M fork than to MP fork.

Color: colored stripes along crossveins and main branches (difficult to detail).

Remark: The wing is partly preserved and the breaking point along the posterior margin likely follows the CuP (even if the vein itself is not preserved and therefore not figured on our drawing).

Discussion

Systematic placement and discussion on the limit of wing venation characters

Following the key to grylloblattodean families of [START_REF] Storozhenko | Sistematika, filogeniya i evolyutsiya grilloblattidovykh nasekomykh (Insecta: Grylloblattida) [Systematics, phylogeny and evolution of the grylloblattids (Insecta: Grylloblattida)[END_REF], the new specimen keys out in the family Geinitziidae because of its wings fully developed; the vein m-cua developed (i.e., strong and oblique crossvein between M and CuA, so-called 'M5'); ScP with anterior branches, ending into C; MP fully developed but medially desclerotized; RA with anterior branches. Additionally, the diagnostic characters of the forewings, as detailed and translated from [START_REF] Storozhenko | Sistematika, filogeniya i evolyutsiya grilloblattidovykh nasekomykh (Insecta: Grylloblattida) [Systematics, phylogeny and evolution of the grylloblattids (Insecta: Grylloblattida)[END_REF], are all present in the new specimen: forewing of medium size, membranous, densely and finely pubescent, with a characteristic pattern of transverse spots and stripes, more or less clearly attached to the crossveins; apex of the wing broadly rounded; costal field moderately wide, usually crossed by simple branches of ScP; RA branches reaching anterior wing margin; fork of R located anteriad wing mid-length; RP forked with branches in same axis as its stem; M divided into MA and MP, with its fork located anteriad wing mid-length; stem of MP (i.e., before its fork) desclerotized; strong crossvein m-cua (M5) present and originating from M; CuA divided into CuA1 and CuA2, both distinct, with three to five branches directed toward posterior wing margin; CuP simple, feebly sclerotized; area between CuA and CuP slightly widened basally, crossed by several crossveins or two-three rows of cells; anal field small, with two-four simple anal veins; crossveins mostly simple, sometimes (exceptionally) forming two-three rows of cells.

To date the delineation of the grylloblattodean families is unsatisfactory, resulting from the lack of apomorphic characters, incompletely documented intraspecific variability, and homoplasy. Nevertheless, the following genera are currently included or considered to belong to the family Geinitziidae: Fletchitzia Riek, 1976, Geinitzia Handlirsch, 1906, Geinitziella Aristov, 2022, Permoshurabia Aristov, 2009, Permovalia Aristov, 2015b, Permuliercula Aristov, 2020, Prosepididontus Handlirsch, 1920, Sauk Aristov, 2018[START_REF] Aristov | New insects (Insecta: Trichoptera, Reculida, Eoblattida) from the Mesozoic of Asia[END_REF], Say [START_REF] Aristov | New insects (Insecta: Trichoptera, Reculida, Eoblattida) from the Mesozoic of Asia[END_REF][START_REF] Aristov | New insects (Insecta: Trichoptera, Reculida, Eoblattida) from the Mesozoic of Asia[END_REF], Shurabia [START_REF] Martynov | Liassic insects from Shurab and Kisyl-Kiya, Part I, Various orders e x c e p t B l a t t o d e a a n d C o l e o p t e r a . A k a d e m i y a N a u k S S S R[END_REF], Sinosepididontus Huang and Nel, 2008, Stegopterum Sharov, 1961, and Sukhonia Aristov, 2013. Cui et al. (2012: 259) proposed a key to the genera of Geinitziidae, and following the latter the new specimen keys out in the genus Shurabia because of the following couplets: M with more than three branches; RA branched; M branching on MA and MP before or near the rise of RP; CuA with three branches; RP irregularly pectinate, forming an anterior or posterior comb of branches. Aristov (2020) further considered that a vein CuA with four branches is typical for the genus Geinitzia, while species of the genus Shurabia have a CuA with only three branches. Therefore, affinities with the genus Shurabia appear likely, but we prefer to investigate additional differences with the other genera currently placed in the Geinitziidae.

The new specimen differs from the genus Fletchitzia because the fork of M is at the same level as the fork of R (vs. fork of M far distad the fork of R in Fletchitzia), the fork of CuA1 is located near wing mid-length (vs. distad), CuA is only composed of dichotomous CuA1 and CuA2 (vs. CuA1 with four branches), and the lack of crossveins between subcostal veinlets (vs. present) (Riek, 1976).

The ScP curved and relatively short (vs. straight and long in Sinosepididontus), and the developed RA with three branches reaching the anterior margin of the wing precludes affinities with the genus Sinosepididontus that only possesses a simple RA [START_REF] Huang | New "Grylloblattida" related to the genus Prosepididontus Handlirsch, 1920 in the Middle Jurassic of China (Insecta: Geinitziidae)[END_REF]. Note that on the drawing on S. chifengensis in Huang and Nel (2008: fig. 2), the CuA2 is in fact the posterior branch of the CuA1, which is dichotomous. Similarly, the CuA3 figured in the original description of S. chifengensis is the CuA2. Cui et al. (2012: fig. 6D) interpreted the CuA of Sinosepididontus with three branches but it is clearly a misinterpretation of the CuA1, with two branches, and the CuA2 simple.

The genus Geinitziella was recently described from the Triassic of Kyrgyzstan and not included in the key of [START_REF] Cui | New and little-known species of Geinitziidae (Insecta: Grylloblattida) from the Middle Jurassic of China, with notes on taxonomy, habitus and habitat of these insects[END_REF]. [START_REF] Aristov | New Geinitziidae (Insecta, Reculida) from the Permian and Triassic of Eurasia[END_REF] proposed the following diagnosis for the genus (vein nomenclature adapted): ScP slightly S-shaped. RP comb-like anteriorly.

CuA1 with three branches, comb-like posteriorly, CuA2 straight. The interradial space is narrowed at its base. The ScP slightly S-shaped is not clearly visible in the photograph or the drawing proposed to illustrate the genus (Aristov, 2022: fig 1e, pl. 4,4). Aristov (2022: 266) suggested that the main difference between the new genus and Geinitzia is the condition of CuA1 (dichotomous in Geinitzia vs. with three branches in Geinitziella). Because Shurabia species have a CuA1 dichotomous (similar condition found in Geinitzia), we consider this difference sufficient to preclude affinities of the new fossil with the genus Geinitziella.

The configuration of the ScP, i.e., short and slightly s-shaped, differs between the new fossil and Permoshurabia. In fact, ScP is nearly straight and not s-shaped in Permoshurabia [START_REF] Aristov | Review of the stratigraphic distribution of Permian Grylloblattida (Insecta), with descriptions of new taxa[END_REF][START_REF] Lara | First records of Geinitziidae (Insecta: Grylloblattida) from the Upper Triassic of Argentina (Mendoza)[END_REF]. Additionally, the new fossil further differs from this genus owing to its RP with five branches (vs. four in Permoshurabia).

Permovalia was created based on specimens from the Belebeevo Formation of Russia (Aristov, 2015b). The genus is considered to be brachypterous (the new fossil is macropterous) and possesses a MP simple (vs. forked) (Aristov, 2015b: fig. 4). The validity of this genus is to be considered as doubtful given the relatively poor preservation of the type material, but we restrain from proposing any taxonomic changes until additional material is found.

The new fossil differs from Permuliercula because of the vein m-cua (M5) connected to CuA (vs. connected to CuA1 in Permuliercula). The particular configuration of m-cua found in Permuliercula is shared with Prosepididontus precluding affinities with the latter genus (Aristov, 2020a). The new specimen further differs from Permuliercula because it possesses a vein RA with three branches (vs. two in Permuliercula) (Aristov, 2020a: fig. 1h).

The genus Prosepididontus was described from the Jurassic of Germany [START_REF] Handlirsch | Palaeontologie[END_REF][START_REF] Ansorge | Identity of Prosepididontus calopteryx Handlirsch 1920 (Insecta: Grylloblattida: Geinitziidae)[END_REF]. The new specimen strongly differs from this genus because of its RA with three branches (vs. simple in Prosepididontus), RP with five branches (vs. three), and MA with four branches (vs. two). Note that m-cua is connected to CuA1 in

Prosepididontus.

The new fossil differs from the genus Sauk because its vein ScP ends near wing midlength (vs. in distal one-third of the wing for the genus Sauk), RP starting before wing midlength (vs. near the middle of the wing), and CuA1 forked (vs. simple) [START_REF] Aristov | New insects (Insecta: Trichoptera, Reculida, Eoblattida) from the Mesozoic of Asia[END_REF].

The genus Say is characterized by a vein RA with four branches (vs. three in the new fossil), MA and MP with more than two branches (vs. with two branches), and CuA1 with three branches (vs. with two branches) [START_REF] Aristov | New insects (Insecta: Trichoptera, Reculida, Eoblattida) from the Mesozoic of Asia[END_REF].

The new fossil cannot be attributed to the genus Stegopterum because it possesses a forewing without a trace of sculpture (vs. forewing with finely tuberculate sculpture in Stegopterum) [START_REF] Sharov | Order Protoblattodea, Paraplecoptera, Plecoptera, in Paleozoiskie n a s e k o m y e k u z n e t s k o g o b a s s e i n a . A k a d e m i y a N a u k S S S R[END_REF]Aristov, 2020a). Additionally, most of the Stegopterum species have long ScP (i.e., terminating distad wing mid-length), with numerous veinlets while the ScP of the new fossil is comparatively shorter and with fewer veinlets [START_REF] Sharov | Order Protoblattodea, Paraplecoptera, Plecoptera, in Paleozoiskie n a s e k o m y e k u z n e t s k o g o b a s s e i n a . A k a d e m i y a N a u k S S S R[END_REF]Aristov, 2020a: fig. 1c,D).

The genus Sukhonia strongly differs from the new fossil because of its thin costal space (vs. broad in the new fossil), RP encompassing wing apex (vs. finishing anteriad wing apex along anterior wing margin), basal branch of MA fused with MP (vs. not fused), and CuA1 with three branches (vs. two) (Aristov, 2013: fig. 73g).

The genus Roemerula [START_REF] Bode | Die Insektenfauna des Ostniedersachsischen Oberen Lias[END_REF] was placed in the Geinitziidae by [START_REF] Bode | Die Insektenfauna des Ostniedersachsischen Oberen Lias[END_REF] but this placement has never been revised. The description of the genus Roemerula is old and the illustration poor (Bode, 1953: pl. 1, fig. 6). The drawing of the wing of Roemerula is not sufficiently accurate to provide a detailed comparison with the new fossil. However, they clearly differ in the length of the vein ScP (short in the new fossil vs. long in Roemerula), the size of the radial field (encompassing about a third of anterior wing margin and distinct from wing apex vs. encompassing wing apex), the patterns of MA and MP (both with two branches vs. MA with three? branches and MP simple), and a different configuration of the cubital field. [START_REF] Cui | New and little-known species of Geinitziidae (Insecta: Grylloblattida) from the Middle Jurassic of China, with notes on taxonomy, habitus and habitat of these insects[END_REF] indicated that 'the combination of a branched RA and the presence of only three branches of the CuA in the forewing is characteristic of Shurabia'. Therefore, the new fossil is a representative of the genus Shurabia because its CuA has three branches.

Therefore, it appears that the new fossil is mostly similar to the genus Shurabia because its vein ScP is relatively short, the costal space broad, the RA has three branches, the RP five branches, MA and MP have both two branches, and CuA has three branches. The delineation of genera within the family Geinitziidae seems to be mainly based on the configuration of the RA, RP, CuA (i.e., number of branches), and the position of the forks of main branches. These characters are extremely plastic in Polyneoptera and have to be deeply investigated before to ensure their robustness and their utilization in genus diagnoses. The plasticity of the wing venation of Polyneoptera could be problematic and has to be investigated (see Plecoptera [START_REF] Béthoux | Wing venation pattern of Plecoptera (Insecta: Neoptera)[END_REF][START_REF] Béthoux | At last, a Pennsylvanian stemstonefly (Plecoptera) discovered[END_REF]. It has notably been demonstrated that the position of the connection between MP and CuA is highly variable in species of the order Plecoptera, and it is likely that similar plasticity existed for the 'Grylloblattodea' or the Reculida. With some restrictions linked to this potential variability, but based one the comparison provided before, we place the new fossil in the genus Shurabia, which is one of the most speciose genera of the family. than that of the posterior one (vs. nearly aligned), the fork of MP is distad that of MA (vs. anteriad in S. izyumica), the forking of RP not anteriorly pectinate (vs. anteriorly pectinate).

The new species differs from S. kapokkraalensis because all its costal veinlets are simple (vs. sometimes forked), its RA is simpler (three vs. six? branches in S. kapokkraalensis), its RP with five branches (vs. three), the fork of MA anteriad the fork of MP (vs. distad MP), MP with two branches (vs. three) [START_REF] Wappler | New orthopteroid insects from the Upper Triassic Molteno Formation, Karoo Basin, southern Africa[END_REF]. Affinities with S. lukashevichae are excluded because the new species has a RP with five branches (vs. three), the fork of MA slightly anteriad the fork of MP (vs. distad), MA with four branches (vs. apparently with fewer branches) [START_REF] Aristov | New and little known Grylloblattida (Insecta) from intertrappean deposits of the Tunguska Basin of Siberia[END_REF].

The new species differs from S. magna because of its smaller size (forewing length 12-13 mm vs. 22-27 mm in S. magna), its RP has five branches (vs. at most four), and the fork of MA is anteriad fork of MP (vs. distad) (Rasnitsyn, 1982: fig. 1f,g). The new species differs from S. minuta because its costal field is broad (vs. thin in S. minuta), its ScP is comparatively short and s-shaped (vs. long and straight), its RP has five branches (vs. four), the fork of MA anteriad the fork of MP (vs. far distad), and its MA has four branches (vs. three) (Rasnitsyn, 1982: fig. 3a). The new species differs from S. minutissima, at least, because of its larger size (ca. 12-13 mm vs. 5 mm in S. minutissima), its broad costal area (vs. thin), its RP with five branches (vs. two?), and its MA with four branches (vs. apparently simple) [START_REF] Aristov | New Geinitziidae (Insecta, Reculida) from the Permian and Triassic of Eurasia[END_REF].

The ovoid wings of S. ovata preclude affinities with the new species [START_REF] Martynov | Liassic insects from Shurab and Kisyl-Kiya, Part I, Various orders e x c e p t B l a t t o d e a a n d C o l e o p t e r a . A k a d e m i y a N a u k S S S R[END_REF] fig. 34). The new species further differs from S. ovata because its ScP is comparatively short (vs. long in S. ovata), its RA has three branches (vs. two), its RP five (vs. three), the fork of R and the fork of M nearly aligned (vs. fork of R far distad fork of M), its MA with four branches (vs. two), and the fork of MA slightly anteriad the fork of MP (vs. fork of MA distad fork of MP) (Martynov, 1937: fig. 34). The new species differs from S. parvula, at least, because its RP has five branches (vs. three in S. parvula), and its MA has four branches (vs. two) (Rasnitsyn, 1982: fig. 2a,d).

The new species differs from S. permiana because of it short and s-shaped ScP (vs. straight and long in S. permiana), RA with three branches (vs. often more), its fork of R and M nearly aligned (vs. fork of M strongly anteriad fork of R), its MA with four branches (vs. at most three), its fork of MA anteriad fork of MP (vs. distad fork of MP) (Aristov, 2013: fig. 73). The species S. postiretis possesses a thin costal space (vs. broad in the new species), a pectinate RP directed toward the anterior margin of the wing (vs. not pectinate), and a fork of MA distad the fork of MP (vs. slightly anteriad) (Huang et al., 1991: fig. 4). The new species differs from S. serrata because its RA has three branches (vs. two), its RP is not pectinate (vs. pectinate), its MA with a different branching pattern, and more space between main branches (Aristov et al., 2009: fig. 1f). The apex of the veins RA and RP of the holotype specimen of S. sogutensis are not preserved but the paratype specimen provides important information on this part of the wing (Rasnitsyn, 1982: fig. 1h,i). The RP of S. sogutensis have three branches (vs. five in the new species), MA seems to possess at most three branches (vs. four), and the fork of MA is distad the fork of MP (vs. anteriad the fork of MP) (Rasnitsyn, 1982: fig. 1h,i).

Finally, the new species differs from the species S. tanga because it lacks the distinctive feature of the latter species i.e., a long s-shaped crossvein in the interradial field (Aristov, 2018: fig. 1d). Both species further differ because the MA is directed toward posterior wing margin in S. tanga while directed toward the anterior margin in the new species.

The temporal ranges of these species are also potential arguments to preclude affinities between each other. For example, it is expected that Cretaceous species or Jurassic species would greatly differ from the Permian ones, an argument already proposed by [START_REF] Rasnitsyn | Triassic and Jurassic Insects of the genus Shurabia (Grylloblattida, Geinitziidae)[END_REF]. This assumption stems from the longevity of fossil species that is often a few million years (sometimes more), and suggesting a 'cyclic decline' of species through time. Additionally, even if the similarity of the wing venation is the main criterion to discriminate Shurabia species, it is also expected that the other parts of the body may have differed between species (e.g., genitalia). Therefore, superficial wing venation similarities cannot be considered unalterable arguments to group specimens from different geographical areas or temporal periods. Biogeographical constraints such as insularism, or the presence of mountains may result in allopatric speciation not necessarily detectible on from wing venation alone.

Conclusion

The description of Shurabia taewani sp. nov. serves to increase the diversity of Geinitziidae documented from Asia during the Triassic period and adds to the morphological diversity of the genus Shurabia of this period. Similarly, it extends the distribution of the genus already known from the Triassic of Australia, Japan, Kyrgyzstan, South Africa and Ukraine.

  Illustrator and Adobe Photoshop CC2019. The specimen of Shurabia taewani sp. nov. is housed in the Gongju National University of Education under the collection number GNUE112003. Venation abbreviations are as follows: C costa; CuA cubitus anterior; CuA1 first branch of cubitus anterior; CuA2 second branch of cubitus anterior; MA media anterior; MP media posterior; m-cua or M5 (supposedly the fifth branch of the media), also interpreted as a strong crossvein between M and R; RA radius anterior; RP radius posterior; ScP subcosta posterior. urn:lsid:zoobank.org:pub: xxxx

Figure 1 .

 1 Figure 1. Shurabia taewani sp. nov., holotype GNUE112003. Photographs with different

Figure 1 .

 1 Figure 1. Shurabia taewani sp. nov., holotype GNUE112003. Photographies avec différentes

Figure 2 .

 2 Figure 2. Shurabia taewani sp. nov., holotype GNUE112003. Interpretative line drawing of

Figure 2 .

 2 Figure 2. Shurabia taewani sp. nov., holotype GNUE112003. Dessin interprétatif de la

Acknowledgment

We are grateful to Prof. Alexandr P. Rasnitsyn and an anonymous reviewer for their insightful comments on the manuscript. We also thank the editor Prof. Didier Néraudeau for managing our article. This work is part of the Ph.D. project of CJ.

Differences with other species of the genus Shurabia

Currently, the genus Shurabia encompasses the following species: S. aleda, S. angustata, S. anomala, S. australis, S. bashkuevi, S. creta, S. elegans, S. ferganensis, S. grandis, S. hissarica, S. inferior, S. izyumica, S. kapokkraalensis, S. lukashevichae, S. magna, S. minuta, S. minutissima, S. ovata, S. parvula, S. permiana, S. postiretis, S. serrata, S. sogutensis, and S. tanga. The new species differs from S. aleda because of its MA with at least four branches (vs. two in S. aleda), and its CuA1 forked (vs. simple) [START_REF] Wappler | New orthopteroid insects from the Upper Triassic Molteno Formation, Karoo Basin, southern Africa[END_REF]. The species S. angustata differs from the new fossil because its RA has fewer branches (at most two vs. three in the new species), and a different branching pattern of RP (Rasnitsyn, 1982: fig. 1b-e).

The new specimen cannot be placed in S. anomala because it lacks the presence of numerous crossveins (vs. present in S. anomala), the completely different branching pattern of the medial field with MA possessing four branches (vs. apparently simple), MP with two branches (vs. with at least three branches), and the space between CuA1 and CuA2 broad (vs. extremely narrow) (Rasnitsyn, 1982: fig. 4e).

The new specimen does not fit inside S. australis because of its wide costal area (vs. thin in S. australis), its ScP s-shaped (vs. straight), its RA with three branches (vs. two), RP with five branches (vs. four), MA apparently with four branches (vs. apparently three) (Rasnitsyn, 1982: fig. 4f). Additionally, the space between main branches appears to be larger in the fossil than in S. australis, notably between the first abscissa of RP and MA, and between the first abscissa of CuA1 and CuA2. The new specimen clearly differs from S. bashkuevi because of its RP with five branches (vs. two in S. bashkuevi), MA with four branches (vs. two), and the wide space between the posterior branch of CuA1 and CuA2 (Aristov, 2011: fig. 2b). The species S. creta is the only species of the genus known from the Cretaceous period. The preservation of the type material is not perfect, but this species differs from our new species because its interradial field is not widened (vs. widened in Shurabia taewani sp. nov.), RP with two(?) (vs. five), MA pectinate backward (vs. encompassing wing apex) (Aristov, 2020b). The temporal range of the two species (Cretaceous vs. Late Triassic) is another argument supporting their placement into two different species. The new specimen differs from S. elegans because of its short ScP (vs. long in S. elegans), its RA with three branches (vs. four), MP with two branches (vs. three), and the forks of R and M nearly aligned (vs. fork of R clearly distad fork M) (Fujiyama, 1973: fig. 16).

Shurabia ferganensis has a long ScP (vs. short in Shurabia taewani sp. nov.), a comparatively thin costal space (vs. broad) the fork of R clearly distad fork M (vs. nearly aligned), MA branches mostly directed toward the posterior wing margin (vs. anterior wing margin and apex) [START_REF] Rasnitsyn | Triassic and Jurassic Insects of the genus Shurabia (Grylloblattida, Geinitziidae)[END_REF] fig. 4a-d). The new species differs from S. grandis because of its ScP short (vs. long), its RP with five branches (vs. three?), and the fork of R nearly aligned with the fork of M (vs. fork of R clearly distad the fork of M) [START_REF] Huang | New "Grylloblattida" related to the genus Prosepididontus Handlirsch, 1920 in the Middle Jurassic of China (Insecta: Geinitziidae)[END_REF] fig. 3). The new species differs from S. hissarica because of its RP with five branches (vs. three in S. hissarica), its MA with four branches (vs. two) [START_REF] Aristov | New and little-known grylloblattids of the family Geinitziidae (Insecta: Grylloblattida) from the Triassic and Jurassic of Europe, Asia, and South Africa[END_REF]. The new species differs from S. inferior because its RA has three branches (vs. two in S. inferior), its RP with five branches (vs. three and of different shape), and its MA with four branches (vs. simple) (Aristov, 2015b: fig. 3d).

The new species differs from S. izyumica by its smaller size (forewing 12-13 mm vs.

Disclosure statement

The authors declare that they have no known competing interests or personal relationships that could have appeared to influence the work reported in this paper.


Shabina Ashraf

Yves Méheust

Jyoti Phirani

Spontaneous imbibition dynamics in two-dimensional porous media: a generalized interacting multi-capillary model

The capillary bundle model, wherein the flow dynamics of a porous medium is predicted from that of a bundle of independent cylindrical tubes/capillaries whose radii are distributed according to the medium's pore size distribution, has been used extensively. But, as it lacks interaction between the flow channels, this model fails at predicting complex flow configuration, including those involving two-phase flow. We propose here to predict spontaneous imbibition in quasi-two-dimensional (quasi-2D) porous media from a model based on a planar bundle of interacting capillaries. The imbibition flow dynamics, and in particular, the breakthrough time, the global wetting fluid saturation at breakthrough, and which capillary carries the leading meniscus, are governed by the distribution of the capillaries' radii and their spatial arrangement. For an interacting capillary system consisting of 20 capillaries, the breakthrough time can be 39% smaller than that predicted by the classic, non-interacting, capillary bundle model of identical capillary radii distribution, depending on the spatial arrangement of the capillaries. We propose a stochastic approach to use this model of interacting capillaries for quantitative predictions. Comparing bundles of interacting capillaries with the same capillary diameter distribution as that of the pore sizes in the target porous medium, and computing the average behavior of a randomly-chosen samples of such interacting capillary bundles with different spatial arrangements, we obtain predictions of the position in time of the bulk saturating front, and of that of the leading visible leading front, that agree well with measurements taken from the literature. This semi-analytical model is very quick to run and could be useful to provide fast predictions on one-dimensional spontaneous imbibition in porous media whose porosity structure can reasonably be considered two-dimensional, e.g., paper, thin porous media in general, or layered aquifers. 49 gle pore/tube was observed to be strongly dependent on the 50 geometries of the capillaries 29-40 . 51 Due to the similarity in the macroscopic laws describing 52 the time evolution of the imbibition length between imbibition 53 in a capillary tube and imbibition in a homogeneous porous 54 medium, the capillary bundle model, considering a bundle of 55 non-interacting capillaries of different radii, is classically con-56

 . In 59 a quasi-two-dimensional (2D) porous medium such as paper, 60 Bico and Quéré 47 showed that there are two imbibing fronts, 61 a leading front in the small pores and a bulk saturating front 62 which lags behind, which is contradictory to the predictions 63 of the classic bundle of (non-interacting) capillaries, where 64 the pores with larger radii have the leading front during imbi-65 bition. 66 The model geometry consisting of interacting capillaries 67 (i.e., a capillary bundle where an opening allowing fluid ex-68 change exists between adjacent capillairies, see e.g. Ref.55) 69 accounts for the effect of the interaction between pores on 70 the pore scale flow dynamics, which in turn affects the Darcy 71 scale flows in porous media 48-55 . In a system of two inter-72 acting capillaries, the imbibition in the capillary of smaller

73

radius is found to be faster than that in the one of larger ra-74 dius, unlike the behavior suggested by Eq. (1). However, a 75 majority of these models were limited to predicting the im-76 bibition dynamics in an ordered arrangement of pores or in 77 two and three interacting capillary systems. For a system con-78 sisting of three interacting non-cylindrical capillaries, Unsal

I. INTRODUCTION

When a wetting fluid is placed in contact with a porous medium, the fluid spontaneously imbibes into the pore space due to capillary suction. Such spontaneous imbibition in the porous matrix is crucial for applications such as oil recovery from reservoirs 1-3 , Paper Analytic Devices (µPADs) [START_REF] Taghizadeh-Behbahani | A paper-based length of stain analytical device for naked eye 1055 (readout-free) detection of cystic fibrosis[END_REF][START_REF] Soda | Equipment-free detection of k+ on mi-1057 crofluidic paper-based analytical devices based on exhaustive replacement 1058 with ionic dye in ion-selective capillary sensors[END_REF] , textiles [START_REF] Dai | Bioinspired janus textile with conical micropores for human body mois-1062 ture and thermal management[END_REF] , inkjet printing [START_REF] Rosello | Dripping-onto-substrate 1064 capillary breakup extensional rheometry of low-viscosity printing inks[END_REF][START_REF] Wang | Fabrication of monolayers of 1067 uniform polymeric particles by inkjet printing of monodisperse emulsions 1068 produced by microfluidics[END_REF] , microfluidics [START_REF] Liu | Microfluidic investigations of crude oil-1070 brine interface elasticity modifications via brine chemistry to enhance oil 1071 recovery[END_REF][START_REF] Gharibshahi | Hybridiza-1073 tion of superparamagnetic fe3o4 nanoparticles with mwcnts and effect of 1074 surface modification on electromagnetic heating process efficiency: A mi-1075 crofluidics enhanced oil recovery study[END_REF][START_REF] Carrell | Beyond the lateral flow assay: A review of paper-1078 based microfluidics[END_REF][START_REF] Schaumburg | Assessing the rapid flow in multilayer 1080 paper-based microfluidic devices[END_REF][START_REF] Rich | Characterization 1083 of glass frit capillary pumps for microfluidic devices[END_REF] , lab-on-chip devices [START_REF] Lin | Performance analysis of 1086 a proton-exchange membrane fuel cell (pemfc) with various hydrophobic 1087 agents in a gas diffusion layer[END_REF][START_REF] Lee | A whole blood sample-to-answer 1089 polymer lab-on-a-chip with superhydrophilic surface toward point-of-care 1090 technology[END_REF] , point-of-care diagnostics [START_REF] Liang | Smartphone-app 1093 based point-of-care testing for myocardial infarction biomarker ctni using 1094 an autonomous capillary microfluidic chip with self-aligned on-chip focus-1095 ing (sof) lenses[END_REF][START_REF] Joung | based multiplexed 1098 vertical flow assay for point-of-care testing[END_REF] , Polymer Electrolyte Membrane Fuel Cell (PEMFC) [START_REF] Xiao | A 1101 novel fractal model for relative permeability of gas diffusion layer in proton 1102 exchange membrane fuel cell with capillary pressure effect[END_REF][START_REF] Carrere | Liquid water in cathode gas diffusion layers of pem 1105 fuel cells: Identification of various pore filling regimes from pore network 1106 simulations[END_REF] , micro heat pipes [START_REF] Singh | Enhanced ther-1109 mal performance of micro heat pipes through optimization of wettability 1110 gradient[END_REF][START_REF] Chernysheva | Simulation of heat and mass transfer in 1112 a cylindrical evaporator of a loop heat pipe[END_REF] , in understanding the motion of blood cells [START_REF] Pozrikidis | Axisymmetric motion of a file of red blood cells through 1115 capillaries[END_REF] and in the design of bio-inspired drainage and ventilation systems [START_REF] Singh | The architectural design of smart ventila-1118 tion and drainage systems in termite nests[END_REF] . Capillary driven imbibition in a homogeneous porous medium follows diffusive dynamics, where the imbibition length is proportional to the square root of time [START_REF] Li | Criteria for applying the 1121 lucas-washburn law[END_REF][START_REF] Gruener | Capillarity-driven oil flow in nanopores: Darcy 1123 scale analysis of lucas-washburn imbibition dynamics[END_REF][START_REF] Cai | Capillary imbibition and flow 1126 of wetting liquid in irregular capillaries: A 100-year review[END_REF] . This kind of dynamics was first characterized by Lucas [START_REF] Lucas | Ueber das zeitgesetz des kapillaren aufstiegs von flüssigkeiten[END_REF] and Washburn [START_REF] Washburn | The dynamics of capillary flow[END_REF] for a horizontal cylindrical capillary tube: during the spontaneous imbibition of a wetting fluid of viscosity µ in a tube of radius r, the imbibition length (which here is simply the longitudinal position of the meniscus along the tube) is given by

l = rσ cos θ w 2µ t, (1)
where σ is the surface tension coefficient and θ w is the wetting angle of the invading fluid on the tube's wall. In Eq. (1), the prefactor of the √ t law is proportional to √ r, which implies that at any given time the meniscus will have advanced a) Electronic mail: jyoti.phirani@strath.ac.uk more along a capillary of larger radius than along one of 48 smaller radius. Later, the phenomenon of imbibition in a sin-et al. [START_REF] Unsal | Co-current and counter-1214 current imbibition in independent tubes of non-axisymmetric geometry[END_REF][START_REF] Unsal | Co-and counter-current 1217 spontaneous imbibition into groups of capillary tubes with lateral connec-1218 tions permitting cross-flow[END_REF][START_REF] Unsal | Bubble snap-off and 1221 capillary-back pressure during counter-current spontaneous imbibition into 1222 model pores[END_REF] showed experimentally that the imbibition speed is fastest in the capillary of least effective radius. On the contrary, Ashraf et al. [START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF] , using a one-dimensional lubrication approximation model and considering a system of three interacting cylindrical capillaries, showed that imbibition is not always fastest in the capillary of smallest radius. Furthermore, both these studies [START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF][START_REF] Unsal | Co-current and counter-1214 current imbibition in independent tubes of non-axisymmetric geometry[END_REF] showed that, for three capillary systems, the random positioning of the capillaries strongly impacts the invasion behaviour. But how the interconnection between capillaries impacts the overall imbibition dynamics is far from being fully understood in the general case of a larger number of tubes. Consequently, interacting capillary systems, despite having a complexity which is intermediate between that of the classical bundle of non-interacting capillaries, have so far not been used to predict the generalized imbibition phenomenon observed in porous media consisting of several pores of irregular sizes and varying connectivity. To this aim, more complex models have been introduced since, based on pore-network geometries inferred from a geometrical analysis of the porous medium in which imbibition is to be investigated [START_REF] Bultreys | Verifying pore network mod-1225 els of imbibition in rocks using time-resolved synchrotron imaging[END_REF][START_REF] Bultreys | Improving the description of two-phase flow in rocks by in-1229 tegrating pore scale models and experiments[END_REF][START_REF] Foroughi | Pore-by-pore modelling, vali-1232 dation and prediction of waterflooding in oil-wet rocks using dynamic syn-1233 chrotron data[END_REF] . We will present here a model of intermediate complexity between those early interacting-capillary models and pore network models. Note that in many practical cases, the detailed porous structure is not known, and only an estimate of the pore size distribution is available; in such cases a pore network model cannot be applied without making assumptions on the unknown structure, whereas the model presented here can be applied directly.

We thus propose a generalized one-dimensional model to predict spontaneous imbibition in a capillary bundle consisting of any number of randomly arranged cylindrical tubes that interact with each other, with any arbitrary distribution of the capillaries' radii. The model generalizes the study by Ashraf et al., [START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF] for systems of two and three interacting capillaries, to an arbitrary number of interacting capillaries. It is meant to model spontaneous imbibition in quasi-2D porous media for which the pore size distribution is known. The model is inspired from a model developed to tackle spontaneous imbibition in stratified geological porous media [START_REF] Ashraf | A generalized model for spontaneous imbibition 1235 in a horizontal, multi-layered porous medium[END_REF] . The two models are formally very similar to each other, but, due to the difference in geometries (flat layers for the stratified geological formation, cylindrical tubes in the present model), the equations are not identical. More importantly, the two studies differ widely in that the relative positioning of the layers in a geological medium is given, whereas, for a quasi-2D porous medium whose pore size distribution is known, the relative positioning of connected capillaries of different diameters within the 2D bundle that can predict the medium's behavior is not known a priori. Here, we explain the underlying physical phenomena causing the menisci to advance at different rates in the different capillaries, and demonstrate that both the spatial arrangement of the interacting capillaries, and, for a given arrangement, the contrasts in the capillaries' radii (i.e., their ratios), are crucial in predicting the imbibition dynamics. In contrast to the standard (non-interacting) capillary bundle, this model provides predictions that are qualitatively consistent with the phenomenology of spontaneous imbibition in real (quasi-)two-dimensional (2D) porous media. In particular, this model correctly predicts that the smaller pores In this study we shall generalize the reduced order model of Ashraf et al., [START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF] to an arbitrary number of capillaries positioned in the same plane and interacting with their neighbours. For a flat bundle of three interacting capillaries, the model of Ashraf et al., [START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF] showed that the distribution of radii and the spatial arrangement of the capillaries impact the imbibition behavior in the capillary system significantly. The meniscus in the capillary of smallest radius does not always move ahead of the others.

In the following sections, we examine the dynamics of menisci during spontaneous imbibition in a flat bundle containing an arbitrary number of interacting capillaries. This generalization of the interacting capillaries' model follows the model development formulations from the study of Ashraf et al., [START_REF] Ashraf | A generalized model for spontaneous imbibition 1235 in a horizontal, multi-layered porous medium[END_REF] for imbibition in stratified porous media. In a stratified porous medium, the contrasts in layer transmissivities and the relative positioning of the layers control the imbibition dynamics, whereas in the present interacting capillaries bundle model, the positioning of the capillaries also plays a crucial role, but the role played by the transmissivities in the stratified medium is played by the product of the capillaries' permeabilities by their cross-sectional areas, both of which are controlled by the contrasts in the capillaries' radii.

We first describe below the one-dimensional model formulation for a system of four interacting capillaries to understand the underlying equations, before generalizing the model to a multiple-interacting capillary system.

B. Model development for four interacting capillaries

To predict the dynamics of spontaneous imbibition in a porous medium using a system of interacting capillaries, we need to take the arrangement of capillaries into account, unlike for the classic capillary bundle (sometimes called bundleof-tubes) model. For a porous medium made of n interacting capillaries, there are n! /2 different arrangements. Fig. 2 shows a bundle of four interacting capillaries that are ordered spatially according to their radii r α > r β > r γ > r δ ; we call this arrangement αβ γδ . The capillary pressure in tube i

(i = α, β , γ, δ
) is given by the Young-Laplace equation as [START_REF] Young | III. an essay on the cohesion of fluids[END_REF][START_REF] De Laplace | Supplément au dixième livre du Traité de mécanique 1243 céleste[END_REF] ,

Pc i = 2σ cos θ w r i , ( 2
)
where σ is the surface tension and θ w the contact angle; hence, Pc α < Pc β < Pc γ < Pc δ . The corresponding imbibition lengths in the tubes at any time t are denoted respectively α β γ δ FIG. 2. Schematic showing the spontaneous imbibition in an ordered system of four interacting capillaries. The imbibition lengths in capillaries α, β , γ, δ of radii r α , r β , r γ , r δ are denoted by z α , z β , z γ , z δ , respectively. The cross section of the system of interacting capillaries is also shown.

by z i (t). We consider the assumptions from Ashraf et al. [START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF] 3. Spontaneous imbibition in a system of four interacting capillaries with a spatial arrangement of β γαδ of the capillaries. The imbibition lengths in capillaries α, β , γ, δ of radii r α , r β , r γ , r δ are z α (t), z β (t), z γ (t), z δ (t), respectively. The schematics of the imbibition phenomenon show the fluid transfer at menisci locations with arrows. Fot this spatial arrangement, depending upon the contrasts in the capillaries' radii, the possible orders in the invasion lengths can be (a)

z α < z β < z γ < z δ , (b) z α < z β < z δ < z γ and (c) z α < z δ < z β < z γ .
The cross section of the system of interacting capillaries is also shown for (a).

z δ (t) relative to z β (t) and z γ (t) depends on the contrasts in the 283 capillaries' radii.

284

The detailed development of the generalized one-285 dimensional model for this system of four interacting capil-286 laries with arrangement β γαδ is described in Appendix A.

287

The pressure drop across each of the sections is determined 288 individually, i.e., for sections (I) 0

< z < z α , (II) z α < z < z β , 289 (III) z β < z < z γ , and (IV) z α < z < z δ . As spontaneous im- 290
bibition is driven by capillary forces, the sum of the pressure 291 drops across all the sections of a capillary is equal to the cap-292 illary pressure of that capillary.

293

Pc i = ∑ j P i ( j) , (3)
where P i ( j) is the pressure drop across the section of index

Pc α = 8µz α (t) r 4 α + r 4 β + r 4 γ + r 4 δ r 2 α dz α dt + r 2 β dz β dt + r 2 γ dz γ dt + r 2 δ dz δ dt , ( 4
) 301 Pc δ -Pc α = 8µ(z δ (t) -z α (t)) r 2 δ dz δ dt , ( 5
) 302 Pc β -Pc α = 8µ(z β (t) -z α (t)) r 4 β + r 4 γ r 2 β dz β dt + r 2 γ dz γ dt , (6)
303

Pc γ -Pc β = 8µ(z γ (t) -z β (t)) r 2 γ dz γ dt . (7)
Eqs. ( 4) to ( 7) are rendered non-dimensional by normaliz-

Z i = z i L , i = α, β , γ, δ and T = Pc α r 2 α 8µL 2 t . ( 8
)
Introducing the contrasts in radii, λ i = r i /r α , and in capillary 308 pressures, ε i = Pc i /Pc α , for i = β , γ, δ , we then obtain the 309 non-dimensional form of Eqs. ( 4) to (7) as

310 1 = Z α 1 + λ 4 β + λ 4 γ + λ 4 δ dZ α dT + λ 2 β dZ β dT + λ 2 γ dZ γ dT + λ 2 δ dZ δ dT , (9)
311 ε δ -1 = Z δ -Z α λ 2 δ dZ δ dT , ( 10
) 312 ε β -1 = Z β -Z α λ 4 β + λ 4 γ λ 2 β dZ β dT + λ 2 γ dZ γ dT . ( 11
) 313 ε γ -ε β = Z γ -Z β λ 2 γ dZ γ dT , (12)
Further assuming that the contact angle θ w is the same in all 314 capillaries, we have ε i = 1/λ i , and upon rearranging the gov-315 erning Eqs. ( 9) to ( 12) and adding them, we obtain,

316 2 1 + ∑ i=β ,γ,δ ε i λ 4 i T = Z 2 α + Z 2 β λ 2 β + Z 2 γ λ 2 γ + Z 2 δ λ 2 δ . (13)
Eq. ( 13) expresses that, in a system of interacting capillaries, C. Generalizing the one-dimensional spontaneous imbibition 329 model in the interacting capillary system 330 Equation ( 13) is readily generalized to a system of n inter-331 acting capillaries, in the form

332 2 n ∑ i=1 ε i λ 4 i T = n ∑ i=1 ψ i Z i (14)
where

ψ i = πr 2 i z i /(πr 2 α L) ( j = 1, 2, ..., n) is the non- 333
dimensional volume imbibed in the capillary of index i.

334

Eq. ( 14) expresses that the sum over all capillaries of 335 the non-dimensional volumes times the corresponding non-336 dimensional imbibition lengths, is proportional to time. This 337 can be compared to the dynamics in a bundle of non-338 interacting capillaries, for which we know that the dynam-339 ics are diffusive, i.e., for each of the capillaries, the imbibed 340 length square is proportional to time.

341

We note from the derivation of Eq. ( 13) for the system con- First, the algorithm searches for the capillary of largest ra-358 dius in the arbitrary arrangement, whose meniscus position is 359 z a at a given time; it is denoted C i in Fig. 4, where the capil-360 laries in the order of arrangement are denoted from C 1 to C n .

361

The pressure drop in the region 0 < z < z a is determined for all 362 the capillaries and the algorithm then considers two regions: We first explore the imbibition of a system of four interact-381 ing capillaries, followed by the imbibition in a system consist-382 ing of 20 capillaries.

592

The capillary radii are identical in the two cases. For non-593 interacting capillaries, by non dimensionalizing the Wash-594 burn's law, z 2 i = (Pc i r 2 i /4µ)t, we obtain 595

Z 2 i = 2ε i λ 2 i T, (15)
where Z i = z i /L is the non-dimensional length imbibed in 596 the capillary of radius r i and L is the total length of the For instance, the leading meniscus for an orderly arranged interacting capillary system is in the smallest radius capillary and we know that the fraction of saturation contributed by the smallest radius capillary is small. For the arrangement 2 shown in Fig. 9, the leading meniscus is in the capillary of radius 100 m. In the capillary bundle model, the crosssection area of the leading capillary (200m) is 13.93% of the total cross-section area, whereas for the ordered arrangement and the arrangement number 2, the respective cross-section area of the leading meniscus capillaries are 0.03% and 3.43%.

= 2Pc eff K µ t , ( 16
)
where K is the medium's permeability and we have assumed 760 that at time t = 0 no wetting fluid has yet invaded the medium.

761

If we choose to non-dimensionalize time by the character-762 iztic time (µL 2 )/(Pc eff K), we obtain from Eq. ( 16) the non-

P(z β (t),t)-P(z α (t),t) = - 8µ(z β (t) -z α (t)) r 2 β v β (t) + ω dq α A β . (A7)
At the meniscus in the capillary β , the capillary pressure jump Similarly, for z > z β (t), the meniscus in the capillary γ trav-922 els with a velocity given by

923 dz γ dt = v γ (t) + (1 -ω) dq α πr 2 γ + dq β πr 2 γ . (A9)
The pressure drop between z = z β (t) and z = z γ (t) in capillary γ is then given by,

P(z γ (t),t) -P(z β (t),t) = - 8µ(z γ (t) -z β (t)) r 2 γ v γ (t) + (1 -ω) dq α πr 2 γ + dq β πr 2 γ . ( A10
)
The pressure jump across the meniscus in each of the 924 capillaries is given by the Young-Laplace equation [START_REF] Young | III. an essay on the cohesion of fluids[END_REF][START_REF] De Laplace | Supplément au dixième livre du Traité de mécanique 1243 céleste[END_REF] , i.e.,

925

Eq. ( 2), from which it follows that 926 P(z i ,t)

-P 0 = -Pc i = - 2σ cos θ w r i , (

138

  carry the leading front, while the larger pores carry the lag-139 ging saturating front responsible for the mass uptake of fluid 140 in the porous medium, as measured in a paper-based porous 141 medium[START_REF] Bico | Precursors of impregnation[END_REF] . Furthermore, we provide a successful quantitative 142 comparison between the measurements of Bico & Quéré on 143 the leading and lagging imbibition fronts to predictions of the 144 model obtained using a stochastic approach: the predicted be-145 havior is the average of those obtained for all possible spatial 146 organizations of the capillaries' diameter distribution. Though 147 less accurate than fully numerical (and much more compli-148 cated) pore network models, this semi-analytical model has 149 the advantage of running within seconds on any computer.150The presentation is organized as follows. We first review 151 the model by Ashraf et al.,[START_REF] Ashraf | Spontaneous imbibition in ran-1211 domly arranged interacting capillaries[END_REF] (section II A). We then pro-152 ceed to extend it to a system consisting of 4 interacting cap-153 illaries (section II B), before presenting the generalized one-154 dimensional model predicting spontaneous imbibition in an 155 interacting multi-capillary system (section II C). We then ex-156 amine the imbibition dynamics in a system of four interacting 157 capillaries (section III A) and in a similar system consisting of 158 20 capillaries (section III B). In the discussion, we first com-159 pare the predictions of our model to those of the classic, non-160 interacting, capillary bundle (section III C 1), and, finally, con-161 front its predictions of the leading and lagging fronts in a real 162 quasi-2D porous medium from the literature to the published 163 experimental measurements (section III C 2). Section IV con-164 tains a summary of the work and conclusive remarks, and dis-165 cusses prospects to this study. 166 II. MODELS 167 A. Capillary imbibition in interacting capillaries 168 Using the capillary system shown in Fig. 1, Ashraf et al., 55 169 used volume of fluid 63 (VOF) two-phase flow simulations to 170 study spontaneous imbibition in a bundle of two or three inter-171 acting capillaries. These CFD (computational fluid dynamics) 172 calculations provided the entire pressure and velocity fields 173 inside the connected capillaries. They showed that (1) the in-174 vading wetting fluid transfers between two adjacent capillar-175 ies from the capillary of larger radius to that of smaller radius, 176 but this transfer occurs only in the immediate vicinity of the 177 (less advanced) meniscus of the capillary of larger radius; (2) 178 that everywhere else (that is, everywhere except in the vicinity 179 of that meniscus), the flow in the capillaries is not perturbed 180 by the transfer of fluid between the capillaries; and (3) that, 181 consequently, the pressure can be considered uniform over all 182 transverse sections of the capillary system where both capil-183 laries are filled with the same fluid, since no flow occurs along 184 the transverse direction (if one neglects the small regions in 185 the vicinity of the less advanced meniscus). These findings 186 (1-3) served as basic assumptions to develop a reduced or-187 der, Washburn-like one dimensional model for a bundle of two 188 and three interacting capillaries that can interact hydrodynam-189 ically with the neighbouring capillaries along their touching 190 sides. The model predicted that in a bundle of two interact-ing capillaries the meniscus in the capillary of smaller radius moves ahead of the other one during the spontaneous imbibition, in consistency with the results of the VOF simulations.

FIG. 1 .

 1 FIG. 1. Spontaneous imbibition in two interacting capillaries, (a) cross-sectional view, (b) lateral view showing the contact angle θ w .

FIG.

  FIG.3. Spontaneous imbibition in a system of four interacting capillaries with a spatial arrangement of β γαδ of the capillaries. The imbibition lengths in capillaries α, β , γ, δ of radii r α , r β , r γ , r δ are z α (t), z β (t), z γ (t), z δ (t), respectively. The schematics of the imbibition phenomenon show the fluid transfer at menisci locations with arrows. Fot this spatial arrangement, depending upon the contrasts in the capillaries' radii, the possible orders in the invasion lengths can be (a)z α < z β < z γ < z δ , (b) z α < z β < z δ < z γ and (c) z α < z δ < z β < z γ .The cross section of the system of interacting capillaries is also shown for (a).

294j

  = (I), (II), (III), (IV ) of the capillary of index i = α, β , γ, δ . 295 By solving the system of equations expressing (i) Darcy's law 296 in each of the capillaries, and (ii) the relations between the 297 meniscii's advancement and the fluid velocities and fluid ex-298 change between the capillaries, we obtain the equations gov-299 erning the flow in the interacting capillaries, which are, 300

304

  ing the positions by the total capillary system's length, L, and 305 time by [8µL 2 /(Pc α r 2 α )], thus defining the non-dimensional 306 positions and times 307

317

  the sum of the squares of the product of the non-dimensional 318 radius with the non-dimensional distance invaded in all the 319 capillaries is proportional to the invasion time T . For dif-320 ferent arrangements of a system of 4 interacting capillaries 321 having the same contrasts in capillary radii, the total capil-322 lary suction of the system remains the same. Therefore, for 323 any of the 4! /2 = 12 possible arrangements, rearranging the 324 equations governing the imbibition process, and adding them, 325 leads to Eq. (13). However, the velocity at which the indi-326 vidual meniscii travels in each of the tubes depends on the 327 particular arrangement of the capillaries.

342FIG. 4 .

 4 FIG. 4. Schematic of spontaneous imbibition in an n-capillary system where the capillaries are positioned randomly. The capillaries in the arrangement are denoted by C 1 ,C 2 , • • •C n . The capillary radii are denoted as r a , r b , r c , • • •, and the corresponding imbibition distances at time t are denoted by z a (t), z b (t), z c (t), • • •.

383A.FIG. 8 .FIG. 9 .Fig. 9 (FIG. 10 .

 89910 FIG. 5. Spontaneous imbibition in a system of four interacting capillaries which are positioned with respect to each other according to the arrangement β γαδ , for three different contrasts in capillary radii. (a), (c), (e) represent the schematics of possible imbibition behavior at a given time t. The distribution of radii predicting the imbibtion phenomenon are indicated in the plots (b), (d) and (f). The non-dimensional times at which the leading meniscus reaches the outlet end of the interacting capillary system (T bt ) for the cases (b), (d) and (f) are 0.43, 0.40 and 0.39, respectively.

597

  FIG. 11. (a) Spatial saturation profile during spontaneous imbibition in a bundle-of-tubes consisting of twenty non-interacting capillaries at T = 0.1, T = 0.3, and T = T bt = 0.5. (b) Average spatial saturation profile for 1000 different spatial arrangements of the system consisting of twenty interacting capillaries of identical radii as in (a), at T = 0.1, T = 0.3, T = T bt .

Consequently, as shown

  by Fig.11(b), the cross-sectional saturation decreases gradually with longitudinal position for the classic capillary bundle model, while in the case of interacting capillaries a steep decrease is observed already at small longitudinal positions. Fig.11(b) also shows that the standard deviation in saturation from the average across the 1000 arrangements at T = 0.1 and T = 0.2, which is as high as 0.2 at Z = 0.59 and 0.69, respectively; whereas for T = T bt , it is 0.18 at Z = 0.76. In real two-dimensional porous media where the spatial arrangement of pores may vary, the interacting capillaries model will be more helpful in predicting the accurate imbibition behaviour than the classic capillary bundle model. The saturation of the porous medium with length and the breakthrough time significantly differ for the classic (noninteracting) capilary bundle and for the different arrangements of the interacting multi-capillary system, although the contrast in the radii of the capillaries is the same.

2 . 2 FIG. 13 .

 2213 FIG. 12. Imbibition in a system of two interacting capillaries having radii r α = 300 m and r β = 170 m. Predictions from our semianalytical model (solid lines) compare well to the data (symbols) of Bico and Quéré 47 . The inset of the figure shows the same comparison for predictions of the classic (non-interacting) capillary bundle model (solid lines), obtained through Eq. 15, which underestimates the advancement Z β of the meniscus in the narrower capillary (red line) by a factor 5.695 696

763 dimensional equation Z 2 =(

 2 2T . Since, according to Bico 764 and Quéré's observation mentioned above, it is the lagging 765 (macroscopic) front that carries most of the interface between 766 the two fluids, Eq. 16, and therefore its non-dimensional coun-767 terpart, can be assumed to describe the behavior of the lag-768 ging front. From the experimental data for the lagging front, 769 Pc eff K)/(µL 2 ) is measured to be 9.7 • 10 -5 s -1 , which we 770 thus use to non-dimensionalize all plots in Fig.13. The de-771 pendence of Z 2 on T for the lagging (macroscopic) front then 772 has a slope 2 (as shown by the dotted gray line in Fig.13),773while that for the leading (microscopic) front exhibits a larger 774 imbibition rate, with a slope 2.67 (as shown by the orange 775 dotted line in Fig.13).776IV. CONCLUSIONS777 In conclusion, we investigated spontaneous imbibition of a 778 wetting fluid in a randomly arranged planar system of inter-779 acting capillaries. This generalized model can predict the im-780 bibition behavior for all the n! /2 possible arrangements of an 781 interacting n-capillary system. It is inspired from a previous 782 work on stratified geological formations, with planar layers 783 instead of cylindrical capillaries. 784 Using an interacting capillary system containing 4 capillar-785 ies, we showed that the imbibition dynamics depends signifi-786 cantly on the arrangement of the capillaries within the capil-787 lary system, for a given distribution of the capillary radii. Sim-788 ilarly, the dynamics are affected by that distribution for a given 789 arrangement of the capillaries. Furthermore, we showed that 790 the arrangement and radii distribution of the capillaries jointly 791 control the relative menisci's locations, the breakthrough time, 792 and which capillary carries the leading meniscus. The cross-793 sectional saturation of the impregnating fluid along the length 794 of the capillary system also changes with a change in the ar-795 rangement of the capillaries. However, the total capillary pres-796 sure driving the flow is identical for all arrangements, there-797 fore, the overall volume fraction occupied by the invading 798 fluid (i.e, the global saturation of the wetting fluid) at a given 799 time remains the same across all arrangements, until break-800 through occurs in one of the arrangements. 801 Similarly, considering 1000 randomly-chosen different ar-802 rangements of an interacting twenty-capillary system having a 803 uniform distribution of radii between 10 m and 200 m, we ob-804 served that, depending on the arrangement of the capillaries, 805 the leading meniscus can be in any of the capillaries whose 806 radii are between 10 m and 100 m, and the non-dimensional 807 breakthrough time lies between T bt = 0.31 and T bt = 0.42. 808 The dynamics of spontaneous imbibition as predicted by 809 this new model is significantly different from that predicted by 810 the classic bundle of non-interactive capillaries (or tubes), for 811 which the leading meniscus is always in the largest radius cap-812 illary. For the interacting multi-capillary system mentioned 813 above, on the contrary, the leading meniscus can be in any 814 of the capillaries having radii between 10 m and 100 m. We 815 observed that the breakthrough occurs earlier than in the clas-816 sic capillary bundle, where it occurs at non-dimensional time 817 T bt = 0.5 for the aforementioned 20-capillary-system, to be 818 compared to the 0.31-0.42 range for the 20-capillary-system 819 mentioned above. Furthermore, for this system the saturation 820 at breakthrough time falls in the range 0.69-0.79, whereas for 821 the classic capillary bundle it is equal to 0.86. The dependence 822 of the saturation as a function of the longitudinal position are 823 also shows a stark contrast between the predictions of the clas-824 sic capillary bundle and the average behavior of the 1000 ar-

917is

  Pc β and some of the impregnating fluid transfers from β 918 to γ, which we assume to correspond to a differential flow 919 rate dq β . The velocity of the meniscus in the capillary β for 920 z > z β (t) is then

A11) 1 .

 1 We initiate the model formulation by finding the largest 941 radius capillary, C i . The pressure field is identical in all 942 capillaries for z < z a (t), and the corresponding pressure 943 gradient is related to the fluid velocity in each capillary 944 by Hagen-Poiseuille's law. Some of the invading fluid 945 from capillary i transfers to other capillaries in the im-946 mediate vicinity of the meniscus position z a (t).

947 2 . 3 .

 23 For z > z a (t), the imbibing fluid in the capillaries C 1 to 948 C (i-1) is separated from the imbibing fluid in the capil-949 laries C (i+1) to C n . We thus classify the capillaries on 950 either sides of the capillary C i in two regions, the cap-951 illaries C 1 to C (i-1) in the first one, the capillaries from 952 C (i+1) to C n in another one. The fluid transfer from the 953 capillary C i is divided among the other capillaries ac-954 cording to their radii. If the fluid transfer to the 'top 955 region' is dq t , the fraction of dq t flowing from capillary 956 C i to a capillary of radius r p would be r 4 p dq t /∑ i-1 q=1 (r 4 q ). 957 Similarly, for the 'bottom region', if dq b is the fluid 958 transfer from C i , the fractional flow in a capillary of 959 radius r r will be r 4 r dq b /∑ n s=i+1 (r 4 s ). This fluid trans-960 fer causes the flow rates to increase in capillaries C 1 to 961 C (i-1) and C (i+1) to C n . 962 The widest capillary among the capillaries C 1 to C (i-1) , 963 C j is now identified. For z a (t) < z < z b (t) the pres-964 sure field in the imbibing fluid is identical in capillaries 965 C 1 to C (i-1) , and is related to the fluid velocity in each 966 capillary by Hagen-Poiseuille's law. In the vicinity of 967 z = z b (t), some of the invading fluid transfers from C j 968 to the capillaries C 1 to C ( j-1) and C ( j+1) to C (i-1) , which 969 increases the flow rate in these capillaries.

970 4 .979 5 .

 45 Similarly, the widest capillary among capilllary C (i+1) 971 to C n , which we denote C k , is chosen. The pressure field 972 is identical in the capillaries C (i+1) to C n for z a (t) < z < 973 z c (t), and the pressure gradient is related to the fluid 974 velocity in each of these capillaries from the Hagen-975 Poiseuille law. A z = z c (t), some of the fluid invad-976 ing C k transfers into the capillaries C (i+1) to C (k-1) and 977 C (k+1) to C n , which increases the flow rate in in these 978 capillaries. The impregnating fluids in the regions encompassing 980 capillaries C 1 to C ( j-1) and C ( j+1) to C (i-1) are sepa-981 rated by displaced fluid in capillary C j for z > z j . Again, 982 the capillary of largest radius among the capillaries C 1 983 to C ( j-1) is identified, as well as the capillary of largest 984 radius among the capillaries C ( j+1) to C (i-1) . The simi-985 lar procedure previously explained for the pressure field 986 and its relation to the fluid velocity is repeated for those 987 two regions.

988 6 . 7 . 8 .

 678 Fig. C.1(a),(b),(f),(g),(i),(j),(k),(l). For the arrangements 1013

FIG. C. 1 .

 1 FIG. C.1. Spontaneous imbibition in a system of four interacting capillaries of radii r α = 80 m, r β = 60 m, r γ = 40 m and r δ = 20 m. The non-dimensional positions of the four menisci are shown as a function of non-dimensional time for all the 12 possible arrangements in (a) to (l). The arrangement, the ordering of the menisci locations, and the breakthrough times for each of the cases (a) to (l) are provided as legends of the plots.

  ,

	239
	240	according to which (1) the pressure equilibrates over the sec-
	241	tions of the capillary system that are entirely filled with the
	242	invading fluid, and (2) fluid transfers from a capillary having
	243	a larger radius to an adjacent capillary having a smaller radius
	244	just before the meniscus, which in the model we assume to oc-
	245	cur at the position of the meniscus. We show this fluid transfer
	246	between adjacent capillaries in the vicinity of the meniscus by
	247	vertical arrows in Fig. 2. We consider the interaction between
	248	the capillaries to be sufficiently low for the Poiseuille flow in
	249	each of the capillaries to be maintained. At any given time
	250	t, the less advanced meniscus (i.e., that for which the imbibi-
	251	tion length is the smallest) will be in the capillary for which
	252	the driving capillary pressure jump across the meniscus is the
		smallest, hence it is will be the meniscus in the α capillary.
	254
	255	illaries. Similarly, the next-less-advanced meniscus is neces-
	256	sarily the β capillary driven by the capillary pressure Pc β , so
	257	at any time t the pressure field is identical in capillaries β , γ
	258	and δ for z α (t) < z < z β (t), and so forth: the pressure field is
	259	identical in the δ and γ capillaries for z β (t) < z < z γ (t).The
	260	imbibition length in capillary δ , z δ (t) is the largest at any time
	261	t.
	262	We now consider one of the random arrangements as shown
	263	in the schematic of Fig. 3, where the order of arrangement
	264	of the capillaries is β γαδ . It was explained by Ashraf et
	265	al., 54 that, for a randomly-arranged interacting capillary sys-
	266	tem, the meniscus in the smallest radius capillary does not
	267	always lead. For this arrangement, depending upon the con-
	268	trasts in the radii, three different positionings of the menisci
	269	are possible as shown in Fig. 3 (a), (b) and (c). At any given
	270	time t, for 0 < z α (t), the pressure field is identical in all capil-
	271	laries, and the pressure drop from the inlet to z α (t) is Pc α . For
	272	z > z α (t), the imbibing fluid is continuous in the capillaries β
	273	and γ, since they are connected. Therefore, the pressure field
	274	is the same in the capillaries β and γ for z α (t) < z < z β (t). As
	275	r β > r γ (meaning that the capillary suction in β is less than
	276	that in γ), during the spontaneous imbibition, z β (t) < z γ (t), at
	277	all times. Although the capillary δ is filled with the imbibing
		phase, the non-wetting fluid in α disconnects it from capillar-

253

For z < z α (t), the pressure field must be identical in all cap-278 ies β , γ for z > z α (t). Therefore, for z > z α (t) the pressure 279 field in δ can be different from that in β , γ. For the arrange-280 ment β γαδ shown in the schematic of Fig. 3, z α < z β < z γ 281 and z α < z δ during the imbibition process and the position of

  (k+1) to C n on the other hand in the 'bottom region'. The

	373
	374	pressure drops are determined in each of the subregions. This
	375	procedure is then performed recursively until the algorithm
	376	has identified the pressure drop in each of the sections for ev-
	377	ery capillary. It can then formulate the governing equations,
		which are consequently solved to obtain the advancement of

363

the 'top region' consisting of the capillaries C 1 to C (i-1) and 364 the 'bottom region' consisting of the capillaries C (i+1) to C n 365 (see Fig.

4

). The largest radius capillaries in each of these 366 two regions are determined and the pressure drop in the re-367 spective regions are determined for sections z a < z < z b and 368 z a < z < z c . Now, each of these two regions is further divided 369 into two subregions each, i.e., containing the capillaires C 1 to 370 C ( j-1) on the one hand and C ( j+1) to C (i-1) on the other hand 371 in the 'top region', and C (i+1) to C (k-1) on the one hand and 372 C 378 all menisci as a function of time. 379 III. RESULTS AND DISCUSSIONS 380

  The capillaries β and γ are on the other side of the capillary 902 α with respect to the capillary δ . As the capillary pressure 903 jump of the capillary β is smaller than that in the capillary γ, 904 the meniscus in β lags behind that in γ. Hence, the imbibing 905 fluid in these capillaries is continuous for z α (t) < z < z β (t).

	the positions of whose axes in a transverse plane would be the
	nodes of a triangular grid.	900	since the pressure in the non-wetting fluid is the atmospheric
			pressure.
		906
		907	Defining ω and (1 -ω) as the fractions of dq α transferred
		908	respectively to β and γ, we can write an equation similar to
		909	Eq. (A4) for both β and α, where ωdq α and (1 -ω)dq α ap-
		910	pear respectively as a differential velocity term arising from
		911	fluid transfer. Considering that the pressure field is the same
		912	in the capillaries β and γ for z α (t) < z < z β (t), we then obtain
		913	in that z range:
		825 826 827	rangements of interacting capillaries. Indeed, the interacting capillary system shows a steep decrease in the saturation with v β (t) + ωdq α πr 2 β πr 2 β = v γ (t) + (1-ω)dq α πr 2 γ γ πr 2 . (A6)
		914	length as compared to the classic capillary bundle. Addition-Combining Eq. (A2) and Eq. (A6), we then obtain the fraction
			ω from the capillaries' radii: ω = r 4 β /(r 4 β +r 4 γ ). Therefore, the
		916
		850
		851	front carried by larger pores, which agrees very well with
		852	the experimental findings. The quantitative predictions for
		853	the positions in time of these two fronts, obtained from av-
		854	eraging over the statistics of randomly-chosen arrangements,
		855	agree well with the measurements.
		856	This generalized model for spontaneous imbibition in a pla-
		857	nar bundle of interacting capillaries, which is semi-analytical
		858	and runs extremely quickly, could be useful for fast as-
		859	sessment of one-dimensional imbibition dynamics in design-
		860	based porous media such as loop heat pipes, diagnostic
			devices and microfluidic devices, or in real porous media
		865
			three-dimensional models by considering parallel capillaries,

828

ally, the interacting multi-capillary system shows that the spa-829 tial arrangement of the capillaries may cause significantly dif-830 ferent saturation values at a given longitudinal position. 831 So, how is this model consisting of a planar bundle of inter-832 acting capillaries to be used to predict spontaneous imbibition 833 in quasi-two-dimensional porous media whose pore size dis-834 tribution is known? We propose to use a stochastic approach, 835 i.e., to consider the average behavior between a large number 836 of randomly-picked spatial arrangements of the capillary di-837 ameters, the distribution of these diameters being equal to the 838 pore size distribution of the real porous medium. We tested 839 that method against data from the literature. Firstly, qualita-840 tive observations relative to which ranges of pore sizes mainly 841 contribute to the leading and lagging fronts of the imbibition 842 interface, and to the longitudinal saturation profile, are consis-843 tent between experiments from the literature and the predic-844 tions of our model. Secondly, to validate the model's quanti-845 tative predictive capacity, we compared its predictions to im-846 bibition measurements in filter paper, performed by Bico and 847 Quéré 47 . The model predicts that the visible leading front is 848 carried by smaller pores and that the bulk saturating front re-849 sponsible for most of the fluid mass invasion is the lagging 861 whose porosity structure can reasonably be considered two-862 dimensional, e.g., paper, thin porous media in general, or lay-863 ered aquifers. 864 Prospects to this work include extending this approach to 901 915 pressure drop in capillaries β and γ for z α (t) < z < z β (t) is

B. Xiao, J. Fan, and F. Ding, "Prediction of relative permeability of unsat-

M. Dong and I. Chatzis, "The imbibition and flow of a wetting liquid along

Accepted to Phys. Fluids 10.1063/5.0123229 Accepted to Phys. Fluids 10.1063/5.0123229 Accepted to Phys. Fluids 10.1063/5.0123229 Accepted to Phys. Fluids 10.1063/5.0123229 Accepted to Phys. Fluids 10.1063/5.0123229

In capillary α, for 0 < z < z α (t), the pressure drop is given by the Hagen-Poiseuille law as,

where µ is the imbibing fluid's viscosity, v α (t) is the instantaneous velocity of the wetting fluid in the capillary α, P 0 is the inlet pressure and P(z α (t),t) is the pressure in the imbibing fluid at z α (t), as shown in Fig. 3. Since the pressure fields are identical in all capillaries for z < z α (t), the pressure gradient is the same in all capillaries , which from Eq. (A1) implies

where the index i (i = α, β , γ, δ ) indicates quantities relative to the capillary of radius r i and v i (t) (i = α, β , γ, δ ) is the velocity of the imbibing fluid for z < z α (t).

The capillary pressure jump through the fluid-fluid interface is Pc α at z α (t), where some of the imbibing fluid transfers from the capillary α to other capillaries. The volumetric fluid transfer from the capillary α to the capillaries β and γ is dq α , whereas the fluid transfer from the capillary α to the capillary δ is dq α . The velocity of the advancing meniscus in capillary α, dz α /dt, is thus given by

For z α (t) < z < z δ (t), the velocity of the fluid in capillary δ is similarly given by

so the pressure drop in the capillary δ between z = z α (t) and

At z = z δ (t), the pressure jump across the meniscus is Pc δ , for i = α, β , γ, δ . Note that the prefactor 2 is controlled by cir-927 cular cross-section of the tube, another geometry (e.g., square 934

Appendix B: Generalization of the model for an arbitrary 938 number of capillaires

939

The following step-by-step procedure must be followed: