Hydrodynamic Modeling of Stratification and Mixing in a Shallow, Tropical Floodplain Lake

Floodplain lakes are widespread and ecologically important throughout tropical river systems, however data are rare that describe how temporal variations in hydrological, meteorological and optical conditions moderate stratification and mixing in these shallow lakes. Using time series measurements of meteorology and water‐column temperatures from 17 several day campaigns spanning two hydrological years in a representative Amazon floodplain lake, we calculated surface energy fluxes and thermal stratification, and applied and evaluated a 3‐dimensional hydrodynamic model. The model successfully simulated diel cycles in thermal structure characterized by buoyancy frequency, depth of the actively mixing layer, and other terms associated with the surface energy budget. Diurnal heating with strong stratification and nocturnal mixing were common; despite considerable heat loss at night, the strong stratification during the day meant that mixing only infrequently extended to the bottom at night. Simulations indicated that the diurnal thermocline up and downwelled creating lake‐wide differences in near‐surface temperatures and mixing depths. Infrequent full mixing creates conditions conducive to anoxia in these shallow lakes given their warm temperatures.

Riverine floodplains constitute important aquatic ecosystems in the tropics (Jardine et al., 2015;Junk, 1997), and seasonal variations in inundated areas and water levels are often large (Hess et al., 2015;Paiva et al., 2013).Shallow, floodplain lakes in the lowland Amazon basin are abundant (Sippel et al., 1992) and ecologically and biogeochemically significant (Melack et al., 2009), and similar floodplain lakes occur elsewhere in the tropics.Characteristic features of these lakes are strong cycles of diel stratification and mixing (Augusto-Silva et al., 2019), with the depth of vertical mixing and persistence of the thermocline related to water depths (MacIntyre & Melack, 1988).Seasonal changes in connectivity to rivers, resulting in changes in water levels that can be over 10 m, and changes in optical properties contribute to altered thermal stratification through the year.
In general, studies that consider seasonal changes in meteorology, latent and sensible heat fluxes and near-surface stratification in relation to the efficacy of mixing due to wind and cooling are rare.Similarly, studies have seldom incorporated the three-dimensionality of the lakes and spatial variations in mixing.Such knowledge is important for assessing the extent of deoxygenation, the fluxes of climate forcing trace gases, algal growth and other ecological conditions.
Our overall objectives are as follows: (a) application and evaluation of a 3D hydrodynamic model for conditions in shallow, tropical lakes based on measurements and computations of key processes over diel and seasonal periods in an Amazon floodplain lake; and (b) analysis and simulation of diel and seasonal changes in thermal structure and 3-dimensional patterns of stratification and mixing in relation to environmental conditions.Time series measurements of meteorology, water-column temperatures and light attenuation spanning 2 years and a wide range of water depths were made.These measurements are used to calculate latent and sensible heat fluxes, shear and convective velocity scales, and buoyancy frequencies (Imberger, 1985).The 3-Dimensional coupled Hydrodynamic-Aquatic Ecosystem Model (AEM3D) (Hodges & Dallimore, 2019) was used to model stratification and mixing.We compared the simulated and measured temperature and density structure over diel periods and other aspects of stratification and mixing.

Study Site
Lake Janauacá (3.38°S, 60.3°W) is a floodplain lake located on the southern side of the Solimões River in the central Amazon basin (Brazil) (Figure 1).The study focused on the northern open lake.Between October 2014 and November 2016, water level in Lake Janauacá fluctuated up to 13 m (Figure 2).Based on the water level change, a hydrological year can be divided into four phases: low water during November and December, rising water from January to May, high water during June and July and falling water from August to October.The range of hydrological, optical, chemical and ecological conditions observed in Lake Janauacá are typical of Amazon floodplain lakes (Junk, 1997;Melack et al., 2021).

Field Measurements
Time-series measurements of water temperature were obtained with thermistors deployed on vertical arrays during field campaigns between 2014 and 2016 (Figure 2).Thermistor arrays had the upper-most loggers suspended below and shaded by a surface float at ∼0.05 m, and deeper loggers were on taut-line moorings with a float ∼0.3 m below the water surface.The thermistors were RBR Solo® thermistors (accuracy 0.002°C), sampled at 0.1 Hz and spaced about 0.25 m in the top meter and about 0.5-1 m below 2 m (depths shown in figures).The temperatures were averaged over 5 min.Wind speed and direction sensors (Onset S-WCA-M003: accuracy ±0.5 m s −1 , resolution 0.2 m s −1 ; direction ±20, given motion of buoy) were deployed 2 m above the water surface on a buoy at the same location as the thermistor arrays.Conductivity was measured using a profiler (resolution 1 µS cm −1 ; Castway, Sontek Inst.Co), sampling at 4 Hz, with data reported at 0.3 m intervals.Vertical profiles of photosynthetically active radiation (400-700 nm, PAR) were obtained using a Licor LI-192 SB underwater sensor, and diffuse attenuation coefficients for those wavelengths (k d ) computed following Beers Law.
At the nearby floating research station (Figure 1a), measurements of downwelling PAR (Onset S-LIA-M003), air temperature and relative humidity at 2 m above water surface (Onset S-THB-M002, precision: ±2.5%, ±0.2°C), rainfall (Onset S-RGB-M002, resolution 0.2 mm) and wind speed and wind direction at 2 and 4 m above water surface (same sensors as used on buoy) were measured.Incoming longwave radiation (Kipp & Zonen CGR 3), and incoming shortwave radiation (Kipp & Zonen CMP 3) were measured for a period in August 2016 near the floating research station.Water levels in Lake Janauacá were manually read daily from a series of leveled stage gauges on the shore (3.42°S, 60.26°W) (Figure 1a) from October 2014 to October 2016.

