# On the relevance of the Godot Engine in the indie game development industry

Julian Holfeld, M. Sc.
University of Kassel
<EMAIL_ADDRESS>
## Abstract

This paper examines the relevance of the Godot Engine in the indie game
industry. The Godot Engine is a relatively new game engine from 2014 and
competes with leading market players. To get to the bottom of its relevance,
two major online sales platforms and the game engines that are commonly used,
Steam and itch.io, are examined. Mainly, these findings are compared with
reference data from 2018. It turns out that the Godot engine has gained
massive relevance in 2020 and now seems to be one of the leading players in
the indie game industry. The exact causes are difficult to determine. However,
this paper provides many clues for further research in this area.

##  1 Introduction

In the realm of video game development, the methodology for creating games has
evolved significantly over time. As early as the 1980s, several companies
began using their own in-house engines. For instance, in the case of Super
Mario Bros (1985), Nintendo used a fast scrolling game engine that was
originally developed for Excitebike (1984) [1]. The concept of licensing game
engines to other companies as a standalone product has been around since the
mid-90’s. At that time, games like Doom (1993) were developed using the id
Tech game engine [2]. This engine was specifically designed for 3D shooters,
although this iteration was not yet a fully-fledged 3D game engine. Over the
years, an increasing number of game engines were developed to speed up the
development of video games.

It should be mentioned that the definition of a game engine differs in
literature. In this paper, game engines are seen as frameworks that accelerate
the development process and provide at least support for 2D or 3D rendering.
It is crucial to recognize that not only does the definition differ
significantly, but so does the use case of a game engine. A prime example is
the id Tech engine, which was developed explicitly for a particular genre of
video games. There are also game engines like Unity, which could be called a
general purpose game engine. Therefore, a comparison between game engines is
not easy. Since they always serve a specific purpose, there is no game engine
that stands above all others. This paper takes a closer look at the purpose of
video game development in the indie game industry, with a particular focus on
the Godot Engine. The decision to investigate this specific game engine is
driven by the insights derived from a few datasets, which will be expounded
upon later in this paper.

##  2 Godot Engine

The Godot Engine was released as open source under the MIT license on GitHub
on February 10 in 2014111The first commit containing the source code can be
found under the SHA:
0b806ee0fc9097fa7bda7ac0109191c9c5e0a1ac [3]. The first stable release was
published in December 2014 [4]. This makes it relatively new in comparison to
popular game engines like Unity and Unreal Engine, which were released in 2005
and 1998, respectively [5, 6]. Considering the different release dates, it
should be clear that due to the young age of the Godot Engine the ecosystem
that formed around this technology is still small in comparison to older
engines.

The Godot Engine supports cross-platform exports and facilitates the
development of 2D and 3D games [7]. Video games can be exported for PC, mobile
devices or web platforms. These games are written with a programming language
specially developed for the Godot Engine called GDScript. GDScript is a
dynamically typed programming language which is syntactically similar to
Python. However, the Godot Engine also officially supports C# and includes
GDNative, a technology for creating bindings to other languages [8]. These
officially include C and C++. Unofficially support for languages such as D,
Kotlin, Nim, Python and Rust have been developed by the community.

It’s noteworthy that the Godot Engine appears to attract not only smaller
companies but also larger players like Microsoft, Epic Games and Meta’s
Reality Labs. Microsoft donated $24,000 to the Godot team in 2017 for
implementing C# as a programming language [9]. In 2020, the Godot Engine
received the Epic MegaGrant of $250,000, which is distributed by Epic Games to
projects that expand the open source 3D graphics ecosystem [10]. This is
noteworthy as Epic Games develops the Unreal Engine, suggesting a potential
interest in supporting other game engines for mutual benefit. The adoption of
an MIT license, as in the case of Godot Engine, may facilitate the
incorporation of new technologies without negative implications for
intellectual property. In December 2020 and 2021 Godot Engine received
donations from Meta’s Reality Labs [11, 12]. The funds will be used
specifically for development in the field of extended reality (XR). There is
no public information available regarding the exact amount of the donation.
This was confirmed in a Reddit comment by project manager Rémi Verschelde
(akien-mga) [13]. Rémi also mentioned that the team actively approaches
companies rather than the other way around.

##  3 Methodology

In the last section it became evident that the Godot Engine has managed to
capture the attention of larger companies, suggesting a degree of relevance in
the video game industry. However, this does not directly address the central
question of how relevant the Godot Engine is for indie game developers. To
effectively answer this question, we must first establish a clear definition
of indie game development. Just as the definition of a game engine is elusive,
so is the definition of indie game development. The definition used here is
from O’Donnell, who states that indie games focus on a few goals and can be
developed on a smaller scale in terms of manpower [14]. An indie game
development team operates without financial backing from a larger company,
distinguishing it from AAA titles developed with substantial budgets and large
teams.

To investigate the relevance of the Godot Engine for indie game developers, a
comprehensive literature review was conducted. It is important to note that
this research was not limited to scholarly articles but also included gray
literature. This is because the Godot Engine is a relatively new game engine,
and there is a lack of scholarly articles on the topic. The gray literature
used in this paper was carefully selected to ensure reliable and trustworthy
information.

For this paper, the digital sales platforms Steam and itch.io were selected
for investigation. Firstly, it aligns with a study conducted by Toftedahl and
Engström that provides comparative data from 2018 [15]. Secondly, reference is
made to a survey named “State of the Game Industry” from 2023 encompassing
2,300 video game developers [16]. Notably, the focus in this paper lies on PC
and not on console or mobile games. This is because the mentioned survey from
2023 showed that 65% of the respondents are actively developing for PC and 57%
plan to do so in the next project as well. In 2019, a survey involving nearly
4,000 video game developers highlighted that, aside from publisher-owned and
individual websites, most games are primarily sold on Steam and itch.io [17].

Identifying the game engine used in games available on Steam presents a
significant challenge due to a lack of detailed information. In the paper by
Toftedahl and Engström mentioned above, a script was used that matches games
with articles on Wikipedia and subsequently assigns game engines. However,
this procedure presents challenges and gaps in the data sets. This is because
not all games provide information about the underlying game engine.

Alternatively, another method involves recognizing the game engines used by
analyzing filenames through pattern matching. The website SteamDB employs this
method to deduce the game engines utilized in development [18]. While this
approach provides valuable insights, it is not without its limitations,
including the potential for gaps or false positives due to ambiguous
filenames. In this paper, the datasets from SteamDB are used because matching
filenames to game engines yield more accurate results than matching the data
to Wikipedia articles. In addition, the different versions of the id Tech
engine were combined as one.

Analyzing the data from itch.io proves to be more straightforward because the
website provides a list of game engines used in development with the number of
games for each engine. This structured presentation of data facilitates ease
of access and analysis. This approach aligns with the methodology employed by
Toftedahl and Engström in their research. However, it is important to
acknowledge that this method captures a static snapshot of the data at a given
point in time. To address this limitation a script was developed that extracts
data from itch.io’s RSS feed in XML format, specifically recording the
publication date of each game [19]. This additional temporal dimension allows
for the creation of a time series, enabling a more in-depth analysis of used
game engines over time.

##  4 Results

When examining the presence of various game engines on Steam for 2018, 2022,
and 2023, it is important to see the significance of these results in the
context of indie game development. The following Table 1 presents the
percentages of game engines on Steam for these years.

Game Engine | 2018 | 2022 | 2023
---|---|---|---
Unity | 25.6% | 61.22% | ${\color[rgb]{1,0,0}\blacktriangledown 60.04\%}$
Unreal | 13.2% | 15.91% | ${\color[rgb]{0,0.88,0}\blacktriangleup 16.42\%}$
Source | 4.0% | 0.27% | ${\color[rgb]{1,0,0}\blacktriangledown 0.23\%}$
CryEngine | 3.5% | 0.23% | ${\color[rgb]{1,0,0}\blacktriangledown 0.19\%}$
GameBryo | 3.2% | - | -
IW | 2.9% | - | -
Anvil | 2.5% | 0.05% | ${\color[rgb]{.5,.5,.5}\bullet 0.05\%}$
id Tech | 1.7% | 0.18% | ${\color[rgb]{1,0,0}\blacktriangledown 0.15\%}$
Essence | 1.1% | - | -
Clausewitz | 1.0% | 0.03% | ${\color[rgb]{1,0,0}\blacktriangledown 0.02\%}$
Other | 48.4% | 22.17% | ${\color[rgb]{0,0.88,0}\blacktriangleup 22.94\%}$

Table 1: Percentage of total games identified on Steam (data collected
2022-09-30 and 2023-10-21)

Game engines that have experienced an increase in their percentage of total
games developed in comparison to the previous year(s) are highlighted in
green, indicating a positive trend, while those that have seen a decrease are
highlighted in red, suggesting a decline in their usage over the same period.
The data from 2018 was collected using the Wikipedia method, while the data
from 2022 and 2023 was collected using the pattern matching method. Color-
coding in Table 1 is missing for the column for 2022. This is due to the
differing data collection methods used in 2018. It is difficult to say whether
the newer method can attribute more game engines or simply more games were
created using a particular game engine. To illustrate, consider Unity, which
has 35.62% more share of all identified games from 2018 to 2022. Additionally,
it is worth noting that among the game engines analyzed, only the Unreal
Engine has consistently exhibited growth over the various years. This trend
holds true when considering both the complete dataset, which includes 2018
data, and the dataset from 2022 and 2023 alone. The Godot Engine is not listed
in Table 1 because the paper from Toftedahl and Engström did not include it in
their analysis. Despite this omission, the data remains valuable in shedding
light on changing market dynamics.

Between 2022 and 2023, it is noteworthy that other game engines have also
gained increased popularity, indicating shifts in the competitive landscape.
At the time of the collected data in 2022, 1.15% of games were created with
the Godot Engine. However, in 2023, this percentage increased to 1.44%.
Compared to other game engines on SteamDB, the Godot Engine ranks 6th out of
65 recognized game engines. In addition to that it is important to understand
that the percentages from Table 1 represent lower bounds. Particularly in the
case of the Godot engine, there could be a lot more games. This is due to the
fact that the Godot Engine allows for the export of executables without a
recognizable filename signature, making them unidentifiable by the SteamDB
algorithm. Consequently, only a subset of games with specific filename
patterns can be matched.

Figure 1:  Number of released Apps on Steam per month

