Quo Vadis Lakes Azuei and Enriquillo: A Future Outlook for Two of the Caribbean Basin’s Largest Lakes

Lakes Azuei (LA) and Enriquillo (LE) on Hispaniola Island started expanding in 2005 and continued to do so until 2016. After inundating large swaths of arable land, submerging a small community, and threatening to swallow a significant trade route between the Dominican Republic and Haiti; worries persisted at how far this seemingly unstoppable expansion would go. The paper outlines the approach to a look forward to answer this question vis-à-vis climate change scenarios developed by the Intergovernmental Panel on Climate Change (IPCC). It uses numerical representations of the two lakes, and it examines how the lakes might evolve, deploying three different forcing mechanisms: that of weather and drift due to climate change, that of extreme events, such as hurricanes, and that of anthropogenic impacts, such as unintended water transfers between adjacent watersheds. Runs are executed Monte Carlo style using 11 different forcing combinations, each with a thousand instances of results generated by varying the numerous parameters that define the numerical models. The results are necessarily not precise and vary significantly as the forecast horizon expands, creating expanding envelopes of outcomes. Although some outcomes suggest a continued rise of the lake levels, most scenarios yield a reduction and recession of the lake waters.


Introduction
When lakes Azuei (LA), on the Haitian, and Enriquillo (LE) on the Dominican Republic side of the border, see Figure 1, started to expand in 2003, nobody knew, or at least anticipated, what would happen during the next 11 years. The initial expansions in the second half of 2003, while coming off relatively low levels, were not different from any previous expansions, typically followed by some shrinkage of the same or similar magnitude. Both lakes experience a constant cycle of expansion and shrinkage, and while, at times, highly asynchronous (LE would expand while LE would shrink and vice versa), the peaks and troughs have been modest through the decades.
Although the small and regular up and down signals, reflecting seasonal changes such as rainy and dry seasons, continued, a general quite strong upward trend established itself, which would last for about 11 years. Over these 11 years, both lakes expanded dramatically, swallowing up an estimated 40,000 acres (LE) and 6000 acres (LA) of farmland on both sides of the border [1]. LA rose 8 m and grew in surface area from about 115 km 2 to 137 km 2 (a 20% increase), while LE rose about 11 m and more than doubled its surface area from 155 km 2 to 350 km 2 . Both lakes, flanked by high mountains in the north (Neiba Mountains and W-extensions) and south (Baoruco Mountains and W-extension to Haiti), can only expand along their W-E axis that features shallow waters and slowly rising topography. In these west and east zones, just a moderate "centimeter rise" would cause another "square kilometers inundation" of farmland. The level rise in LA flooded communities in Lunette and Letant, inundated most of the customs buildings at the border crossing in Jimani, and for LE eventually forced the evacuation and abandonment of Boca de Cachon re-settlement During the years 2013 to 2015, when researchers of The City College of New York (USA) and the Instituto Tecnologicó de Santo Domingo (the Dominican Republic) carried out several data collection campaigns [2] they also carried out numerous interviews with residents, as well as representatives from the Instituto Nacional de Recursos Hidráulicos (INDRHI), the Oficina Nacional de Metorologicá (ONAMET), both in the Dominican Republic, and the Observatoire Nacional de Environnement at Vulnérabilité (ONEV), on the Haitian side. Although these anecdotal data were very helpful in understanding some of the lakes' histories, i.e., pre-Landsat-imagery-times, meaning before 1974, it was also clear that there was much concern. The lakes had been rising for 8-9 years straight, and there was no end in sight (the lakes eventually stopped growing in 2014 and have receded somewhat since then), making people wonder what would be happening in the future [3]. The concerns focused on the potentially detrimental impact on continuing lake level rise and what that would mean for the communities around the lakes. Questions asked were: will there be more community re-settlement situations (such as for people from Boca de Cachon)? How much arable land would be inundated by the lake's salty waters and, consequently, render the soil unusable for years to come? What impact would this have on the, albeit limited at the Haitian side and underdeveloped in the Dominican Republic, tourism of the lakes? What impact would this have on the fauna and flora of the lakes and their immediate surroundings? And lastly, to what extent would the trade route on the Jimani border crossing be further affected, potentially completely disrupting trade because the road could not be further elevated to keep up with rising water levels? During the years 2013 to 2015, when researchers of The City College of New York (USA) and the Instituto Tecnologicó de Santo Domingo (the Dominican Republic) carried out several data collection campaigns [2] they also carried out numerous interviews with residents, as well as representatives from the Instituto Nacional de Recursos Hidráulicos (INDRHI), the Oficina Nacional de Metorologicá (ONAMET), both in the Dominican Republic, and the Observatoire Nacional de Environnement at Vulnérabilité (ONEV), on the Haitian side. Although these anecdotal data were very helpful in understanding some of the lakes' histories, i.e., pre-Landsat-imagery-times, meaning before 1974, it was also clear that there was much concern. The lakes had been rising for 8-9 years straight, and there was no end in sight (the lakes eventually stopped growing in 2014 and have receded somewhat since then), making people wonder what would be happening in the future [3]. The concerns focused on the potentially detrimental impact on continuing lake level rise and what that would mean for the communities around the lakes. Questions asked were: will there be more community re-settlement situations (such as for people from Boca de Cachon)? How much arable land would be inundated by the lake's salty waters and, consequently, render the soil unusable for years to come? What impact would this have on the, albeit limited at the Haitian side and underdeveloped in the Dominican Republic, tourism of the lakes? What impact would this have on the fauna and flora of the lakes and their immediate surroundings? And lastly, to what extent would the trade route on the Jimani border crossing be further affected, potentially completely disrupting trade because the road could not be further elevated to keep up with rising water levels?
The answer to all these questions depends on the future climate of the Caribbean, and, hence, the lakes' region. The Caribbean region is considered susceptible to climatological changes, and studies demonstrate that even transient changes in climatological patterns significantly affect the Caribbean environment [4]. Recent studies have already shown evidence of change in climate patterns in the Caribbean region, i.e., an increasing trend in temperature and changes in precipitation patterns and water budgets [5][6][7][8][9][10]. Expectations are that the region's average temperature increases by approximately 1-3 • C by the end of the 21st century, accompanied by increased frequency and intensity of extreme rainfall events [11][12][13][14]. In case of regular precipitation events, the Providing Regional Climates for Impact Studies (PRECIS, [15][16][17][18]) model predicts up to 25% rainfall reduction in the Caribbean by the 2080s (medium-high emissions scenario), with more drying in the southern regions (below 18 • N, which includes Hispaniola) than in the north [8]. Other climate models, such as the ensemble models of CMIP3 and CMIP 5, also confirm a drier south Caribbean in the future [4].
North Atlantic cyclone activities have exhibited a significant impact directly on LE's water balance and indirectly on LA's [19]. Hence, future projections on tropical cyclone (TC) occurrence and intensity in the Atlantic basin must also be investigated. Future projections of cyclone characteristics have been challenging, yielding conflicting results, as have studies regarding their past changes [20][21][22]. Most studies agree with a reduction in cyclone frequency in the late 21st century. However, some studies have shown an increase in the number of storm events in the future [23,24]. In both cases, storm-based precipitation is projected to increase between 5% to 25% by 2100 [25][26][27][28][29][30][31].
How the two factors above combine and affect the lakes in question is the objective of this paper. The paper introduces the conditions and constraints of the model input parameters, followed by a selection of run time scenarios to expand the forecast horizon for the lakes' levels. Finally, it presents and discusses results and provides a summary along with an outlook.

