Júlia Jacoby
email: julia.jacoby.s@gmail.com

Patricia Luciano Mancini

Márcio Sophie Lanco Bertrand

Amorim Efe

Leandro Bugoni

Guilherme Tavares Nunes

Márcio Amorim Efe

Guilherme Tavares Nunes

Biogeographic variation on dietary aspects of a widely distributed seabird

Keywords: Brown booby, Local adaptation, Marine biogeography, Seabirds, Stable isotopes, Trophic plasticity

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Trophic relationships of widely distributed species may not be homogeneous along their home range, as local conditions influence the availability of food resources [START_REF] Garvey | Trophic ecology[END_REF]. This can allow the identification of dietary ariations between groups, which can be used as tools to understand local adaptation [START_REF] Sanford | Local selection and latitudinal variation in a marine predator-prey interaction[END_REF]Barnagaud t al. 2019). In addition, understanding prey composition could potentially illustrate and refine biogeographic boundaies at regional [START_REF] Lozano | Biogeographical patterns in the diet of the wildcat, Felis silvestris Schreber, in Eurasia: factors affecting the trophic diversity[END_REF][START_REF] Piatt | Biogeography of pelagic food webs in the North Pacific[END_REF] or lobal scales [START_REF] Duffy | Global trophic ecology of yellowfin, bigeye, and albacore tunas: understanding predation on micronekton communities at ocean-basin scales[END_REF][START_REF] Romano | Global biogeographical patterns in the diet of a cosmopolitan avian predator[END_REF].

Trophic plasticity, i.e. the capacity of an organism to adapt its diet according to the variability of food resources vailable in time and space [START_REF] Larkin | Interspecific competition and population control in freshwater fish[END_REF][START_REF] Gerking | Feeding ecology of fish[END_REF], is an advantageous strategy for organisms inhabiting places subject to seasonal fluctuations in environmental conditions [START_REF] Gerking | Feeding ecology of fish[END_REF][START_REF] Abelha | Plasticidade trófica em peixes de água doce[END_REF]. For instance, seabirds are able to adjust their foraging behavior according to oceanic features, whether static (depth, slope) or ephemeral (chlorophyll, sea surface temperature, sea-ice distribution) [START_REF] Widmann | Habitat use and sex-specific foraging behavior of Adélie penguins throughout the breeding season in Adélie Land, East Antartica[END_REF][START_REF] Gilmour | Plasticity of foraging behaviors in response to diverse environmental conditions[END_REF]. Such plasticty can also be advantageous for organisms that depend on resources with unequal and ephemeral distributions in space [START_REF] Harding | Prey density and the behavioral flexibility of a marine predator: the common murre (Uria aalge)[END_REF], especially for species with high mobility and wide geographic distribution such as sharks [START_REF] Drymon | Trophic plasticity in the Atlantic sharpnose shark (Rhizoprionodon terraenovae) from the north central Gulf of Mexico[END_REF], marine mammals [START_REF] Spitz | Dietary plasticity of the oceanic striped dolphin, Stenella coeruleoalba, in the neritic waters of the Bay of Biscay[END_REF], and migratory birds [START_REF] Parrish | Behavioral, energetic and conservation: implications of foraging plasticity during migration[END_REF][START_REF] Briggs | Feeding ecology of the common tern Sterna hirundo in a wintering area in southern Brazil[END_REF]).

In the long term, intraspecific variation in diet could lead to local adaptation and even influence genetic structure [START_REF] Pilot | Dietary differentiation and the evolution of population genetic structure in a highly mobile carnivore[END_REF].

Diet heterogeneity can contribute to local adaptation, intensifying processes that lead to gene flow disruption. Isolation by Environment/Ecology (IBE) is a pattern in which population differentiation is increased according to environmental differences [START_REF] Wang | Isolation by environment[END_REF], and is the most common model to explain population structure in animals [START_REF] Sexton | Genetic isolation by environment or distance: which pattern of gene flow is most common?[END_REF], including seabirds such as prions and petrels [START_REF] Quillfeldt | Cool, cold or colder? Spatial segregation of prions and blue petrels is explained by differences in preferred sea surface temperatures[END_REF], and brown boobies Sula leucogaster (Nunes and Bugoni 2018). According to IBE, gene flow is higher between environmentally similar sites, while in distinct environments, local adaptation intensifies selection against immigrants, leading to isolation [START_REF] Richardson | Microgeographic adaptation and the spatial scale of evolution[END_REF][START_REF] Sexton | Genetic isolation by environment or distance: which pattern of gene flow is most common?[END_REF]) and consequently genetic differentiation between subpopulations (i.e. population structure) [START_REF] Hartl | Principles of population genetics[END_REF]. In this context, variability in the diet of a widely distributed species may act as a complementary cause to population structure because variations in food availability and composition can contribute to the development of intraspecific genetic [START_REF] Sanford | Local selection and latitudinal variation in a marine predator-prey interaction[END_REF], morphological [START_REF] Grant | Unpredictable evolution in a 30-year study of Darwin's finches[END_REF], physiological [START_REF] Fohringer | Integrating omics to characterize eco-physiological adaptations: how moose diet and metabolism differ across biogeographic zones[END_REF], and behavioral variations [START_REF] Moseley | Foraging ecology and ecophysiology of Cape gannets from colonies in contrasting feeding environments[END_REF].

Information on the distribution of organisms has been historically used to establish biogeographic zonings in the oceans [START_REF] Forbes | The natural history of the European seas[END_REF][START_REF] Briggs | Feeding ecology of the common tern Sterna hirundo in a wintering area in southern Brazil[END_REF] or to refine and illustrate patterns previously established [START_REF] Floeter | Atlantic reef fish biogeography and evolution[END_REF][START_REF] Menni | Marine biodiversity at the community level: zoogeography of sharks, skates, rays and chimaeras in the southwestern Atlantic[END_REF][START_REF] Pereira | Parasites of Urophycis brasiliensis (Gadiformes: Phycidae) as indicators of marine ecoregions in coastal areas of the South American Atlantic[END_REF]. Alternatively, additional proposals used environmental parameters such as ocean currents, bathymetry, and primary productivity as input for marine zonings [START_REF] Sherman | Sustainability, biomass yields, and health of coastal ecosystems: an ecological perspective[END_REF][START_REF] Longhurst | An estimate of global primary production in the ocean from satellite radiometer data[END_REF]. Fish is the most common vertebrate group used to identify biogeographic boundaries in the marine realm, from sedentary species associated with coral reefs [START_REF] Pinheiro | South-western Atlantic reef fishes: zoogeographical patterns and ecological drivers reveal a secondary biodiversity centre in the Atlantic Ocean[END_REF] to highly mobile species such as tunas and swordfishes [START_REF] Reygondeau | Biogeography of tuna and billfish communities[END_REF]). However, an inherent challenge in using these groups as indicators of marine life distribution is associated with the need to sample large areas in the ocean, whose sampling is hampered by logistical and methodological constraints [START_REF] Robinette | Spatial patterns of recruitment in a demersal fish as revealed by seabird diet[END_REF][START_REF] Duffy | Global trophic ecology of yellowfin, bigeye, and albacore tunas: understanding predation on micronekton communities at ocean-basin scales[END_REF][START_REF] Piatt | Biogeography of pelagic food webs in the North Pacific[END_REF]. Thus, assessing the prey composition of widely distributed marine organisms, such as seabirds, sampled along their distribution, can provide relevant information to refine biogeographic patterns without the biases of fish sampling methods applied to different locations.

Seabirds depend on the marine environment for at least part of their lives [START_REF] Votier | Seabirds[END_REF]) and present consistent philopatry (i.e. tendency to return to the natal site to breed), despite their high mobility (Schreiber and Burger 2002;[START_REF] Friesen | Speciation in seabirds: why are there so many species and why aren't there more?[END_REF]. Due to the adaptation and specialization to conditions around colonies, seabirds can present population differentiation even between nearby colonies (Nunes and Bugoni 2018;[START_REF] Danckwerts | Extreme philopatry and genetic diversification at unprecedented scales in a seabird[END_REF]. Therefore, the use of food resources available around breeding areas makes local adaptation a potential mechanism of population differentiation in the group [START_REF] Friesen | Mechanisms of population differentiation in seabirds[END_REF][START_REF] Friesen | Speciation in seabirds: why are there so many species and why aren't there more?[END_REF], although historical aspects may also influence population structure [START_REF] Lombal | Identifying mechanisms of genetic differentiation among populations in vagile species: historical factors dominate genetic differentiation in seabirds[END_REF]. Additionally, such processes make seabirds indicators of spatial and temporal variability of food resources, so that dietary information and its variations among colonies can contribute to understanding the distribution patterns of prey species [START_REF] Montevecchi | Dietary changes of seabirds indicate shifts in pelagic food webs[END_REF][START_REF] Piatt | Biogeography of pelagic food webs in the North Pacific[END_REF].

The brown booby is a strictly marine bird widely distributed in all ocean basins in tropical and subtropical regions [START_REF] Nelson | The Sulidae: gannets and boobies[END_REF]. This is a piscivorous species that feed on prey around colonies (i.e. central-place foragers) in an opportunistic way [START_REF] Castillo-Guerrero | Environmentally-mediated flexible foraging strategies in brown boobies in the Gulf of California[END_REF]Nunes et al. 2018) in addition to consuming discards from fisheries [START_REF] Schreiber | Brown booby (Sula leucogaster)[END_REF]. Brown boobies hold a marked phylogeographic structure, so basal populations inhabit the northwest Pacific Ocean, from where dispersal historically occurred towards the Indian and Atlantic Oceans [START_REF] Morris-Pocock | Mechanisms of global diversification in the brown booby (Sula leucogaster) revealed by uniting statistical phylogeographic and multilocus phylogenetic methods[END_REF]). In the southwestern Atlantic Ocean, brown booby colonies are distributed under a heterogeneous seascape [START_REF] Nunes | When Bergmann's rule fails: evidences of environmental selection pressures shaping phenotypic diversification in a widespread seabird[END_REF]) and colony landscapes (Nunes et al. 2018) which are suggested to influence population structure. Low-latitude colonies are located around the Equator in oceanic zones, while the southern breeding areas are on the continental shelf and close to the mainland, ranging from 27°50′S to 0°55′N, comprising 28° of latitude.