Examining the monthly app releases depicted in Figure 1, a clear trend
emerges: in 2020, there was a significant increase in the number of released
games. This increase in game releases prompts further investigation to
understand the underlying reasons. Prior to conducting a more in-depth
analysis of the underlying reasons, we first draw a parallel by comparing the
data from itch.io. In contrast to Steam, the digital sales platform itch.io
predominantly showcases indie games, aligning with the specific focus of this
paper. Unlike Steam, itch.io possesses an official statistics website that
provides information on game engines [20]. However, it’s essential to
acknowledge that this data might be incomplete due to its reliance on self-
reporting. Therefore, there may be gaps or inaccuracies in the information, a
limitation inherent to self-reported data.

Game Engine | 2018 | 2022 | 2023
---|---|---|---
Unity | 47.3% | ${\color[rgb]{0,0.88,0}\blacktriangleup 49.75\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 46.33\%}$
Construct | 12.3% | ${\color[rgb]{0,0.88,0}\blacktriangleup 13.12\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 13.82\%}$
GameMaker | 11.0% | ${\color[rgb]{1,0,0}\blacktriangledown 7.32\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 6.95\%}$
Twine | 6.2% | ${\color[rgb]{1,0,0}\blacktriangledown 5.35\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 6.03\%}$
RPG Maker | 3.9% | ${\color[rgb]{1,0,0}\blacktriangledown 2.74\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 2.76\%}$
Bitsy | 3.3% | ${\color[rgb]{1,0,0}\blacktriangledown 3.11\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 3.18\%}$
PICO-8 | 2.9% | ${\color[rgb]{1,0,0}\blacktriangledown 2.68\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 2.60\%}$
Unreal | 2.8% | ${\color[rgb]{0,0.88,0}\blacktriangleup 2.92\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 3.01\%}$
Godot | 2.5% | ${\color[rgb]{0,0.88,0}\blacktriangleup 5.55\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 7.51\%}$
Ren’Py | 2.0% | ${\color[rgb]{1,0,0}\blacktriangledown 1.93\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 2.07\%}$
Other | 5.9% | ${\color[rgb]{1,0,0}\blacktriangledown 5.55\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 5.75\%}$

Table 2: Percentage of total games identified on itch.io (data collected
2022-09-30 and 2023-10-21)

Similar to Steam, the data from itch.io can be compared to the 2018 reference
paper using the most current data available from 2022 and 2023. The consistent
data collection method allows for direct comparisons, as illustrated in Table
2. Upon comparing the data presented in Table 1 and Table 2, Unity stands out
as a prevalent game engine. In both 2022 and 2023, Unity consistently
constitutes over half of all games on Steam and nearly half on itch.io. This
dominance is particularly striking when considering the second-ranking game
engines, which represent approximately 16% of games on Steam and roughly 14%
on itch.io in 2023. A similar trend can be observed for 2022.

It can also be seen that only three game engines have consistently gained
popularity on itch.io. These three game engines are Construct, Unreal and the
Godot Engine. Construct is a game engine primarily aimed at non-programmers.
With the help of visual programming, it is intended to make it easier for
novice programmers to get started. The notable increase in its percentage of
games may be attributed to the appeal of game development to a broader
audience with limited coding expertise. However, this is only an assumption
that needs to be investigated further.

Looking at Table 2, it becomes evident that the Godot Engine has experienced
remarkable growth in popularity on itch.io. To delve deeper into this trend,
the previously mentioned script was used to compile games from itch.io along
with their release dates. It should be mentioned that not all games could be
fetched with this script. Despite the missing data sets of up to about 10% of
the games, it is still possible to create a trend plot for the game engines
listed in Table 2.

Figure 2:  Cumulative sum of the number of projects on itch.io

This graph is shown in Figure 2. Notably, Unity, due to its large number of
projects, has been excluded from the graph to ensure better visibility and
clarity for the other game engines. It is evident that the Godot Engine has
gained a significant surge in popularity in 2020, surpassing both Bitsy and
Twine. This analysis has its limitations, as it relies on the quantity of
games published. Quality and quantity of games by individual publishers were
not considered. Determining the exact percentage of indie game studios is
challenging as not every publisher was individually analyzed.

Game Engine | 2020 | 2021 | 2022 | 2023
---|---|---|---|---
Unity | 62.2% | ${\color[rgb]{1,0,0}\blacktriangledown 61.6\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 61.1\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 59\%}$
Godot | 12.2% | ${\color[rgb]{0,0.88,0}\blacktriangleup 13.1\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 15.6\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 19\%}$
GameMaker | 10.9% | ${\color[rgb]{1,0,0}\blacktriangledown 8.9\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 6.1\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 5\%}$
Unreal | - | 4.2% | ${\color[rgb]{0,0.88,0}\blacktriangleup 4.8\%}$ | -
Construct | 1.7% | ${\color[rgb]{0,0.88,0}\blacktriangleup 2.4\%}$ | ${\color[rgb]{1,0,0}\blacktriangledown 1.7\%}$ | -
Stencyl | 0.1% | 0.1% | - | -
Other | 12.9% | ${\color[rgb]{1,0,0}\blacktriangledown 6.5\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 7.0\%}$ | ${\color[rgb]{0,0.88,0}\blacktriangleup 17\%}$
No engine | - | 3.1% | ${\color[rgb]{0,0.88,0}\blacktriangleup 3.7\%}$ | -

Table 3: Used game engines by GMTK Game Jam participants over the years [21]

An interesting distinction on itch.io, which is not as prominent on Steam, is
the prevalence of games created during game jams, which are defined as
follows:

> A game jam is an accelerated opportunistic game creation event where a game
> is created in a relatively short timeframe exploring given design
> constraint(s) and end results are shared publically [22].

With around 18,000 - 23,000 participants and around 6,000 - 7,000 submissions,
GMTK Game Jam is the largest past game jam event on itch.io [23, 24, 25].
Table 3 shows the game engines used by participants in the events from 2020 to
2023. This table is clearly showing that the Godot Engine is the only game
engine from the table to outperform the usage of the previous year three years
in a row.

To explore more about the Godot Engine, the Godot Community Poll from 2022 and
2023 was examined [26, 27]. In 2022, there were 5,315 responses to the survey,
while in 2023 the number increased to 7,671. These surveys contain answers to
a variety of questions. Regarding previous experiences with other game
engines, respondents in the 2022 survey mentioned Unity, other third-party
engines, and GameMaker. In the 2023 survey, the trends remained consistent.
The surveys show that most people heard about the Godot Engine for the first
time between 2018 and 2020. However, most did not start development with the
engine until the following years 2019 to 2022. The surveys also clearly show
that the game engine is mainly used for 2D game development.

The survey from 2022 sheds light on further facts that are relevant to indie
game development. A substantial 84% of participants stated that they use the
Godot Engine as a hobby, while 9% identified themselves as full-time indie
game developers. According to the definition outlined in this paper, hobby
developers fall under the category of indie game developers, as they share
common characteristics such as focusing on specific goals, lacking financial
support from larger companies, and working in small teams. This is further
supported by the finding that 97.9% work on their project alone or in a team
of up to five people. Furthermore, the survey reveals that 83.7% of
respondents do not earn any income from their games, emphasizing their indie
status. Additionally, the two most prevalent digital sales platforms among
respondents are itch.io (7.1%) and Steam (6%).

##  5 Discussion

When examining various data sets, it becomes clear that the popularity of the
Godot Engine has increased. Starting with the digital sales platform Steam, it
is noteworthy that the Godot Engine ranks 6th out of 65, which indicates its
importance even at lower bounds, which arise due to data collection. For a
more meaningful discussion, Unity and Unreal Engine will also be included in
this section. According to Table 1 Unity lost a few percent from 2022 to 2023
while Unreal Engine gained a few percent. A possible reason for this could be
that AAA titles are also published on Steam and Unreal Engine is known to be
used for developing such games [28].

On the other hand, this argument explains why growth on a platform like
itch.io for Unreal Engine is rather slow compared to the Godot Engine. While
the Unreal Engine prominently appears as the second most used game engine on
Steam, its presence on itch.io has seen only marginal growth of roughly 0.2%
over the past four to five years. Furthermore, Table 2 shows that during the
same period of time, the Godot Engine has gained a share of roughly 5%. This
is particularly interesting as it suggests that the Unreal Engine has limited
relevance among indie developers, while the Godot Engine has become
increasingly relevant in the indie game sector over time. This observation is
further supported by the fact that the Godot Engine ranks third on itch.io.

Although Unity is the biggest player in the indie games sector, Table 2
clearly shows that other game engines are gaining popularity for indie game
development. Another game engine under competitive pressure is GameMaker,
which has lost roughly 4% of its share over the years. This observation can
also be reproduced in the context of game jams, where Unity is gradually being
displaced by other competitors, while GameMaker has been rapidly losing ground
in recent years. In contrast, the Godot Engine has seen a significant increase
in popularity in game jams, where it has become the second most used game
engine in the GMTK Game Jam. It is important to consider that the overall
popularity of a game engine could lead to its increased utilization in game
jams. Otherwise, it should also be noted that some games published in the GMTK
Game Jam end up on itch.io and in this way also increase the number of games
developed with the Godot Engine. Further insight would require a precise
analysis of this interaction.

Unity and GameMaker will persist as references in the Godot Engine Community
Polls, as they are mentioned in relation to participants who have experience
with these game engines. Giving the declining shares of Unity and GameMaker,
it is possible to assume that certain developers may have switched entirely to
the Godot Engine. There could be many reasons for this. On one hand, the
license of the Godot Engine offers great freedom and on the other hand,
certain features could be favored. However, in order to draw a direct
comparison, a deeper analysis would also have to be carried out. It should
also be noted that not every game engine can be compared with every other, as
mentioned at the beginning. Each game engine has its own purpose and therefore
a comparison is not easy. If a comparison is nevertheless desired, one
possible useful comparison would be between Unity and the Godot Engine, as
both can be classified as general-purpose game engines that cover a wide range
of game genres.

The results of this study make it clear that the Godot Engine definitely has
relevance in the indie game industry. In all three data sets analyzed (Steam,
itch.io and GMTK Game Jams), the Godot Engine ranks better than most game
engines. In addition, it can also be said that various larger companies that
have shown an interest in the Godot Engine through financial support
indicating its relevance in the industry. The only question that remains is
why this is the case. Both Figure 1 and Figure 2 clearly show 2020 as a
successful year for the Godot Engine. If you combine this finding with the
answers from the Godot Community Poll, you can see that people had already
heard of the game engine before 2020, but only developed with it later. It is
also quite possible that events before 2020 were decisive for the rapidly
growing popularity in the indie game industry. However, this data reveals
little about the reasons for this rise. Even after extensive research, it is
difficult to pinpoint the exact reason for the growth. This is probably due to
the fact that many factors are responsible for this growth.