Simulation Scenarios
IPCC has developed emission scenarios (SRES) based on possible future demographic, technological, and economic developments to provide input for evaluating climatic and environmental consequences of future greenhouse gas emissions and assessing alternative mitigation and adaptation strategies [32,33]. The scenarios consist of six sets, including A1FI, A1B, A1T, A2, B1, and B2, corresponding to low and high emission situations, among which A1B (middle-impact scenario) [34,35], A2 (medium-high) [36][37][38][39][40][41], and B2 (mediumlow) [38,39,42] are the most used scenarios in the assessment of the impact of climate change in hydrological studies. Out of the three scenarios, scenarios A2 and B2 were deployed in this study.
Changes in the region, such as land cover and land use, were not considered in the simulation scenarios, primarily due to insufficient data. Instead, emission scenarios were combined with predictions for North Atlantic cyclone alterations subject to climate change scenarios. Combining scenarios A2 and B2 with predicted cyclone changes yielded four base simulation scenarios: Scenarios A2 and B2 both project a decreasing trend in precipitation for Hispaniola and an increasing trend in evaporation [43]. Based on the results of PRECIS (ECHAM4 and HadCM3.5 models), Hispaniola's precipitation decreases by 22.66% (in 2050) and 24.17% (in 2100) under the A2 scenario, and 14.23% (2050) and 18.86% (2100) under the B2 scenario. Potential evaporation, on the other hand, increases by 6.34% (2050) and 15.04% (2100) under the A2 scenario, and 6.34% (2050) and 6.34% (2100) under the B2 scenario, with respect to the 1961-1990 baseline [43]. These rates were adopted for the lakes' future projections on the 2050 and 2100 horizons.
Considering a warmer climate, most models predict fewer storm occurrences [25][26][27]. However, simulations of TC formation are challenging, especially in the Atlantic basin due to the inadequate horizontal resolution of the climate models [25]. The predicted range of cyclone activity decrease differs between studies due to the variations in climate models and emission scenarios applied to the simulations. In 2010, Knutson et al. [20] suggested changes between +6% to −34% globally in TC frequency. In another study, based on CMIP5 models, simulations (RCP4.5 and RCP8.5 emission scenarios) showed a decrease between 21.9% to 60.6% in the total number of North Atlantic storms [21]. A recent study reported a decrease of 9.4% under the RCP4.5 scenario; however, the predicted value of decrease was not statistically significant in the North Atlantic basin. A few studies are showing an increase or no change in the frequency of the TCs. The global and regional synthesis of projections was assumed to be either a 50% increase or decreased TC activity for the North Atlantic basin [24]. Another projection, using the RCP2.6 emission scenario, yielded an increase in TCs in the North Atlantic by 14.9% and a 2.5% decrease using RCP8.5 [23].
Although there is no consensus on the projections for the frequency of TCs, there is a consensus on the increase in storm intensity and, consequently, in precipitation rates. Knutson et al. [22] suggest that the precipitation rate in the North Atlantic Basin is estimated to increase between 9.4% (Hurricanes Category 4-5) to 20.5% (Hurricanes Category 1-5) with an average of 17.3% for all TCs, by the late 21st century using emission scenario RCP4.5 [22].
Since the rate of change regarding North Atlantic cyclone frequency predictions varies, for this research, an increase and a decrease of ±20% and ±50% (forecast horizon is 2100) were introduced to examine the lakes' response to these imposed changes. Past studies agree on an increase in cyclone intensity and precipitation rates by the end of the 21st century. This study adopted a constant rate of +17.3% for TC precipitation for the simulations scenarios.
Combining climate change impact on precipitation, evaporation, and North Atlantic cyclone frequency and intensity, the runtime design yielded eleven simulation scenarios as follows: • B2 scenario without a change in cyclones; • B2 scenario and +50% increase in cyclone frequency; • B2 scenario and +20% increase in cyclone frequency; • B2 scenario and −50% decrease in cyclone frequency; • B2 scenario and −20% decrease in cyclone frequency; • A2 scenario without change in cyclones; • A2 scenario and +50% increase in cyclone frequency; • A2 scenario and +20% increase in cyclone frequency; • A2 scenario and −50% decrease in cyclone frequency; • A2 scenario and −20% decrease in cyclone frequency; • No climate change in future.
One additional (baseline) scenario was added for comparison purposes, which assumed no future climate change, thus no altered future precipitation and evaporation, and no change in cyclone characteristics.