Differences in prey availability around brown booby colonies may expose them to distinct selective pressures associated with the use of food resources, which may represent a potential explanation or consequence for the population structure observed in the southwest Atlantic Ocean. In this region, colonies located on the continental shelf share similar environmental conditions, and individuals exhibit gene flow between them, while oceanic colonies maintain higher gene flow between each other (Nunes and Bugoni 2018). However, the colony located in São Pedro e São Paulo (hereafter 'SPSP'), a remote small oceanic archipelago, is genetically isolated from the remaining colonies, which has been associated with local adaptation promoted by intense intraspecific competition for nesting areas (Nunes et al. 2018). Considering that prey distribution may be associated with oceanographic variables such as chlorophyll-α and sea surface temperature, the spatial diet patterns associated with population structure can contribute to demonstrating that IBE is an important mechanism in seabird populations. Additionally, patterns of variation in the diet between breeding sites of the species may represent useful information for the refinement of biogeographic boundaries in the marine realm since it is a proxy of important biotic and environmental variables in the seascape. In this context, the use of complementary techniques to assess the diet of an organism provides accurate information on variations in the use of food resources in space and time [START_REF] Garvey | Trophic ecology[END_REF].

The acquisition of dietary information from seabirds is facilitated because they spontaneously regurgitate the ingested food in stressful situations (Mallet-Rodrigues 2010). From this, prey can be identified at the species level, even if partially digested. Nevertheless, this method is limited because of the short time window represented (i.e. last meal) (Schreiber and Burger 2002). Thus, the use of stable isotopes is complementary to the analysis of regurgitated material as they can act as natural markers of trophic position (i.e. δ 15 N) and feeding environment (i.e. δ 13 C) [START_REF] Fry | Stable isotope ecology[END_REF], such as between coastal and oceanic zones [START_REF] Magozzi | Using ocean models to predict spatial and temporal variation in marine carbon isotopes[END_REF]. Isotopic measurements can be obtained from different tissues, which have distinct turnover rates and therefore provide information from wider temporal windows than information obtained from regurgitates [START_REF] Dalerum | Resolving temporal variation in vertebrate diets using naturally occurring stable isotopes[END_REF]. Additionally, stable isotope values can also be used as coordinates to an area (δ-space) defined as 'isotopic niche' [START_REF] Newsome | A niche for isotope ecology[END_REF]). Therefore, the use of complementary techniques makes it possible to understand the use of food resources more accurately, due to the greater representativeness of the diet in the combined datasets.

In the present study, we aimed to assess the variation in the brown booby diet in the southwest Atlantic Ocean and to test its association with pre-established biogeographic zonings and population structure, as well as their potential use as biological samplers of the fish communities. For this, we reviewed and compiled data on regurgitate material (hereafter 'regurgitates') and stable isotopes of carbon and nitrogen from previous publications and unpublished databases. We expect the diet to vary spatially according to environmental heterogeneity, fitting into biogeographic patterns and the population structure already described for the species in the southwestern Atlantic Ocean.

Materials and methods

Study area

The study area comprises the entire distribution of brown booby colonies in the southwestern Atlantic Ocean, from Moleques do Sul islands (hereafter 'Moleques'; 27°50′42″S, 48°25′49″W) to SPSP (0°55′00″N, 29º20′45″W) [START_REF] Sick | Abundance and spatial-temporal distribution of the shrimp Xiphopenaeus kroyeri (Decapoda: Penaeidae): an exploited species in southeast Brazil[END_REF]) (Fig. 1a). Coastal colonies (Moleques, Currais, Cagarras, Cabo Frio, Santana, and Abrolhos) are located on the continental shelf with distance from mainland ranging from less than 1 to a maximum of 70 km and are influenced by shelf waters, such as Tropical Water and Subtropical Shelf Water, in addition to river discharges [START_REF] Piola | Subtropical shelf front off eastern South America[END_REF] and upwelling processes (e.g. South Atlantic Central Water in Cabo Frio) [START_REF] Valentin | The Cabo Frio upwelling system, Brazil[END_REF]. Offshore the continental shelf, the oceanic colonies are located at the Fernando de Noronha archipelago ('Noronha'), Rocas atoll ('Rocas'), and SPSP, ranging from 200 to a maximum of 1000 km from the mainland and under the influence of the South Equatorial Current [START_REF] Richardson | Mapping climatological seasonal variations of surface currents in the tropical Atlantic using ship drifts[END_REF]. In addition, SPSP is also influenced by the Equatorial Undercurrent [START_REF] Araujo | Modelagem matemática da circulação oceânica na região equatorial[END_REF].

Data collection and analysis

Regurgitates

We gathered data on regurgitates from an active and non-systematic search in articles, book chapters, technical reports, and conference abstracts. We obtained six sources of datasets published between 2004 and 2023, two Ph.D. theses, one conference abstract, one technical report, and one unpublished database from Noronha, totaling nine breeding sites from Moleques to SPSP (Table 1, and Table S1), which were classified into at least two distinct biogeographical zones (Fig. 1b). The prey taxonomy was updated according to Eschmeyer's Catalog of Fishes [START_REF] Fricke | Eschmeyer's catalog of fishes: genera, species, references[END_REF] for fish and the World Register of Marine Species (WoRMS 2022) for other taxa (e.g. squids). Then, we extracted information on the number of regurgitates (i.e. samples), and prey abundance at the family and species levels per colony from each data Table 1 Data sources and their respective sample sizes obtained by non-systematic search that described the diet of Sula leucogaster from regurgitated material and stable isotopes in colonies along the southwestern Atlantic Ocean source using the 'vegan' [START_REF] Oksanen | Vegan: community ecology package. R package version 2.5-7[END_REF]) and 'plyr' [START_REF] Wickham | The split-apply-combine strategy for data analysis[END_REF] packages in the R software (R Core Team 2021). For descriptive metrics of prey species and families, we calculated the frequency of occurrence (FO%), relative abundance (RA%), richness (R), and α-diversity (Shannon Index) [START_REF] Magurran | Measuring biological diversity[END_REF] from each data source separately. Prey items not identified at species and family levels were disregarded from species and family analysis, even if it was the only representative of a taxonomic lineage. Items derived from cannibalism were also removed from the analysis (as showed in A. Kohlrausch unpubl. data).

To assess differences between prey communities, we used complementary statistical approaches. Due to the differences in sampling efforts between studies and databases, the dataset was converted to a binary matrix of presence/ absence. Paired distances between colonies were estimated with the beta-sim index (β sim ) with subsequent generation of a dendrogram based on Euclidean distances using the 'betapart' [START_REF] Baselga | Betapart: partitioning beta diversity into turnover and nestedness components[END_REF]) and 'ggdendro' (Vries and Ripley 2020) packages. The β sim was calculated according to the following equation [START_REF] Kreft | A framework for delineating biogeographical regions based on species distributions: global quantitative biogeographical regionalizations[END_REF]:

where "a" = total number of shared species between colonies, and "b,c" = number of exclusive species of each colony. β sim varies between 0 (low dissimilarity, identical taxa list) to 1 (high dissimilarity, no shared taxa). To visually identify similarities between colonies and location (Fig. 1a), genetic (as described below), and biogeographic criteria (Fig. 1b), we performed a non-metric multidimensional scaling (NMDS) with the Bray-Curtis index. Additionally, we ran

sim = 1 - a min(b,c) +a ,
a permutational analysis of variance (PERMANOVA) with 999 permutations to test differences between the prey dataset and the criteria, assuming no dietary differences as the null hypothesis [START_REF] Legendre | Numerical ecology[END_REF]. For graphing, we used the 'ggplot2' package [START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF]).

The criteria used for the analyses were categorical and considered geographic, genetic, and biogeographical aspects. Colony location concerning the continental shelf was assigned as 'coastal' or 'oceanic' whether the colony was on or off the shelf, respectively. The genetic criteria followed the previously published clustering based on microsatellites (Nunes and Bugoni 2018), where the K = 2 arrangement splits the colonies into two groups (SPSP; and Others), and the K = 3 suggests isolation of SPSP, in addition to a weaker gene flow between "Coastal" (Moleques, Cagarras, and Abrolhos) and "FN-RO" (Noronha and Rocas) groups. For standardization purposes, we included the coastal colonies not studied by Nunes and Bugoni (2018) in the "Coastal" group. Finally, we considered distinct biogeographical classifications of the marine environment: the Biogeochemical Provinces of Longhurst (hereafter "Longhurst's Provinces") [START_REF] Longhurst | An estimate of global primary production in the ocean from satellite radiometer data[END_REF]; the Large Marine Ecosystems, (hereafter "LME") [START_REF] Sherman | Sustainability, biomass yields, and health of coastal ecosystems: an ecological perspective[END_REF]; and the Provinces and Ecoregions established by [START_REF] Spalding | Marine ecoregions of the world: a bioregionalization of coastal and shelf areas[END_REF] (hereafter "Spalding's Provinces" and "Spalding's Ecoregions", respectively). The Longhurst's Provinces were outlined based on chlorophyll-α concentration, which is a proxy of primary productivity for coastal and oceanic regions. LME were delimited from the characterization of the margins of ocean basins based on bathymetric, hydrographic, and depth data; in addition to fish distributions (Sherman 1994). Spalding's Provinces and Ecoregions were developed from a review of pre-existing classifications which combined benthic and pelagic taxa distribution with multiple oceanographic and geomorphologic factors (i.e. islands, current dynamics, upwellings, and salinity), in addition to adjustments according to geopolitical limits. In this classification, the provinces are identified by names and subdivided into ecoregions, represented by numbers [START_REF] Spalding | Marine ecoregions of the world: a bioregionalization of coastal and shelf areas[END_REF].

