Nitta et al. Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-020-00383-7(2020) 7:68

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

Development of Integrated Land Simulator
Tomoko Nitta1, Takashi Arakawa2, Misako Hatono3,4, Akira Takeshima1 and Kei Yoshimura1*

Abstract
Accurate simulations of land processes are crucial for many purposes, such as climate simulation, weather, flood,and drought prediction, and climate change impact assessment studies. In this paper, we present a new land
simulator called the Integrated Land Simulator (ILS). The ILS consists of multiple models that represent processes
related to land (hereafter, referred to as “land models”). They are coupled by a general-purpose coupler, Jcup, andexecuted using the Multiple Program Multiple Data approach. Currently, ILS includes a physical land surface model,
the Minimal Advanced Treatments of Surface Interaction and Runoff model, and a hydrodynamic model, theCatchment-based Macro-scale Floodplain model, and the inclusion of additional land models is planned. We
conducted several test simulations to evaluate the computational speed and scalability and the basic physicalperformance of the ILS. The results will become a benchmark for further development.
Keywords: Land model, Earth system modeling, General-purpose coupler, Energy and water cycle, Climate
simulation

1 Introduction
Humans live on land, which is one of the most influential elements of the Earth system in the phenomenon ofclimate change. Therefore, it is important to model land
processes as precisely and as detailed as possible to estimate the impact on and from the land. Land surface andhydrological modeling communities have exerted considerable effort in this direction (e.g., van den Hurk et al.
2011; Archfield et al. 2015; Decharme et al. 2019). However, progress in this direction has not been sufficient interms of preciseness (i.e., more precisely represent geological, physical, chemical, and biological processes) ordetailedness (i.e., representing land processes in higher
spatiotemporal resolution). Some breakthroughs in landmodel development have been expected to fill the gap
between societal expectations and reality (Wood et al.
2011). For example, there is a need to establish adaptation strategies for local municipalities (e.g., coastal fisheries, forestry, and water resources), and numericalexperimental results of climate models are often used as
reference. In such cases, the more precise and detailed* Correspondence: kei@iis.u-tokyo.ac.jp1
Institute of Industrial Science, the University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8574, Japan
Full list of author information is available at the end of the article

simulations are for particular land regions, the more useful and beneficial the adaptation strategies are for research as well as for society as a whole.
Furthermore, there are common long-lasting biases inclimate models. One such bias is the hot and dry summer over the mid-latitudinal semi-arid regions (Mueller
and Seneviratne 2014; Hirota et al. 2016). The soil andsnow hydrology elements of climate models seem to exhibit some inadequacies, generating surface fluxes that
may trigger this bias; this issue has not yet been fully resolved. Therefore, improvements in land simulation are
essential.
Climate sensitivity is one of the seven grand challengesof the World Climate Research Programme, and theland surface mechanism also impacts clouds and atmospheric circulation (Bony and Stevens 2012). Similarly,
the land surface mechanism is known as a key process
in other grand challenges, such as carbon feedbacks inthe climate system, weather and climate extremes, andwater which helps supply the food baskets of the world.
Why has land surface and hydrological model development not advanced sufficiently to help solve these issues? It is a complicated problem, as partly reported by
Blöschl et al. (2019) with 23 unsolved problems in hydrology. However, particularly for the improvement of© The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in this article are included in the article's Creative Commonslicence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commonslicence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain
permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

climate modeling, one technical element that prevents
the smooth evolution of the land model is being one ofthe physical parameterizations of an atmospheric model.
This affects both preciseness and detailedness, as explained below.
The preciseness of the land model is realized by betterand more appropriate expressions for all types of land
phenomena. Although there are independent modelingstudies in hydrology, biology, geology, and other relevant
fields, these studies typically require long periods for implementation in climate models or are never incorporated. This is because, historically, the land surface
component in climate models was developed as one ofthe physical parameterizations of the atmospheric general circulation model to produce lower boundary conditions for the atmosphere and because the typicalcoordination of a land surface model is significantly different from those of the more precise independent
models. Furthermore, climate models tend to prioritize
globally impacted processes. In contrast, most land processes tend to impact the local scale; thus, they tend not
to be prioritized in climate modeling.
The detailedness of the land model is even more technical than the preciseness. Although land has inherent
strong heterogeneity, it cannot have a high spatial resolution because the land model is technically attached asa boundary condition of the atmospheric column, which
involves a rather lower resolution. This problem essentially remains in higher-resolution atmospheric models
because the heterogeneity of land is always greater thanthat of the atmosphere. Mosaic or tiled representation ofland-type classification is a simple solution to partially
mimic the behavior at a higher resolution. It should be
noted that the spatial coordinates of the geophysicalmodel are very important for expressing the target behavior. For example, ocean and atmospheric models usually use different spatial coordinates. For land models,there is also a spatial coordinate system that is more appropriate than simply using coordinates more optimized
for atmospheric models, as most stand-alone landmodels (hydrological models, in particular) adopt basinunit coordinates.
In this study, we create a new framework to improve
the Model for Interdisciplinary Research on Climate(MIROC) land surface component. This framework iscalled the Integrated Land Simulator (ILS). The Minimal
Advanced Treatments of Surface Interaction and Runoff(MATSIRO) model was first developed in 2003 (Takata
et al. 2003) as the MIROC land component. It was originally developed as one of the physical parameterizations of the atmospheric general circulation model
previously known as the Center for Climate System Research Atmospheric General Circulation Model, currently known as MIROC, with a vertical one-

Page 2 of 14

dimensional representation of the canopy, snow, and soil
layers. ILS enables the improvement of preciseness in
two ways. First, the MATSIRO code was rewritten in amodern way to enhance readability and userfriendliness. This makes it easier to implement new parameterizations. Second, new models can be coupledwithin the ILS. The extent of benefit from these improvements depend on the type of application; theformer may be more beneficial for snow or soil parameterizations, while the latter may be more beneficial forindependently developed vegetation models. The detailedness can be achieved by using appropriate resolution and grid systems. As the first step of thisframework, we adopt the CaMa-Flood (Catchment-based
Macro-scale Floodplain) hydrodynamic model (Yamazaki
et al. 2011) and coupled these two models (MATSIRO
and CaMa-Flood) in the ILS. The second section presents the ILS framework and components, followed by
the experimental settings of the three experiments. The
third section provides the results of the experiments. Finally, the last section provides the summary andconclusions.