Model Setup
Volatile dynamics of LE and the more stable behavior of LA is formulated through a conceptual water balance equation as follows [19]. where L, W, H, G, and I are lake, watershed, hydrological interactions of LA and LE, groundwater, and inter-basin water transfer components, respectively. This equation can be expanded as: where P L is the direct precipitation over the lake, E L is the direct evaporation depth from the lake, A L is the lake's surface area, r c is the run-off coefficient, P w is the precipitation depth over the watershed, A w is the watershed surface area excluding the lake, inf c is the average infiltration rate, f is the discharge factor, and Q is the inter-basin water transfer discharge to the lakes. In Equation (2), P D is equal to P w for LE and P w − P w for LA, in which P w is the equilibrium modifier of the LA water balance. In Equation (2), both H and G terms from Equation (1) are embedded in (in f c P D A w ), accounting for the connectivity of the lakes and their groundwater interactions. When (in f c P D A w ) is negative, outgoing seepage from the lake occurs. In this manner, the outgoing seepage computed for the LA model is added to the LE model and at the same time the elevation of LE is incorporated into the calculation of P w in the LA model (for more information, refer to [19]). To account for the impact on LA due to LE's fluctuations, volume-elevation rating curves, extracted from the lakes' bathymetry [44], were coupled with Equation (2).
For forecast simulations, infc was kept constant in the model (estimated value can be found in [19]). For the computation of lake and watershed surface area (A L and A w ), the rating curves of volume-area [44] for LE and LA were coupled with the lakes' models. Precipitation values (P L and P w ) were generated stochastically using long-term precipitation statistics in the region with a monotonically decreasing trend to account for A2 and B2 projections. The values for evaporation estimation (E L ) were gradually increased based on A2 and B2 projections.
It is known that cyclone activity influences runoff generation and the inter-basin water transfer in the Enriquillo basin [19,45]. Run-off coefficient (r c ) was kept constant in the simulations for both LA and LE assuming no land cover change in the region. However, the possibility of a storm incident with significant precipitation in a specific year changed the selection of r c value for LE (for more information, refer to [19]). The chances of TC activity in the vicinity of the lakes in a single year were around 23%. To account for TC future projections, ±20% and ±50% changes in storm frequency were applied to the annual TC occurrence probability. A gradual +17.3% increase in the intensity of the storm precipitationl within 100 km of the TC center was also considered in the simulations.
Although, in general, a significant transport of water occurs towards LE via the Cristobal Canal, this discharge is not present every year [19]. Based on [19] findings, interbasin water transfer is triggered mainly when consecutive storms hit the island, causing a high water production in neighboring watersheds, eventually draining into LE through the region's canal system. The only exception was through unexpected one-time events, such as the 2007 failure of the control structure on the canal system, resulting in uncontrolled water transfer to LE. Although the probability of such one-time occurrences was very low and they were not included in future simulations, the possibility of consecutive storms had to be considered. During years of storm activity close to the lakes, the probabilities for a single TC occurrence and consecutive incidents were approximately 76% and 24%, respectively. For future projections, storm distribution in any given year remained unaltered since the projection of future storm tracks was impossible. The average amount of Cristobal Canal discharge (Q) was kept constant for the years with occurrences of consecutive storms (f = 1). On the Haitian side, the inter-basin water transfer was set to zero (f = 0) due to its insignificant impact on the lake's water budget and the difficulty in predicting its occurrences in the future.
Both lake models were coupled together to reflect their connectivity when calculating the water balances. Using the described constraints in conjunction with the simulation scenarios provided the base for the future modeling effort with projection horizons of 2050 and 2100.