Stable isotopes

We used carbon (δ 13 C) and nitrogen (δ 15 N) isotopic ratios from brown booby whole blood samples obtained from scientific articles, technical reports, and unpublished datasets (Table 1). All data were collected from breeding individuals and without distinction between males and females. Sample processing of the unpublished datasets consisted in place a few drops of the whole blood on tubes, which were air-dried. Blood samples were freeze-dried, ground, and homogenized.

Subsamples of 1 mg were weighed into tin cups and analyzed in a mass spectrometer (Mancini et al. 2014;Nunes et al. 2018). Standards used were Vienna Pee Dee belemnite and atmospheric air for carbon and nitrogen, respectively. Measurement precision of both δ 15 N and δ 13 C is described in each reference as well as secondary isotopic reference materials used by each laboratory (Table 1). Lipids were not extracted from seabird blood due to their expected low concentration [START_REF] Bearhop | Influence of lipid and uric acid on δ 13 C and δ 15 N values of avian blood: implications for trophic studies[END_REF]. We calculated the mean and standard deviation of each colony using the 'FSA' package [START_REF] Ogle | FSA: fisheries stock analysis. R package version 0[END_REF]. Then, we tested univariate differences between colonies using the Kruskal-Wallis test, and Mann-Whitney U test as post-hoc, using False Discovery Rate for adjusting P values [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. Additionally, we calculated the two-dimensional isotopic niche breadth based on Bayesian ellipses as implemented in the 'SIBER' package [START_REF] Jackson | Comparing isotopic niche widths among and within communities: SIBER-stable isotope Bayesian ellipses in R[END_REF].

Results

Regurgitates

In total, 877 regurgitates composed the dataset, with 4270 prey belonging to at least 105 species, 43 families, and 17 orders (Table S2), showing an increasing accumulation curve of prey species relative to the number of datasets (Fig. S1). Different sample sizes were reported among studies, with Rocas (n = 1128) presenting the highest prey numbers and Noronha (n = 85) the lowest ones (Fig. S2).

For the analysis at the species-level only, we considered 3727 prey belonging to 99 species. Moleques exhibited the highest species richness (R = 39) and Noronha the lowest (R = 8) (Fig. 2a). Considering all colonies, the most frequent species were Sardinella aurita and Harengula clupeola (FO = 55.5%), followed by Micropogonias furnieri, Ctenosciaena gracilicirrhus, Pellona harroweri, Hirundichthys affinis, Cetengraulis edentulus and Opisthonema oglinum (FO = 44.4%). Additionally, C. edentulus also showed the highest relative abundance in Moleques (RA = 24.9%) and Cagarras (RA = 39.6%), and S. aurita in Currais (RA = 15.7%) (Table S3).

For the family-level analysis, 43 families and 4200 prey were included. The lowest number of families was recorded in Noronha (n = 4) and the highest in Moleques (n = 24) (Fig. 2b). The highest frequency of occurrence was observed for Clupeidae (FO = 88.9%) followed by Scombridae, Sciaenidae, Hemiramphidae, Engraulidae, and Carangidae (FO = 66.7%). Exocoetidae and Hemiramphidae occurred only in Abrolhos, among sites over the continental shelf, and in the oceanic colonies and together accounted for over 90% of the prey consumed by SPSP boobies. Among colonies south of Abrolhos, the most representative prey families were Sciaenidae, Engraulidae, and Clupeidae (Table S4 and Fig. S3).

Diversity values ranged from 2.68 in Currais to 0.29 in Noronha (Fig. S4a) and oceanic colonies showed lower diversity in comparison to coastal colonies (Fig. S4b). Considering the genetic criteria, SPSP had lower diversity compared to the remaining colonies for K = 2 (Fig. S4c); while for K = 3, the coastal group had higher diversity compared to the FN-RO group and SPSP (Fig. S4d). For the LME and Spalding's Provinces, the highest diversity was observed, respectively, in the South Brazilian Shelf (SBS) and Warm Temperate Southwestern Atlantic (WTSA) regions, which correspond to the southernmost portions of the boundaries established for these classifications (Fig. 3a,c). The Longhurst's Provinces showed decreasing diversity values from the Brazil Current Coastal Province (BRAZ) to the Western Tropical Atlantic Province (WTRA) (Fig. 3b). Spalding's Ecoregions in the Tropical Southwestern Atlantic (TSA, named TSA74 and TSA73), which correspond to the FN-RO and SPSP colonies, respectively, had lower diversity (Fig. 3d).

Colonies were clustered into two large groups following the β sim index: the oceanic Noronha, Rocas, and SPSP; and the colonies on the continental shelf (Fig. 1c, Table S5). PERMANOVA grouping the prey community with the predefined criteria generated significant results (P ≤ 0.01) for geographic, genetic, and biogeographic, except for the clustering K = 2 (P = 0.02). NMDS showed that colonies were grouped according to all criteria (Fig. 4 and Fig. S5), except for the delineation of Spalding's Provinces (Fig. 4c).

Stable isotopes

The analyzed dataset comprised 194 blood samples from seven colonies (Table 1). The mean carbon and nitrogen isotopic ratios obtained from the data sources ranged from -17.7‰ in Cagarras to -16.5‰ in Noronha for δ 13 C, and from 9.9‰ in Noronha to 14.3‰ in Santana for δ 15 N (Table S6). Univariate differences between colonies were significant both for δ 13 C and δ 15 N values (P ≤ 0.01). Pairdifferences between colonies resulted in significant differences for δ 15 N between all colonies, except between Noronha and Rocas (P = 0.42). For δ 13 significant pairwise differences were identified among all colonies, except between Moleques and Abrolhos (P = 0.27), Moleques and Rocas (P = 0.03), Moleques and SPSP (P = 0.03), and Abrolhos and Rocas (P = 0.40). Differences occurred between genetic and biogeographic criteria for δ 15 N, except between Spalding's Ecoregions WTS180 and TSA76 (P = 0.02), which correspond to the Moleques-Cagarras and Santana-Abrolhos, respectively. Differences in δ 15 N between coastal and oceanic colonies were not significant (P = 0.86). Considering the genetic criteria = 3, paired significant differences were identified for δ 13 C between FN-RO and the coastal group, and also between FN-RO and SPSP (P ≤ 0.01). Additionally, the biogeographic criteria for δ 13 C resulted in significant differences between Longhurst's Provinces, LME, Spalding's Provinces, and most Ecoregions (P ≤ 0.01), except between TSA73 and TSA76 (P = 0.66), which corresponds to SPSP and Santana-Abrolhos, respectively. Finally, we also observed similarities between ellipse areas among the brown booby colonies, with Noronha and Rocas holding similar areas, with lower δ 15 N values and Cagarras, Santana, and Moleques showing higher δ 15 N values and sharing similar δ 13 C ranges. Population from Abrolhos showed the widest ellipses among the colonies (ellipse area = 9.58) (Fig. 5).

Discussion

In the present study, we evidenced the diversity and spatial patterns in the diet of brown boobies, which forage over a hundred prey species in the southwest Atlantic Ocean. Diet variation was associated with colony location, genetic clustering, and biogeographical zonings. Findings suggest that the diet of brown boobies can be directly influenced by the prey available around breeding sites, considering the number of demersal prey species associated with trawling in coastal colonies and tropical prey species in oceanic colonies. These results reinforce the role of local adaptation in shaping population structure, considering the association of the diet (i.e. prey community and pairwise differences of δ 15 N between colonies) with the known brown booby population structure in the southwest Atlantic Ocean (Nunes and Bugoni 2018), especially between colonies on and off the continental shelf. Additionally, the results demonstrate the usefulness of dietary aspects of widely distributed seabirds in informing biogeographic boundaries, considering that these species can be used as samplers of local biodiversity.

The diet of the brown booby varies according to the prey availability around the colonies throughout its global distribution, which can even make such interactions predictable. For example, differences in foraging behavior and diet can be identified according to variations in the abundance of their main prey [START_REF] Mellink | Diet of eastern Pacific brown boobies Sula leucogaster brewsteri on Isla San Jorge, north-eastern Gulf of California, and an April comparison with diets on the middle Gulf of California[END_REF], to oceanographic dynamics in the Gulf of California [START_REF] Castillo-Guerrero | Environmentally-mediated flexible foraging strategies in brown boobies in the Gulf of California[END_REF], and to the risk of kleptoparasitism in Australia [START_REF] Miller | Sexual segregation in tropical seabirds: drivers of sex specific foraging in the brown booby Sula leucogaster[END_REF]. Similarly, brown boobies breeding in the north [START_REF] Harrison | Hawaiian seabird feeding ecology[END_REF]) and south [START_REF] Harrison | The diet of the brown booby Sula leucogaster and masked booby Sula dactylatra on Rose Atoll, Samoa[END_REF] of the Pacific Ocean take advantage of any prey available around the colonies. Over a hundred prey species were reported in the diet of brown booby in the southwest Atlantic Ocean compiled here, with marked differences associated with previously proposed biogeographic patterns, which in turn are based on oceanographic parameters (i.e. biological, physical, and chemical). Therefore, the analysis of environmental covariates around the colonies and their temporal variations can predict aspects of the brown booby diet, which may represent an interesting model to predict potential impacts of climate change scenarios in trophic interactions of seabirds [START_REF] Barbraud | Effects of climate change and fisheries bycatch on Southern Ocean seabirds: a review[END_REF][START_REF] Ramírez | Stable isotopes in seabirds reflect changes in marine productivity patterns[END_REF][START_REF] Watanuki | Seabird reproductive responses to changing climate and prey communities are mediated by prey packaging[END_REF].