Another insight is that most games developed with the Godot Engine are mainly
2D games. This is also in line with the statements from the Godot community
survey. On one hand, the Godot Engine might offer more robust features for 2D
game development compared to its 3D capabilities. On the other hand, the
amount of work involved in developing 3D games is usually higher than that of
2D games, which can be a challenge for indie developers, causing them to favor
2D development.

The methods of analysis used here have raised many new questions. These should
continue to be investigated. The market relevance in general or specifically
for indie game studios could be further analyzed. Possible suggestions include
interviews with indie game studios that use the Godot Engine. Additionally, in
the case of Unity and GameMaker, it would be interesting to do further
research to see if the Godot Engine is slowly pushing them out of the market.

One suggestion is to interview indie game studios that use the Godot Engine to
gain further insights. Furthermore, direct comparisons, where appropriate,
would be useful to compare technical aspects between the Godot Engine and
other game engines. Methods such as a feature matrix, benchmarks or
alternative methodologies would be useful for this. Although the reasons for
Godot Engines growth in popularity are unclear, it can be concluded that the
Godot engine is highly relevant in the indie game industry and will probably
continue to grow over the next few years.

## References

  * [1] Andrew Williams.  History of Digital Games: Developments in Art, Design and Interaction, pages 152–154.  CRC Press, Taylor & Francis Group, 2017.
  * [2] Henry Lowood.  Game Engines and Game History.  In History of Games International Conference Proceedings, 2014\.
  * [3] GitHub - Godot.  https://github.com/godotengine/godot.  2023-12-23.
  * [4] Godot Docs - Release policy - Release support timeline.  https://docs.godotengine.org/en/stable/about/release_policy.html#release-support-timeline.  2023-12-23.
  * [5] Unity Forum - Unity 1.0 is shipping.  https://forum.unity.com/threads/unity-1-0-is-shipping.56/.  2023-12-23.
  * [6] Wired - The Imagination Engine: Why Next-Gen Videogames Will Rock Your World.  https://www.wired.com/2012/05/ff-unreal4/.  2023-12-23.
  * [7] Godot Engine - Homepage.  https://godotengine.org/.  2023-12-23.
  * [8] Godot Docs - What is GDNative.  https://docs.godotengine.org/en/stable/tutorials/scripting/gdnative/what_is_gdnative.html.  2023-12-23.
  * [9] Godot Engine - Introducing C# in Godot.  https://godotengine.org/article/introducing-csharp-godot.  2023-12-23.
  * [10] Godot Engine - Godot Engine was awarded an Epic MegaGrant.  https://godotengine.org/article/godot-engine-was-awarded-epic-megagrant.  2023-12-23.
  * [11] Godot Engine - Godot Engine receiving support funded by Facebook Reality Labs.  https://godotengine.org/article/godot-engine-receiving-support-funded-facebook-reality-labs.  2023-12-23.
  * [12] Godot Engine - Godot Engine receiving a new grant from Meta’s Reality Labs.  https://godotengine.org/article/godot-engine-receiving-new-grant-meta-reality-labs.  2023-12-23.
  * [13] Reddit - Godot Engine receiving a new grant from Meta’s Reality Labs.  https://www.reddit.com/r/godot/comments/rm8d6s/comment/hpog9o1.  2023-12-23.
  * [14] Casey O’Donnell.  The Video Game Industry - Formation, Present State, and Future, chapter 5 The North American Game Industry, pages 99–115.  Routledge, 2012.
  * [15] Marcus Toftedahl and Henrik Engström.  A Taxonomy of Game Engines and the Tools that Drive the Industry.  In Digital Games Research Association, editor, DiGRA 2019, The 12th Digital Games Research Association Conference, 2019.
  * [16] UBM - Game Developers Conference.  State of the Game Industry, 2023.
  * [17] UBM - Game Developers Conference.  State of the Game Industry, 2019.
  * [18] SteamDB - Tech.  https://steamdb.info/tech/.  2023-12-23.
  * [19] GitHub - Morphclue/trend-itch.  https://github.com/Morphclue/trend-itch.  2023-12-23.
  * [20] itch.io - Top Engines.  https://itch.io/game-development/engines.  2023-12-23.
  * [21] Twitter - Game Maker’s Toolkit.  https://twitter.com/gamemakerstk/status/1552183527847989248.  2023-12-23.
  * [22] Annakaisa Kultima.  Defining Game Jam.  In Proceedings of the 10th International Conference on the Foundations of Digital Games, 2015.
  * [23] itch.io - GMTK Game Jam 2021.  https://itch.io/jam/gmtk-2021.  2023-12-23.
  * [24] itch.io - GMTK Game Jam 2022.  https://itch.io/jam/gmtk-jam-2022.  2023-12-23.
  * [25] itch.io - GMTK Game Jam 2023.  https://itch.io/jam/gmtk-2023.  2023-12-23.
  * [26] Google Forms - Godot Community Poll 2022.  https://docs.google.com/forms/d/e/1FAIpQLSe-OIpxXqou9cDnPXEAjxzpICbf8_YZB3jUizdECXRydtB8cA/viewanalytics.  2023-12-23.
  * [27] Google Forms - Godot Community Poll 2023.  https://docs.google.com/forms/d/e/1FAIpQLSeXRE1nF64PUilO6fA7Pevh2lWukJtpdBvc2_A3fGfuciy-gQ/viewanalytics.  2023-12-23.
  * [28] Reece A. Boyd and Salvador E. Barbosa.  Reinforcement Learning for All: An Implementation Using Unreal Engine Blueprint.  In International Conference on Computational Science and Computational Intelligence (CSCI), pages 787–792. IEEE, 2017.


# Transport scaling in porous media convection

Xiaojue Zhu1<EMAIL_ADDRESS>Yifeng Fu 1  Marco De Paoli2,3
<EMAIL_ADDRESS>1Max Planck Institute for Solar System Research, 37077
Göttingen, Germany 2Physics of Fluids Group and Max Planck Center for Complex
Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics, University of
Twente, P.O. Box 217, 7500AE Enschede, The Netherlands 3Institute of Fluid
Mechanics and Heat Transfer, TU Wien, 1060 Vienna, Austria

###### Abstract

We present a theory to describe the Nusselt number
($\operatorname{\mathit{Nu}}$), corresponding to the heat or mass flux, as a
function of the Rayleigh–Darcy number ($\operatorname{\mathit{Ra}}$), the
ratio of buoyant driving force over diffusive dissipation, in convective
porous media flows. First, we derive exact relationships within the system for
the kinetic energy and the thermal dissipation rate. Second, by segregating
the thermal dissipation rate into contributions from the boundary layer and
the bulk, which is inspired by the ideas of the Grossmann and Lohse theory
(_J. Fluid Mech._ , vol. 407, 2000; _Phys. Rev. Lett._ , vol. 86, 2001), we
derive the scaling relation for $\operatorname{\mathit{Nu}}$ as a function of
$\operatorname{\mathit{Ra}}$ and provide a robust theoretical explanation to
the empirical relations proposed in previous studies. Specifically, by
incorporating the length scale of the flow structure into the theory, we
demonstrate why heat or mass transport differs between two-dimensional and
three-dimensional porous media convection. Our model is in excellent agreement
with the data obtained from numerical simulations, affirming its validity and
predictive capabilities.

###### keywords:

porous media , convection , heat transport, solute transport, mixing

##  1 Introduction

Carbon dioxide (CO2) sequestration is a process aimed at long-term storage of
large volumes of CO2 (Schrag, 2007), primarily to mitigate climate change and
support energy transition. One of the most promising sequestration strategies
involves natural underground reservoirs. In this case, liquid CO2 is injected
in saline aquifers, geological porous formations located hundreds of meters
beneath the Earth’s surface and confined by horizontal impermeable layers
(Huppert & Neufeld, 2014; De Paoli, 2021). Saline aquifers are filled with
brine, highly salted water denser than CO2. Because of this density
difference, the injected volume of CO2 will sit on top of the brine, and a
critical configuration takes place: in case of fractures in the top confining
layer of the formation, CO2 would migrate upwards and eventually reach the
upper strata up to the atmosphere (Hidalgo et al., 2015). However, CO2 is
partially soluble in brine and the resulting mixture, which is heavier than
both starting fluids, sinks downward through the porous rocks and makes CO2
safely trapped (Emami-Meybodi et al., 2015; Letelier et al., 2023). To
determine the optimal CO2 injection rate and predict the long-time behavior of
the injected CO2, it is therefore imperative to conduct a meticulous
assessment of the flow dynamics and the associated mixing laws (MacMinn et
al., 2012; Guo et al., 2021). An idealized representation of this complex
system consists of a porous domain fully saturated with fluid and confined
between a heated bottom plate and a cooled top plate (Hewitt et al., 2012; Wen
et al., 2018). The top-to-bottom temperature difference induces a density
gradient that drives the flow. This configuration, which we label here as
Rayleigh–Darcy (RD) convection, serves as a fundamental model of the
aforementioned process, where the CO2 concentration field, responsible for the
density increase in the case of geological carbon sequestration, is replaced
by a temperature field. Indeed, thermal and solutal convective porous media
flows can be considered equivalent and controlled by the same governing
equations provided that: (i) in the temperature-driven flow, the solid phase
is locally in thermal equilibrium with the liquid phase, and (ii) in the
corresponding concentration-driven system, the solid is impermeable to the
solute. Additional factors to be accounted for a proper comparison between
these systems are the dependency of the viscosity and the fluid density with
respect to the value of the scalar. While viscosity is generally weakly
sensitive to temperature variations, it may be considerably affected by the
local value of solute concentration. However, it has been previously shown by
Hidalgo et al. (2012) that, in convective porous media flows, concentration-
induced viscosity variations do not significantly affect the global transport
properties of the system. In contrast, the shape of the density-concentration
(or density-temperature) curves is shown to be crucial. For a general
introduction to the RD convection, we refer to the reviews by Hewitt (2020)
and De Paoli (2023). In this work, we will refer to fluid characterized by a
constant viscosity and a linear dependency of density with the transported
scalar (temperature).