Data Analysis
To understand controls on the temporal and spatial variability of thermal structure, we calculated surface energy budgets from meteorological data and near-surface temperatures following MacIntyre et al. (2002MacIntyre et al. ( , 2014) ) and calculated effective heat fluxes and buoyancy fluxes into the actively mixing layer (AML); Imberger (1985) provides the conceptual and empirical basis for the AML.Depth of the AML (z AML ) was calculated as the first depth with temperature 0.04°C lower than the near-surface temperature.The criterion was based on visual inspection that z AML captures shoaling during heating and deepening with increased wind speed or cooling.Variables associated with the surface energy budget and dynamics of stratification and mixing are listed in Table 1.Bulk transfer coefficients for momentum and heat were computed taking into account atmospheric stability (Hicks, 1975;Imberger & Patterson, 1990), followed by the computation of sensible and latent heat fluxes and water friction velocity, u *w .Incoming shortwave radiation was computed from PAR multiplied by 2.15, based on comparison of measurements of incident PAR and shortwave radiation at the site.Net longwave radiation (LWnet) was determined based on mean shortwave radiation and daily averaged LWnet during August 2016.On the cloudiest day, LWnet had a minimum of −50 W m −2 whereas on other days it averaged −80 W m −2 .Using the time series of rainfall and shortwave radiation, we identified the appropriate daily value at other times.Net shortwave radiation in the AML was estimated following three steps: (a) reflected shortwave radiation at the surface was computed with albedo calculated following Fresnel's Law (Neumann & Pierson, 1966); (b) penetration of shortwave radiation in the AML was computed following Beer's Law for seven-wavelength bands suggested by Jellison and Melack (1993) and k d (Figure 2), with angle of refraction computed following Snell's Law; (c) effective heat flux into the AML was computed as the sum of sensible heat flux, latent heat flux, net longwave radiation and net shortwave radiation retained in the AML (Imberger, 1985).Density was computed following Chen and Millero (1986) with salinity scaled from specific conductance based on analyses of major solute concentrations and specific conductances (Lesack, 1988).

Model Description
Numerical simulations were done with the hydrodynamic module of the Aquatic Ecosystem Model 3D (AEM3D), a coupled 3-D hydrodynamic-ecosystem model (Hodges & Dallimore, 2019).The hydrodynamic module uses a z-coordinate Cartesian grid allowing non-uniform spacing in three directions and solves the unsteady Reynolds-averaged Navier-Stokes equations for momentum under the Boussinesq approximation, the continuity equations for incompressible fluids, and the transport equation for scalars (temperature, salinity and passive tracer).The numerical scheme is a modification of TRIM-3D (Casulli & Cheng, 1992), a semi-implicit finite difference method solving 3-D shallow-water momentum equations for stratified flows.Temperature is modeled with the ULTIMATE and the QUICKEST methods for scalar transport (Leonard, 1991).Turbulence in the horizontal is characterized indirectly in a way that turbulent shear stress is linearly proportional to the gradients of   1∕3 This term is also known as the turbulent velocity scale due to heat loss.
F q m 3 s −3  Mechanical energy flux for mixing, computed as  = 0.5 Table 1 Derived Variables (From Imberger, 1985) Based on Field Data mean velocities in which the proportionality factor is the eddy viscosity (Hodges & Dallimore, 2019), assuming zero off-diagonal eddy viscosity and zero molecular viscosity.Turbulence in the vertical is characterized by the mixed-layer model, applied to the whole water column, originally developed by Imberger and Patterson (1981) and Spigel et al. (1986) in which the turbulent kinetic energy for mixing is induced by wind and convection at the surface and bottom friction contributes to inducing the shear that drives the mixing from the surface to bottom.

Model Setup for AEM3D
For the simulations in this study, AEM3D was not calibrated, and the model parameters were assigned values derived from laboratory and field experiments and theory.To evaluate AEM3D's capability to simulate key hydrodynamic features, we conducted 17 simulations during the periods with measured water temperatures between November 2014 and September 2016 (Figure 2).We evaluated the model based on comparisons of the measured and simulated results.
AEM3D was configured with a domain of 100 m × 100 m mesh grids in the horizontal and 0.1 m layers in the vertical; the time step was 30 s.The horizontal resolution was the same as the digital elevation model available, and the vertical resolution was higher than the measurements to ensure proper characterization of the near-surface thermal structure.The time step was chosen through trials in which the Courant-Friedrichs-Lewy condition and the stable simulation of turbulent kinetic energy dissipation were met.
The bathymetry was calculated using a combination of field measurements with a GPS-equipped echo sounder and acoustic Doppler current profiler plus a vegetation-corrected digital elevation model (Pinel et al., 2015).The southeast boundary is considered an open boundary with changes in water level applied but no exchanges of water or constituents.The bathymetry is shown with depths on 22 June 2015 when water level was the highest during the period of study (Figure 1b).
The model parameters were configured with non-site-specific literature values as suggested by Hodges and Dallimore (2019) (Table S1 in Supporting Information S1) with two exceptions: (a) the diffuse attenuation coefficient varied among simulations based on the measured PAR profiles, and the attenuation coefficients for five wavelength bands of near-infrared radiation were adopted from Jellison and Melack (1993); (b) the two coefficients, used to quantify the extent of mixing driven by potential energy, were evaluated and selected based on comparisons of measurements and modeling.The coefficients consider the efficiency of mixing and its rapidity, and are called C C and C CT , respectively.We evaluated C CT of 0.5 or 1 and C C varying from 0.2 to 0.8 in 4 steps (see Text S2 and Figure S1 in Supporting Information S1).Results improved with a value C CT of 0.5 and with a value C C of 0.6.Comparisons of simulations with and without Coriolis force had similar thermal structure; hence, Coriolis force was not included.
Initial values of the bulk transfer coefficients for heat and momentum (Table S1 in Supporting Information S1) are those for neutral atmospheric stability.The bulk transfer coefficients are updated at each time step of the simulation following the algorithms in Imberger and Patterson (1990), taking into account the effects of non-neutral atmospheric stability.
The AEM3D validation simulations were forced with wind speed and wind direction data measured on platforms floating near thermistor arrays.Other meteorological data were obtained at the floating research station, as described above; during campaign 17 in late August 2016, the 2-m meteorological station was deployed on a buoy at the measurement site.

