MAESPA-Style Synthesis​#

Let's combine the earlier tutorials together to calculate CO₂ uptake and water loss by leaves, and their interactions with soil water availability. The example is inspired by MAESPA, which couples photosynthesis, transpiration and soil–plant water relations. Here, a small stand contains two species and five leaves. Leaf and canopy fluxes and soil water are calculated hourly. Within each hour, a controller repeatedly runs the leaf models to find a consistent canopy air temperature and humidity.

The original MAESPA model does not include carbon allocation. We deliberately add daily carbon allocation and a separate leaf area index (LAI) model to show how PlantSimEngine lets you extend a simulation by connecting additional models, each with its own time step. Allocation distributes the carbon assimilated by each plant among leaf, wood and reserve pools; the daily LAI calculation sums leaf areas per unit ground area. Leaf areas stay fixed in this example, so allocation does not change LAI.

This is an uncalibrated teaching example inspired by MAESPA's process structure. It is not a validated MAESPA implementation. Leaf illumination is uniform with complete absorption, the canopy has one air layer, and the soil water model uses prescribed withdrawal fractions and bounds rather than a complete soil hydraulic balance. The example tests coupling and carbon accounting; its results should not be used to predict a real stand's behaviour.

This is an advanced example; start with the individual tutorials if these steps are unfamiliar. The table in How the pieces compose links each part of this example to the tutorial that explains it.

Run the reference case​#

The complete, tested source is in examples/maespa_model_example.jl. Run it with 73 hourly weather records to show three daily cycles. The daily models run on the first step and again every 24 hours, giving four calls:

using PlantSimEngine, DataFrames

include(joinpath(
    pkgdir(PlantSimEngine),
    "examples",
    "maespa_model_example.jl",
))

result = run_maespa_example(; nhours=73, check=true)
simulation = result.simulation
model = result.model
nothing

The model contains one object for the whole scene, one soil object, and two plants with the same model connections. Each plant has its own template, species parameters, and number of leaves:

(
    instances=DataFrame(Diagnostics.explain_instances(model)),
    plants=length(model_objects(model; scale=:Plant)),
    leaves=length(model_objects(model; scale=:Leaf)),
    species_A=length(model_objects(model; scale=:Leaf, species=:A)),
    species_B=length(model_objects(model; scale=:Leaf, species=:B)),
)

Plot dynamics at three scales​#

Use the same simulation to compare scene transpiration, carbon allocated by each plant, and water content in the two soil layers. CairoMakie is already part of the documentation environment. Separate panels keep the different units readable, while a shared time axis shows the hourly and daily changes.

using CairoMakie

function history(object, variable)
    sort!(collect_outputs(simulation, object, variable; sink=DataFrame), :timestep)
end

transpiration = history(:model, :scene_transpiration)
plant_carbon = [history(id, :accounted_carbon) for id in (:plant_A, :plant_B)]
soil_water = [history(:soil, variable) for variable in (:theta1, :theta2)]

figure = Figure(size=(840, 780), fontsize=16)
scene_axis = Axis(figure[1, 1];
    title="Scene · transpiration",
    ylabel="Water loss\n(mm per hourly interval)",
    xticks=[1, 25, 49, 73],
)
plant_axis = Axis(figure[2, 1];
    title="Plants · cumulative allocated carbon",
    ylabel="Allocated carbon\n(g C per plant)",
    xticks=[1, 25, 49, 73],
)
soil_axis = Axis(figure[3, 1];
    title="Soil · water content",
    xlabel="Hourly simulation step",
    ylabel="Water content\n(m³ m⁻³)",
    xticks=[1, 25, 49, 73],
)

lines!(scene_axis, transpiration.timestep, Float64.(transpiration.value);
    color=:steelblue, linewidth=2.5)