The single control parameter in RD convection is the Rayleigh–Darcy number
$\operatorname{\mathit{Ra}}$, which indicates the relative strength between
driving forces (convection) and dissipative effects (diffusion and viscosity),
while the major response parameter of the system is the Nusselt number
$\operatorname{\mathit{Nu}}$, a dimensionless measure of the amount of heat
(or solute) exchanged. Similar to the Rayleigh-Bénard (RB) convection (i.e., a
fluid heated from below and cooled from the top, in the absence of any porous
medium), in recent years, major efforts have been put into understanding the
scaling relation between $\operatorname{\mathit{Nu}}$ and
$\operatorname{\mathit{Ra}}$, where $\operatorname{\mathit{Ra}}$ is intended
as the thermal Rayleigh number (Ahlers et al., 2009). The classical theory
(Priestley, 1954; Malkus, 1954; Howard, 1966) posits that at significantly
high $\operatorname{\mathit{Ra}}$, the buoyancy flux should be independent of
the layer’s height ($L$). In the high-$\operatorname{\mathit{Ra}}$ regime
within a porous medium, this argument predicts a linear scaling of
$\operatorname{\mathit{Nu}}\sim\operatorname{\mathit{Ra}}$. It has also been
rigorously demonstrated that the linear scaling serves as an upper bound
(Doering & Constantin, 1998; Otero et al., 2004; Wen et al., 2012; Hassanzadeh
et al., 2014). Interestingly, such scaling also means that the dimensional
flux is independent of thermal diffusivity and, as a result, a realization of
the scaling indicates the system reaches the so-called ultimate regime (Hewitt
et al., 2012; Pirozzoli et al., 2021). In comparison, in RB convection, a
similar argument leads to
$\operatorname{\mathit{Nu}}\sim\operatorname{\mathit{Ra}}^{1/3}$ (Priestley,
1954; Malkus, 1954; Howard, 1966), different from the diffusion-free ultimate
scaling $\operatorname{\mathit{Nu}}\sim\operatorname{\mathit{Ra}}^{1/2}$
proposed by Kraichnan (1962) and Spiegel (1963). A detailed introduction to
the scalings in RB is provided by Ahlers et al. (2009), Chillà & Schumacher
(2012) and Xia et al. (2023).

Direct numerical simulations (DNS) have been conducted in both two and three
dimensions for RD convection to investigate heat transfer scaling at finite
$\operatorname{\mathit{Ra}}$. The two-dimensional DNS at high
$\operatorname{\mathit{Ra}}$ ($10^{3}\leq\operatorname{\mathit{Ra}}\leq
10^{4}$) by Otero et al. (2004) suggested a slightly sub-linear
$\operatorname{\mathit{Nu}}(\operatorname{\mathit{Ra}})$ scaling. Subsequent
DNS, as reported by Hewitt et al. (2012), extended up to
$\operatorname{\mathit{Ra}}=4\times 10^{4}$, and indicated that the scaling
$\operatorname{\mathit{Nu}}\sim\operatorname{\mathit{Ra}}$ is asymptotically
attained, albeit with a correction to the linear scaling. A simple fit,
$\operatorname{\mathit{Nu}}=0.0069\operatorname{\mathit{Ra}}+2.75$, was
proposed to accommodate the data within this range. For three-dimensional RD
convection, DNS conducted by Pirozzoli et al. (2021) and De Paoli et al.
(2022) reached up to $\operatorname{\mathit{Ra}}=8\times 10^{4}$ and suggested
that the appropriate scaling for $\operatorname{\mathit{Nu}}$ is given by
$\operatorname{\mathit{Nu}}=0.0081\operatorname{\mathit{Ra}}+0.067\operatorname{\mathit{Ra}}^{0.61}$.
This contrasts with an alternative fit proposed by Hewitt et al. (2014), where
$\operatorname{\mathit{Nu}}=0.0096\operatorname{\mathit{Ra}}+4.6$ was
considered for the data within the range up to
$\operatorname{\mathit{Ra}}=2\times 10^{4}$.

It is crucial to underline that the previously mentioned corrections to the
linear scaling are purely empirical in nature. This leads to the fundamental
question: can we provide an explanation for these corrections and quantify
them? To address this, we turn our attention to the Grossmann-Lohse (GL)
theory (Grossmann & Lohse, 2000, 2001), a key tool for comprehending the
effective scaling of Nusselt and Reynolds numbers in relation to
$\operatorname{\mathit{Ra}}$ in turbulent RB convection. The central premise
of the GL theory can be summarized as follows: Firstly, it establishes a
connection between $\operatorname{\mathit{Nu}}$ and
$\operatorname{\mathit{Ra}}$ by considering their relationship with the
kinetic energy dissipation rate and thermal dissipation rate through exact
relations. Secondly, the theory dissects these dissipation rates into
contributions from the boundary layer and the bulk flow. In this work, we will
derive corresponding exact relations for RD convection. By applying the
principles of the GL theory, we can deduce the boundary layer and bulk
contributions to the thermal dissipation rate, shedding light on the origins
and expressions of the corrections to the linear scaling.

##  2 Governing equations

We consider a fluid-saturated porous domain heated from below and cooled from
above, as sketched in Figure 1(a-i). Although we discuss here the problem of a
thermally-driven flow in a porous medium that is locally in thermal
equilibrium with the fluid, the same conclusions apply when the scalar is a
solute, provided that the governing parameter ($\operatorname{\mathit{Ra}}$)
is matched. The size of the domain considered in $W$ is the wall-parallel
directions $x,y$, and $L$ in wall-normal direction $z$, along which gravity
($\mathbf{g}$) acts. The flow is visualized in terms of dimensionless
temperature $T$. For sufficiently high $\operatorname{\mathit{Ra}}$, a
columnar flow structure develops both in 2D and in 3D, as one can observe from
the cross section relative to the domain mid height in Figure 1(a-ii), while
the near-wall region (Figure 1a-iii) is populated by thin filamentary plumes.
This structure differs significantly from the classical RB turbulence,
reported in Figure 1(b), which is controlled in the bulk by large-scale rolls.
In RB convection, large-scale structures span the entire domain, with typical
length scales comparable to the height of the system. In comparison, in RD
convection, columnar-like structures prevail, and are characterized by length
scales distinct from the system’s height. This difference will be taken into
account when we build up a theory for RD convection in the following sections.

Figure 1: Instantaneous dimensionless temperature field $T$ for convection in
(a) porous media for $\operatorname{\mathit{Ra}}=10^{4}$, $W/L=1$ (Pirozzoli
et al., 2021), labeled as Rayleigh-Darcy (RD) convection, and (b) in classical
Rayleigh-Bénard turbulence (thermal Rayleigh number $10^{9}$, Prandtl number
$1$, $W/L=1$), labeled as RB. The temperature distribution is shown over the
entire volume (a-i, b-i), at the centerline $z=1/2$ (a-ii, b-ii) and near the
upper wall (a-iii, b-iii).

Before presenting the dimensionless equations, we introduce the scaling
quantities employed. The equations are made dimensionless with respect to
convective flow scales (Pirozzoli et al., 2021), namely, velocities are scaled
with $\mathscr{U}=\alpha g\Delta K/\nu$, where $\alpha$ is the thermal
expansion rate, $\Delta$ the temperature difference between the bottom and the
top plate, $g$ the acceleration due to gravity, $\nu$ the kinematic viscosity
and $K$ the permeability of the porous medium, which we assume to be
homogeneous and isotropic. Lengths are scaled with $L$ and time with $\phi
L/\mathscr{U}$. The dimensionless temperature is obtained as
$T=(T^{*}-T^{*}_{\text{top}})/\Delta$, being $T^{*}$ and $T^{*}_{\text{top}}$
the dimensional temperature field and the temperature value at the top
boundary, respectively. Finally, pressure is scaled by
$gL(\rho^{*}_{\text{top}}-\rho^{*}_{\text{bot}})$, where the top-to-bottom
fluid density difference ($\rho^{*}_{\text{top}}-\rho^{*}_{\text{bot}}$) is
used. In an incompressible RD system, the heat transport is controlled by the
dimensionless advection-diffusion equation (Pirozzoli et al., 2021):

$\displaystyle\frac{\partial T}{\partial t}+\mathbf{u}\cdot\nabla T$
$\displaystyle=$
$\displaystyle\frac{1}{\operatorname{\mathit{Ra}}}\nabla^{2}T\text{ ,}$ (1)

where $\mathbf{u}$ and $T$ are the velocity and temperature fields,
respectively, $t$ is time and $\operatorname{\mathit{Ra}}$ is the Rayleigh-
Darcy number defined as

$\operatorname{\mathit{Ra}}=\frac{\alpha g\Delta KL}{\kappa\nu}.$ (2)

In this parameter, the medium ($K$), domain ($g,L$) and fluid
($\alpha,\Delta,\nu,\kappa$) properties are included, where $\kappa$ is the
thermal diffusivity. The momentum transport and the flow incompressibility are
accounted by the Darcy law and continuity equations, respectively:

$\mathbf{u}=-(\nabla p-T\mathbf{k})$ (3) $\nabla\cdot\mathbf{u}=0,$ (4)

where $p$ is the reduced pressure field and $\mathbf{k}$ is the unit vector
aligned with $z$.

At the horizontal boundaries, we consider the temperature constant and equal
to $T=1$ at the bottom plate and $T=0$ at the top, so that an unstable
configuration is achieved and the flow is driven by convection. No-penetration
boundary conditions are assumed at both plates for the velocity, while the
sides are considered periodic. Equations (1), (3) and (4) together with these
boundary conditions determine the flow dynamics, which is controlled by two
dimensionless parameters, namely $\operatorname{\mathit{Ra}}$ and the
horizontal domain width $W/L$. The latter does not appear explicitly in the
equations, but may play a significant role in determining the flow structure,
especially at low $\operatorname{\mathit{Ra}}$.

##  3 Nusselt number and exact conservation equations

First, we will correlate the thermal (Nusselt number) and the kinetic (Péclet
number) response parameters to the control parameter (Rayleigh-Darcy number),
and then the thermal dissipation will be linked to the Nusselt number.

The temporal and horizontal average of (1) can be written as

$\frac{\partial}{\partial
z}\left(\operatorname{\mathit{Ra}}\overline{\left\langle
u_{z}T\right\rangle}_{A}-\overline{\left\langle\frac{\partial T}{\partial
z}\right\rangle}_{A}\right)=0.$ (5)

The Nusselt number then reads (Letelier et al., 2019; Ulloa & Letelier, 2022)

$\operatorname{\mathit{Nu}}=\operatorname{\mathit{Ra}}\overline{\left<u_{z}T\right>}_{A}-\overline{\left<\frac{\partial
T}{\partial z}\right>}_{A}.$ (6)