Isotopic niche variation among the colonies may be associated with behavioral aspects driven by food availability, but also with differences in isotopic baselines. Nitrogen values could be influenced by the feeding ecology of the main prey species at each colony. In southern sites, there is a substantial contribution of carnivorous species, such as sciaenids [START_REF] Lucena | Seasonal variation in the diets of bluefish, Pomatomus saltatrix (Pomatomidae) and striped weakfish, Cynoscion guatucupa (Sciaenidae) in southern Brazil: implications of food partitioning[END_REF][START_REF] Sealey | Setting geographic priorities for marine conservation in Latin America and the Caribbean. The Nature Conservancy[END_REF], whereas the predominant prey in Noronha and Rocas are planktivorous, such as clupeids (Whitehead 1985) and flying fish [START_REF] Van Noord | Flyingfish feeding ecology in the eastern Pacific: prey partitioning within a speciose epipelagic community[END_REF]. Boobies from Abrolhos forage on both planktivorous and carnivorous species (e.g. Carangidae), the latter being usually a bycatch of shrimp trawling in this region [START_REF] Santos | Avaliação quali-quantitativa da ictiofauna acompanhante na pesca do camarão setebarbas, Xiphopenaeus kroyeri (Heller, 1862) no município de Caravelas (Bahia-Brasil)[END_REF]. Considering all colonies, Abrolhos showed the widest and largest isotopic niche, which may be associated with the variety of prey consumed at this site, both pelagic and demersal. This wider isotopic niche in Abrolhos could indicate that seabirds are feeding on a greater variety of prey, while a narrow niche, as observed in Noronha, would suggest a lower variety of prey ingested (Mancini et al. 2014), as corroborated by regurgitates. Different isotopic values may also be associated with spatiotemporal variations in isotopic baselines among colonies, which may be influenced by fluctuations of environmental conditions in foraging areas [START_REF] Bond | A practical introduction to stable-isotope analysis for seabird biologists: approaches, cautions and caveats[END_REF], such as the influence of the upwelling in Cabo Frio [START_REF] Soares | Seasonal variability of δ 13 C and δ 15 N of fish and squid in the Cabo Frio upwelling system of the southwestern Atlantic[END_REF]) that can impact northern areas, e.g. Santana [START_REF] Bauer | Marine fishes (Elasmobranchii and Teleostei) from the Santana Archipelago, a marine protected area in the southwestern Atlantic[END_REF]. However, the variations in isotopic niche breadth observed in the present study suggest this is closely associated with the resources available around the colonies, which represents additional evidence of adaptation to local conditions.

Because brown boobies use fishery discards as a food source, the variation in prey composition and predictability may also be associated with the spatial distribution of fisheries. Demersal prey (e.g. Sciaenidae and Batrachoididae) discarded from shrimp trawling were frequent in the diet of boobies from southern coastal colonies [START_REF] Vianna | Bony fish bycatch in the southern Brazil pink shrimp (Farfantepenaeus brasiliensis and F. paulensis) fishery[END_REF][START_REF] Branco | Análise quali-quantitativa da ictiofauna acompanhante na pesca do camarão sete-barbas, na Armação do Itapocoroy, Penha, Santa Catarina[END_REF][START_REF] Cattani | Avaliação da ictiofauna da fauna acompanhante da pesca do camarão sete-barbas no município de Pontal do Paraná, Brasil[END_REF][START_REF] Paranhos | Spatial patterns of distribution and the influence of seasonal and abiotic factors on demersal ichthyofauna in an estuarine tropical bay[END_REF]) and absent in oceanic colonies, suggesting an influence of such fishery in the dietary composition of brown boobies. However, prey usually associated with discards were scarce in Cagarras, despite the occurrence of bottom trawling around the archipelago [START_REF] Moraes | A pesca no entorno do Monumento Natural das Ilhas Cagarras[END_REF][START_REF] Amorim | Marine protected area and the spatial distribution of the gill net fishery in Copacabana, Rio de Janeiro, RJ, Brazil[END_REF]. It is possible that the dataset from Cagarras is biased by insufficient sampling (i.e. one breeding season) (Torres et al. unpubl. data), but the availability of engraulids and clupeids around the colony [START_REF] Jablonski | Fisheries and conflicts in Guanabara Bay, Rio de Janeiro, Brazil[END_REF][START_REF] Paranhos | Spatial patterns of distribution and the influence of seasonal and abiotic factors on demersal ichthyofauna in an estuarine tropical bay[END_REF]) may favor "natural prey" over discards from trawling, as observed in other Sulidae species [START_REF] Votier | Individual responses of seabirds to commercial fisheries revealed using GPS tracking, stable isotopes and vessel monitoring systems[END_REF][START_REF] Moseley | Foraging ecology and ecophysiology of Cape gannets from colonies in contrasting feeding environments[END_REF]). Another possible influence may be from interspecific interactions such as kleptoparasitism of frigatebirds [START_REF] Cunha | Aves do Monumento Natural das Ilhas Cagarras[END_REF]. At oceanic colonies, fisheries targeting large pelagic fishes, such as tunas and sharks, are carried out by using live bait, such as clupeids in Noronha [START_REF] Sazima | Occupational hazards: brown boobies (Sula leucogaster) as a nuisance to fishermen at Fernando de Noronha Island, with comments on injuries inflicted to the birds[END_REF] and flying fish in SPSP [START_REF] Viana | Fisheries in the Saint Peter and Saint Paul archipelago: 13 years of monitoring[END_REF]. Surplus baits are discarded after fishing and then used by brown boobies as a food source around Noronha and SPSP (pers. obs. GT Nunes), which could explain the high contribution of these prey to the diet of boobies, in addition to the ones that are naturally predated. In this context, the variety of fisheries throughout the study area could play an important role in the prey composition of boobies and therefore decisions targeting the management of fisheries and fishery discards should consider their effects on the persistence of seabird populations.

Prey community and isotopic niche variations among colonies were associated with the population structure of brown boobies in the southwest Atlantic Ocean, reinforcing the idea of local adaptation as a relevant mechanism for the promotion/disruption of gene flow. Genetic clustering among coastal colonies, as well as between Noronha and Rocas (Nunes and Bugoni 2018) was mirrored in dietary aspects explored in this study, demonstrated by the clustering generated from the β sim index and the associations between the genetic clusters and isotopic ratios. This suggests that population structure could be influenced by local adaptation in obtaining available food resources. Dispersal ability could be more likely among colonies with similar environmental conditions (e.g. sea surface temperature, chlorophyll-α concentration, fisheries), representing similar prey communities and selective pressures, as it was observed among coastal colonies. Additionally, Abrolhos could represent a steppingstone for dispersal between coastal and oceanic colonies, due to the occurrence of prey species found both in the northern (i.e. Hemiramphidae and Exocoetidae) and southern colonies (i.e. fishery discards), which is supported in the genetic clustering model (Nunes and Bugoni 2018). Finally, SPSP is genetically isolated from the remaining colonies, which is suggested to be associated with selective pressures in the nesting area (Nunes et al. 2018). However, the high availability and consumption of flying fish in SPSP (RA ~ 63%) can also influence population differentiation through local adaptation, even when compared to adjacent colonies, such as Noronha and Rocas. Therefore, spatial patterns of diet variation associated with population structure provide support to the IBE model [START_REF] Wang | Isolation by environment[END_REF] and local adaptation as an important mechanism of population structuring, even in highly mobile organisms such as seabirds.

Dietary aspects investigated from regurgitates and isotopic values (especially δ 15 N) were demonstrated to be relevant tools to inform biogeographic patterns, even considering differences between boundaries of each zoning. Comparing the boundaries of Spalding's Provinces with LME, Santana is placed differently, grouped with the northern colonies in the first and with the southern colonies in the second one. [START_REF] Spalding | Marine ecoregions of the world: a bioregionalization of coastal and shelf areas[END_REF] proposed delimitations based on the zoning of [START_REF] Sealey | Setting geographic priorities for marine conservation in Latin America and the Caribbean. The Nature Conservancy[END_REF] in South America, which assumes the northern limit of the Falkland Current as the boundary between WTSA and TSA provinces, in Cabo Frio. However, upwelling of the South Atlantic Central Water (SACW) in this same location can extend environmental conditions (e.g. lower sea surface temperature) northward, and make Santana more similar to the southern colonies, regarding environmental conditions and species composition [START_REF] Soares | Seasonal variability of δ 13 C and δ 15 N of fish and squid in the Cabo Frio upwelling system of the southwestern Atlantic[END_REF][START_REF] Bauer | Marine fishes (Elasmobranchii and Teleostei) from the Santana Archipelago, a marine protected area in the southwestern Atlantic[END_REF]. Therefore, the biogeographic models which place Santana in the same zoning as the southern colonies are better supported by dietary variation, given that primary productivity distribution is the criteria used by [START_REF] Longhurst | An estimate of global primary production in the ocean from satellite radiometer data[END_REF], while LME are based on bathymetry, ocean currents, and fish community [START_REF] Sherman | Sustainability, biomass yields, and health of coastal ecosystems: an ecological perspective[END_REF]. Considering that diet composition may represent local biodiversity, these findings can contribute to the refinement of distinct marine zonings, helping to identify inconsistencies and adjust boundaries.