for (series, label, color) in zip(plant_carbon, ("Plant A", "Plant B"),
                                (:seagreen, :darkorange))
    stairs!(plant_axis, series.timestep, Float64.(series.value);
        step=:post, color, label, linewidth=2.5)
    scatter!(plant_axis, series.timestep, Float64.(series.value);
        color, markersize=8)
end
for (series, label, color) in zip(soil_water, ("Upper layer", "Lower layer"),
                                (:sienna, :mediumpurple))
    lines!(soil_axis, series.timestep, Float64.(series.value);
        color, label, linewidth=2.5)
end
axislegend(plant_axis; position=:lt, framevisible=false)
axislegend(soil_axis; position=:rt, framevisible=false)
linkxaxes!(scene_axis, plant_axis, soil_axis)
hidexdecorations!(scene_axis; grid=false)
hidexdecorations!(plant_axis; grid=false)
xlims!(scene_axis, 0, 74)
figure

Transpiration follows the daily weather cycle. The plant curves change only when allocation runs, at steps 1, 25, 49 and 73; markers show those saved values, and the steps show the value held between calls. This is allocated elemental carbon, not biomass. The first allocation covers only the starting hour. Soil water decreases as transpiration draws on the two layers; this example supplies no rain or infiltration.

The horizontal positions come from the saved timestep column, including for the daily models. These are global simulation steps, not row numbers within each output series. The plot uses hourly steps because this example's environment provider does not attach calendar dates to the saved outputs.

Check when models run and where their inputs come from​#

Check when each calculation runs. The scene controller decides when to call leaf energy balance and soil water; these models do not run independently. Allocation and LAI run once per day. In the table, root_scheduled means a model runs directly from the simulation schedule, manual_call_only means another model calls it, and dt_steps gives the interval in hourly steps:

schedule = DataFrame(Diagnostics.explain_schedule(result.compiled))
select(
    filter(
        row -> row.application_id in (
            :scene_eb,
            :soil_water,
            :lai_dynamic,
            :plant_A__energy_balance,
            :plant_A__allocation,
            :plant_B__allocation,
        ),
        schedule,
    ),
    :application_id,
    :root_scheduled,
    :manual_call_only,
    :dt_steps,
)

The scene energy-balance model calls all five leaf models with Many and the single soil model with One. These are hard calls: the calling model controls when the other models run. The table lists the models and objects called from :scene_eb:

calls = DataFrame(Diagnostics.explain_calls(result.compiled))
select(
    filter(row -> row.application_id == :scene_eb, calls),
    :call,
    :callee_application_ids,
    :callee_object_ids,
    :publication_policy,
)

Each plant's allocation model reads carbon values from its own leaves. The scene model reads areas from all leaves and water potential from the soil. Here AllocA and AllocB use the same fixed-fraction allocation equation with different parameters. For two different allocation rules running on two plants together, see Instantiate Several Plants. Inspect source_ids below to check where each input comes from. The values are shared by reference, meaning the reader sees the current source values without copying them:

bindings = DataFrame(Diagnostics.explain_bindings(result.compiled))
select(
    filter(
        row -> (
            row.application_id in (
                :plant_A__allocation,
                :plant_B__allocation,
            ) && row.input == :leaf_carbon
        ) || (
            row.application_id == :scene_eb &&
            row.input in (:leaf_areas, :psi_soil)
        ),
        bindings,
    ),
    :application_id,
    :input,
    :source_ids,
    :carrier_kind,
    :copy_semantics,
)

Follow trial canopy air to its accepted state​#

The scene controller starts from the above-canopy weather, named :forcing. It tries different canopy air conditions and runs the leaf models for each trial with publish=false, so these intermediate results are not saved in the output history. Once the calculation converges, it commits the accepted air conditions with sink=:canopy and runs the leaves once more to publish their accepted results.

The leaf models then read the accepted conditions from :canopy. The table below shows which environmental variables each model reads or changes. Its handle column contains an internal identifier used to retrieve those conditions; you do not need to interpret that identifier:

environment_bindings = DataFrame(
    Diagnostics.explain_environment_bindings(result.environment),
)
select(
    filter(
        row -> row.application_id in (
            :scene_eb,
            :plant_A__energy_balance,
            :plant_B__energy_balance,
        ),
        environment_bindings,
    ),
    :application_id,
    :object_id,
    :handle,
    :required_inputs,
    :produced_outputs,
)

MaespaSingleLayerEnvironment supplies one set of air conditions for the whole canopy. It keeps the above-canopy weather separate from the canopy conditions that the controller changes. You could replace it with an environment that represents several layers or 3D cells, while keeping the same variables available to the process models. See the two-cell example in Understand Environments and the implementation guide in Environment Backend Extensions.

Check units and carbon accounting​#

Check the units and conversions as values pass between models:

Every leaf receives the same above-canopy irradiance here. There is no shading, scattering, or leaf-angle calculation; a radiation model would need to replace that assumption for a realistic stand.

Use final_state to read the latest values, whether or not you saved their history:

scene = final_state(simulation, :model)
soil = final_state(simulation, :soil)
plants = final_state(simulation, Many(scale=:Plant))

(
    lai=scene.lai,
    canopy_temperature_C=scene.canopy_tair,
    hourly_transpiration_mm=scene.scene_transpiration,
    soil_water_potential_MPa=soil.psi_soil,
    allocation_interval_g_C=Dict(
        id => state.daily_growth
        for (id, state) in plants
    ),
)

Allocation reads cumulative leaf C without resetting it. Each plant stores the cumulative amount it has already accounted for and allocates only the difference at the next daily call. This prevents the same carbon being allocated again on later days. Carbon gained or lost since the last allocation remains pending until the next one:

DataFrame([
    (
        plant=id,
        cumulative_net_C_g=sum(state.leaf_carbon),
        accounted_C_g=state.accounted_carbon,
        pools_C_g=state.leaf_pool + state.wood_pool + state.reserve_pool,
        pending_C_g=sum(state.leaf_carbon) - state.accounted_carbon,
    )
    for (id, state) in plants
])

At each allocation, the three pools sum to accounted_carbon. Adding pending C recovers cumulative net assimilation. These accounts can increase or decrease: negative net assimilation reduces them. The example does not model initial biomass, construction respiration, dry-matter conversion, or limits on withdrawing reserves, so the pools are not predictions of organ mass.

Count the saved values to check how often the models ran. Hourly scene and leaf variables have 73 samples; daily LAI and allocation variables have four:

output_summary = DataFrame(Diagnostics.explain_outputs(simulation))
select(
    filter(
        row -> (
            row.object_id == :model &&
            row.variable in (:scene_transpiration, :lai)
        ) || (
            row.object_id in (:plant_A, :plant_B) &&
            row.variable == :daily_growth
        ) || (
            row.object_id == :plant_A_leaf_1 &&
            row.variable == :λE
        ),
        output_summary,
    ),
    :application_id,
    :object_id,
    :variable,
    :nsamples,
)

How the pieces compose​#

Plant structure stays fixed during this 73-hour example. To add or remove organs with a growth model, follow Modify Plant Structure.

What the tests check​#

The automated tests check that the parts work together as intended:

• the stand contains two plants, with five leaves assigned to the right plant;

• each plant's allocation model reads only its own leaves;

• the scene controller calls every leaf and the shared soil model;

• hourly and daily output counts match their scheduled intervals;

• accepting new canopy conditions does not replace the above-canopy weather;

• leaf fluxes are finite and their totals agree with the scene results;

• hourly transpiration agrees with latent heat loss and the change in soil water storage;

• PAR energy is bounded by shortwave energy and converted to photon units;

• a longer run with 73 hourly samples covers three daily intervals plus the initial call, with no carbon added by rejected trial calculations;

• each plant allocates every C increment once, preserves cumulative leaf C, and conserves the sum of its leaf, wood, reserve, and pending C accounts.