Results
The lake models (paired-lake-catchment model) were run using the eleven simulation scenarios mentioned earlier. For each scenario, the model was executed 1000 times to capture all the possible future lake' volumes and to generate a broad spectrum of lake responses. The random nature of storm occurrence yielded an envelope of possible results for each lake where smaller (narrower) envelopes reflect a higher degree of certainty in the lake responses. Figure 2 shows the model results for both lakes. Only three scenarios, i.e., "No climate change", "B2 scenario without a change in cyclones", and "A2 scenario without a change in cyclones", are presented for clarity's sake. The 90% confidence interval (90% CI), which is defined by the difference between the 5th and 95th percentile values of the cumulative distribution of the output variable, was used to quantify the uncertainty of the lakes' future at the projected horizons (shaded in grey in Figure 2). This region is the uncertainty of the model responses due to the variability in cyclone occurrence. scenarios provided the base for the future modeling effort with projection horizons of 2050 and 2100.

Results
The lake models (paired-lake-catchment model) were run using the eleven simulation scenarios mentioned earlier. For each scenario, the model was executed 1000 times to capture all the possible future lake' volumes and to generate a broad spectrum of lake responses. The random nature of storm occurrence yielded an envelope of possible results for each lake where smaller (narrower) envelopes reflect a higher degree of certainty in the lake responses. Figure 2 shows the model results for both lakes. Only three scenarios, i.e., "No climate change", "B2 scenario without a change in cyclones", and "A2 scenario without a change in cyclones", are presented for clarity's sake. The 90% confidence interval (90% CI), which is defined by the difference between the 5th and 95th percentile values of the cumulative distribution of the output variable, was used to quantify the uncertainty of the lakes' future at the projected horizons (shaded in grey in Figure 2). This region is the uncertainty of the model responses due to the variability in cyclone occurrence. However, when using the "no climate change" scenario, both lakes' variations are more expansive, suggesting they may grow by the 2050 and 2100 horizons. The chance of growth is higher for LA. Although the probability is low for the 2100 horizon, the 2050 horizon offers a considerably higher chance for LA to maintain its 2017 volume or grow beyond it (based on 90% CI). This is due to the interconnectedness of the two lakes that causes LA to rise when LE reaches or surpasses a certain level, impeding LA drainage. The downward trend of LE volume for no climate change is because of the assumption of controlled Cristobal Canal discharge, which was assumed only to occur when consecutive storms strike the island.
The results of all scenarios also indicated that uncertainty decreases with time (uncertainty is less in 2100 in comparison to 2050), as the envelope of possible results narrows However, when using the "no climate change" scenario, both lakes' variations are more expansive, suggesting they may grow by the 2050 and 2100 horizons. The chance of growth is higher for LA. Although the probability is low for the 2100 horizon, the 2050 horizon offers a considerably higher chance for LA to maintain its 2017 volume or grow beyond it (based on 90% CI). This is due to the interconnectedness of the two lakes that causes LA to rise when LE reaches or surpasses a certain level, impeding LA drainage. The downward trend of LE volume for no climate change is because of the assumption of controlled Cristobal Canal discharge, which was assumed only to occur when consecutive storms strike the island.
The results of all scenarios also indicated that uncertainty decreases with time (uncertainty is less in 2100 in comparison to 2050), as the envelope of possible results narrows down. The exception is for the "no climate change" scenario, which displays high uncertainty in both 2050 and 2100 horizons for LE and LA. All model responses have skewed distributions with all outliers on one side. These outliers are related to the North Atlantic cyclone activities, which cause a sudden growth in LE and subsequently in LA.
down. The exception is for the "no climate change" scenario, which displays high uncertainty in both 2050 and 2100 horizons for LE and LA. All model responses have skewed distributions with all outliers on one side. These outliers are related to the North Atlantic cyclone activities, which cause a sudden growth in LE and subsequently in LA. The difference between the result of A2 and B2 scenarios is statistically significant in 2050 and 2100 (1% significance level, p-value < 0.01). The A2 and B2 scenarios drift apart more when approaching the end of 21st century. This drift is more distinct for LA for the 2100 horizon. The difference between B2 sub-scenarios and A2 sub-scenarios also increases from 2050 to 2100. The comparison of A2 and B2 sub-scenarios shows the sensitivity of the lakes' future to changes in TC activities. For A2 sub-scenarios of LE, the projections are statistically different for 2050 and 2100 horizons. The only exception is between the two A2 sub-scenarios of 20% and 50% increase in TC activity, for which their mean and variance are not significantly different in 2050 (p-value > 0.1). This means that the 20% increase in cyclone frequency does not differ much from those with 50% increase. For B2, three sub-scenarios with no change, 20% increase, and 20% decrease in TC activity show similar projections in 2050. Other B2 sub-scenarios have generated significantly different results. As expected, LA is less sensitive to the changes in TC activities when comparing the difference within B2 sub-scenarios and A2 sub-scenarios. In staying consistent with LE results, the difference between "no change in cyclone", "20% increase in TC activity", and "20% decrease in TC activity" B2 sub-scenario results for LA are insignificant in 2050. Comparing A2 sub-scenarios also reveals the similarity of projections between 20% and 50% increases in TC activity and between no change, 20% decrease, and 50% decrease in TC activity in 2050. The values of mean, minimum, and maximum volume for each simulation scenario are summarized in Table 1. These values were derived from integrating model iterations for 2050 and 2100 horizons. The difference between the result of A2 and B2 scenarios is statistically significant in 2050 and 2100 (1% significance level, p-value < 0.01). The A2 and B2 scenarios drift apart more when approaching the end of 21st century. This drift is more distinct for LA for the 2100 horizon. The difference between B2 sub-scenarios and A2 sub-scenarios also increases from 2050 to 2100. The comparison of A2 and B2 sub-scenarios shows the sensitivity of the lakes' future to changes in TC activities. For A2 sub-scenarios of LE, the projections are statistically different for 2050 and 2100 horizons. The only exception is between the two A2 sub-scenarios of 20% and 50% increase in TC activity, for which their mean and variance are not significantly different in 2050 (p-value > 0.1). This means that the 20% increase in cyclone frequency does not differ much from those with 50% increase. For B2, three sub-scenarios with no change, 20% increase, and 20% decrease in TC activity show similar projections in 2050. Other B2 sub-scenarios have generated significantly different results. As expected, LA is less sensitive to the changes in TC activities when comparing the difference within B2 sub-scenarios and A2 sub-scenarios. In staying consistent with LE results, the difference between "no change in cyclone", "20% increase in TC activity", and "20% decrease in TC activity" B2 sub-scenario results for LA are insignificant in 2050. Comparing A2 sub-scenarios also reveals the similarity of projections between 20% and 50% increases in TC activity and between no change, 20% decrease, and 50% decrease in TC activity in 2050. The values of mean, minimum, and maximum volume for each simulation scenario are summarized in Table 1. These values were derived from integrating model iterations for 2050 and 2100 horizons.
The mid-response results suggest that LA shows a smaller shrinkage rate ranging between 8% and 20% (2050) and 19% and 43% (2100) in comparison to LE, which decreases by 84% for the "A2 scenario with a 50% decrease in cyclone frequency". The reference volume value was 3.22 km 3 for LE (2017) and 1.96 km 3 for LA (2017). Moreover, the model mid-response predicts that the volume of LE will be smaller than that of LA in 2050 and 2100 irrespective of the simulation scenario. The minimum and maximum possible volume of LE across all simulation scenarios are 0.3 km 3 and 4 km 3 in 2050, and 0.2 km 3 and 3.4 km 3 in 2100, which is a broader range than projections for LA. In 2050, LA minimum and maximum possible volume values are 0.4 km 3 and 2.3 km 3 , respectively, and 1.1 km 3 and 2.4 km 3 in 2100. This proves LE's higher vulnerability to climate change than LA. Changes to the lakes' volume can be translated to corresponding water levels, which can be computed from volume-elevation rating curves. The results show that the difference between the elevation estimation of mid-response within B2 and A2 scenarios was less than few centimeters (<10 cm for LE and <20 cm for LA). Therefore, the results were regrouped into "No Climate Change", "Climate Change under IPCC B2 SRES", and "Climate Change under IPCC A2 SRES" scenarios for each time horizon. Table 2  In the absence of climate change impacts, the model results indicate the possibility of lake expansion and shrinkage with about the same magnitudes of fluctuations as it has been historically. When considering the impact of climate change, the change in precipitation and evaporation rates have more influence on the lakes' dynamics than imposing a 50% increase in cyclone frequency (along with a 17.3% increase in the cyclone rainfall rate) which cannot compensate for the amount of water lost through evaporation. However, this does not guarantee a constant future shrinkage of the lake throughout time. Although the mid-response of the model iterations demonstrates a prolonged future shrinkage pattern of the lakes, the actual pattern of their fluctuations would vary depending on the time and intensity of the North Atlantic cyclones and consequent inter-basin water transfer. Analyzing single responses highlights that both lakes will still experience sudden growth (resulting from storm incidents and inter-basin water transfer) and lengthened periods of shrinkage and drying in between (Figure 4).  In the absence of climate change impacts, the model results indicate the possibility of lake expansion and shrinkage with about the same magnitudes of fluctuations as it has been historically. When considering the impact of climate change, the change in precipitation and evaporation rates have more influence on the lakes' dynamics than imposing a 50% increase in cyclone frequency (along with a 17.3% increase in the cyclone rainfall rate) which cannot compensate for the amount of water lost through evaporation. However, this does not guarantee a constant future shrinkage of the lake throughout time. Although the mid-response of the model iterations demonstrates a prolonged future shrinkage pattern of the lakes, the actual pattern of their fluctuations would vary depending on the time and intensity of the North Atlantic cyclones and consequent inter-basin water transfer. Analyzing single responses highlights that both lakes will still experience sudden growth (resulting from storm incidents and inter-basin water transfer) and lengthened periods of shrinkage and drying in between ( Figure 4).