Evaluation of Model Performance
Hipsey et al. (2020) identified four aspects of model performance: level 0-conceptual validation; level 1-state validation, comparison of simulated variables with observations; level 2-process validation, comparisons with measured and calculated processes, and level 3-system validation, assessment of system-level emergent properties.To evaluate AEM3D performance we considered levels 0, 1, and 2. Conceptual validation focused on the value of incorporating 3-dimensional, sub-daily features to represent the system.Statistical metrics and graphical comparisons are used for level 1 evaluations of simulated temperatures and terms derived from temperature.Level 2 performance is illustrated by graphical comparisons and use of statistical metrics for variables derived from the surface energy budget and dynamics of stratification and mixing (Table 1).
Root mean square error (RMSE) and percent relative error (PRE) were computed to examine the absolute and relative deviation of the simulated results to measurements of temperature and terms derived from temperature and meteorological data.Correlation coefficients (R) were computed to examine how well the variation of simulated results match with the measurements.The three statistics (Bennett et al., 2013;Bruce et al., 2018) were computed as follows: Root mean square error (RMSE) and   are the "ith" simulated and observed data,   and   the mean simulated and observed data, and   the total number of data points.
The model evaluation was done for both state validation and process validation as categorized by Hipsey et al. (2020).For the state validation, two groups of the three statistics were computed with measured and simulated temperatures, one with the time series of temperatures at the depths of thermistor sensors (T(z)) and the other with each vertical temperature profile at each timestep of the measured data (T(t)).To ensure the same depth-time coordinates of the simulated data and measured data, the simulated temperatures were averaged to 5-min bins and resampled at the depths nearest to the depths of measurements.
For the process validation, the three statistics were computed with times series of z AML and H * (Table 1) derived from meteorological data and measured and simulated temperature data.We chose these variables for the process validation because they quantify the energy fluxes and the mixing dynamics in the water column which are the main interests of the study.

Water Depths, Underwater Light Attenuation and Meteorology
As a result of changing stage in the Solimões River and local inputs, water depths in L. Janauacá varied from 1.6 to 14.7 m during the period from October 2014 to November 2016 (Figure 2).Surface area in the modeled domain varied from 0.1 to 51.4 km 2 , and during lowest water the lake was restricted to a channel.Surface area increased from 2% to 68% of the maximum area when the maximum depth increased from 2 to 4 m.The rate of change per day (ΔZ d , cm day −1 ) varied among the four hydrological phases and included periods of rising and falling water in all phases: During low water (November and December 2014 and 2015), ΔZ d varied from −20.1 to 12.1 cm day −1 ; during rising water (January-May 2015 and 2016), ΔZ d varied from −5.1 to 23.2 cm day −1 ; during high water (June and July 2015 and 2016), ΔZ d varied from −7.1 to 3.0 cm day −1 ; during falling water (August-October 2015 and 2016), ΔZ d varied from −26.3 to 8.1 cm day −1 .Water depth during falling water changed nearly twice as fast as rising water.
The diffuse attenuation coefficient varied with water depth (Figure 2).Highest values occurred when water depth was shallowest (e.g., 8.2 m −1 in January 2016), and values declined with rising water, such that k d ranged from 4 to 1 m −1 over the remainder of the year.
Meteorological data during the 17 field campaigns illustrate both diel and seasonal patterns (Figures S2 and S3 in Supporting Information S1).The local rainy season occurred from January to June, and less rain usually occurred from July to December (called the dry season).The frequency and intensity of rainfall events was higher during rising water than other months (Table S2 in Supporting Information S1).Rainfall varied from tenths of mm hr −1 with occasional events of 15 mm hr −1 during campaigns.The rainy season had lower insolation due to considerable cloud cover.Air temperatures varied by ∼8 o C over the course of a day with air temperature exceeding surface water temperature in the afternoon.The dry season had increased solar insolation, relative to the rainy season.Air temperatures had a similar range as in the rainy season.During the dry season of 2016, air temperatures remained cooler than surface water temperatures throughout the day.Afternoon air temperatures were highest at low water.Relative humidity was highest at night and decreased to 50 or 70% in the day with the highest variability at low water and a tendency for reduced variability in the rainy season.
Wind velocity varied over the course of a day.Afternoon winds were highest during the dry season, with maxima often exceeding 4 m s −1 .Highest winds occurred during local squalls which occurred in all seasons.Variations of wind speed were accompanied by changes in wind direction, particularly when shifting from the calm winds in morning to higher winds beginning around mid-day.Wind directions showed seasonal variations, with a tendency for more southerly winds during rising water and northerly winds during falling water.