2 Methods
In this section, we describe the concept of ILS, its physical processes, the input/output (I/O) component, coupling system, regridding table calculator, boundary
condition creator, and benchmarking. Then, we explain
three experiments to (1) examine computational speed
and scalability, (2) show the impact of coupling concerning computational reproducibility, and (3) provide an example of benchmarking.
2.1 Framework and components
2.1.1 Concept of ILS

Figure 1 shows a conceptual diagram of ILS. Multiple
terrestrial models are coupled using a general-purpose
coupler Jcup (Arakawa et al. 2020). Here, the ILS consists of two land models: a physical land surface model,
MATSIRO (Takata et al. 2003; Nitta et al. 2014), and ahydrodynamic model, CaMa-Flood (Yamazaki et al.
2011). An I/O component, used to input and output
files, is developed as an independent executable. The ILS
employs the Multiple Program Multiple Data (MPMD)
approach and executes these models and the I/O component in parallel. Each model can run in its preferredspatial coordinate and resolution. The physical land surface model, MATSIRO, is managed within the ILS because it models essential land processes. However, there
are independent groups that have developed other landmodels. In such a case, the coupling procedure is as follows: (1) port the latest version of the component model,(2) employ the Jcup Application Programming Interface(API) to be coupled with ILS, and (3) prepare regridding

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

Page 3 of 14

Fig. 1 Concept of Integrated Land Simulator

