The Sacramento River is a vital source of water for landowners and irrigation water suppliers in the Colusa Subbasin in addition to providing 80 percent of inflow to the Sacramento-San Joaquin Delta . The Glenn-Colusa Canal originates from the Sacramento River and flows through the Colusa Subbasin into Colusa County via the local canal system, providing irrigation water for the Glenn-Colusa Irrigation District, which is the largest irrigation district in the Sacramento Valley . The only principal aquifer in the Colusa Subbasin is contained in the freshwater-bearing sediments and stream channel deposits of the Modesto Formation, Tehama Formation, Riverbank Formation and Tuscan Formation, with the majority of the fresh groundwater in the project area contained within the Tehama Formation . Soils conducive to recharge, based on the SAGBI rating, range from excellent to very poor suitability, with the most suitable soils in the project area falling within the good to poor suitability range .Hydrologic conditions of interest in the project area include a natural cone of depression in groundwater elevation levels in the middle to southern half of the parent model domain . Based on recent groundwater monitoring data from the California Department of Water Resources , seasonal variations in groundwater elevation contours show a cone of depression forming in the center of the parent model domain during the Fall 2020 season due to stresses to the aquifer, such as excessive groundwater pumping for agriculture during the summer. Groundwater levels usually recover in the springtime, depending on how much precipitation occurred during the winter, but Spring 2021 contours still show lower groundwater levels in the same area with steeper groundwater head gradients,plastic garden pot as opposed to other areas in the parent model domain. Groundwater typically flows away from topographically high areas like the foothills of the Coast Ranges in the west, away from the Sacramento River in the east, and flows towards the topographically low areas where the cone of depression forms.
This area is heavily influenced by agriculture, and high amounts of pumping for agricultural use may be the cause of groundwater level depletion and the subsequent formation of the cone of depression and areas of low groundwater elevation.In 2019 and 2020, TNC recruited several landowners into their fall season bird habitat and ag-MAR incentive program. A total of 8 field sites were enrolled in TNC’s incentive program, which are all located in Colusa County, California . The nearest major cities to the ag-MAR sites are Williams, Colusa, and Meridian. The field sites modeled in this study are located near the city of Meridian, along the western bank of the Sacramento River . The fields are part of the Davis Ranches property. This area was chosen for the child model domain because Field 2 was the only field site that participated in flooding during consecutive years . The crop type grown on these fields is rice, but the fields were fallow during the period of flooding during the fall. The SAGBI recharge rating for the area in and around Field 2 and Field 15 are good to moderately good suitability for deep percolation in the soils .The data collected and used for the two groundwater models were primarily extracted from the Sacramento Valley Groundwater-Surface Water Simulation Model , which is similar to the California Central Valley Groundwater-Surface Water Simulation Model . Both SVSim and C2VSim were developed by the California Department of Water Resources, but SVSim is more specific to the northern Sacramento Valley, while C2VSim encompasses the entire Central Valley. Other sources of data include the California DWR, the United States Geological Survey , and Davids Engineering, Inc., the consulting company hired by TNC to conduct the field methods and technical analysis of the recharge program at each field site during 2019 and 2020. The data extracted from these sources were used in the MODFLOW model development process for parameterizing the packages for the parent and child models. More information on the development and calibration of C2VSim and SVSim can be found in Brush et al., 2013, and California Department of Water Resources, 2022, respectively. Layer stratigraphy data was extracted from SVSim, as were the layer and aquifer parameters, such as specific storage , specific yield , and vertical and horizontal hydraulic conductivity .
The layer stratigraphy and parameter data were used to define the aquifer parameters of all 9 groundwater model layers in both the parent and child models in their discretization packages. The uppermost layers in both models are relatively thin, with an average thickness ranging from 10 to 50 m. The lower layers are the thickest, ranging from about 50 to 200 m thick. The total thickness of all layers combined in both models is approximately 600 m . Each field site was flooded with a depth of 4 inches of water maintained for 30 continuous days during the fall, and 60 days of deep percolation data was collected , with deep percolation rates calculated for each site . These calculations were made using a mass balance approach – in which the inflows to the system are equal to the outflows – along with groundwater level measurements that were taken either periodically or continuously throughout the flooding period at groundwater monitoring wells adjacent to the field sites . Recharge took place at Field 2 during two periods, 9/18/2019-11/16/2019 and 9/17/2020-11/16/2020. Recharge took place at Field 15 during one period, 9/18/2019-11/16/2019. The average recharge rates at the three field sites ranged from 0.011 to 0.023 m/day . Deep percolation rates for all other days, or model stress periods, in the remainder of the domain are rates averaged by month from SVSim. Davids Engineering also monitored groundwater levels in four wells adjacent to the field sites . We used the groundwater level data provided by Davids Engineering for our child model’s head observation package. The average depth to groundwater between all four of the monitoring wells ranged from 5.3 to 7.7 m of depth . The software used for the models constructed for this study is the MODFLOW-2005 Three-Dimensional Finite-Difference Ground-Water Model developed by the U.S. Geological Survey . MODFLOW-2005, hereafter referred to as MODFLOW, uses packages that simulate the effects on groundwater flow processes of wells, rivers, lakes, and other relevant aspects and boundary conditions. A list of the packages used in both models is included in Table 4. MODFLOW can also vertically simulate aquifer systems with different geologic layers that may be specified as confined or unconfined. MODFLOW’s mathematical solution for simulating groundwater flow through the center of each cell in a model follows the partial-differential groundwater flow equation .The larger of the two models developed for this study was a steady state model that encompassed most of the area of Colusa County . As is the case with steady state models, all model parameters remained constant over a specified amount of time,draining pots so all input data that varied over time were averaged to one value for input.
The purpose for making the large model steady state was to have all inputs to the system equal to the outputs, stabilizing the system as a whole. This provided a more general simulation of the regional hydrologic processes , which allowed for faster model runs and processing times. The results of this large steady state model provided information on flow magnitude and direction that was used to constrain the boundaries of the child model. The domain for the large model spans an area of approximately 160,000 square meters and includes 9 layers, modeled after the layers in SVSim.The parent model has 140 rows and 113 columns, with a cell size of 400 m × 400 m. The number of stress periods, or days, is 1 since the parent model is run under steady state conditions. The boundary conditions used to define the model domain include 2 general head boundaries , 2 constant head boundaries , and the river package to simulate the effects of rivers, streams, and canals that flow through the model domain . The smaller of the two models developed for this study was a transient model that was located within the area of the parent model domain, including the ag-MAR field sites . Since the child model is transient, all data input for each package varied over time. The starting date of the child model was defined to be 8/18/2019, one month before the flooding occured, and the end date of the model is 11/16/2020, the last day of the flooding period. Like the parent model, the coordinate system used to define the child model was NAD83 California Albers. The child model has 58 rows and 58 columns, and each cell in the model has a size of 100 m × 100 m. The number of stress periods, or days, included in the child model is 456. Groundwater fluxes and directions of flow were extracted from the results of the parent model in order to inform the boundary conditions of the child model . The well package was used to simulate those fluxes. We defined the flux from each well in all model cells on the child model’s boundaries to represent the general flow that occurred in the same area of the parent model. In more specific terms, cell-by-cell flow data was extracted for each cell in every layer from the parent model, and the cell-by-cell flows were used to inform the groundwater fluxes at the boundaries of the child model. We also used the River package to represent the small part of the Sacramento River that acts as the eastern boundary of the child model. However, the RIV package is only present in the first layer of the model, so we used the GHB package for the eight layers directly below the river cells. The recharge package was used to represent the 2019 and 2020 ag-MAR field sites, with deep percolation data from Davids Engineering.Three different scenarios were run using the child model to test different parameters and their influence on the hydrologic response seen in the model outputs. The three scenarios were designed to assess whether the timing, frequency, and amount of recharge significantly influences the hydrologic results of the child model. Scenario 1 was the baseline model run, using the original and unaltered deep percolation data from Davids Engineering. No alterations or calibrations to the model inputs were made for this scenario. For Scenario 2, we increased deep percolation rates from Davids Engineering by one order of magnitude to ascertain whether we could see a significant hydrologic response. Only the deep percolation data was altered for Scenario 2; all other inputs remained the same as Scenario 1. It should also be clarified that a tenfold increase in recharge rates is not realistic. However, we wanted to test this hypothetical extreme solely to see a significant response. For Scenario 3, we increased the duration of recharge to span ten years.
We chose to use the 2019 deep percolation data, and repeated that year ten times to see how increasing the duration of recharge would affect the child model results. The 2019 deep percolation data input for Scenario 3 was unaltered; only the duration of recharge and the time discretization of the child model was increased to span ten years. Scenario 3 was designed to represent ten consecutive years of implementing ag-MAR on these field sites.The results analyzed from the steady-state parent model include model performance plots and model outputs. We assessed the performance of the parent model with a one-to-one plot of the observed heads that were input into the model, with the simulated equivalent heads that the model calculated at each groundwater monitoring well in the parent model domain . The different colors on the plot represent the binned data of the residual heads, which were calculated as the difference between the observed and simulated equivalent heads. The different bins represent different groups of residual data ranges. The gray line represents the line of equal values, or one-to-one line, between the observed and simulated equivalent heads.The same binned residual head data were plotted on a map of the parent model domain, showing the location of the groundwater monitoring wells and the residual head value . Like the residual head bins plotted in the one-to-one plot in Figure 12, the bins for the residual map were calculated from the difference between the observed and simulated equivalent head values at each groundwater monitoring well on the map.