Finding spatial distribution patterns and associations with biogeographic models in highly mobile species of wide home range is challenging, but meta-analysis approaches can enable representative sampling in time and space, despite the limitations they may offer in analytical and interpretation terms. Publications and databases used in the present study comprised variable sampling efforts, periods, sex ratios, and researchers, which could make results biased regarding prey abundance and richness. In addition, the comparison of isotopic values requires caution, as the study covers colonies distributed over a 28° latitudinal gradient and therefore include potential differences in isotopic baselines. Temporal variation in sampling the colonies could lead to variations in isotopic information given the dynamics of primary productivity and prey availability around the colonies (Mancini et al. 2014). However, the heterogeneity of the sample for each colony also has an important positive aspect, since the representativeness of the dietary variations in the colony is greater compared to a point-in-time sampling. In this context, the framework applied in this study was successful in integrating information about the prey species consumed by the brown booby in a vast study area and was sufficient to detect spatial patterns of the brown booby diet along the southwest Atlantic Ocean. This allowed testing for fitting into pre-established marine biogeographic zonings and demonstrated the suitability of using top predators as samplers of fish fauna at large biogeographical scales. Additionally, as central-place foragers, seabirds use extensive areas around colonies to obtain food and breed in aggregations, which facilitates sampling. Spontaneous regurgitation increases the potential of seabirds as a source of information about the pelagic and demersal species community because regurgitates provide accurate taxonomic information about prey. Finally, the use of complementary approaches, such as regurgitates and stable isotopes, can take a more accurate picture of the dietary aspects and therefore of the prey community. Thus, studying the seabird diet is an important and useful tool to define and refine marine biogeographic zoning, especially from the perspective of widely distributed and hard-tosample organisms, such as flying fish.

Fig. 1 a

 1 Fig. 1 a Distribution of brown booby Sula leucogaster colonies sampled along the southwest Atlantic Ocean. SPSP = São Pedro e São Paulo. b Biogeographic boundaries of Large Marine Ecosystems, following Sherman (1994) (bottom left): EBS = East Brazilian Shelf, and SBS = South Brazilian Shelf; Spalding's Provinces (initials) and Ecoregions (numbers), following Spalding et al. (2007) (bottom middle): TSA = Tropical Southwestern Atlantic, WTSA = Warm Temperate Southwestern Atlantic; Longhurst's Provinces, following Lon-

n

  = total samples by each study; m/y = month/year of sampling; years separated by bars indicate isolated years of sampling and years separated by a dash indicate continuous sampling over the years. SPSP = São Pedro e São Paulo *References of 'unpbl data' are in