Thermal Structure, Heat Budgets and Validation of Simulations
We illustrate changes in thermal structure and buoyancy frequency over several diel cycles and compare simulations with measurements, by focusing on campaign 17 (22-27 August 2016), as it has the longest data set and occurred at moderate water levels (Figure 3), and with discussion of simulations from low and rising water (Figures 4 and 5).Results for all simulations are summarized in Figures S4 and S5 in Supporting Information S1 with descriptions in Text S3 in Supporting Information S1; further comparisons of measurements and simulations are illustrated for the four hydrological phases (Text S4; Figures S8-S15 in Supporting Information S1).Results of calculations of surface energy budget and of within-lake mixing dynamics from field data and simulations are contrasted in Figure 6 for campaign 17 and for all the simulations in Figures S6 and S7 in Supporting Information S1.Solar radiation during the day combined with nocturnal cooling created conditions favorable for the development of strong stratification in the day and mixing at night (Figures 3-5).Daily maxima of incoming shortwave radiation approached 1,000 W m −2 .Heating within the AML, including the diurnal thermocline, occurred when H * was positive, from shortly after sunrise to late afternoon.Heat losses at night varied with wind speed with typical losses ranging from −100 to −400 W m −2 .Air temperatures were less than surface water temperatures such that the stratification in the atmosphere would be unstable except briefly during late afternoon in late August 2016 and during rising water in January 2015.Winds were low in early morning and on some days rose around mid-day to 4 or 5 m s −1 .During squalls, as captured in August 2016, winds briefly reached 12 m s −1 .The frequent shifts in wind direction were independent of time of day.
Diurnal thermoclines formed each morning in the upper meter (Figures 3d, 4d, and 5d).Their thickness and temperature varied with the magnitude of wind near sunrise and over the next few hours and with cloud cover.They began to deepen in the afternoon with slight changes in wind velocity and sometimes despite continued heating.Once cooling began, the diurnal thermocline deepened further.During low water, mixing extended to the bottom, but as the lake deepened to 5 m and when winds were light at night, mixing was incomplete (Figure 5).Observations during low water and in August 2016 indicate mixing to the bottom on some nights but not others.By the following morning, water cooler than measured the previous night occurred near the bottom.The arrival of cooler water could be a result of density driven flows from greater cooling in shallow regions.However, the up and down motions associated with the cool water over the course of the day occur with changes in wind velocity implying they are the result of the diurnal thermocline tilting along the axis of the lake.Such tilting at night followed with a relaxation in the morning would bring cooler water to the measurement site.This inference is supported by the simulations presented in Section 3.3 on spatial variability.Warmer water also flowed to the measurement site in the upper water column after periods of cooling.Such was evident from mid to late afternoon on 30 November 2011 (Figure 4c), the nights of January 29 and 31 (Figure 5c) and early morning 23 August 2016 (Figure 3c).These advective flows also occurred with changes in wind speed or direction with the first a response to an initial upwelling of cool water and its subsequent downwelling as wind speeds declined.
Buoyancy frequencies (N) in the diurnal thermocline during morning heating were high with the magnitudes dependent on solar heating, wind speed, and k d (Figures 3e, 4e, and 5e) Thus, the values during falling water when The simulated temperatures captured the predominant features, including the high near-surface temperatures and buoyancy frequencies during morning heating and downwelling of the diurnal thermocline (Figures 3f,3g,4f,4g,5f,and 5g).While the observational data showed the diurnal thermocline to begin to downwell in response to changes in wind speed and direction in the afternoon, the response was often delayed by an hour or two in the simulations.This discrepancy sometimes led to thickening of the downwelling diurnal thermocline, as during rising water and falling water (Figures 3d-3g and 5d-5g).The simulations also captured the arrival of cool water near the bottom in early morning (Figure 3d), the intrusion of warm water in the upper water column at low water (late afternoon and extending into the night on 30 November 2014, Figure 4d), and the intrusion of warm water during rising water (early morning, 31 January 2015, Figure 5d).Simulations did not capture the warm water intrusion late on 29 January 2015 (Figure 5d).While simulated temperatures were similar to measured ones throughout the water column at low water, subtle differences occurred at other times.For example, the model underestimated near-surface temperatures and temperatures in the lower water column by ∼0.5°C during falling water (Figure 3).Despite these differences, the model correctly simulated the high values of N near the surface and below.During falling water, stratifica tion more often persisted mid-water column in the simulations, especially when wind speeds were low.Overall, simulations captured the patterns of diel stratification and mixing with occasional underestimates of the depth of nocturnal mixing during rising and falling water.
The comparisons of measurements and simulations for the other campaigns have a strong correspondence with respect to near-surface temperatures and depths of the diurnal thermocline, timing of descent of the diurnal thermocline, extent of cooling, near bottom upwelling at night or in the early morning, and magnitude of the buoyancy frequency (Text S3 and S4; Figures S4, S5, S8, S10, S12, and S14 in Supporting Information S1).The largest discrepancies are at higher water.Further evaluations of model performance are provided in Sections 3.3 and 3.4.

The Surface Energy Budget and Metrics of Mixing
Accuracy in modeling depends on appropriately quantifying the surface energy budget.Comparisons of key terms for simulation 17 and for all the campaigns are presented in Figure 6 and Figures S6 and S7 in Supporting Depths of the AML, computed using thermistor data and simulated temperatures, were similar during morning heating, but diverged with changes in wind velocity, typically being delayed by an hour, indicating the model was simulating a slower downwelling of the diurnal thermocline in response to changes in wind speed or direction.When simulations generate a shallower z AML at night, calculated values of w * are smaller than when based on observations, but the differences in the calculated F q are minor given the shallow z AML .Results of the comparisons for the other simulations also indicate minor discrepancies in the calculated variables (Figures S6 and S7 in Supporting Information S1).
While the turbulent kinetic energy for mixing can be modeled based on u * w and w * (Imberger, 1985), mixing by convection in AEM3D depends on the actual heat loss at the air-water interface and the resulting near-surface density instability.As illustrated in Figure 6 and in Figures S6 and S7 in Supporting Information S1, the calculated heat losses are similar for observations and simulations.Results of the surface energy budget calculations for the other 16 field campaigns and AEM3D simulations are shown in Figures S6 and S7 in Supporting Information S1.In these other simulations, z AML was essentially the same for the measurements and simulations, so the variables which depend on it had only minor differences between those calculated and observed (Figures S6c and S7c in Supporting Information S1).

Metrics of Model Performance
The statistics RMSE, PRE and R were computed with the time series of temperatures at each sampling depth (T(z)) and the profiles at each sampling timestep (T(t)) for the state validation (level 1), and with the time series of H * and z AML for process validation (level 2) (terms as used in Hipsey et al., 2020) for each simulated period.
The RMSE of T(z) (Figure 7a 1 ) shows that the absolute deviations in temperature decreased from the surface downwards.Deviations tended to be less than 1°C near surface and close to zero near the bottom, with 95% of the RMSE values below 1.0°C and less than 0.5°C below 2 m.The PRE of T(z) (Figure 7a 2 ) ranged from −2% to 4% near surface and decreased with depth such that it was close to zero near the bottom.The correlation between measured and simulated temperature, as shown by R (Figure 7a 3 ), was in the range from 0.5 to 1 near the surface with values closer to 1 when water depths were least.Deeper in the water column, values ranged from −0.5 to 0.5; the weaker correlation at deeper depths was attributed to the small temperature variability in deeper water and indicate that R values are misleading in this situation (Figures 7a 1 and 7a 2 ).AEM3D simulates well the temporal variation of temperature at the sampling depths.
The absolute and relative deviations in temperature profiles at the sampling times were similar to the deviations in time series of temperature at the sampling depths.RMSE of T(t) was mostly between 0 and 1.5°C (Figure 7b 1 ) and PRE of T(t) was between −2% and 4% (Figure 7b 2 ).R of T(t) was generally above 0.75 when the water depth was greater than 4 m and intermittently lower than 0.5 when the water column was shallower than 4 m; temperature variability tended to be low at the same time.Overall, AEM3D was able to simulate the vertical variation of the temperature in the water column at the sampling sites.Key hydrodynamic processes are depths of the AML and energy transfer in the near-surface water.RMSE of  AML (Figure 7c 1 ) during the day varied between 0.3 and 2 m.PRE of  AML (Figure 7c 2 ) during the day was usually less than 30%; simulation 1 was 54% and simulation 5 was 81%.R of  AML (Figure 7c 3 ) was usually greater than 0.7.Simulation 13 is anomalous because the water was only about 1 m deep.
As shown by the RMSE of H * (Figure 7d 1 ), AEM3D simulated net energy fluxes in the AML during the day deviated from those measured by 45-115 W m 2 with no clear pattern or correlation with water depth.These discrepancies occur during heating and result from differences in the calculated values of  AML from the observations and simulated data.During heating, H * varies from 0 to 600 W m −2 (Figures 3-5), hence the magnitude of the errors is not large.As evidence for that interpretation, PRE of H * (Figure 7d 2 ) are primarily between −20% and 20%, and the variation of simulated H * was well correlated to the measured because R of H * was ∼1.At night, z AML does not influence H * .RMSE of z AML (Figure 7e 1 ) during the night varied between 0.2 and 2.6 m.PRE of z AML (Figure 7e 2 ) during the night was usually less than 20%; simulations 3, 6, and 16 are outliers.R of z AML (Figure 7e 3 ) was usually greater than 0.5.The variability is indicative of observed and modeled differences in depth of nocturnal mixing.