Discussion
The development of two distinctively different models for each of the lakes and the introduction of climate change scenarios (providing a wide range of climate forcings) served as the computational background for future predictions of the two lakes. Three main simulation scenarios: (1) No climate change, (2) IPCC B2 SRES, and (3) IPCC A2 SRES were considered and then expanded to include projected changes in the North Atlantic cyclone activities in 2100, yielding eight more scenarios for a total of 11. The model was conditioned with optimized parameters derived from the hydrologic model design. Future projections for precipitation, evaporation, cyclone frequency, and cyclone intensity were varied in a Monte Carlo simulation to yield thousands of result datasets.
The key finding of this study answers the original question posed at the beginning of this paper: will the two lakes continue to increase in size and volume, thus aggravating the situation around them for the inhabitants and commerce? Although not 100% in the affirmative, the answer is clear: all climate change scenarios (3) yield the majority of parameter combinations that result in lake size decrease, suggesting that the lakes are very likely transitioning back to "normal" levels. Observations since 2017, the year when the

Discussion
The development of two distinctively different models for each of the lakes and the introduction of climate change scenarios (providing a wide range of climate forcings) served as the computational background for future predictions of the two lakes. Three main simulation scenarios: (1) No climate change, (2) IPCC B2 SRES, and (3) IPCC A2 SRES were considered and then expanded to include projected changes in the North Atlantic cyclone activities in 2100, yielding eight more scenarios for a total of 11. The model was conditioned with optimized parameters derived from the hydrologic model design. Future projections for precipitation, evaporation, cyclone frequency, and cyclone intensity were varied in a Monte Carlo simulation to yield thousands of result datasets.
The key finding of this study answers the original question posed at the beginning of this paper: will the two lakes continue to increase in size and volume, thus aggravating the situation around them for the inhabitants and commerce? Although not 100% in the affirmative, the answer is clear: all climate change scenarios (3) yield the majority of parameter combinations that result in lake size decrease, suggesting that the lakes are very likely transitioning back to "normal" levels. Observations since 2017, the year when the study stopped, have confirmed the recession of the lakes' surface areas that started a few years earlier in 2015.
Of all the factors and mechanisms that could potentially influence the outcome of any given run, changes in precipitation and evaporation rates for the climate change scenarios have a higher impact on the lakes' fate than changes in cyclone characteristics. This is relative, however. A 20% decrease in cyclone activities for LA does not cause a significant change in the lake's response. However, an increase of 50% in cyclone frequency yields a significant impact when considering scenarios B2 and A2. For LE, 50% decreases and increases yield significant impacts on the lakes' behavior, i.e., significant expansion and shrinking patterns which matter because of an ever-evolving climate change. In 2100, all responses of the two models are significantly different from each other. Although the average response of the models shows a smooth decreasing trend, some responses of the models suggest that both lakes are subject to sudden future rises and falls. Thus, despite a future with smaller lakes, the occurrence of sudden changes persists, causing similar impacts as before. Therefore, strategies to prevent adverse effects must be centered around proactive watershed management plans for the lakes and their neighboring watersheds.
The images in Figure 2 show a significant spread of the outcomes the further the simulation stretches, to the point where the results (Lake volume estimates) vary by 400 to 500%. This fact gives rise to question about the usefulness of the information presented and indeed introduces significant uncertainty. Uncertainty arises from three primary sources: climate predictions, model setup, and internal variability of model responses. Furthermore, uncertainties in climate predictions due to the scenario uncertainty, the GCM model uncertainty, and downscaling of the climate models' results add to the scope of uncertainty. There is an additional source of uncertainty called "internal variability of model responses" that express themselves as an internal variability in the hydrologic models. These are due to the stochastic behavior of the North Atlantic cyclones, where variability in intensity and distance causes changes in the lakes' states. The simulations are carried out assuming no land cover change, constant equilibrium index for the Lake Azuei watershed, no contribution of the Desaguas Canal discharge, and a regulated discharge of Cristobal Canal. These constraints, while necessary, limit the outcome horizon. All uncertainties add up and propagate through the model outputs, resulting in higher levels of prediction uncertainty and a large scatter in projected values. Yet, the presented results constitute a best-effort-possible given the dearth of data encountered and the complexity of the various levels of conceptualizations that needed to be accounted for and included in the models.

Future Outlook
As a way forward to obtain less uncertain results (the envelope of the outcomes would be much narrower), other emission scenarios, other climate change models (perhaps ensemble runs), and, more importantly, the spatial distribution of their results on Hispaniola should be incorporated. This could include more elaborate data collection campaigns or more focused data assimilation and sophisticated data imputation efforts. There is also a possibility of establishing competing model setups that see the introduction for geospatially distributed models that utilize a higher degree of time-variant hydro-climatic input datasets. In this vein, key hydrologic processes that have a high impact on model results, such as evaporation, should be better represented through an improved data spectrum needed for its calculation. Another alternative would be to examine hydrologic features such as the Canal contributions one by one, rather than the somewhat simplified approaches used in this study.

Data Availability Statement:
The data underlying the analyses in this paper are available from the lead author Mahrokh Moknatian. The interested reader is encouraged to contact her using the email given in the header of this manuscript. Some data are also available from the CUNY Digital Library at: Moknatian M., Piasecki M., "Report: Observational Time Series for Lakes Azuei and Enriquillo: Surface Area, Volume, and Elevation", New York, USA: Academic Works Digital Library, January 2019, https://academicworks.cuny.edu/cc_pubs/625/.