Fig. 2 a

 2 Fig. 2 a Species richness and b number of prey families consumed by Sula leucogaster in colonies along the southwest Atlantic Ocean. The colonies are ordered in a decreasing latitudinal gradient (from left

Fig. 3

 3 Fig. 3 Shannon diversity index of prey species consumed by Sula leucogaster in the southwest Atlantic Ocean considering colony location with a Large Marine Ecosystems (LME); b Longhurst's Provinces; c Spalding's Provinces; and d Spalding's Ecoregions. SBS15 = South Brazilian Shelf; EBS16 = East Brazilian Shelf;

Fig. 4

 4 Fig. 4 Two-dimensional NMDS plots considering a binary matrix of the prey community consumed by Sula leucogaster in relation to the framework of each colony in the following biogeographic zonings: a Large Marine Ecosystems; b Longhurst's Provinces; c Spalding's Provinces; d Spalding's Ecoregions. EBS16 = East Brazilian Shelf;

Fig. 5

 5 Fig. 5 Bayesian ellipses generated from carbon (δ 13 C) and nitrogen (δ 15 N) isotopic ratios from blood samples of Sula leucogaster obtained in colonies along the southwest Atlantic Ocean. Ellipses comprise 95% of the data and resulted in the following areas: Abrolhos = 9.58; Cagarras = 2.69; SPSP = 2.68; Santana = 1.94; Rocas = 1.5; Moleques = 1.28; Noronha = 0.89. SPSP = São Pedro e São Paulo

Table S1

 S1

	Colony	Regurgitates		Stable isotopes
		n	Data sources	Sampling years	n	Data sources	m/y
	Moleques	257	Kohlrausch, AB (unpbl data*)	1997/2002	18	Nunes and Bugoni (2018)	Feb/2014
		46	Branco et al. (2005)	2002-2004
	Currais	43	Krul (2004)	1995-1996	0	-	-
	Cagarras	25	Torres et al. (unpbl data*)	2012	27	Nunes and Bugoni (2018)	Dec/2014
	Cabo Frio	81	Coelho et al. (2004)	1990	0	-	-
	Santana	54	Mancini et al. (2023)	2017-2018	40	Bighetti et al. (2021)	Feb-Nov/2018
					62	Mancini et al. (2023)	Jun/2017 to May/2018
	Abrolhos	17	Alves et al. (2004)	1990-91/1994-96	31	Mancini et al. (2014)	Feb/2011
		19	Mancini, PL (unpbl data*)	2011	12	Nunes et al. (2022)	Feb/2019
		5	Rede Rio Doce Mar (unpbl data*)	2018-2019
	Noronha	15	Mancini, PL (unpbl data*)	2011	27	Mancini et al. (2014)	Mar/2011
		5	Bertrand, SL (unpbl data*)	2019
	Rocas	92	Kohlrausch, AB (unbpl data*)	2000-2002	28	Mancini et al. (2014)	Sep/2010
		20	Mancini, PL (unpbl data*)	2010
	SPSP	93	Kohlrausch, AB (unpbl data*)	2000-2002	29	Mancini et al. (2014)	Aug/2010 and Aug/2011
		33	Mancini and Bugoni (2014)	2010-2011	98	Nunes et al. (2018)	Jul/2015
		72	Nunes et al. (2018)	2014-2015

Acknowledgements

We are grateful to many researchers who studied the diet of brown boobies during the last 30 years in the Atlantic Ocean. We also thank E.E.M. Valim, T.P.X. Nascimento, L.R.M. Porto, G.P. Bighetti, F.P. Marques, C. Campolina, C. Barbraud, and K. Delord for their support in fieldwork. F.L. Rodrigues provided important contribution in prey identification and classification. We also thank the two anonymous reviewers for the comments and suggestions that contribute to improving the manuscript.

Author contributions GTN and JJ conceptualized the idea and designed the study. GTN, PLM, SLB and LB collected and shared data. GTN and JJ analyzed the data and wrote the manuscript and PLM, SLB, MAE, and LB contributed critically to the final revision.

Funding

The unpublished datasets were obtained with funding from the IRD Young Associated Team Program (JEAI TABASCO), the IRD mixed international laboratory program (LMI TAPIOCA), CPER Celimer (France), Fundação Boticário (Brazil), the TRIATLAS project (European Union's Horizon 2020 research and innovation program, Grant No. 817578), the PADDLE project (European Union's Horizon 2020 RISE program (Grant No. 734271), and the "Projeto Costões Rochosos", which is an environmental offset measure established through a Consent Decree/Conduct Adjustment Agreement between Petrorio and the Brazilian Ministry for the Environment, with the Brazilian Biodiversity Fund (FUNBIO) as an implementer. PLM received a PNPD scholarship from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES) Finance Code 001. LB had been funded by National Research Council (CNPq; Grant No. 405497/2012-1) and is a fellowship from CNPq (No. 311409/2018-0). GTN is funded by CNPq (Grant No 443328/2019-6).

Data availability All data and the R script are publicly available through the GitHub repository: https:// github. com/ Seabi rdEco logyU FRGS/ DietB iogeo graphy

Declarations

Conflict of interest The authors declare no conflict of interest.

Ethical approval

The sampling of the unpublished datasets was approved by SISBIO 64234-6 and by the ethics committee at UFRGS.


Didier Lasseux

Francisco J Valdés-Parada

Upscaled dynamic capillary pressure for two-phase flow in porous media

Keywords:

A closed expression for the average pressure difference (often called the macroscopic dynamic capillary pressure in the literature) is proposed for two-phase, Newtonian, incompressible, isothermal and creeping flow in homogeneous porous media. This upscaled equation complements the average equations for mass and momentum transport derived in a previous article. Consistently with this work, the expression is derived employing a simplified version of the volume-averaging method that makes use of elements of the adjoint method and Green's formula. The resulting equation for the average pressure difference is novel, as it shows that this quantity is controlled by the pressure gradient (and body forces) in each phase, as well as interfacial effects, and is applicable to situations in which the fluid-fluid interface is not necessarily at its steady position. The effective-medium quantities associated with the sources are all obtained from the solution of a single adjoint (or closure) problem to be solved on a (periodic) unit cell representative of the process. The average pressure difference predicted by the derived expression is validated through excellent comparisons with direct numerical simulations performed in a model porous structure.

Introduction

Modelling two-phase flow in porous media has been a long-lasting subject of research, and remains a challenging task. Over the past two decades, developments in image acquisition and processing (see e.g. [START_REF] Singh | Capillary-dominated fluid displacement in porous media[END_REF], and references therein) as well as progress in direct numerical simulations (DNS) at the pore scale (e.g. [START_REF] Konangi | Comparison of pore-scale capillary pressure to macroscale capillary pressure using direct numerical simulations of drainage under dynamic and quasi-static conditions[END_REF] have been proved to be useful in providing valuable pore-scale information that confirms, and sometimes questions, existing theoretical macroscopic models. Nevertheless, mathematical modelling of two-phase flow in porous media remains of major importance.

A traditional approach for slow, incompressible and immiscible two-phase flow, originating from the pioneering work of [START_REF] Muskat | The flow of heterogeneous fluids through porous media[END_REF], consists in using the so-called empirical generalized Darcy's equations for momentum transport, together with the macroscopic mass conservation equations. However, to close the macroscopic model, an additional relationship is required [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF]. On the basis of the assumption that the flow is slow enough to be considered as a succession of equilibria, this closure relationship is taken as the macroscopic pressure difference, often called the macroscopic capillary pressure. This terminology was adopted due to the quasi-equilibrium assumption, so that the pore-scale pressure jump at the fluid-fluid interface is compensated only by capillary effects together with the fact that the pressure in each phase is quasi-uniform and is thus assumed to be equal to that at the fluid-fluid interface. In addition, since the macroscopic flow description is given in terms of the local (wetting-phase) saturation, the macroscopic capillary pressure has been considered as a function of this unique variable.

The shortcomings of such a heuristic description have been widely discussed in the literature over the past four decades. First, the fact that, due to dynamic effects, thedifference between thevolume-averaged pressures should be distinguished from the interfacial average of the equilibrium pressure difference has been pointed out. In particular, [START_REF] Gray | Non-hysteretic functional form of capillary pressure in porous media[END_REF] emphasized that, under dynamic conditions, the macroscopic capillary pressure should be evaluated in terms of the interfacial average of the pressure difference between the two fluids, although, at equilibrium, it corresponds to the volume-average pressure difference. This point of view was shared by [START_REF] Starnoni | On the concept of macroscopic capillary pressure in two-phase porous media flow[END_REF], who derived a macroscale Laplace-type equation using volume averaging. Their result is similar to that anticipated by [START_REF] Whitaker | The closure problem for two-phase flow in homogeneous porous media[END_REF]. Second, the unique dependence of the (equilibrium) macroscopic capillary pressure on the saturation has also been questioned.

On a thermodynamic basis, [START_REF] Hassanizadeh | Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries[END_REF][START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF] developed a relationship for the macroscopic pressure difference that differs from the macroscopic capillary pressure by a dynamic term. Assuming a linear theory, this term was proposed to be proportional to the time rate of variation of the saturation as anticipated by [START_REF] Marle | From the pore scale to the macroscopic scale: equations governing multiphase fluid flow through porous media[END_REF]. However, a closure relationship to estimate the corresponding proportionality coefficient was not provided. Furthermore, these authors pinpointed the dependence of the (equilibrium) macroscopic capillary pressure on many parameters, which, under classical flow conditions, can be restricted to the fluid-fluid interfacial area and the saturation. This allowed them to conclude that the equilibrium macroscopic capillary pressure should be viewed as a unique surface parametrized by the interfacial specific area and saturation,and that the hysteretic dependence of the capillary pressure on the saturation, which has been widely reported, results from a projection of this surface onto the saturation axis.

Experimental measurements in transparent micromodels of the equilibrium macroscopic pressure were carried out by [START_REF] Pyrak-Nolte | Relating capillary pressure to interfacial areas[END_REF] to investigate functionalities proposed in the above two references. Their results highlighted the fact that the terms of the capillary pressure related to the interfacial area were nearly constant for any equilibrium situation, and, consequently, that terms associated with the phase-free energy variations with respect to the saturation and related to saturation gradients play a very important role. Under dynamic conditions, Joekar-Niasar, [START_REF] Joekar-Niasar | Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling[END_REF] and Joekar-Niasar & Hassanizadeh (2012) used a three-dimensional pore-network modelling on a regular lattice and confirmed that a classical capillary pressure function of the saturation was inappropriate. They concluded that a unique surface relating the capillary pressure to the saturation and interfacial area was in fair agreement with the observations, although the dependence upon the capillary number was not fully captured.

Using the thermodynamically constrained averaging theory, [START_REF] Gray | TCAT analysis of capillary pressure in non-equilibrium, two-fluid-phase, porous medium systems[END_REF] proposed an expression for the capillary pressure evolution to equilibrium that incorporates the rate of change of the specific fluid-fluid interfacial area. Since no closure scheme was provided to determine the weighting coefficients of both the rate of change of saturation and interfacial specific area, their determination was left to experimental investigation.

In a recent development proposed by Lasseux & Valdés-Parada (2022), a formal derivation of the macroscopic momentum and mass conservation equations for isothermal, creeping Newtonian flow of two immiscible and incompressible phases, β and γ , in homogeneous porous media, was obtained, which reads (see the definitions of the different quantities in § 2)

v α α = - K * αα μ α • ∇P α - K * ακ μ κ • ∇P κ + 1 μ α V A βγ (2γ Hn βγ + ∇ s γ ) • D αα dA, (1.1a) ∂ε α ∂t + ∇ • v α α = 0, (1.1b)
with α, κ = β, γ and κ / = α. In the above momentum equation, the (dominant and coupling) permeability tensors are defined as K * ακ = D ακ α , where D ακ are closure variables that solve two ancillary closure problems (see (5.14) and (5.15) in Lasseux & Valdés-Parada (2022)).

For completeness, it is of major interest to derive a closed upscaled equation for the average pressure difference between the two fluid phases (which, for the sake of semantic exactness, is not referred to as the dynamic capillary pressure). This is the purpose of the present companion work developed in the next sections as follows. In § 2, the pore-scale flow problem is formulated for a two-phase, Newtonian, incompressible, isothermal and creeping flow. Assuming the existence of an elementary volume representative of the process, the formulation is provided in the corresponding periodic unit cell. Using a simplified version of the volume-averaging method and elements of the adjoint method, an adjoint problem is proposed that is further combined with the physical pore-scale problem with the aid of Green's formula. The procedure, which is free of any simplifying assumption and yields a closed equation for the average pressure difference, is detailed in § 3. Section 4 is dedicated to first evaluate the effective-medium quantities resulting from the adjoint (or closure) problem solution and, second, to a validation through comparison between DNS of the pore-scale flow problem and predictions from the upscaled equation in a model two-dimensional configuration. Conclusions are drawn in § 5.

Pore-scale flow model

The development starts by considering the incompressible, Newtonian, isothermal and creeping flow of two immiscible fluid phases β and γ (indifferently denoted the α-phase in the following) through a rigid and homogeneous porous medium, the solid skeleton of which is the σ -phase. No mass transport is assumed to take place at the fluid-fluid, A βγ , and solid-fluid, A ασ , interfaces. Both fluid phases are assumed to occupy a connected region of the pore space, and the analysis is focused on a periodic unit cell, which is assumed to be representative of the process and implies that both A βγ and A ασ are periodic. The flow is driven by volume forces per unit mass, b α , in each phase and by macroscopic pressure gradients applied at the outer boundaries of the unit cell, ∇ p α α , all being considered constant in the remainder of this work. The pore-scale velocities can therefore be considered as periodic fields. The flow problem, in dimension N of space, can be written in a unit cell as (α = β, γ )

∇ • v α = 0, in V α ,
(2.1a)

0 = ∇ • T p α + ρ α b α , in V α , (2.1b) v β = v γ , at A βγ , (2.1c
)

n βγ • (T p β -T p γ ) = 2Hγ n βγ + ∇ s γ, at A βγ , (2.1d) v α = 0, at A ασ , (2.1e) v α S - αi = v α S + αi , i = 1, . . . , N, (2.1f ) (n α • T p α ) S - αi = -(n α • T p α ) S + αi -l i • ∇ p α α n α S + αi , i = 1, . . . , N, (2.1g) p α = 0, at r 0 α .
(2.1h)

In the above equations, V α is the space (of measure V α ) occupied by the α-phase within the unit cell of volume V. In addition, ρ α , v α and T p α = -p α I + μ α (∇v α + ∇v T α ), respectively, represent the fluid density, velocity vector and total stress tensor in each phase, with p α and μ α , respectively, denoting the pore-scale α-phase pressure and dynamic (constant) viscosity, whereas I is the identity tensor. In (2.1d), the stress jump is due to capillary effects, including a possible surface gradient, ∇ s γ , of the interfacial tension, γ , along A βγ at which the mean curvature and unit normal vector, directed from the β-phase towards the γ -phase, are respectively denoted by H and n βγ . Note that (2.1h) is required for the problem to be well posed. Moreover, the flow is supposed to be free of acceleration, yielding a steady version of the momentum equation (2.1b). Nevertheless, flow unsteadiness can arise due to the displacement of the fluid-fluid interfaces, with or without any change of the fluids' volume fraction, which, in the former case, justifies the unsteady character of the macroscopic mass balance equation (1.1b) and, in the latter, yields a flow that is macroscopically steady.

In (2.1f ) and (2.1g), S - αi and S + αi represent the unit cell outer boundaries, so that

r α S + αi = r α S - αi + l i , where l i is the unit cell lattice vector in the ith direction (α = β, γ ).
In the rest of this work, r α is used to locate a point within the α-phase with respect to a fixed coordinate system. This vector can be decomposed as r α = x α + z α , with x α locating the α-phase barycentre and z α locating the same point as r α , albeit with respect to x α . Furthermore, the intrinsic average, ψ α α , of a pore-scale quantity, ψ α , locates the resulting average at the phase barycentre, i.e.

1 V α V α ψ α dV = ψ α α x α , α = β, γ . (2.2)
Note that x α = r α α and, consequently, z α α = 0. In addition, the superficial average, ψ α α , of ψ α is defined as ψ α α = ε α ψ α α , where ε α = V α /V is the α-phase volume fraction. Following [START_REF] Gray | A derivation of the equations for multi-phase transport[END_REF], ψ α can be decomposed in terms of its intrinsic average and spatial deviations, ψα , according to

ψ α r α = ψα r α + ψ α α r α , α = β, γ .
(2.3) Note that ψα is l i -periodic. Performing a Taylor-series expansion about x α of the above decomposition and averaging the result yields (α = β, γ )

ψ α α x α = ψ α α x α + z α • ∇ ψ α α x α + 1 2 z α z α : ∇∇ ψ α α x α + • • • α x α + ψα α x α . (2.4) Since ψ α α x α
and its successive gradients are constant within the unit cell, they can be taken out of the first average on the right-hand side of this expression, and making use of the order-of-magnitude estimate 1 2 z α z α : ∇∇ ψ α

α x α = O(r 2 0 /L 2 ψ α α )
, where r 0 and L are the characteristic sizes of the averaging and macroscopic domains, respectively, it can be concluded from the above expression that ψα α 0, as long as r 2 0 L 2 . This average constraint serves the purpose of bounding the fields of ψ and is necessary to make the problem well-posed, as shown below.

The decomposition defined in (2.3) can be applied to the pore-scale pressure in both phases, so that (2.1b), (2.1d), (2.1g) and (2.1h), respectively, can be replaced by (α = β, γ )

0 = ∇ • T pα -∇P α , in V α ,
(2.5a)

n βγ • T pβ = n βγ • T pγ + n βγ ( p β β -p γ γ ) r βγ + 2Hγ n βγ + ∇ s γ, at A βγ , (2.5b) (n α • T pα ) S - αi = -(n α • T pα ) S + αi , i = 1, . . . , N, (2.5c) pα α = 0. (2.5d)
Here, T pα = -p α I + μ α (∇v α + ∇v T α ) and ∇P α = ∇ p α α -ρ α b α , yielding the formulation used in Lasseux & Valdés-Parada (2022) in which b α was taken as the gravitational acceleration. In (2.5b), r βγ ≡ r β = r γ is used to locate points at A βγ . Note that the average constraint expressed in (2.5d) replaces (2.1h). In the following, ∇P α is treated as a constant, and a formalism similar to that reported in the above-cited reference is employed to derive an expression for the average pressure difference between the two fluid phases.

Derivation of the average pressure difference

At this point, it is convenient to introduce a new pair of closure variables in each phase, f α and f α , which solve the following adjoint (or closure) boundary-value problem in a periodic unit cell (α = β, γ ):

∇ • f β = 1 V β , in V β , ∇ • f γ = - 1 V γ , in V γ , (3.1a,b) 0 = ∇ • T f α , in V α , (3.1c) f β = f γ , n βγ • T f β = n βγ • T f γ , at A βγ , (3.1d,e) f α = 0, at A ασ , (3.1f ) f α S - αi = f α S + αi , (n α • T f α ) S - αi = -(n α • T f α ) S + αi , i = 1, . . . , N, (3.1g,h) f α = 0, at r 0 α . (3.1i)
In these equations,

T f α = -f α I + μ α (∇ f α + ∇ f T α )
is a stress-like second-order tensor, whereas, in (3.1i), r 0 α locates an arbitrary point in the α-phase (α = β, γ ). This new closure problem can be combined with the flow problem in a unit cell (namely (2.1a), (2.1c), (2.5), (2.1e) and (2.1f )) by means of a Green's formula. Taking two arbitrary scalar fields, a and q, and two arbitrary vector fields, a and q, all of them defined in the α-phase and having sufficient regularity, this Green's formula can be written as (see the derivation in Appendix A of Sánchez-Vargas, Valdés-Parada & Lasseux (2022))

V α [a • (∇ • T q ) -(∇ • T a ) • q -q∇ • a + a∇ • q] dV = A α [a • (n • T q ) -n • T a • q] dA.
(3.2)

Here A α represents the enclosing surfaces of V α , whereas T a ≡ -aI + μ α (∇a + ∇a T ) and T q ≡ -qI + μ α (∇q + ∇q T ) (α = β, γ ). a ≡ pβ , a ≡ v β , q ≡ f β , q ≡ f β and μ α ≡ μ β in the above equation, and substituting the corresponding differential equations and boundary conditions, together with the constraint (2.5d), leads to the following expression:

- V β f β dV • ∇P β = A βγ [v γ • (n βγ • T f γ ) -(n βγ • T pγ ) • f γ ] dA - A βγ ( p β β -p γ γ ) r βγ n βγ • f β dA - A βγ (2Hγ n βγ + ∇ s γ ) • f β dA. (3.3)
Note that ∇P β is taken out of the integral on the left-hand side, as it is assumed to be constant. However, this is not the case for the average pressure difference inside the second area integral on the right-hand side of this equation, as both average pressures are considered at points r βγ at A βγ .

To progress towards an expression for the average pressure difference, a Taylor-series expansion of p α α | r βγ can be performed about x α . Maintaining the first two terms of this expansion (the remaining ones are indeed zero if ∇ p α α is constant), and recalling that ∇P β = ∇ p β β -ρ β b β , the above equation takes the following form:

- V β f β dV • (∇ p β β -ρ β b β ) = A βγ [v γ • (n βγ • T f γ ) -(n βγ • T pγ ) • f γ ] dA -( p β β x β -p γ γ x γ ) A βγ n βγ • f β dA - A βγ n βγ • f β z β dA • ∇ p β β + A βγ n βγ • f β z γ dA • ∇ p γ γ - A βγ (2Hγ n βγ + ∇ s γ ) • f β dA. (3.4)
Note that in the area integrals on the right-hand side of this expression, z α = r βγx α (α = β, γ ).

Dynamic capillary pressure in porous media two-phase flow

To simplify the above equation, it is pertinent to note that

A βγ n βγ • f β dA = V β ∇ • f β dV = 1.
In addition, the following nomenclature is proposed for the effective-medium coefficients (α = β, γ and α / = κ)

ϕ α = V α f α dV, ψ α = ϕ α - A βγ n ακ • f α z α dA, (3.5a,b) s βγ = A βγ (2Hγ n βγ + ∇ s γ ) • f β dA, (3.5c)
and for the average pressure difference

P = p γ γ x γ -p β β x β . (3.6)
After rearranging, (3.4) can be written in the following more compact form:

P = - A βγ [v γ • (n βγ • T f γ ) -n βγ • T pγ • f γ ] dA -ψ β • ∇ p β β + ρ β b β • ϕ β -(ψ γ -ϕ γ ) • ∇ p γ γ + s βγ . (3.7)
The first area integral in (3.7) can be expressed in an alternative form by considering again Green's formula given in (3.2) with a ≡ pγ , a ≡ v γ , q ≡ f γ , q ≡ f γ and μ α ≡ μ γ .

After substitution of the corresponding differential equations and boundary conditions, the resulting expression is

A βγ [v γ • (n βγ • T f γ ) -(n βγ • T pγ ) • f γ ] dA = ϕ γ • ∇P γ . (3.8)
Substituting this result into (3.7) yields the following expression for the average pressure difference:

P = -ψ β • ∇ p β β + ρ β b β • ϕ β -ψ γ • ∇ p γ γ + ρ γ b γ • ϕ γ + s βγ . (3.9a)
An alternative equivalent form can be derived by combining the adjoint problem defined in (3.1) with the flow problem given in (2.1) by means of Green's formula (3.2). Employing the same procedure as above results in the following expression:

P = N i=1 S + βi n β • f β dA l i • ∇ p β β + S + γ i n γ • f γ dA l i • ∇ p γ γ + ρ β b β • ϕ β + ρ γ b γ • ϕ γ + s βγ . (3.9b) Using the identity V α ∇ • ( f α z α ) dV = A α n α • f α z α dA = V α ∇ • f α z α dV + V α f α dV
(note that ∇z α = I is employed here), taking into account (3.1a,b) and boundary conditions for f α , it can be readily shown that

N i=1 S + αi n α • f α dA l i = -ψ α (α = β, γ
), proving that (3.9b) coincides with (3.9a). Note that the average pressure difference can also be expressed at the centroid, x, of the unit cell by using the relationship

p α α x = p α α x α -y α α • ∇ p α α
, with the definition r α = x + y α . Indeed, equality of the two pressure gradients in both phases is a compatibility requirement with the assumption of a periodic unit cell representative of the process (see details in Appendix B in Lasseux & Valdés-Parada (2022)). The effective-medium coefficient associated with the pressure gradient in that case can be shown to be -ϕ β -ϕ γ + x γx β .

The expressions for P given in (3.9), which are the important result of this work, only rely on the assumptions made for the pore-scale flow problem together with scale separation and the possibility of a local flow description in a representative periodic unit cell. These assumptions may appear to be overly severe and were adopted with the intention of deriving a closure problem for the computation of the effective-medium quantities involved in (3.9). Nevertheless, the upscaled expressions may be used even in situations when these assumptions are not completely met. Note that the derivation does not require A βγ to be stationary (i.e. motionless fluid-fluid interfaces), or a time-independent saturation, as explained in § 2. The first four terms on the right-hand side of (3.9a) represent the influence of the hydrodynamic forcing in each phase, whereas the last term represents the influence of interfacial tension. The latter results from a complex interfacial average of the Laplace effects (i.e. 2Hγ ) and interfacial tension gradient effects at the pore scale, and involves f β (= f γ ) at A βγ as a weighting function.

Obviously, P depends not only on the wetting-phase saturation, but also on the configuration of A βγ as well as on the dynamic effects, as anticipated, for instance,by [START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF] and later by [START_REF] Gray | TCAT analysis of capillary pressure in non-equilibrium, two-fluid-phase, porous medium systems[END_REF], but yet in an unclosed form. This is clearly different from the Laplace equation typically found in the literature [START_REF] Whitaker | The closure problem for two-phase flow in homogeneous porous media[END_REF]). This result is only recovered in the absence of any macroscopic forcing. Indeed, in that case, 2H must be constant, and, assuming that γ is also constant, (3.9a) reduces to the classical macroscopic Laplace relation, i.e. P ≡ P cL = 2Hγ . The vectors ψ α and ϕ α can be conceived as effective characteristic lengths for, respectively, ∇ p α α and ρ α b α .

The effective-medium coefficients, ψ α , ϕ α and s βγ , can be predicted from the solution of the associated closure problem defined in (3.1). The equation for P is compatible with an expression proposed in the literature (see (43) in [START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF]). It has the advantage of being closed and it allows identifying and evaluating the role of each source. An analysis of the prediction from this upscaled equation is proposed in the next section.

Results

In this section, numerical simulation results are reported for a two-dimensional (N = 2) two-phase flow configuration schematically represented in figure 1 and already employed by Lasseux & Valdés-Parada (2022) to investigate the macroscopic momentum equations. Results are presented for b α = 0, ∇ s γ = 0 and a prescribed pressure gradient in each phase given by ∂ p α α /∂x =-h and ∂ p α α /∂y = 0(α = β, γ ). Moreover, the flow is considered to be fully steady, i.e. with A βγ at its stationary position for the sake of brevity in presentation.

Flow and closure problem solutions

The first step in the analysis is the solution of the pore-scale flow problem, i.e. DNS and the computation of the average pressure difference. The boundary integral element method (BIEM) employed for DNS and results on the configuration under concern are detailed in Lasseux & Valdés-Parada (2022) (see Appendix C for the BIEM and § 6 for the results). Accordingly, the same dimensionless quantities as in this reference are employed, i.e.

r * α = r α / , ∇ * = ∇, v * α = μ γ v α /(h 2 ) and p * α = μ γ p α /(μ α h ) (α = β, γ )
, where is the x and y unit cell size (see figure 1). The solution depends on the porosity,

ε = (V β + V γ )/V, wetting-phase saturation, S β = ε β /ε = V β /(V β + V γ )
, capillary number, Ca = h 2 /γ , and viscosity ratio, μ * = μ β /μ γ . Once the steady shape of A βγ is determined, the pressure fields are computed (see (C10) in Lasseux & Valdés-Parada ( 2022)), together with their intrinsic average, p α α x α , and their difference, yielding the DNS result for this quantity.

In the second step, the closure problem is solved taking the shape of A βγ determined from DNS in the unit cell. Using the following dimensionless changes of variables, , from the upscaled equation given by (3.9b). Under the conditions employed here, it reduces to

d * α = ( 3 /μ α )f α , d * β = 2 ( f β -r β /NV β ) = f * β -r * β /Nε β and d * γ = 2 ( f γ + r γ /NV γ ) = f * γ + r * γ /
P * = - S + γ 1 f * γ x dA * - S + β1 f * βx dA * + s * βγ , (4.1)
where

s * βγ = (1/Ca) A βγ 2H * n βγ • f * β dA * , f * αx is the x component of f * α (α = β, γ ) and H * = H.
Numerical simulations are performed for a range of S β with ε = 0.8, Ca = 0.1 and 10, and μ * = 0.1 and 10. The numerical results obtained with the BIEM were verified to match those obtained using the finite-element method with Comsol Multiphysics 6.0.

Validation with numerical simulations

Before focusing on the predictions of P * , it is pertinent to examine the dimensionless pressure and closure variable pore-scale fields, the boundary values of the latter providing the effective-medium quantities involved in (4.1). An illustration of p * α as well as of the x and y components of f * f * α , are also evident. Note that, due to the non-symmetry of A βγ with respect to the y-axis, no y symmetry is observed on these fields.

A comparison of P * obtained from DNS and predicted from the upscaled equation is represented versus S β in figure 3. In all cases, an excellent agreement between the two approaches is obtained over the whole range of S β . Indeed, the relative error, taking the DNS result as the reference, is less than 0.8 % among all cases. This confirms the pertinence of the expression for P derived here. In figure 3, the classical macroscopic Laplace capillary pressure is also reported. The dimensionless form of this quantity is given by P * cL = 2 H * βγ /Ca, where H * βγ = (1/A * βγ ) A βγ H * dA * , and A * βγ is the dimensionless measure of A βγ . From these results, it must be emphasized that P * cL cannot represent a reasonable approximation of P * , as it is always much smaller, by at least five orders of magnitude, than P * . The dominant contribution to P * is from s * βγ , which, however, overpredicts it. Indeed, as can be observed from figure 3, the two hydrodynamic terms remain significant, except in some specific circumstances, i.e. for S β 0.6 for μ * = 0.1 and Ca = 10 (figure 3a) or S β 0.4 for μ * = 0.1 and Ca = 0.1 (figure 3b).

It was also numerically verified that the upscaled equation prediction of P * perfectly matches the DNS results for other shapes of A βγ , different from the stationary one that is characterized by a zero normal interfacial velocity. Moreover, although the flow under consideration here is macroscopically steady (∂S β /∂t = 0), it clearly appears that P does not reduce to the equilibrium capillary pressure remaining as the unique term if the relationship proposed by [START_REF] Hassanizadeh | Thermodynamic basis of capillary pressure in porous media[END_REF] is considered. This was suggested by [START_REF] Gray | TCAT analysis of capillary pressure in non-equilibrium, two-fluid-phase, porous medium systems[END_REF], who stated that this relationship 'does not account properly for disequilibrium due to the dynamics of interfacial rearrangement at constant saturation'; the upscaled equation derived here accounts for this mechanism.

Conclusions

In this work, an upscaled equation was derived to predict the average pressure difference, P, between two immiscible fluid phases in homogeneous porous media using a simplified version of the volume-averaging method, combined with elements of the adjoint method and Green's formula. The resulting expression shows that P includes contributions from surface tension effects, and also from macroscopic pressure gradients and volume forces acting in each phase. The derivation of the expression for P is not restricted to situations in which the fluid-fluid interface is at its stationary position. The effective-medium quantities involved in this expression are predicted from the solution of a single associated closure (adjoint) problem in a representative periodic unit cell. Predictions of P from its expression derived here were validated through comparisons with DNS, showing excellent agreement. It is shown that, for the particular porous structure and phases distribution considered here, the term related to capillary effects is dominant in the prediction of the average pressure difference (especially as the wetting-phase saturation is closer to unity), albeit the influence of the macroscopic pressure gradients cannot be disregarded.

Imaging techniques may be envisaged to obtain fluid-phase distributions in a representative domain as an input for the closure problem solution and to further investigate real configurations. Alternatively, an exhaustive numerical analysis of porous skeleton topologies and fluid-phase distributions, which may be seen as representative of a real specific configuration of interest, could be performed to propose correlations that could be implemented in standard two-phase flow simulators of common use. This would avoid the dependence on the solution of closure problems and allow one to express the upscaled equation for P in terms of macroscopic parameters only.

The average pressure difference derived in this work completes the macroscopic mass and momentum equations previously reported by Lasseux & Valdés-Parada (2022). Finally, the analysis provided here opens the way to study more complex flow conditions that will be addressed in future works.

Figure 1 .

 1 Figure 1. Two-dimensional two-phase flow configuration considered for the numerical simulations, made of a square pattern of parallel cylinders of circular cross-section (a) and corresponding periodic unit cell (b). The flow results from a constant horizontal macroscopic pressure gradient in each phase. The β-phase flows as a wetting layer attached to each cylinder row, whereas the non-wetting γ -phase flows in the intermediate layers.

  Nε γ , and noticing that ∇ • r α = N and ∇r α = I, d * α and d * α satisfy a homogeneous incompressible Stokes-like problem along with the corresponding boundary conditions. Its solution is obtained with the same BIEM as for DNS. The solution on d * α (and hence on f * α ) at the boundaries is then used to predict the dimensionless average pressure difference, P * = p *

Figure 2 .

 2 Figure 2. (a) Dimensionless pressure field resulting from DNS at steady state; (b) x component and (c) y component of the dimensionless closure variable f * α (α = β, γ ). Corresponding isolines are represented in white. Black lines materialize A βσ and A βγ . Here ε = 0.8, Ca = 10, μ * = 10 and S β = 0.5.

Figure 3 .

 3 Figure 3. Dimensionless average pressure difference, P * , obtained from DNS and predicted from the upscaled equation (UE), along with the interfacial term contribution, s * βγ (see (4.1)), on the left-hand vertical axes, and classical Laplace capillary pressure, P * cL (see text), on the right-hand vertical axes, versus the wetting-phase saturation, S β , for: ε = 0.8; and (a) μ * = 0.1, Ca = 10; (b) μ * = 0.1, Ca = 0.1; (c) μ * = 10, Ca = 10; (d) μ * = 10, Ca = 0.1.

Funding. Part of this work was performed while both authors were visiting professors at the University of Genoa (UniGe) in collaboration with Professor A. Bottaro. Financial support from UniGe is sincerely appreciated. Data of the numerical results presented in this work are available from the authors upon request.

Declaration of interests. The authors report no conflict of interest.

Author ORCIDs.

Didier Lasseux https://orcid.org/0000-0002-6080-8226; Francisco J. Valdés-Parada https://orcid.org/0000-0003-4208-1075.