tables for spatial interpolation. If there is an updated
version of that model, steps (1) and (2) are repeated.
This approach allows the latest versions of each landmodel to be included in the ILS with minimal code
modifications. We can use the same procedure for a newdifferent land model, if necessary. As Jcup has been used
to couple an atmospheric general circulation model withan oceanic general circulation model (Miyakawa et al.
2017), it is relatively easy to couple the ILS to the atmospheric and oceanic general circulation models. Previously, offline land simulations were conducted bymanually modifying the code. However, ILS enables siteexperiments, global offline experiments, and coupled climate experiments in a single framework without any
modification of the code. Our future goal is to contribute to improving climate model simulations by incorporating more sophisticated land models with rapid
development and validation cycles. The details of eachmodel of the ILS are described in the following subsections in this section.
To speed up the development and validation cycles,
we decided to employ a software management environment. We manage the codes through version-control
software and use GitLab (https://about.gitlab.com/) toshare information via Wiki, manage issues, and test thecode automatically when a merge is requested. Continuous testing is essential for model development. Currently, a 2-day simulation of the global offline
experiment and a site experiment are performed, and weare planning to add benchmarking tools as well. The
benchmark system is explained in a separate subsectionof this section.
2.1.2 Physical processes

Currently, ILS includes two physical models: MATSIRO
and CaMa-Flood. As explained above, MATSIRO has
been adopted as the MIROC land surface component. It
consists of a single-layer canopy and multiple userdefined layers of snow and soil and employs a tilescheme. MIROC participated in several phases of the
Coupled Model Intercomparison Project. As one of themulti-models, many studies have evaluated its simulated
land-related variables, such as surface air temperature
and/or precipitation over land (Mueller and Seneviratne
2014; Miao et al. 2014), land-atmosphere interactions(Ukkola et al. 2018), evapotranspiration (Lian et al.
2018), surface albedo (Wang et al. 2016), snow cover(Thackeray et al. 2015), and permafrost (Koven et al.

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

2013). It has also been used in studies assessing water
resources and flood impact, which combine offline experiments with observation-based meteorological data as
inputs (Yoshimura et al. 2018).
First, we rewrote the MIROC5 version of MATSIRO(Watanabe et al. 2010; hereafter referred to as MATS
IRO5) to improve user-friendliness such that new
schemes are easier to implement, and thus, the preciseness can be improved. We decoupled MATSIRO from
MIROC5, removed all unnecessary parts produced by
the decoupling, and separated the rest into MATSIRO’s
library and driver. The library includes only calculationsof land processes, whereas the driver includes calendar
management, parameter settings, reading and writing restart files, and subroutines for coupling with other landmodels. Users can modify the driver to conduct different
types of simulations; ILS can handle any combination ofland cover tiles, soil type tiles, and elevation tiles bymodifying the driver. Users will choose which tiles touse and how to use them depending on their purposes.
The original MATSIRO is a vertical one-dimensional

Page 4 of 14

model and has an IJ loop, a two-dimensional loop in
space, in the innermost part of the code. We replaced
this IJ loop with an expandable universal loop, which
can be applied for multiple purposes, not just IJ but also
ensemble members or land cover tiles. This simplifies
the code, improves readability, and is advantageous forvector processors.
To confirm the reproducibility of the rewritten code,
we performed two experiments with the same driving
data provided by Plumbing of Land Surface Model(PALS; Best et al. 2015; Abramowitz 2012) using the ILS
and MATSIRO5. The experimental settings are described in the experimental section. We examined sevenmain output variables: soil temperature and moisture in
the first layer, snow water equivalent, canopy water,
sensible heat flux, latent heat flux, and runoff for 20 sites
with a 30-min temporal resolution. The results show bitidentical reproducibility in the single-point floatingpoint data type in most cases, and the differences are
small enough compared to their standard deviations(Fig. 2).

Fig. 2 a The normalized standard deviation of 7 variables simulated by the MIROC framework. The vertical axis shows soil moisture andtemperature in the first layer, snow water equivalent, canopy water, latent heat flux, sensible heat flux, and runoff. The horizontal axis shows 20
Fluxnet sites. The standard deviation values are normalized by their mean. b As in a, but for the difference between MIROC framework and ILS
normalized by mean

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

We implemented iterative calculations in the flux calculation to solve the well-known instability of the skintemperature and surface fluxes. The surface heat balance
equation can be written as follows:Δg ¼ H g þ Rg þ E g þ E s þ E sn − F g − F sn
Δc ¼ H c þ Rc þ E i þ E is þ E t
where Δg and Δc are the energy divergences at the
ground surface and canopy, respectively; Hg and Hc are
the respective sensible heat fluxes; Rg and Rc are the respective net radiations; Eg, Es, Esn, Ei, Eis, and Et are thebare ground evaporation flux, bare ground sublimation
flux, snow sublimation flux, latent heat flux from canopy
interception, sublimation flux from canopy snow interception, and transpiration flux, respectively; and Fg and
Fsn are the heat conduction fluxes to the soil and snow
surface, respectively. If there is no melting of the groundsurface, the above equations are solved by linearizing theskin temperatures Tg and Tc of the bare ground and the
canopy, respectively. If the change in Tg or Tc from theprevious time step is large, the error of the solution becomes large. In the initial verification of this study, there
were some sites where the average daily variation fluctuated. We considered that one of the causes of the calculation instability might be that the skin temperatureassumed an unrealistic value owing to an error in theheat balance calculation. Therefore, we implemented the
Newton–Raphson method in MATSIRO.
CaMa-Flood (Yamazaki et al. 2011, 2013) is a globalhydrodynamic model that uses daily river runoff simulated by land surface models as an input. It is currently
one of the most widely used hydrodynamic models,owing to its unique method of representing flood-plain
inundation using subgrid topographical parameters. In
CaMa-Flood, an adaptive timestep based on the
Courant-Friedrichs-Lewy condition is used. Aside fromthe implementation of modules necessary for the coupling system, the only modifications to the original
CaMa-Flood model are related to acquiring daily river
runoff data and outputting the simulated results. To input the river runoff data, Jcup is used, and only unit
conversion is conducted in the modified model. Outputting the simulated results is conducted by the I/O component, rather than the original modules.
2.1.3 I/O component

In developing the ILS, we decided to develop an I/Ocomponent for reading and writing files independently.
The I/O component is included in the MPMD programas an independent executable. Models can avoid therelatively time-consuming file input/output, such thatthe overall execution time can be reduced. The I/Ocomponent includes a user interface for temporal

Page 5 of 14

interpolations that users can implement freely. The gridinterpolation function of the Jcup application programming interface (API) permits the use of input data withdifferent spatial resolutions from that of the experiment.
For example, we can use meteorological data with different spatial resolutions for an experiment. It is also possible to combine simulated results for only land gridswith missing values for ocean grids before outputting
the results. The function allows single-process execution
and two options with multiple-process (MP) executions.
The first MP option uses the Message Passing Interface(MPI) I/O function. In the second option, the root processor collects data from other processors and outputs
them. The I/O component supports the binary and Network Common Data Form (NetCDF) formats. It is designed to be extendable to other output formats such as
gtool. For NetCDF, metadata are prepared in JavaScript
Object Notation (JSON) format. File names and metadata follow the Assistance for Land Surface Modeling
Activities (ALMA) and Climate and Forecast (CF)conventions.

2.1.4 Coupling system

As described above, the ILS is designed to perform simulations by coupling a plurality of components. This section describes the coupling system built into the ILS.
The coupling software used in ILS is Jcup (Arakawa
et al. 2020). Jcup can couple multiple components combining both parallel and serial implementations. In
addition, it has high versatility because there is no restriction on applicable grid coordinates, and users canfreely implement interpolation schemes. However, because of its high flexibility, the interface is complicated,making it difficult for beginners to use. Therefore, an
API with a more compact interface customized to the
ILS was developed. The most distinguishable feature of
Jcup is that it requires a correspondence between thesending-side grid point index, receiving-side grid pointindex, and interpolation coefficients in the interpolation
calculation as input information (hereafter, referred to as“regridding table”). Figure 3 shows a typical use of the
Jcup API subroutines. Before the time-integration loop,
there are three types of subroutines (and moj_put_data),such as the initialization of the coupler, setting of grid
information, and setting of the regridding table. In thetime-integration loop, moj_set_time is called near the
beginning of the loop. Data exchange and interpolation
calculations were performed inside this subroutine. Subroutine moj_get_data obtains the data from the sending
component, and moj_put_data is a subroutine thattransfers the data to be sent to the receiving component
to the coupler. These two subroutines can be called anywhere in the time-integration loop.

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

Page 6 of 14

Fig. 3 Example of typical usage of the Jcup API subroutines. The blue squares represent the subroutines that can be called multiple times. The
light blue squares indicate the subroutines that can be called only once. The red arrows indicate that communication is performed in thesubroutines. The cream squares represent operations inside moj_set_time. Data exchange and interpolation calculations were performed
inside moj_set_time

As described above, Jcup and its API require a regridding table as input information; therefore, the user needs
to calculate a regridding table in advance. In previous
cases, the calculation tool was developed individually forthe models to be coupled. For example, to couple theicosahedral atmospheric model NICAM (Satoh et al.
2014) and the global ocean model COCO (Hasumi2015), a program that searches for grid points includedin each polygon of the opposite model grid and calculates the area of the overlapped region has been developed (Arakawa et al. 2014). For the seismic model
Seism3D (Furumura 2005) and the structure model
FrontISTR++ (Okuda 2019) coupling, the regridding
table could be calculated easily because the grid structure of the seismic model was a simple rectangular gridwith a constant grid spacing (Matsumoto et al. 2015).
Unlike these cases, the ILS has the possibility of coupling various models that have their own grid coordinates,
and the combination of coupling should be flexible;
therefore, a regridding table calculation tool applicable
to a wide range of grid coordinates is required. Hence, ageneral-purpose regridding table calculation tool called

Spheroidal Coordinates Regridding Interpolation Table
Generator (SPRING) has been developed. The next section details this software.
2.1.5 Regridding table calculator

SPRING (Takeshima 2020) has been developed tocouple models with any type of grid coordinate system,
such as the latitude–longitude grid, polygonal grid, or
grid defined by an aggregation of tiny latitude–longitude
cells (Fig. 4). The basin-shaped coordinate system of the
CaMa-Flood is of the last type. This tool provides regridding tables, a dataset of overlapping grid indices, andtheir interpolation weights for these coordinates. In conservative data communication in the ILS, the data value
XIJ in the receiving-side grid point IJ is calculated by thefollowing equation:
X IJ ¼

X
ij

0

cij xij ; cij ¼

sij
S IJ

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

Page 7 of 14

Fig. 4 Examples of grids to which SPRING is applicable. a Latitude–longitude grid, b NICAM grid (Satoh et al. 2008), a type of convex polygon
grid, c concave polygon grid, and d basin-shaped grid, a type of grid defined by the aggregation of tiny cells

where ij is the sending-side grid point index that overlaps with the grid IJ, cij is the interpolation weight, xij is0
the data value, sij is the overlapping area between grid
IJ and grid ij, and SIJ is the area of the grid IJ. In thedefault configuration, MATSIRO and CaMa-Flood
have basin-shaped coordinates, whose grids consist of
tiny latitude–longitude cells, and the forcing data
used in this study are expressed on a latitude–longi0
tude grid. Therefore, sij and SIJ are calculated by thefollowing equation regarding Earth as a perfect spherewith radius a

Fig. 5 Concept of boundary condition creatorS¼

πa2
ð sinφ1 − sinφ2 Þðθ2 − θ1 Þ
180

where S is the area of a rectangular shape, θ1 and θ2
are the longitudes of the western/eastern boundaries,
respectively, and φ1and φ2 are the latitudes of thenorthern/southern boundaries, respectively. The detailed evaluation of regridding tables will be presented
in a separate paper and is beyond the scope of thisstudy.

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

2.1.6 Boundary condition creator

The ILS includes a tool to create the boundary conditions of MATSIRO. They were manually created for different experiments. However, the information was not
sufficiently shared, and sometimes the data included human errors. To avoid these issues, we manage the version of the codes to create boundary conditions. Thus,
we can use the common parts repeatedly. Figure 5 showsa conceptual diagram of the boundary creation tool. By
combining the original datasets and regridding tables described in the previous section, we can create newboundary condition data for any horizontal resolution,
which include land cover type, soil type, soil albedo,
standard deviation and gradient of elevation, and leaf
area index (LAI). First, the land cover categories in thehigh-resolution original data are reclassified into categories to be used in MATSIRO. The regridding table isthen used to calculate the dominant categories. For soil
types, the proportions of sand, silt, and clay are interpolated using a regridding table. The soil type is then determined using the triangle diagram of the soil
properties. The soil albedo is simply interpolated. The
standard deviation and slope of elevation are calculated
from high-resolution topographic data. The LAI is derived by a simple interpolation without the tilingscheme. If a land cover tiling scheme is used, the LAI isfirst classified by land use data and then averaged. The
boundary tool creator allows the implementation of aPage 8 of 14

new method, as the method varies depending on thesource data.
2.1.7 Benchmarking

ILS incorporates benchmarking as part of the model development cycle. It is important to validate the model
against a variety of terrestrial observations during eachmodel development cycle and to evaluate how new parameterizations or models for preciseness and changes inresolution for detailedness can improve the simulated results. In this paper, we present the results of an experiment using PALS, which includes 20 Fluxnet sites (Fig. 6).
The experiment is suitable for benchmarking because thesites cover the locations with various conditions, and its
computational cost is low. The experimental settings areshown in the next section.
2.2 Specification of three experiments

First, we conducted global offline land simulations andevaluated the computational performance. The ILS wasforced using three-hourly meteorological datasets prepared by the Global Soil Wetness Project Phase 3(GSWP3; Kim et al. 2019). GSWP3 forcing data is basedon downscaling of the Twentieth-Century Reanalysis(Compo et al. 2011), including rainfall, snowfall,temperature, wind speed, specific humidity, shortwave
and longwave radiation, and surface pressure. Cloud
cover from the same downscaling experiment was also

Fig. 6 Locations of 20 Fluxnet sites in PALS land cover types in 0.5° grids are shown on the map. Circles on the map show the locations withland cover types in color. There are 11 land cover types: 1: continental ice, 2: broadleaf evergreen forest, 3: broadleaf deciduous forest andwoodland, 4: mixed coniferous and broadleaf deciduous forest and woodland, 5: coniferous forest and woodland, 6: high-latitude deciduous
forest and woodland, 7: wooded c4 grassland, 8: shrubs and bare ground, 9: tundra, 10: cultivation, 11: desert

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

used to partition shortwave radiation into visible andnear-infrared, and direct and diffusive components. In
this experiment, ILS consists of MATSIRO, CaMaFlood, and I/O components. The temporal time step of
MATSIRO was 1 h. The soil consists of six layers withdepths of 5, 20, 75, 100, 200, and 1000 cm, and the snow
consists of a maximum of three layers. The number ofsnow layers changes depending on the amount of snow.
Figure 7 shows the computational flow of the global
offline experiment. The I/O component reads the meteorological forcing and boundary data, interpolates
them temporally, and sends them to MATSIRO. Then,
MATSIRO calculates land processes and sends simulated runoff to the CaMa-Flood routine, where the riverprocesses are calculated. The output data from MATS
IRO and CaMa-Flood are sent to the I/O component tobe written to files. The experiment was run on a computer cluster with Intel Xeon processors at the University of Tokyo. We changed the number of cores for acomponent from 1 to 128 while fixing the number ofcores (16 cores) for the other two components.
Second, we examined the difference between
uncoupled (stand-alone) and coupled CaMa-Flood simulations to check whether coupling, including the regridding tables created by SPRING, works as expected. Notethat we did not consider feedback from CaMa-Flood to
MATSIRO in the coupled simulation and take theuncoupled simulation as the truth. We conducted astand-alone river discharge simulation using the original

Page 9 of 14

CaMa-Flood with runoff data from MATSIRO in ILS(uncoupled simulation). We compared it with the simulation mentioned above (coupled simulation). This experiment intends to confirm if the coupled simulation is
as similar to the uncoupled simulation by CaMa-Flood,in which its original regridding scheme is implemented.
Finally, we conducted site experiments for 20 Fluxnet
sites from PALS, as shown in Fig. 6, to (1) validate the rewritten code and (2) show the results of the first benchmark. The sites include tropical, arid, temperate, and cold
regions, as assigned by the Koppen climate classification(Peel et al. 2007), and forest and non-forest regions. The
dataset includes meteorological forcing data: downward
shortwave and longwave radiation, temperature, precipitation, specific humidity, wind speed, and surface pressure.
In addition to these data, cloud cover from ERA5 (C3S
2017) was used. PALS also includes surface flux data fromthe observations and estimated by one-, two-, andthree-variable regression equations for validation andbenchmarking. The time steps of the input data andsimulation were 30 min without any missing data. The
simulation periods range from 2 to 10 years, dependingon the site. We used the MIROC5 version of the ILS
for the verification shown above. Then, we implemented the iteration of the flux calculation to MATS
IRO as described above and evaluated the PALS experiment. Further observation stations, validation of global
experiments, and increasing the amount of observation
data are clearly in our future work.

Fig. 7 The computational flow of a global experiment with MATSIRO, CaMa-Flood, and IO components, executed with the MPMD approach

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

3 Results and discussion
In this section, we show the results from the three experiments explained in the previous section.
3.1 Computational speed and scalability

Figure 8a shows the total time required to perform global simulations with different numbers of cores. The 1month simulation with 16 cores for MATSIRO, CaMaFlood, and I/O, respectively (48 cores in total), takes approximately 590 s. When we changed the number ofcores for MATSIRO (blue) while fixing the number ofcores for CaMa-Flood and the I/O component, thesimulation with one core takes approximately 9430 s,
and the computational time decreases with the numberof cores to approximately 100 s with 128 cores. In contrast, the differences in computational time are small
when changing the number of cores for CaMa-Flood(orange) and the I/O component (green) while fixing thenumber of cores for the other two components. This isbecause MATSIRO has the highest computational costand is a bottleneck in these simulations.
Then, we examined the total time spent on each component of MATSIRO: initial setup, main MATSIRO calculation, communication by the coupler, receiving the
input data from the coupler, and sending the output
data to the coupler. In this analysis, we only used the experiments with fixed 16 cores for the CaMa-Flood and I/
O components, respectively, and increasing the numberof cores for MATSIRO from 1 to 128 cores (Fig. 8b).
The initial setup with 64 cores took longer than with 32
cores; the initial setup with multiple nodes seems to take

Page 10 of 14

more time. The coupler communication cost is relatively
fixed. In all settings, the MATSIRO main calculation is
the most expensive component; it accounts for approximately 95% of the MATSIRO computation with 16
cores. However, it shows close to perfect scalability suchthat we would expect ILS to become faster with morecores. Although the computational performance depends
on machines, we believe that high-resolution computation with a high-performance computer is promising.
The validation of the global offline experiment is beyond
the scope of this study.
3.2 Impact of coupling with respect to computational
reproducibility

The impact of coupling CaMa-Flood to ILS with respect
to computational reproducibility is displayed in Fig. 9.
Figure 9a shows the standard deviation of river discharge, flood-plain discharge, river storage, flood-plain
storage, and flooded area fraction from stand-alone
CaMa-Flood simulation normalized by the mean foreach 10-degree band from 60° S to 90° N. The maximum
value is approximately 11, found in flood-plain discharge
at 30–40° N. The minimum value is approximately 1,
found in the flooded area fraction at 60–50°S. Figure 9b
is the same as Fig. 9a, but for the difference between
coupled and uncoupled simulations, also normalized by
the mean. The maximum difference is found in floodplain discharge at 50–60° N, which is less than 0.001% of
the standard deviation. The differences were due toregridding tables, I/O, and simulation framework andare negligible compared to the standard deviation.

Fig. 8 a Computational speed and scalability. Relationship between number of cores for each component (x axis) and computational time (yaxis). We only changed the number of cores from 1 to 128 for one component and fixed the number of cores for the other two components to
16. b Computational speed of MATSIRO’s component in relation to the number of cores

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

Page 11 of 14

Fig. 9 a The standard deviation of uncoupled simulation normalized by mean. The vertical axis shows river discharge, flood-plain discharge, riverstorage, flood-plain storage, and flooded area fraction. The horizontal axis shows latitudes. b As in a but for the difference normalized by their mean

Fig. 10 Taylor diagram for half-hourly a latent and b sensible heat fluxes. Blue shows ILS results and gray shows two-variable regression.
Standard deviation and root-mean-square difference are standardized by the standard deviation of observations

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

3.3 First result of benchmarking

Figure 10 summarizes the performance of the halfhourly latent and sensible heat fluxes using a Taylor
plot. Blue indicates ILS simulations, and each plot showseach site. For reference, the results of linear regressionsof the downward shortwave radiation and surface airtemperature against surface fluxes from the PALS project are shown in gray. The statistics were normalized by
the standard deviation of the reference data.
Figure 10a shows the results of the latent heat flux.
The median correlation coefficients were 0.78 and 0.83,
and the median root-mean-square differences (RMSDs)
were 0.66 and 0.59 for the ILS and empirical methods,
respectively. For these two metrics, the empiricalmethod showed better reproducibility. However, the median standard deviations were 0.93 and 0.74 for the ILS
and the empirical methods, respectively. ILS results arecloser to unity, which is the normalized standard deviation of the reference data. The empirical approach
underestimated the standard deviation in most places.
Figure 10b shows the results of the sensible heat flux.
In the ILS and empirical methods, the median correlation coefficients were 0.84 and 0.89, the median RMSDs
were 0.60 and 0.50, and the median standard deviations
were 1.14 and 0.90, respectively. The empirical methodtended to be better for all three metrics.
Comparing the results of latent and sensible heat
fluxes, the reproducibility of sensible heat fluxes tends tobe better than that of latent heat fluxes. The result ismore scattered in the sensible heat flux than the latent
heat flux. The detailed analysis is beyond the scope ofPage 12 of 14

this study and will be shown in a separate validation
paper.
The empirical method had a problem with the reproducibility of extreme values. For example, when the
99.9th percentile values of the latent heat and sensible
heat flux at each site were examined for 30 min, the ILS
tended to produce better results (Fig. 11). The mean absolute errors of latent heat flux were 77.1 W/m2 and203.7 W/m2 for the ILS and empirical equation, respectively, and the mean absolute errors of sensible heat flux
were 76.8 W/m2 and 149.9 W/m2 for the ILS and empirical equation, respectively. The ILS also showed a better
tendency for the 99th and 99.99th percentile values.

4 Summary and conclusions
The improvement of the land surface model is expected
not only by the climate community but also by society;however, there are two characteristics that are necessary
for such a feat—preciseness and detailedness. Preciseness can be achieved by including geological, physical,
chemical, and biological processes and their interactions
into the comprehensive model. Detailedness can beachieved by computing a land model simulation withhigh spatiotemporal resolution, considering the inherent
heterogeneous conditions of the land. A barrier to realizing these two targets is that the land model was developed as one of the physical parameterizations of theatmospheric model. Therefore, in this study, we developed a new framework to enhance the speed of the improvement of our land surface model. We aim toFig. 11 Absolute biases of 99.9th percentile values of latent and sensible heat fluxes from PALS experiments. Blue shows ILS results and grayshows two-variable regression

Nitta et al. Progress in Earth and Planetary Science(2020) 7:68

develop a new framework to facilitate the realization of
improved preciseness and detailedness in a land model
simulation, a framework named ILS. Note that thispaper only addresses a way to construct the framework.
Further application of the framework, for example, adding new land models to the ILS for better preciseness or
executing the ILS in a much higher spatiotemporal resolution using heterogeneous boundary conditions, will be
conducted in the future.
First, we removed the land surface model, MATSIRO,
from the climate model, MIROC. The vertical onedimensional core of MATSIRO is used as a component
in the ILS. Then, the code was verified. The new hydrodynamic model, CaMa-Flood, was implemented in ILS.
One of the advantages of the ILS is to port a new standalone model almost as it is, using the new generalpurpose coupler, Jcup. Furthermore, we developed ancillary but important tools, such as a component to handle
input and output files, modules to create necessary
boundary conditions automatically, and modules to generate regridding tables to interpolate the geospatial information. Next, we conducted some numericalexperiments to confirm two aspects of performance
evaluation—computational efficiency and validity. Computational efficiency was evaluated by conducting global
simulations using different combinations of CPUs and
ILS components, that is, MATSIRO, CaMa-Flood, and I/
O, and good scalability was demonstrated for simulations at very high resolutions with a many-core architecture. The validation of the land model by comparing itwith reality is still preliminary; however, for selectedsites used in the PALS project, the ILS was satisfactory.
As mentioned above, this paper presents the first description of the new framework. Applications of the ILS
have already been planned in emerging fields, such as
human activities (irrigation/pumping/damming; Hanasaki et al. 2008), sediment dynamics (Hatono and Yoshimura 2020), groundwater dynamics (Miura andYoshimura 2020), water temperature and water quality
in rivers (Tokuda et al. 2019), vegetation dynamics, subgrid representability of hillslope processes, and isotopic
representations (Yoshimura et al. 2006). We believe thatthe basic development conducted in this study is needed
to further advance land/climate models and improve theunderstanding and prediction of the earth system.
Acknowledgements
We thank Dr. Hisashi Yashiro for the discussion. We thank Dr. Gab
Abramowitz for providing the PALS data.
Authors’ contributions
TN designed and coded the new framework and carried out the ILS
experiments. TA developed the Jcup API and I/O. MH coupled the CaMa-Floodmodel and carried out a stand-alone CaMa-Flood simulation. AT developed theregridding table tool. KY proposed the topic and conceived and designed thestudy. The authors have read and approved the final manuscript.

Page 13 of 14

Funding
This work was supported by the Integrated Research Program for Advancing
Climate Models (TOUGOU) Grant Number JPMXD0717935457, ArCS II(Program Grant Number JPMXD1420318865), and DIAS from the Ministry of
Education, Culture, Sports, Science and Technology (MEXT), Japan, and theJapan Society for the Promotion of Science (JSPS) via Grants-in-Aid
18H03794, 18J00585, and 16H06291, JRPs-LEAD with DFG .
Availability of data and materials
The source code used in the current study is available in the repository at
http://isotope.iis.u-tokyo.ac.jp:8000/gitlab/ils-dev/ILS. Please contact thecorresponding author to request access to the repository. The datasets
analyzed during the current study are available from the corresponding
author upon reasonable request.
Competing interests
The authors declare that they have no competing interests.
Author details
1
Institute of Industrial Science, the University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8574, Japan. 2Research Organization for Information
Science and Technology (RIST), 1-18-16 Hamamatsu-cho, Minato-ku, Tokyo
105-0013, Japan. 3School of Engineering, Tohoku University, 6-6-06 Aoba,
Aramaki, Aoba-ku, Sendai 980-8579, Japan. 4Graduate School of Advanced
Science and Engineering, Hiroshima University, 1-4-1 Kagamiyama,
Higashi-Hiroshima City, Hiroshima 739-8527, Japan.
Received: 11 March 2020 Accepted: 23 October 2020

References
Abramowitz G (2012) Towards a public, standardized, diagnostic benchmarking
system for land surface models. Geosci Model Dev 5:819–827
Arakawa T, Inoue T, Sato M (2014) Performance evaluation and case study of acoupling Software ppOpen-MATH/MP. Int Conference Comput Sci 29:924–
935. https://doi.org/10.1016/j.procs.2014.05.083
Arakawa T, Inoue T, Yashiro H, Satoh M (2020) Coupling library Jcup3: its
philosophy and application. Prog Earth Planet Sci 7:6. https://doi.org/10.1186/
s40645-019-0320-z
Archfield SA, Clark M, Arheimer B, Hay LE, McMillan H, Kiang JE, Seibert J, Hakala
K, Bock A, Wagener T, Farmer WH, Andréassian V, Attinger S, Viglione A,
Knight R, Markstrom S, Over T (2015) Accelerating advances in continental
domain hydrologic modeling. Water Resour Res 51:10078–10091. https://doi.
org/10.1002/2015WR017498
Best MJ, Abramowitz G, Johnson HR, Pitman AJ, Balsamo G, Boone A, Cuntz M,
Decharme B, Dirmeyer PA, Dong J, Ek M, Guo Z, Haverd V, van den Hurk BJJ,
Nearing GS, Pak B, Peters-Lidard C, Santanello JA Jr, Stevens L, Vuichard N(2015) The plumbing of land surface models: benchmarking model
performance. J Hydrometeor 16:1425–1442. https://doi.org/10.1175/JHM-D14-0158.1
Blöschl G, Bierkens MFP, Chambel A, Cudennec C, Destouni G, Fiori A, Kirchner
JW, McDonnell JJ, Savenije HHG, Sivapalan M, Stumpp C, Toth E, Volpi E, Carr
G, Lupton C, Salinas J, Széles B, Viglione A, Aksoy H, Allen ST, Amin A,
Andréassian V, Arheimer B, Aryal SK, Baker V, Bardsley E, Barendrecht MH,
Bartosova A, Batelaan O, Berghuijs WR et al (2019) Twenty-three unsolved
problems in hydrology (UPH) – a community perspective. Hydrological
Sciences Journal 64:1141–1158. https://doi.org/10.1080/02626667.2019.
1620507
Bony S, Stevens B (2012) Clouds, circulation and climate sensitivity: how theinteractions between clouds, greenhouse gases and aerosols affecttemperature and precipitation in a changing climate
Compo GP, Whitaker JS, Sardeshmukh PD, Matsui N, Allan RJ, Yin X, Gleason BE,
Vose RS, Rutledge G, Bessemoulin P, Brönnimann S, Brunet M, Crouthamel RI,
Grant AN, Groisman PY, Jones PD, Kruk MC, Kruger AC, Marshall GJ, Maugeri
M, Mok HY, Nordli Ø, Ross TF, Trigo RM, Wang XL, Woodruff SD, Worley SJ(2011) The Twentieth Century Reanalysis Project. QJR Meteorol Soc 137:1–28.
https://doi.org/10.1002/qj.776
Copernicus Climate Change Service (C3S) (2017) ERA5: Fifth generation of ECMW
F atmospheric reanalyses of the global climate, Copernicus Climate Change
Service Climate Data Store (CDS), vol 2019. https://cds.climate.copernicus.eu/
cdsapp#!/homeNitta et al. Progress in Earth and Planetary Science(2020) 7:68

Decharme B, Delire C, Minvielle M, Colin J, Vergnes J-P, Alias A, Saint-Martin D,
Séférian R, Sénési S, Voldoire A (2019) Recent changes in the ISBA-CTRIP land
surface system for use in the CNRM-CM6 climate model and in global offline hydrological applications. J Adv Model Earth Syst 11:1207–1252. https://doi.org/10.1029/2018MS001545
Furumura T (2005) Large-scale parallel simulation of seismic wave propagationand strong ground motions for the past and future earthquakes in Japan. J
Earth Simulator 3:29–38. https://doi.org/10.32131/jes.3.29
Hanasaki N, Kanae S, Oki T, Masuda K, Motoya K, Shirakawa N, Shen Y, Tanaka K(2008) An integrated model for the assessment of global water resources –Part 1: model description and input meteorological forcing. HESS 12:1007–1025
Hasumi H (2015) CCSR Ocean Component Model (COCO) version 4.0. In: Centerfor climate system research report. https://ccsr.aori.u-tokyo.ac.jp/~hasumi/
COCO/coco4.pdf. Accessed 5 Mar 2020
Hatono M, Yoshimura K (2020) Development of a global sediment dynamics
model. Prog Earth Planet Sci. https://doi.org/10.1186/s40645-020-00368-6
Hirota N, Takayabu YN, Hamada A (2016) Reproducibility of summer precipitation
over Northern Eurasia in CMIP5 multiclimate models. J Climate 29:3317–3337.
https://doi.org/10.1175/JCLI-D-15-0480.1
Kim H, Watanabe S, Chang E-C, Yoshimura K, Compo GP, Hirabayashi Y,
Famiglietti JS, Oki T (2019) A century-long global surface meteorology foroffline terrestrial simulations. http://hydro.iis.u-tokyo.ac.jp/GSWP3/. Accessed
5 Mar 2020
Koven CD, Riley WJ, Stern A (2013) Analysis of pthermal dynamics and response
to climate change in the CMIP5 Earth system models. J Climate 26:1877–
1900. https://doi.org/10.1175/JCLI-D-12-00228.1
Lian X, Piao S, Huntingford C, Li Y, Zeng Z, Wang X, Ciais P, McVicar TR, Peng S,
Ottlé C, Yang H, Yang Y, Zhang Y, Wang T (2018) Partitioning global landevapotranspiration using CMIP5 models constrained by observations. Nature
Clim Change 8:640–646
Matsumoto M, Arakawa T, Kitayama T, Mori F, Okuda H, Furumura T, Nakajima K(2015) Multi-scale coupling simulation of seismic waves and building
vibrations using ppOpen-HPC. Procedia Comput Sci 51:514–1523. https://doi.
org/10.1016/j.procs.2015.05.341
Miao C, Duan Q, Sun Q, Huang Y, Kong D, Yang T, Ye A, Di Z, Gong W (2014)
Assessment of CMIP5 climate models and projected temperature changes
over Northern Eurasia. Environ Res Lett 9:055007. https://doi.org/10.1088/
1748-9326/9/5/055007
Miura Y, Yoshimura K (2020) Development and verification of a threedimensional variably saturated flow model for assessment of future global
water resources. J Adv Model Earth Sy. https://doi.org/10.1029/
2020MS002093
Miyakawa T, Yashiro H, Suzuki T, Tatebe H, Satoh M (2017) A Madden-Julian
Oscillation event remotely accelerates ocean upwelling to abruptly terminate
the 1997/1998 super El Niño. Geophys Res Lett 44:9489–9495. https://doi.
org/10.1002/2017GL074683
Mueller B, Seneviratne SI (2014) Systematic land climate and evapotranspirationbiases in CMIP5 simulations. Geophys Res Lett 41:128–134. https://doi.org/10.
1002/2013GL058055
Nitta T, Yoshimura K, Takata K, O’ishi R, Sueyoshi T, Kanae S, Oki T, Abe-Ouchi A,
Liston GE (2014) Representing variability in subgrid snow cover and snow
depth in a global land model: offline validation. J Climate 27:3318–3330.
https://doi.org/10.1175/JCLI-D-13-00310.1
Okuda H (2019) Nonlinear structual analysis open software FrontISTR
Peel MC, Finlayson BL, McMahon TA (2007) Updated world map of the KöppenGeiger climate classification. HESS 11:1633–1644
Satoh M, Tomita H, Yashiro H, Miura H, Kodama C, Seiki T, Noda AT, Yamada Y,
Goto D, Sawada M, Miyoshi T, Niwa Y, Hara M, Ohno T, Iga S-I, Arakawa T,
Inoue T, Kubokawa H (2014) The non-hydrostatic icosahedral atmosphericmodel: description and development. Prog Earth Planet Sci 1:18
Takata K, Emori S, Watanabe T (2003) Development of the minimal advanced
treatments of surface interaction and runoff. Glob Planet Change 38:209–222
Takeshima A (2020) SPheroidal Coordinates Regridding INterpolation table
Generator (SPRING). http://hydro.iis.u-tokyo.ac.jp/~akira/page/Works/
contents/SPRING. Accessed 2 Jul 2020.
Thackeray CW, Fletcher CG, Derksen C (2015) Quantifying the skill of CMIP5
models in simulating seasonal albedo and snow cover evolution. J Geophys
Res Atmos 120:5831–5849
Tokuda D, Kim H, Yamazaki D, Oki T (2019) Development of a global river watertemperature model considering fluvial dynamics and seasonal freeze-thaw

Page 14 of 14

cycle. Water Resour Res 55:1366–1383. https://doi.org/10.1029/
2018WR023083
Ukkola AM, Pitman AJ, Donat MG, De Kauwe MG, Angélil O (2018) Evaluating thecontribution of land-atmosphere coupling to heat extremes in CMIP5
models. Geophys Res Lett 45:9003–9012. https://doi.org/10.1029/
2018GL079102
van den Hurk B, Best M, Dirmeyer P, Pitman A, Polcher J, Santanello J (2011)
Acceleration of land surface model development over a decade of glass. Bull
Amer Meteor Soc 92:1593–1600. https://doi.org/10.1175/BAMS-D-11-00007.1
Wang L, Cole JNS, Bartlett P, Verseghy D, Derksen C, Brown R, Salzen von K(2016) Investigating the spread in surface albedo for snow-covered forests in
CMIP5 models. J Geophys Res Atmos 121:1104–1119. doi: https://doi.org/10.
1002/2015JD023824
Watanabe M, Suzuki T, O’ishi R, Komuro Y, Watanabe S, Emori S, Takemura T,
Chikira M, Ogura T, Sekiguchi M, Takata K, Yamazaki D, Yokohata T, Nozawa T,
Hasumi H, Tatebe H, Kimoto M (2010) Improved climate simulation by
MIROC5: mean states, variability, and climate sensitivity. Journal of Climate
23:6312–6335. https://doi.org/10.1175/2010JCLI3679.1
Wood EF, Roundy JK, Troy TJ, van Beek LPH, Bierkens MFP, Blyth E, de Roo A, Döll
P, Ek M, Famiglietti J, Gochis D, van de Giesen N, Houser P, Jaffé PR, Kollet S,
Lehner B, Lettenmaier DP, Peters-Lidard C, Sivapalan M, Sheffield J, Wade A,
Whitehead P (2011) Hyperresolution global land surface modeling: meeting a
grand challenge for monitoring Earth’s terrestrial water. Water Resour Res 47:
54. https://doi.org/10.1029/2010WR010090
Yamazaki D, de Almeida GAM, Bates PD (2013) Improving computational
efficiency in global river models by implementing the local inertial flow
equation and a vector-based river network map. Water Resour Res 49:7221–
7235. https://doi.org/10.1002/wrcr.20552
Yamazaki D, Kanae S, Kim H, Oki T (2011) A physically based description offloodplain inundation dynamics in a global river routing model. Water
Resour Res 47:174. https://doi.org/10.1029/2010WR009726
Yoshimura K, Miyazaki S, Kanae S, Oki T (2006) Iso-MATSIRO, a land surface modelthat incorporates stable water isotopes. Global Planet Change 51:90–107
Yoshimura K, Nitta T, Ishitsuka Y, Tada M, Suzuki K, Takemura T (2018) Impact ofshort-lived climate pollutants on terrestrial water circulation. J Japan Soc Civil
Eng Ser B1 74:I_217–I_222. https://doi.org/10.2208/jscejhe.74.I_217

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.