Here the following notations are used for different averaging procedures.
Overbars $\overline{\cdots}$ correspond to the time average of a dimensionless
value $f$, while an average over the horizontal surface and an average over
the whole volume domain are denoted by $\left<\cdots\right>_{A}$ and
$\left<\cdots\right>$, respectively:

$\displaystyle\overline{f}$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}\int_{t_{0}}^{t_{0}+\tau}f\mathrm{d}t$ (7)
$\displaystyle\left\langle f\right\rangle_{A}$ $\displaystyle=$
$\displaystyle\frac{1}{A}\int_{0}^{W/L}\int_{0}^{W/L}f\mathrm{d}x\mathrm{d}y$
(8) $\displaystyle\left\langle f\right\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{V}\int_{0}^{W/L}\int_{0}^{W/L}\int_{0}^{1}f\mathrm{d}x\mathrm{d}y\mathrm{d}z,$
(9)

where $A=(W/L)^{2}$ is the dimensionless horizontal surface area and
$V=(W/L)^{2}$ is the dimensionless volume of the whole domain based on our
characteristic length scale $L$. Two exact relations exist in our system and
can be derived from the governing equations. Using the dimensionless velocity
$\mathbf{u}$ to dot product both sides of (3) and combining with the
incompressible continuity equation (4), we get:

$\left|\mathbf{u}\right|^{2}=-\nabla\cdot(p\mathbf{u})+Tu_{z}.$ (10)

The volume and time average of (10) reads

$\overline{\left<\left|\mathbf{u}\right|^{2}\right>}=-\frac{1}{V}\iint_{\Sigma}(p\mathbf{u})\cdot\hat{\mathbf{n}}\mathrm{d}S+\overline{\left<Tu_{z}\right>},$
(11)

where $\Sigma$ is the boundary surface of the domain, $\mathrm{d}S$ denotes
the surface element on the boundary, and $\hat{\mathbf{n}}$ is the normal unit
vector for the surface elements. The mean power given by pressure gradient
vanishes due to the non-penetration boundary condition:

$\frac{1}{V}\iint_{\Sigma}(p\mathbf{u})\cdot\hat{\mathbf{n}}\mathrm{d}S=-\frac{1}{V}\iint_{\Sigma(z=0)}pu_{z}\mathrm{d}S+\frac{1}{V}\iint_{\Sigma(z=1)}pu_{z}\mathrm{d}S=0.$
(12)

The mean buoyancy power in (11) can be written as

$\overline{\left<Tu_{z}\right>}=\frac{1}{V}\int_{0}^{1}\overline{\left<Tu_{z}\right>}_{A}A\mathrm{d}z=\frac{1}{\operatorname{\mathit{Ra}}}\int_{0}^{1}\left(\operatorname{\mathit{Nu}}+\overline{\left<\frac{\partial
T}{\partial z}\right>}_{A}\right)\mathrm{d}z.$ (13)

Here the last equivalence comes from the $\operatorname{\mathit{Nu}}$
definition (6). Since we use $L$ as our length scale, $z\in[0,1]$. The last
term in the above equation reads:

$\frac{1}{\operatorname{\mathit{Ra}}}\int_{0}^{1}\overline{\left<\frac{\partial
T}{\partial
z}\right>}_{A}\mathrm{d}z=\frac{1}{\operatorname{\mathit{Ra}}}\int_{0}^{1}\frac{\partial\overline{\left<T\right>}_{A}}{\partial
z}\mathrm{d}z=\frac{1}{\operatorname{\mathit{Ra}}}\left(\overline{\left<T\right>}_{A}|_{z=1}-\overline{\left<T\right>}_{A}|_{z=0}\right)=-\frac{1}{\operatorname{\mathit{Ra}}}.$
(14)

Reintroducing (14) back into (13), after some algebraic manipulations we get:

$\overline{\left<Tu_{z}\right>}=\frac{1}{\operatorname{\mathit{Ra}}}\int_{0}^{1}\left(\operatorname{\mathit{Nu}}+\overline{\left<\frac{\partial
T}{\partial
z}\right>}_{A}\right)\mathrm{d}z=\frac{1}{\operatorname{\mathit{Ra}}}(\operatorname{\mathit{Nu}}-1).$
(15)

Combining (11), (12) and (15), we obtain an expression for the mean
dimensionless velocity square:

$\overline{\left<\left|\mathbf{u}\right|^{2}\right>}=\frac{1}{\operatorname{\mathit{Ra}}}(\operatorname{\mathit{Nu}}-1).$
(16)

We introduce the velocity scale

$\mathscr{V}=\mathscr{U}\sqrt{\overline{\left<\left|\mathbf{u}\right|^{2}\right>}},$
(17)

with $\mathscr{U}=\alpha g\Delta K/\nu$, and one finally obtains an exact
relation:

$\operatorname{\mathit{Pe}}^{2}=(\operatorname{\mathit{Nu}}-1)\operatorname{\mathit{Ra}}$
(18)

with

$\operatorname{\mathit{Pe}}=\frac{\mathscr{V}L}{\kappa}=\operatorname{\mathit{Ra}}\frac{\mathscr{V}}{\mathscr{U}},$
(19)

where $\operatorname{\mathit{Pe}}$ is the Péclet number. To the best of our
knowledge, this is the first time that a relation among
$\operatorname{\mathit{Pe}}$, $\operatorname{\mathit{Ra}}$, and
$\operatorname{\mathit{Nu}}$ has been presented for RD convection. Note that
for RB convection, the analogous exact relation is
$\epsilon_{u}=\nu^{3}/L^{4}(\operatorname{\mathit{Nu}}-1)\operatorname{\mathit{Ra}}\operatorname{\mathit{Pr}}^{-2}$,
where $\epsilon_{u}$ is the kinetic energy dissipation rate and $Pr$ is the
Prandtl number (Ahlers et al., 2009). To assess the validity of Eq. (18), we
consider the numerical measurements available in literature. For 2D flows,
$\operatorname{\mathit{Pe}}$ is measured by De Paoli et al. (2024) using (17)
and (19). The velocity Root Mean Square (rms) at the centerline (2D and 3D) is
computed by Hewitt et al. (2012, 2014) and reported in Figure 2(a). Since in
all directions no mean flow exists, we have that $\mathscr{V}$ is obtained
from the rms of the velocity components ($u_{i}$), namely:

$\mathscr{V}=\mathscr{U}\left(\overline{\left<[\text{rms}(u_{x})]^{2}+[\text{rms}(u_{y})]^{2}+[\text{rms}(u_{z})]^{2}\right>}\right)^{1/2}.$
(20)

Assuming that the centerline flow is representative of the kinetic energy of
the system, we have that:

$\displaystyle\mathscr{V}\approx\mathscr{U}\left(\overline{\left<[\text{rms}(u_{x})]^{2}+[\text{rms}(u_{y})]^{2}+[\text{rms}(u_{z})]^{2}\right>_{\Sigma(z=1/2)}}\right)^{1/2}.$
(21)

We use this approximation to compute $\mathscr{V}$ and verify the validity of
(18) for the data of Hewitt et al. (2012, 2014). We finally observe in Figure
2(b) that Eq. (18) (dashed line) is in excellent agreement with the
measurements obtained from the exact definition of
$\operatorname{\mathit{Pe}}$ (2D and $\mathscr{V}$ computed with Eq. (17) from
De Paoli et al., 2024) and also with measurements obtained from the
approximated definition of $\operatorname{\mathit{Pe}}$ (2D and 3D,
$\mathscr{V}$ computed with Eq. (21) from Hewitt et al., 2012, 2014).

Figure 2:  (a) Root Mean Square (rms) of velocity used to estimate the Péclet
number, $\operatorname{\mathit{Pe}}$. Horizontal component
($\textnormal{rms}(u)$), vertical component ($\textnormal{rms}(w)$) and total
(defined as $\sqrt{[\textnormal{rms}(u)]^{2}+[\textnormal{rms}(w)]^{2}}$ for
2D, and $\sqrt{2[\textnormal{rms}(u)]^{2}+[\textnormal{rms}(w)]^{2}}$ for 3D)
are reported. Data from Hewitt et al. (2012, 2014). (b) Comparison of
$\operatorname{\mathit{Pe}}$ compensated with
$(\operatorname{\mathit{Nu}}-1)\operatorname{\mathit{Ra}}$ as from (18).
$\operatorname{\mathit{Pe}}$ is computed with $\mathscr{V}$ defined as in (21)
for the data from Hewitt et al. (2012, 2014), and as in Eq. (17) for data from
De Paoli et al. (2024). (c) Ratio of the Nusselt number to the dimensionless
thermal dissipation, derived in (29). Data from Pirozzoli et al. (2021) and De
Paoli et al. (2024).

We will now derive an equation to correlate the mean thermal dissipation to
the Nusselt number. Multiplying the dimensionless thermal advection-diffusion
equation (1) by $T$, integrating over the whole domain and time-averaging, we
have:

$\overline{\left<\frac{\partial}{\partial
t}\left(\frac{T^{2}}{2}\right)\right>}=-\overline{\left<T\mathbf{u}\cdot\nabla
T\right>}+\frac{1}{\operatorname{\mathit{Ra}}}\overline{\left<T\nabla^{2}T\right>}.$
(22)

Under the assumption of statistically steady state:

$\overline{\left<\frac{\partial}{\partial
t}\left(\frac{T^{2}}{2}\right)\right>}=\frac{\partial}{\partial
t}\overline{\left<\frac{T^{2}}{2}\right>}=0.$ (23)

The first term on the right hand side of (22) can be further written as

$-\overline{\left<T\mathbf{u}\cdot\nabla
T\right>}=\overline{\left<\frac{T^{2}}{2}\nabla\cdot\mathbf{u}\right>}-\overline{\left<\nabla\cdot\left(\mathbf{u}\frac{T^{2}}{2}\right)\right>}=-\frac{1}{V}\iint_{\Sigma}\left(\mathbf{u}\frac{T^{2}}{2}\right)\cdot\hat{\mathbf{n}}\mathrm{d}S=0.$
(24)

Here we employed again continuity (4) and the no-penetration boundary
condition. The second term on the right hand side of (22) reads

$\frac{1}{\operatorname{\mathit{Ra}}}\overline{\left<T\nabla^{2}T\right>}=\frac{1}{\operatorname{\mathit{Ra}}}\left[\overline{\left<\nabla^{2}\left(\frac{T^{2}}{2}\right)\right>}-\overline{\left<\left\lvert\nabla
T\right\lvert^{2}\right>}\right].$ (25)