Limitations of Modeling
AEM3D model achieved satisfying model performance by simulating key hydrodynamic features including diel cycles of stratification and mixing, without calibration and by using model parameters and coefficients derived from hydrodynamic theory and field experiments.However, improving the estimated depth of convective mixing required evaluation of parameter coefficients.Three coefficients are used to evaluate mixing during convection in AEM3D: C C is the fraction of available potential energy being consumed for mixing, C CT is the timescale for mixing due to convective overturn, and C TT is the timescale for mixing due to turbulence (for further explanation, see Hodges & Dallimore, 2019).Coefficients C C and C CT were varied and values selected that matched patterns in the August 2016 data; varying C TT did not affect the results (Text S2; Figure S1 in Supporting Information S1).The depth of mixing at night compared well with the field data except at the onset of falling water in 2015 and2016 (Figures S4 andS5 in Supporting Information S1).While the model correctly captured the formation of the diurnal thermocline and the extent of heating, differences between measurements and simulations during cooling at higher water levels may be a result of deficiencies in the input data, including longwave radiation and wind speeds and directions (Text S1 in Supporting Information S1).Measurements of wind speeds and directions from several locations would likely improve the simulations.
Rainfall can cause near surface cooling and moderate the depth of deepening.The extent of cooling from rainfall depends on the temperature of the rain relative to near-surface water temperature.For example, during a rain event in August 2016 heat loss caused by the rain was estimated as −1685 W m −2 while maximum heat loss, as the sum of SE, LE, and LWnet was −900 W m −2 with an average for the first hour of −550 W m −2 (Text S5 in Supporting Information S1).Temperature measurements indicated more rapid deepening than the AEM3D simulations (Figure 3).Hence, adding an algorithm for the temperature of the rain may lead to enhanced cooling and improved estimates of the rate of deepening of the diurnal thermocline and mixing during rain events.
Whether the effects of the Earth's rotation are necessary to include seemed unlikely because of the equatorial location of L. Janauacá.However, the Rossby number was estimated to vary around the critical value 1 with simulated horizontal velocities indicating that it was possible that the Earth's rotation may influence simulations.A simulation with Coriolis force included was done for August 2016 (simulation 17), and indicated that including the Coriolis force led to insignificant changes in thermal structure.

Spatial Variability in Thermal Structure
The simulations illustrate variability in thermal structure which cannot be visualized with profiles at one location.The variability results from differences in sheltering from wind, wind direction, and water depth.For example, warmest temperatures were simulated at Col. 4 in the sheltered arm (Figures 1 and 8).The shallow diurnal thermocline at Col. 1 relative to Col. 3 during heating on August 23 and 24 is indicative of the northerly transport of water when winds were southerly.Similarly, the warmer water found to a depth of 6 m at Col. 1 simulated from early afternoon August 25 to early morning August 26 is indicative of downwelling from the southward transport of water resulting from northeasterly winds (Figures 1b and 8a 1 -8e 1 ).While the strength of stratification, is the highest during heating with light winds, its magnitude in mid-water column varies with the up and downwelling induced by changes in wind velocity (Figures 8a 2 -8e 2 ).For example, with the southerly winds on August 23 and through mid-day on August 24, N was ∼20 cph in the diurnal thermocline as it downwelled from ∼0.5 to ∼5 m at northerly Col. 3.
The modeling also provides spatial and cross-sectional views of the dynamics of the thermal structure (Figures 9 and 10, Movie S1).For example, with westerly winds mid-afternoon on 25 Aug, warmed water flowed to the  eastern side of the basin near the surface (Figure 9b) causing downwelling evident at 2 m (Figure 9c) and at the southern ends of the curtains showing cross sections of thermal structure (Figures 9d, 9e, and 9g).In contrast, the diurnal thermocline remained level along curtain 4 with its north to south orientation.
The modeling also shows that prior differences in heating and wind velocity influenced subsequent conditions.For example, winds were southerly to southeasterly and fluctuated from near 0-1 m s −1 for ∼2 hr prior to the simulation shown for 1,344 hr on 27 August 2016 (Figure 10) As a consequence of the light winds and net heating, temperatures were near uniform over the surface of the lake with the heating constrained to a thin layer near the surface.Winds had been south-southwesterly with a magnitude of 4 m s −1 from sunrise until ∼0930 hr.The downwelling to the north and east, visible as warmer water to 3-6 m over the first several km of each curtain (Figures 10d-10g), is a result of the previous northeasterly transport of warmed water and some continued transport even with the later light winds.