Combining results from (22)-(25), one obtains

$\overline{\left<\left\lvert\nabla
T\right\lvert^{2}\right>}=\overline{\left<\nabla^{2}\left(\frac{T^{2}}{2}\right)\right>}.$
(26)

We can use the following procedure to further simplify the right hand side of
(26)

$\displaystyle\overline{\left<\nabla^{2}\left(\frac{T^{2}}{2}\right)\right>}$
$\displaystyle=$
$\displaystyle\frac{1}{V}\iint_{\Sigma}\left(\overline{T\nabla
T}\right)\cdot\hat{\mathbf{n}}\mathrm{d}S$ (27) $\displaystyle=$
$\displaystyle\frac{A}{V}\left(\overline{\left<T\frac{\partial T}{\partial
z}\right>}_{\Sigma(z=1)}-\overline{\left<T\frac{\partial T}{\partial
z}\right>}_{\Sigma(z=0)}\right)=\operatorname{\mathit{Nu}}.$ (28)

Here we considered that $V=A=(W/L)^{2}$, applied the boundary conditions for
$T$ and $u_{z}$, as well as the $\operatorname{\mathit{Nu}}$ definition (6).
Combining (26) and (28), we get

$\overline{\left<\left\lvert\nabla
T\right\lvert^{2}\right>}=\operatorname{\mathit{Nu}},$ (29)

where $\overline{\left<\left\lvert\nabla T\right\lvert^{2}\right>}$ represents
the dimensionless mean thermal dissipation. This relation, which holds also
for RB convection (Ahlers et al., 2009), has been also presented before for RD
flows (Otero et al., 2004; Hidalgo et al., 2012; De Paoli, 2023) and for Hele-
Shaw convection the limit of infinitely thin domains (Letelier et al., 2019;
Ulloa & Letelier, 2022). The ratio of the Nusselt number to the dimensionless
thermal dissipation is compared in Figure 2(c) for the numerical results of
Pirozzoli et al. (2021) and De Paoli et al. (2024). We observe that, also in
this case, the agreement between theory and simulations is good, confirming
the validity of (29). Finally, one obtains the dimensional thermal dissipation
rate:

$\epsilon=\kappa\frac{\Delta^{2}}{L^{2}}\overline{\left<\left\lvert\nabla
T\right\lvert^{2}\right>}=\kappa\frac{\Delta^{2}}{L^{2}}\operatorname{\mathit{Nu}}.$
(30)

Equations (18) and (30) represent the two exact relations we derived in our
system.

##  4 Application of GL theory and scaling relation for the Nusselt number

Figure 3:  (a) Horizontal- and time-averaged profiles of temperature are shown
for different $\operatorname{\mathit{Ra}}$ (De Paoli et al., 2022). The near-
wall region is magnified in (b), where the wall-normal coordinate is rescaled
with $\operatorname{\mathit{Nu}}$. The profiles are self-similar, and in the
boundary layer follow a linear behavior (solid line) with unitary slope. (c)
The thickness of the boundary layer is determined as the distance from the
wall of the intersection between the linear profile fitting the bulk ($0.4\leq
z\leq 0.6$) and the near-wall regions ($z\leq
0.1\operatorname{\mathit{Nu}}/\operatorname{\mathit{Ra}}$). The case relative
to $\operatorname{\mathit{Ra}}=10^{4}$ is reported. A close-up view of the
near-wall region is shown in the inset. (d) Thickness of the thermal boundary
layer $\lambda/L$ as a function of the Rayleigh number. 3D measurements
computed as discussed (bullets) are very well fitted by the correlation
$50/\operatorname{\mathit{Ra}}$. The correlation proposed for 2D flows by
Otero et al. (2004) ($\lambda/L\sim 15/\operatorname{\mathit{Ra}}$) is also
reported, as well as the value obtained from the Nusselt number
$\lambda/L=1/(2\operatorname{\mathit{Nu}})$. Data from De Paoli et al. (2022).

With the two exact relations derived in the previous section, we can now apply
the main ideas of the GL theory to RD convection. The key idea of the GL
theory (Grossmann & Lohse, 2000, 2001) is to split the kinetic and thermal
dissipation rates into contributions from the boundary layers (BL) and bulk.
In RD convection, the procedure becomes even simpler than in RB, as only the
thermal dissipation rate appears in the exact relations. We separate the mean
thermal dissipation as

$\epsilon=\epsilon_{BL}+\epsilon_{bulk},$ (31)

and apply the respective scaling relations for $\epsilon_{BL}$ and
$\epsilon_{bulk}$, based on the boundary layer theory and fully developed flow
in the bulk. The horizontal- and time-averaged profiles of temperature, shown
in Figure 3(a), confirm that the flow can be split into two distinct regions:
a well-mixed bulk with nearly-uniform properties, and a thin boundary layer
characterized by a linear temperature profile, the slope of which is unitary
when $z$ is rescaled with $\operatorname{\mathit{Nu}}$ (Figure 3b). The
thickness of this boundary layer, $\lambda/L$, can be determined as the
distance from the wall at which the linear function fitting the temperature
profile in the bulk ($0.4\leq z\leq 0.6$) intersects the near-wall temperature
fit. The measurement procedure is illustrated in Figure 3(c), where the
intersection is marked by the bullet. The Nusselt number sets the thickness of
the boundary layer $\lambda/L=1/(2\operatorname{\mathit{Nu}})$ (De Paoli et
al., 2022). In RD convection, it has been proposed by Otero et al. (2004) that
the thermal boundary layer thickness scales as
$\lambda/L\sim\operatorname{\mathit{Ra}}^{-1}$ (consistent with $Nu\sim Ra$,
from the classical theory (Priestley, 1954; Malkus, 1954; Howard, 1966) and
the upper bound scaling derived by Doering & Constantin (1998)). As
illustrated in Figure 3(d), this approximation is verified. Although both in
2D and 3D the thermal boundary layer thickness follows
$\lambda/L\sim\operatorname{\mathit{Ra}}^{-1}$, the prefactor differs (see
Figure 3d). This discrepancy arises from the distinct flow structures in 2D
and 3D. In 3D, owing to the additional degree of freedom compared to the 2D
case, plumes can freely move and reorganize towards the most efficient
configuration, resulting in a different value of $\operatorname{\mathit{Nu}}$.
Consequently, the boundary layer thickness also varies. This phenomenon is
analogous to RB convection, where 2D and 3D flows exhibit different boundary
layer thicknesses due to variations in $\operatorname{\mathit{Nu}}$ (Van Der
Poel et al., 2013).

The profiles of dimensionless thermal dissipation $\overline{\langle|\nabla
T|^{2}\rangle_{A}}$ obtained from De Paoli et al. (2022) are shown in Figure
4(a) for different values of the Rayleigh number. In the inset, the
dissipation is rescaled by the respective Nusselt number, and shown up to 1.
We observe that the boundary layer contribution to the dissipation is more
pronounced as $\operatorname{\mathit{Ra}}$ is increased. A more quantitative
description is provided in the following. The thermal dissipation rate in the
boundary layer scales as $\sim\kappa(\Delta/\lambda)^{2}$. Therefore, taking
into account the layer extension
($\lambda/L\sim\operatorname{\mathit{Ra}}^{-1}$), the boundary layer
contribution towards the total thermal dissipation rate reads

$\epsilon_{BL}\sim\kappa\frac{\Delta^{2}}{\lambda^{2}}\frac{\lambda}{L}\sim\kappa\frac{\Delta^{2}}{L^{2}}\operatorname{\mathit{Ra}}.$
(32)

Assuming the flow in the bulk is well mixed (Grossmann & Lohse, 2000; Bader &
Zhu, 2023; Song, Shishkina & Zhu, 2024), we get

$\epsilon_{bulk}\sim\frac{\mathscr{V}\theta^{2}}{\ell}.$ (33)

Here, $\mathscr{V}$ and $\theta$ are the typical velocity and temperature
scales, respectively. The characteristic length scale is $\ell$, defined as
the wavelength associated with the power-averaged mean wave number at the mid
height $(\overline{k})$, i.e. $\ell/L=2\pi/\overline{k}$. In numerical
simulations, $\overline{k}$ is obtained from the time-averaged spectrum of the
temperature field at $z=1/2$ (Hewitt et al., 2014). The importance of
$\epsilon_{BL}$ and $\epsilon_{bulk}$ relative to the total dissipation
$\epsilon$, is reported in Figure 4(b) for 3D RD convection (De Paoli et al.,
2022). Here, $\epsilon_{BL}$ and $\epsilon_{bulk}$ are obtained from the
profiles, and represent the mean value within the respective regions. In RB,
to derive the scalings, the length scale is assumed to be the height of the
domain $L$, which makes sense as there exist large-scale rolls. However, here
in RD, the typical flow structures are columnar-like, making the length scale
$\ell$ different from $L$. This difference is clearly visible in Figures
1(a-i) and 1(b-i). From the definition of $\operatorname{\mathit{Nu}}$
(equation 6) and assuming that in the bulk,
$\theta\sim\operatorname{\mathit{Nu}}(\kappa\Delta)/(\mathscr{V}L)$, and when
$\operatorname{\mathit{Nu}}$ only comes from the fluctuation in the bulk,
$\theta\sim(\epsilon_{bulk}/\mathscr{V})(L/\Delta)$, we get

$\epsilon_{bulk}\sim\kappa\frac{\Delta^{2}}{L^{2}}{\operatorname{\mathit{Pe}}}\frac{\ell}{L}.$
(34)

The same bulk scaling has also been reported for rapidly rotating convection
(Song, Shishkina & Zhu, 2024) and magnetoconvection with strong vertical
magnetic fields (Bader & Zhu, 2023). In all these three systems, in the bulk
there exists a new horizontal dominant length scale that is different from the
height of the domain. In each of these three systems, a new dominant length
scale emerges, distinct from the domain’s height. This disparity constitutes a
significant deviation from the original GL theory.

Figure 4:  (a) Horizontal- and time-averaged profiles of dimensionless thermal
dissipation, $\overline{\langle|\nabla T|^{2}\rangle_{A}}$, are shown for
different $\operatorname{\mathit{Ra}}$ (De Paoli et al., 2022). In the inset,
the dissipation is rescaled by the respective Nusselt number, and shown up to
$\overline{\langle|\nabla T|^{2}\rangle_{A}}/\operatorname{\mathit{Nu}}=1$.
The boundary layer contribution to the dissipation is more pronounced as
$\operatorname{\mathit{Ra}}$ is increased. (b) Contributions to the mean
scalar dissipation $\epsilon$ within the boundary layer ($\epsilon_{BL}$) and
in the bulk region ($\epsilon_{bulk}$). Data are from De Paoli et al. (2022).