Seasonal Variations in Temperature and Thermal Structure
The modest seasonal variability in the surface energy budget combined with the large changes in water depth contribute to determining whether the lake mixes fully or whether it is seasonally stratified.With the diel cycles of heating and cooling, the lake stratified on most of the days.At low water, near isothermy was often found at night, but even when the water column was deeper, water cooler than the previous night sometimes was found near the bottom in the early morning.As the water depth increased, mixing only intermittently extended to the bottom (Figures S4 and S5 in Supporting Information S1).Temperatures throughout the water column were coolest during rising water.Near-surface water temperatures were highest during the dry season with its reduced cloud cover and at low water.Buoyancy frequencies also varied seasonally with this difference linked to the higher k d at lower water levels.Diurnal thermoclines in the upper 0.5-2 m had buoyancy frequencies reaching 100 cph in April 2015 and during low water in 2016 (Figures S4d and S5d in Supporting Information S1).Near-surface values of N decreased after mid-day with downward movement of the diurnal thermocline or the onset of cooling.Downward movement of the diurnal thermocline rarely extended below 6 m; hence, when water depths exceeded 6 m, stratification below 6 m was weaker than in the overlying water (<10 cph).Hence, stratification near the surface and below 6 m varied seasonally.

Discussion
Tropical floodplain lakes differ from other types of lakes given the large seasonal variations in water depths.That said, water depths are shallow so the diel cycles of stratification and mixing observed in other tropical lakes are prevalent.Hence, the lakes are classified as polymictic with stratification due to diurnal, as opposed to seasonal, thermoclines.Modeling such lakes thus requires quantifying changes in the upper mixed layer with its AML, diurnal thermocline, and subsurface layer (Imberger, 1985).Few modeling studies have focused on these dynamics on diel time scales.Here, we have demonstrated that a three-dimensional hydrodynamic model can accurately simulate the diel changes in temperature and stratification in a representative floodplain lake.The model was validated using observational data obtained on 17 field campaigns over two hydrological cycles as water level rose and fell by ∼15 m.Additional validation followed guidelines in Hipsey et al. (2020) and included assessing accuracy in the model's simulating the depth of the AML and heat fluxes within it.Several striking features emerged that extended our understanding of tropical floodplain lakes, including extensive up and downwelling of the diurnal thermocline despite weak temperature stratification and weak winds.In the following, we provide further perspectives on these features and on future studies that can be conducted given the successful modeling demonstrated here.
Three-dimensionality of thermal structure in shallow, tropical floodplain lakes: The simulations indicated the dynamic nature of shallow floodplain lakes.Rather than the dominant feature being one-dimensional heating with subsequent vertical mixing during cooling, the modeling showed the formation of shallow diurnal thermoclines in the morning with light winds.Only subtle changes in wind speed and direction were required for them to tilt in the direction of prevailing winds, further shallowing the diurnal thermocline upwind and causing the AML to deepen downwind with an intensification of stratification in the diurnal thermocline at its base (Figures 8-10 and Movie S1).These motions lead to horizontal variability of temperature at each depth and along cross-sections.The simulations provide new understandings of physical processes in polymictic lakes with numerous ecological consequences.
While prior observations have suggested that shallow tropical lakes often mix to the bottom at night (Talling & Lemoalle, 1998), observational studies in tropical lakes using high-resolution thermistor arrays have found incomplete mixing in shallow tropical lakes.Using thermistor arrays separated by several kilometers in a floodplain lake Augusto-Silva et al. ( 2019) found that internal wave motions led to up and downwelling similar to our observations.MacIntyre et al. (2002) showed the importance of density driven flows from morphologically induced horizontal differences in rates of heating and cooling.The simulations in L. Janauacá indicated that internal wave motions led to the arrival of cool bottom waters in early morning and intrusions of warmer water in the afternoon.These motions led to either persistence or intensification of stratification that restricted complete mixing at night.
Our simulations enable discrimination of several processes retaining stratification.At locations where the diurnal thermocline upwelled in the day, the relaxation of the wind or a change in its direction would allow warmer water to flow back over cooler water at night and intensify stratification (Movie S1).The up and downwelling also led to spatial variability of stratification with N at inshore sites of order 10 to 20 cph (N 2 = 0.0003 and 0.0012 s −2 , respectively) and ∼60 cph (N 2 = 0.011 s −2 ) offshore where downwelling intensified the stratification at mid-water column during the day (Figure 8).The rate of deepening depends on β * , the effective buoyancy flux into AML (Table 1), and N.That is, dh/dt = β * /(hN 2 ) where h is the depth of the AML (Pernica et al., 2014).Heat losses at night, −100 W m −2 to −450 W m −2 , were appreciable with corresponding effective buoyancy fluxes ranging from −7 × 10 −8 to 3 × 10 −7 m 2 s −3 .Inshore for the two values of N, letting β * = 2 × 10 −7 m 2 s −3 , and starting with a mixing layer depth of 2 m, dh/dt over 12 hours would be 7 and 3.6 m.Thus, mixing could extend to the bottom frequently inshore.Offshore, assuming h = 2 m and similar β * , the diurnal thermocline would only descend an additional 0.4 m at water column 5. Thus, while stratification could erode quickly at inshore locations (e.g., water column 1) where stratification was also weakened by up and downwelling of the diurnal thermocline (Figure 8), it would persist over the 12 hr of darkness further offshore (e.g., water column 5).With the resulting horizontal differences in temperature due to differences in the depth of mixing at night, density driven flows develop to restore stratification in the vertical, that is, horizontal advection occurs.The simulations extend our understanding that horizontal flows contribute to persistent vertical stratification, despite weak temperature differences and despite considerable nocturnal cooling in shallow, tropical lakes.
The simulations show that with weak winds under heating, the diurnal thermocline extends across the lake.Even when winds begin to increase, heat remains trapped in a thin layer.Strong stratification near the surface can lead to enhanced gas transfer velocities and more rapid flux of climate forcing trace gases across the air-water interface (MacIntyre et al., 2021).Additionally, the currents that cause the downwelling transport solutes.
Brief convective rain storms with rainfall rates from a few mm hr −1 to 30 mm hr −1 occur in both the rainy and dry seasons (data not shown).As discussed above and in Text S5 in Supporting Information S1, cooling during these events is larger than computed with the surface energy budget and has the potential to cause mixing to the bottom.Further effort on the influence of rainfall on mixing in tropical floodplain lakes is warranted.
Seasonality: The successful modeling of L. Janauacá allows extending understanding of processes moderating seasonality of mixing dynamics in tropical lakes and warm-water lakes elsewhere.L. Januaucá's warm temperatures, and meteorological and optical conditions are shared with many other tropical lakes and lakes in temperate zones during summer.Due to muted seasonal variations in solar radiation and air temperature in tropical regions, seasonal variations in thermal structure in tropical lakes have been attributed to variations in wind speed, latent heat fluxes, and cloud cover (Talling & Lemoalle, 1998).In eastern African lakes movement of the Intertropical Convergence Zone (ITCZ) leads to alternation of rainy and dry seasons with increases in wind speed in the dry season (MacIntyre et al., 2014).In the Amazon basin, northward ITCZ movement coincides with rising water (Paiva et al., 2013), increased cloud cover and therefore more intermittent heating and cooling.Reduced cloud cover, warmer air temperatures, and increased winds occur around high water with clearer skies and warmer air temperatures during falling water.These seasonal differences combined with changing depths lead to frequent mixing to the bottom at lowest water, lower near-surface water temperatures and more frequent mixing events that reduce stratification as water rises during the rainy season, and warmest near-surface water temperatures and stronger and more persistent stratification during the dry season.Thus, the success of the simulations implies that extended simulations can be done at floodplain lakes to gain understanding of seasonal changes of stratification and ecological consequences.
Comparison with other 3-D models: AEM3D simulations were evaluated by examining statistics computed with measured and simulated quantities that can be categorized as state validations and process validations (Section 3) (Hipsey et al., 2020).Performance statistics of AEM3D for L. Janauacá were compared with simulations by other 3-D models (Si3D-L, Delft3D and MIKE3) and by AEM3D in other lakes.No published examples of process validation for other models or sites are available, and few 3-D hydrodynamic simulations have been done in shallow warm waters.Based on comparisons of simulations of temperatures (i.e., state validations), AEM3D's performance for L. Janauacá was similar to or better than performances of other 3-D model applications.
Among the comparisons, Poyang Lake (China) and Clear Lake (California, USA) are similar to Amazon floodplain lakes.Li et al.'s (2018) simulations of Poyang Lake with MIKE3 over 1 year had RMSEs, computed from daily near-surface temperatures, from 1.4 to 1.7°C.Rueda and Schladow (2003) applied Si3D-L to Clear Lake, a warm, polymictic lake with diurnal patterns of heating and wind events.Reported RMSE values based on temperature profiles ranged from 0.41 to 0.61°C.Woodward et al. (2017) and Zhou et al. (2021) applied AEM3D to Lake Argyle (Australia), a tropical reservoir; computed RMSEs were 0.5 and 0.96°C, respectively.RMSEs at L. Janauacá (Figures 6a 1 and 6b 1 ) are similar to these examples.While metrics such as RMSE are commonly applied, from the simulations it is clear that there can be large spatial and temporal variability in near-surface temperatures at different locations.For shallow lakes with low wind speeds, the extent of error, up to ∼1°C, can also be contrasted with cross basin temperatures differences of 2°C in the morning and even larger in the afternoon.From this perspective, as well as that of the percent relative error and the correlation coefficient, the modeling errors are relatively small.