Determining the dominant wavelength in RD convection is a challenging task.
The reason is linked to the complex way in which the dynamic near-wall flow
structures interact with the stationary, large-scale columnar plumes
controlling the bulk. A detailed review is provided by Hewitt (2020), which we
summarize here with additional details including later results (De Paoli et
al., 2022). In 3D, using asymptotic stability theory, Hewitt & Lister (2017)
derived that

$\ell/L\sim\operatorname{\mathit{Ra}}^{-1/2}\quad(3D).$ (35)

Numerical simulations by Hewitt et al. (2014) and De Paoli et al. (2022), best
fitted by scaling exponents of $-0.52$ and $-0.49$, respectively, agree well
with this prediction (see also Figure 5a). Therefore, we employ this scaling
relation ($\ell/L\sim\operatorname{\mathit{Ra}}^{-1/2}$) for the centerline in
3D. The situation is more complex in 2D. Hewitt et al. (2013b) derived
analytically the scaling relation

$\ell/L\sim\operatorname{\mathit{Ra}}^{-5/14}\quad(2D).$ (36)

Wen et al. (2015) have shown that for $\operatorname{\mathit{Ra}}\leq 19976$
the centerline dominant length scale is well approximated by
$\ell/L\sim\operatorname{\mathit{Ra}}^{-0.40}$. However, one can observe in
Figure 5(b) that when $\operatorname{\mathit{Ra}}\geq 39716$, the inter-plume
spacing measured by Wen et al. (2015) is not unique. The conclusion is that in
2D RD convection, at $\operatorname{\mathit{Ra}}\geq 39716$, a precise scaling
remains to be established by running simulations in very wide domains and for
very long times. In view of this, we consider the scaling proposed by Hewitt
et al. (2013b), which represents the best theoretical prediction available, to
be valid. Therefore, combining (34) with (18), (35) and (36), we get

$\epsilon_{bulk}\sim\kappa\frac{\Delta^{2}}{L^{2}}\operatorname{\mathit{Nu}}\quad(3D).$
(37)

and

$\epsilon_{bulk}\sim\kappa\frac{\Delta^{2}}{L^{2}}{\operatorname{\mathit{Nu}}\operatorname{\mathit{Ra}}^{1/7}}\quad(2D).$
(38)

Figure 5:  Dominant length scales $\ell/L$ in the bulk ($z=1/2$) for (a) 3D
(Hewitt et al., 2014; De Paoli et al., 2022) and (b) 2D (Hewitt et al., 2012;
Wen et al., 2015) simulations. Theoretical scaling relations defined in Eq.
(35) and (36), corresponding to 3D and 2D flows respectively, are also
indicated (dashed lines).

Finally, combining (30), (31), (32), (37) and (38), we reach an expression for
$\operatorname{\mathit{Nu}}$ as a function of $\operatorname{\mathit{Ra}}$ for
the 3D and the 2D cases:

$\operatorname{\mathit{Nu}}=A_{3}\operatorname{\mathit{Ra}}+B_{3}\operatorname{\mathit{Nu}}^{1/2}\quad(3D),$
(39)
$\operatorname{\mathit{Nu}}=A_{2}\operatorname{\mathit{Ra}}+B_{2}\operatorname{\mathit{Nu}}^{1/2}\operatorname{\mathit{Ra}}^{1/7}\quad(2D).$
(40)

As reported in Figure 6, these scaling relations fit very well the data
obtained from numerical simulations, both in the 2D and in the 3D cases. The
values of the coefficients $A_{2}$, $A_{3}$, $B_{2}$, $B_{3}$, indicated in
Figure 6, are obtained as best fitting from the data shown, representing the
numerical results available and with $\operatorname{\mathit{Ra}}>2\times
10^{3}$. The choice of considering values larger than this threshold is
motivated by the flow topology: at low $\operatorname{\mathit{Ra}}$ the bulk
flow structure is not columnar, as it is dominated by large-scale convective
rolls (Graham & Steen, 1994), and therefore our theory does not apply. The
expressions of $\operatorname{\mathit{Nu}}$ derived in (39) and (40) take the
form of a linear scaling with a sublinear correction. The linear scaling was
previously proposed for 2D (Hewitt et al., 2012) and 3D (Hewitt et al., 2014)
flows. The scaling proposed here provides similar results to the linear
scaling with sublinear corrections proposed for the 3D case by Pirozzoli et
al. (2021). However, in this case one fitting parameters less is used, i.e.,
all scaling exponents are known. The good agreement of the present scaling
relations with the numerical measurements available in literature suggests
that the theory proposed is indeed valid and promising for higher
$\operatorname{\mathit{Ra}}$. The RD system is completely defined by two
parameters, namely the Rayleigh number $\operatorname{\mathit{Ra}}$ and the
aspect ratio $W/L$. In 3D flows at high-$\operatorname{\mathit{Ra}}$, it has
been observed that all major flow statistics converge for an aspect ratio of
$W/L=1$ (De Paoli et al., 2022). This differs in 2D systems: also at high
$\operatorname{\mathit{Ra}}$ and due to the additional lateral confinement,
the aspect ratio may have an influence on $\ell$ (Wen et al., 2015).
Therefore, it may be required to take $W/L$ into account in the present theory
to better describe the transport properties in 2D RD flows. To this aim and
also to assess the physics of the scaling prefactors, additional simulations
in large domains and at larger $\operatorname{\mathit{Ra}}$ are needed.

Figure 6:  Compensated Nusselt number as a function of Rayleigh number.
Results are for 2D (Hewitt et al., 2012; Wen et al., 2015; De Paoli et al.,
2016; Pirozzoli et al., 2021; De Paoli et al., 2024) and 3D (Pirozzoli et al.,
2021; Hewitt et al., 2014) simulations. New fitting curves are obtained
considering all data with $\operatorname{\mathit{Ra}}>2\times 10^{3}$, with
numerical values of the coefficients being $A_{2}=(6.386\pm 0.007)\times
10^{-3}$, $B_{2}=0.2533\pm 0.0050$ and $A_{3}=(8.592\pm 0.033)\times 10^{-3}$,
$B_{3}=1.376\pm 0.059$. The coefficients of determination ($R^{2}$) for the
best fitting curves provided are 0.991 and 0.999 in the 2D and in the 3D
cases, respectively.

##  5 Conclusions

In summary, we have established two exact relationships, one pertaining to the
Péclet number and the other to the thermal dissipation rate, in the context of
Rayleigh-Darcy convection - a fundamental system for heat and mass transport
in porous media. To the best of our knowledge, this is the first time that a
relation for Péclet number has been proposed for RD convection. Inspired by
previous models developed for Rayleigh-Bénard convection (Grossmann & Lohse,
2000, 2001), we have formulated a scaling theory for heat transfer in 2D and
3D Rayleigh-Darcy flows, where the Nusselt number is expressed as a function
of the Rayleigh number as described by (39) and (40). This theory enables us
to provide a theoretical explanation to the sublinear empirical corrections
proposed in prior studies (Pirozzoli et al., 2021). Our investigations,
supported by both 2D and 3D literature results, confirm the validity of the
proposed theory. Moreover, by taking the length scale of the flow structures
into account, we also shed new light on the physical origins of the
disparities in scaling relations between 2D and 3D Rayleigh-Darcy convection.

Our findings are relevant to convective flows in homogeneous and isotropic
porous media where the top-to-bottom density difference is defined. However,
these hypotheses represent idealized conditions not taking into account
additional flow features that occur in realistic processes, such as
hydrodynamic dispersion (Liang et al., 2018; Tsinober et al., 2022, 2023),
medium heterogeneity (Dentz et al., 2023; Simmons et al., 2001), anisotropy
(De Paoli et al., 2016; Ennis-King et al., 2005) and alternative flow
configurations (Hidalgo et al., 2012; De Paoli et al., 2017; Letelier et al.,
2023). Nonetheless, present findings represent a crucial step required to
develop a robust and physically-sound theory for convection in porous media
flows, which can be further expanded to include the presence of the different
variations mentioned above. We consider for instance the sequestration of
carbon dioxide in saline aquifers. Such a system is usually modelled as a
rectangular domain initially filled with brine and confined by two horizontal
low-permeability layers, and therefore it is assumed to be impermeable at the
bottom boundary (no-flux) (Huppert & Neufeld, 2014). Here, the solute enters
from the top, in correspondence of which the concentration of CO2 is constant.
This flow configuration, defined as “one-sided” (Hewitt et al., 2013a) or
“semi-infinite” (Ennis-King et al., 2005), is subject to a transient
behaviour: the average CO2 concentration within the system will increase over
time, and it will be progressively harder to keep dissolving solute. In
quantitative terms this means that, after a short initial phase in which
dissolution increases due to the formation and growth of the fingers, the flux
of solute through the top boundary will later reduce as a result of the
saturation of the domain. The dynamics of such a system has been thoroughly
characterized (Slim, 2014), and it is shown to be quantitatively related to
the dynamics observed in RD convection (Hewitt et al., 2013a; De Paoli et al.,
2017). In order to describe the evolution of the one-sided system with a
simple box-model, accurate predictions of the transport scaling in RD
convection are essential. In presence of high-permeability formations as the
Utsira Sand reservoir at Sleipner (Bickle et al., 2007), the Rayleigh-Darcy
number may be as high as $6\times 10^{5}$ (Hewitt et al., 2013a), which is
beyond current numerical capabilities. As a result, the theoretical results
provided in our work will play a crucial role as a tool to determine the long-
term evolution of flows in semi-infinite domains.

[Funding]We gratefully acknowledge the financial support from the Max Planck
Society, the German Research Foundation through grants 521319293 and
540422505, the Alexander von Humboldt Foundation, and the Daimler Benz
Foundation. This project has also received funding from the European Union’s
Horizon Europe research and innovation programme under the Marie Sklodowska-
Curie grant agreement No. 101062123.

[Declaration of interests]The authors report no conflict of interest. [Author
ORCID]
Xiaojue Zhu, https://orcid.org/0000-0002-7878-0655;
Yifeng Fu, https://orcid.org/0009-0006-6332-045X;
Marco De Paoli, https://orcid.org/0000-0002-4709-4185;