Ecological and biogeochemical implications:
The process-based analysis of the seasonal variations in stratification and mixing and their modeling are relevant to ecological and biogeochemical studies in shallow warm, water conditions that are common in the tropics and occur during summer in temperate locations.For example, Amaral et al. (2018) suggested that variations in stratification and mixing linked to spatial and seasonal changes in water depths influence carbon dioxide dynamics and that sites with more wind exposure accumulate less dissolved carbon dioxide than wind protected areas.
Biogeochemical models of the production and atmospheric exchange of gases, such as methane, carbon dioxide or oxygen, typically use 1-D parameterizations of vertical mixing (Guo et al., 2023;Stepanenko et al., 2016;Tan et al., 2015Tan et al., , 2017)).Though computationally more demanding than 1-D formations, 3-D models enable improved representation of lateral and vertical exchanges.Applying AEM3D to an embayment, where differences in temperature between sheltered shallow areas and deeper water occur, revealed the importance of lateral exchanges driven by differential heating and cooling to CO 2 concentrations and fluxes (Amaral et al., 2021).
As we have shown, the duration of stratification and depth of mixing varies considerably over diel periods and over a range of depths of water and light attenuation.Hence, our results are pertinent to a wide variety of shallow lakes in different regions and trophic states.These variable periods of stratification have important implications for depletion of dissolved oxygen and accumulation and exchanges of dissolved nutrients, methane and carbon dioxide.Examples include Barbosa et al. (2020) for methane, Amaral et al. (2021) for carbon dioxide, and Melack and Fisher (1983) for dissolved oxygen and nutrients.These papers illustrate large vertical variations in concentrations of these gases in tropical lakes despite reduced stratification relative to results from other latitudes.Cortés et al. (2021)combined measurements and modeling to show the development of hypoxia and its alleviation with intermittent mixing in a temperate polymictic lake.MacIntyre et al. (2019MacIntyre et al. ( , 2021) demonstrated enhanced gas transfer velocities during periods with strong near-surface stratification in Amazon inland waters.
Future work: To extend our modeling results to other Amazon floodplain lakes and over whole river reaches would require information about regional meteorological and hydrological conditions.Meteorological measurements needed to drive hydrodynamic models are sparse for Amazon floodplains.Alternatives for regionalization may be weather forecasting or reanalysis results.Yang and Dominguez (2019) provide an example for the Amazon region of output from the Weather Research and Forecasting Model with outputs four times per day at 20 km resolution.Remote sensing using synthetic aperture radar or optical sensors can provide estimates of the variable area of the lakes and of associated vegetated areas (Fleischmann et al., 2022;Hess et al., 2015).Calculations of flooding and water movements produced by hydraulic models, such as those by Rudorff et al. (2014aRudorff et al. ( , 2014b) ) or Ji et al. (2019), would complement remote sensing results and improve temporal and spatial resolution.Large-scale and multiyear estimates of river discharge and inundation are also available (e.g., Paiva et al., 2013;Yamazaki et al., 2011).
Attempts to forecast and evaluate how climate change scenarios might influence stratification in lakes and related ecological processes are usually based on 1-D models of grid-averaged lake areas and depths and projected climatic conditions (e.g., Jansen et al., 2022;Woolway & Merchant, 2019).While instructive for identification of broad trends in seasonal patterns, as we have shown, the approach is not well suited to shallow floodplain lakes with their large variations in depth, area, transparency and strong diel mixing cycles.Moreover, simulations of climate-change induced variations in inundation suggest large regional differences across the Amazon basin (Sorribas et al., 2016) with further perturbations expected from construction of large hydroelectric dams (Flecker et al., 2022;Forsberg et al., 2017).Assessing the complex spatial and temporal variations expected among lakes, in general, and among tropical floodplain lakes, in particular, will benefit from the increasing availability of in situ and remote sensing systems to complement improved applications of one, two and three dimensional hydrodynamic and ecological models.

Conclusions
We quantified varying meteorology, water depths and optical conditions and the resulting diel and seasonal patterns in stratification and mixing in a representative Amazon floodplain lake.The field data were collected during 17 campaigns spanning two hydrological years.These data were used to apply and evaluate the three-dimensional hydrodynamic model AEM3D.The model successfully simulated the diel cycles of the variations of the thermal structure during the campaign periods, characterized by buoyancy frequency, depth of the AML, and other terms associated with the surface energy budget.Performance of AEM3D was evaluated by examining statistics computed with measured and simulated quantities during the campaign periods and the evaluation can be categorized as state validations and process validations.Based on comparisons of simulations of temperatures (i.e., state validations), AEM3D's performance for L. Janauacá was similar to or better than performances of other 3-D model applications.Mixing dynamics in the AML were controlled by diel cycles of heat fluxes, with wind shear and convection both contributing.Morning heating led to formation of a thin diurnal thermocline that up and downwelled with changes in wind velocity and led to differences in the degree of stratification.Despite considerable heat loss at night, the variable and often strong stratification meant that mixing only infrequently extended to the bottom at night.The infrequent full mixing creates conditions conducive to anoxia in these shallow lakes given their warm temperatures.The demonstration of successful modeling enables simulations over longer time periods to evaluate near-surface stratification and duration of stratification as these physical conditions modulate ecological and biogeochemical processes.

Data Availability Statement
A dataset has been generated and deposited into the Knowledge Network for Biocomplexity archive at the National Center for Ecological Analysis and Synthesis via https://doi.org/10.5063/F11R6P0Wwith the Creative Commons Attribution 4.0 International License (Zhou et al., 2022).The dataset contains meteorological data, measured and simulated temperature and buoyancy frequency data reported above.
The AEM3D model used for the hydrodynamic simulations in this paper has been obtained from HydroNumerics (Dallimore, 2019).The demo version of the model is available at https://www.hydronumerics.com.au/software/aquatic-ecosystem-model-3d.The version with full capacity can be purchased from HydroNumerics under a single-user license.

Figure 1 .
Figure 1.(a) Lake Janauacá and surroundings.Red line marks the boundaries of domain used for AEM3D simulations; yellow diamond marks the location of the research station (RB); white cross marks the location of the water level gauges.(b) Bathymetry of modeled region.LW marks the location of measurements of temperature profiles at low water, and HW marks the approximate locations of measurements of temperature profiles most other times including high water.The white dots mark the locations where vertical profiles were simulated (called water columns), and the black lines mark the horizontal extent of cross basin simulations called curtains.Curtain 1 also marks the open boundary.

Figure 2 .
Figure 2. Maximum depths of the lake (blue curve) from October 2014 to November 2016, and diffuse attenuation coefficient for photosynthetically active radiation (k d ) (blue diamonds).Blue dashed lines mark beginning of each field campaign (numbered from 1 to 17) and when AEM3D simulations were conducted.Number of days and initial depths of the AEM3D simulations are denoted by red and black numbers, respectively.
gravitational acceleration in m s −2 , ρ is density (kg m −3 ) at depth z and ρ 0 is the mean density of the profile SE W m −2 Sensible heat flux k d was least were typically 30 to 60 cycles per hour (cph) and reached 80 cph when cloud cover and winds were low.They were highest during low water, with N > 100 cph with k d > 7 m −1 (Figure 1).Typical values of N in the diurnal thermocline as it downwelled exceeding 10 cph whereas during low waters values exceeded 40 cph.When brief periods of near-isothermy occurred, values decreased below 5 cph.

Figure 3 .
Figure 3. Meteorological conditions and thermal structure between 22 and 27 August 2016 (simulation 17; HW site, Figure 1b): (a) incoming shortwave radiation (black) and effective heat flux in actively mixing layer (red); (b) air temperature (black), surface water temperature (blue dashed) and relative humidity (red); (c) wind speed (black) and wind direction (red); (d, e) measured and simulated temperatures; (f, g) buoyancy frequency derived from measured and simulated temperatures.White dots in panels d and f mark the depths of thermistors.Black lines in panels d and e are isotherms with 0.5°C intervals.30°C isotherm is labeled.

Figure 7 .
Figure 7. Root mean square error, percent relative error and correlation coefficient (R) computed with T(z) (a 1 -a 3 ), T(t) (b 1 -b 3 ), z AML (c 1 -c 3 ) and H * (d 1 -d 3 ) during the periods of validation simulations.All results were shaded with grayscale proportional to the initial depth of each simulation.

Figure 8 .
Figure 8. Contours of temperatures (a 1 -e 1 ) and buoyancy frequencies (a 2 -e 2 ) at five locations (Figure 1b) for simulation 17 (August 2016).Vertical scale varies among panels.WD is easterly August 22, southerly mid-day August 23 to early am August 25 with westerly winds for a brief period on the afternoon of August 24, westerly late afternoon August 25, northeasterly winds from the evening August 25 to morning August 27, and southerly winds the rest of the day.Depths of diurnal thermoclines are indicated by depths with highest buoyancy frequencies.

Figure 10 .
Figure 10.The same as Figure 9 but for early afternoon of 27 August 2016 (simulation 17).
with ρ a the air density, ρ w the density of water at the surface, C D the bulk transfer coefficient for momentum at instrument height, here 2 m, U z the wind speed measured at instrument height.This term is also known as the turbulent velocity scale due to wind.