## References

  * Ahlers et al. (2009) Ahlers, G., Grossmann, S. & Lohse, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Rev. Mod. Phys. 81 (2), 503.
  * Bader & Zhu (2023) Bader, S. H. & Zhu, X. 2023 Scaling relations in quasi-static magnetoconvection with a strong vertical magnetic field. J. Fluid Mech. 976, A4.
  * Bickle et al. (2007) Bickle, M, Chadwick, A, Huppert, H E, Hallworth, M & Lyle, S 2007 Modelling carbon dioxide accumulation at Sleipner: implications for underground carbon storage. Earth Planet. Sci. Lett. 255 (1), 164–176.
  * Chillà & Schumacher (2012) Chillà, F. & Schumacher, J. 2012 New perspectives in turbulent Rayleigh-Bénard convection. Eur. Phys. J. E 35, 1–25.
  * De Paoli (2021) De Paoli, M. 2021 Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide. Phys. Fluids 33 (1), 016602.
  * De Paoli (2023) De Paoli, M. 2023 Convective mixing in porous media: a review of Darcy, pore-scale and Hele-Shaw studies. Eur. Phys. J. E 46 (12), 129.
  * De Paoli et al. (2022) De Paoli, M., Pirozzoli, S., Zonta, F. & Soldati, A. 2022 Strong Rayleigh–Darcy convection regime in three-dimensional porous media. J. Fluid Mech. 943, A51.
  * De Paoli et al. (2024) De Paoli, M., Yerragolam, G. S., Lohse, D. & Verzicco, R. 2024 A finite difference solver for numerical simulations of convective porous media flows. (in preparation) .
  * De Paoli et al. (2016) De Paoli, M, Zonta, F & Soldati, A 2016 Influence of anisotropic permeability on convection in porous media: implications for geological CO2 sequestration. Phys. Fluids (1994-present) 28 (5), 056601.
  * De Paoli et al. (2017) De Paoli, M, Zonta, F & Soldati, A 2017 Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown. Phys Fluids (1994-present) 29 (2), 026601\.
  * Dentz et al. (2023) Dentz, M, Hidalgo, JJ & Lester, D 2023 Mixing in porous media: concepts and approaches across scales. Transp. Porous Media 146 (1-2), 5–53.
  * Doering & Constantin (1998) Doering, CR & Constantin, P 1998 Bounds for heat transport in a porous layer. J. Fluid Mech. 376, 263–296.
  * Emami-Meybodi et al. (2015) Emami-Meybodi, H., Hassanzadeh, H., Green, C. P. & Ennis-King, J. 2015 Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments. Int. J. Greenh. Gas Con. 40, 238–266.
  * Ennis-King et al. (2005) Ennis-King, J, Preston, I & Paterson, L 2005 Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids (1994-present) 17 (8), 084107\.
  * Graham & Steen (1994) Graham, MD & Steen, PH 1994 Plume formation and resonant bifurcations in porous-media convection. J. Fluid Mech. 272, 67–90.
  * Grossmann & Lohse (2000) Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. J. Fluid Mech. 407, 27–56.
  * Grossmann & Lohse (2001) Grossmann, S. & Lohse, D. 2001 Thermal convection for large Prandtl numbers. Phys. Rev. Lett. 86 (15), 3316\.
  * Guo et al. (2021) Guo, R, Sun, H, Zhao, Q, Li, Z, Liu, Y & Chen, C 2021 A novel experimental study on density-driven instability and convective dissolution in porous media. Geophys. Res. Lett. 48 (23), e2021GL095619.
  * Hassanzadeh et al. (2014) Hassanzadeh, P., Chini, G. P. & Doering, C. R. 2014 Wall to wall optimal transport. J. Fluid Mech. 751, 627–662.
  * Hewitt (2020) Hewitt, DR 2020 Vigorous convection in porous media. Proc. Math. Phys. Eng. Sci. 476 (2239), 20200111.
  * Hewitt & Lister (2017) Hewitt, DR & Lister, JR 2017 Stability of three-dimensional columnar convection in a porous medium. J. Fluid Mech. 829, 89–111.
  * Hewitt et al. (2012) Hewitt, DR, Neufeld, JA & Lister, JR 2012 Ultimate regime of high Rayleigh number convection in a porous medium. Phys. Rev. Lett. 108 (22), 224503.
  * Hewitt et al. (2013a) Hewitt, DR, Neufeld, JA & Lister, JR 2013a Convective shutdown in a porous medium at high Rayleigh number. J. Fluid Mech. 719, 551–586.
  * Hewitt et al. (2013b) Hewitt, DR, Neufeld, JA & Lister, JR 2013b Stability of columnar convection in a porous medium. J. Fluid Mech. 737, 205–231.
  * Hewitt et al. (2014) Hewitt, DR, Neufeld, JA & Lister, JR 2014 High Rayleigh number convection in a three-dimensional porous medium. J. Fluid Mech. 748, 879–895.
  * Hidalgo et al. (2015) Hidalgo, JJ, Dentz, M., Cabeza, Y. & Carrera, J. 2015 Dissolution patterns and mixing dynamics in unstable reactive flow. Geophys. Res. Lett. 42 (15), 6357–6364.
  * Hidalgo et al. (2012) Hidalgo, JJ, Fe, J, Cueto-Felgueroso, L & Juanes, R 2012 Scaling of convective mixing in porous media. Phys. Rev. Lett. 109 (26), 264503.
  * Howard (1966) Howard, L. N. 1966 Convection at high Rayleigh number. In Proceedings of the 11th International Congress of Applied Mechanics, pp. 1109–1115. Springer.
  * Huppert & Neufeld (2014) Huppert, H. E. & Neufeld, J. A. 2014 The fluid mechanics of carbon dioxide sequestration. Annu. Rev. Fluid Mech. 46, 255–272.
  * Kraichnan (1962) Kraichnan, R. H. 1962 Turbulent thermal convection at arbitrary Prandtl number. Phys. Fluids 5 (11), 1374–1389.
  * Letelier et al. (2019) Letelier, JA, Mujica, N. & Ortega, JH 2019 Perturbative corrections for the scaling of heat transport in a Hele-Shaw geometry and its application to geological vertical fractures. J. Fluid Mech. 864, 746–767.
  * Letelier et al. (2023) Letelier, JA, Ulloa, HN, Leyrer, J & Ortega, JH 2023 Scaling CO2–brine mixing in permeable media via analogue models. J. Fluid Mech. 962, A8.
  * Liang et al. (2018) Liang, Y, Wen, B, Hesse, MA. & DiCarlo, D 2018 Effect of dispersion on solutal convection in porous media. Geophys. Res. Lett. 45 (18), 9690–9698.
  * MacMinn et al. (2012) MacMinn, CW, Neufeld, JA, Hesse, MA & Huppert, HE 2012 Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers. Water Resour. Res. 48 (11).
  * Malkus (1954) Malkus, W. V. R. 1954 Discrete transitions in turbulent convection. Proc. R. Soc. A 225 (1161), 185–195.
  * Otero et al. (2004) Otero, J, Dontcheva, LA, Johnston, H, Worthing, RA, Kurganov, A, Petrova, G & Doering, CR 2004 High-Rayleigh number convection in a fluid-saturated porous layer. J. Fluid Mech. 500, 263–281.
  * Pirozzoli et al. (2021) Pirozzoli, S., De Paoli, M., Zonta, F. & Soldati, A. 2021 Towards the ultimate regime in Rayleigh–Darcy convection. J. Fluid Mech. 911, R4.
  * Priestley (1954) Priestley, C. H. B. 1954 Convection from a large horizontal surface. Aust. J. Phys. 7 (1), 176–201.
  * Schrag (2007) Schrag, DP 2007 Preparing to capture carbon. science 315 (5813), 812–813.
  * Simmons et al. (2001) Simmons, C. T., Fenstemaker, T. R. & Sharp Jr, J. M. 2001 Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges. J. Contam. Hydrol. 52 (1-4), 245–275.
  * Slim (2014) Slim, A. C. 2014 Solutal-convection regimes in a two-dimensional porous medium. J. Fluid Mech. 741, 461–491.
  * Song et al. (2024) Song, J., Shishkina, O. & Zhu, X. 2024 Scaling regimes in rapidly rotating thermal convection at extreme Rayleigh numbers. J. Fluid Mech. 984, A45.
  * Spiegel (1963) Spiegel, E. A. 1963 A generalization of the mixing-length theory of turbulent convection. Astrophys. J. 138, 216.
  * Tsinober et al. (2022) Tsinober, A., Rosenzweig, R., Class, H., Helmig, R. & Shavit, U. 2022 The Role of Mixed Convection and Hydrodynamic Dispersion During CO2 Dissolution in Saline Aquifers: A Numerical Study. Water Resour. Res. 58 (3).
  * Tsinober et al. (2023) Tsinober, A, Shavit, U & Rosenzweig, R 2023 Numerical investigation of the influence of hydrodynamic dispersion on solutal natural convection. Water Resourc. Res. p. e2023WR034475.
  * Ulloa & Letelier (2022) Ulloa, HN & Letelier, JA 2022 Energetics and mixing of thermally driven flows in Hele-Shaw cells. J. Fluid Mech. 930, A16.
  * Van Der Poel et al. (2013) Van Der Poel, E. P., Stevens, R. J. A. M. & Lohse, D. 2013 Comparison between two-and three-dimensional Rayleigh–Bénard convection. J. Fluid mech. 736, 177–194.
  * Wen et al. (2018) Wen, B, Chang, KW & Hesse, MA 2018 Rayleigh-Darcy convection with hydrodynamic dispersion. Phys. Rev. Fluids 3 (12), 123801.
  * Wen et al. (2015) Wen, B, Corson, LT & Chini, GP 2015 Structure and stability of steady porous medium convection at large Rayleigh number. J. Fluid Mech. 772, 197–224.
  * Wen et al. (2012) Wen, B., Dianati, N., Lunasin, E., Chini, G. P. & Doering, C. R. 2012 New upper bounds and reduced dynamical modeling for Rayleigh–Bénard convection in a fluid saturated porous layer. Commun. Nonlinear Sci. Numer. 17 (5), 2191–2199.
  * Xia et al. (2023) Xia, K.-Q., Huang, S.-D., Xie, Y.-C. & Zhang, Lu 2023 Tuning heat transport via coherent structure manipulation: Recent advances in thermal turbulence. Natl. Sci. Rev. p. nwad012.