Modelica-Based Modeling and Simulation of District Cooling Systems: A Case Study

While equation-based object-oriented modeling language Modelica can evaluate practical energy improvements for district cooling systems, few have adopted Modelica for this type of large-scale thermo-fluid system. Further, to our best knowledge, district cooling modeling studies have yet to include hydraulics in piping networks alongside plant models featuring realistic mechanical systems and controls. These are critical details to include when looking to make energy and control improvements in many physical system installations. To fill these gaps, this study released new open-source district cooling models at the Modelica Buildings Library and leveraged these models for a real-world case study at the University of Colorado Boulder. The site includes six buildings connected to a central chiller plant featuring a waterside economizer. Several energy saving strategies are pursued based on the validated model, including control setpoint optimization, equipment modification, and pump setpoint adjustments. Results indicate that a combination of the studied measures can save the campus annually 84.6 MWh of energy, 8.9% of electricity costs, 58.0 metric tons of carbon dioxide emissions, while the waterside economizer cuts down chillers’ run times by 201 days/year, reducing maintenance costs and extending chiller life.


Introduction
Aggregating cooling equipment at the district scale provides important opportunities for cities and organizations to serve the cooling needs of their communities through financially viable means. Space cooling continues to grow faster than any other building end use [1], largely due to the improving standards of living and warming climate conditions across the globe. To meet climate action targets and reduce fossil fuel consumption, many experts recommend a shift from treating buildings as isolated entities to operating communities of buildings as interconnected systems. District cooling (DC) is one such communityscale technology that provides benefits both financially (e.g., economy of scale, centralized maintenance) and environmentally (e.g., more efficient equipment, higher renewable energy integration).
Computer modeling and simulation is an effective means to further improve the energy utilization of DC systems. Several tools have been adopted throughout literature for DC applications. Using TRNSYS, Anderson et al. [2] adjusted the operating time of chillers (compressor and absorption), cooling towers, and thermal storages with a lumped campus thermal load to minimize exergy destruction. Chow et al. [3] evaluated several water-cooled chiller plant configurations (including thermal storage, several heat rejection methods, and various pumping configurations) with TRNSYS for the plant/district simu-lation and HEVACOMP/DOE2 for the building thermal load models. Using MATLAB, Oppelt [4] implemented dynamic thermo-hydraulic models tailored specifically for DC networks, with buildings and the central plant modeled as single heat exchangers. Some works did not present their modeling framework nor their simulation tools, yet still ought to be mentioned. For example, Gang et al. [5] evaluated chiller plants (containing electric chillers, double-effect chillers for ice production, and an absorption chiller) with ice and hydro pumped storage systems to reduce peak cooling loads. In addition, Matsouka and Hil [6] optimized the number of cooling towers and condenser water (CW) flow rate for a central plant containing turbo-chillers for an online process simulation.
However, there are limitations among these modeling tools and existing methodologies. First, many of these programs use imperative programming languages that tightly intertwine the model equations with the numerical solution method, such as with C/C++ (e.g., EnergyPlus) or Fortran (e.g., TRNSYS). Further, some traditional languages are causal (e.g., MAT-LAB/Simulink) and do not support acausal modeling. While these traditional languages can be characteristically slow for development of heterogeneous cyber-physical models [7], including DC systems, they also are less flexible to support new use cases [8]. Further, most traditional tools have predefined numerical solvers; however, different solvers may be required to solve the differential equation problem for different use cases [9].
As a result, the past literature generally did not consider both the hydraulics in the distribution network and realistic plant configurations and controls, and had to make simplifying assumptions in one of the two categories instead. Yet, hydraulics in the distribution network can significantly impact the pump energy consumption [10], can cause controllability problems [11], and affect the effectiveness of the DC system. Meanwhile, as central plants improve their energy utilization through free cooling economizers and more advanced controls, it becomes increasingly critical to model these plant details to identify realistic energy solutions [12]. The prevalence of economizers and their associated control will only increase too, as the International Energy Conservation Code now requires air or waterside economizers for chilled water systems above certain capacities in most climate zones (all except extremely/very hot and humid 0A and 1A) [13].
Modelica [14] is a modeling language, governed by an open standard, that can effectively address these limitations. First, as an equation-based and object-oriented language that supports both causal and acausal modeling, Modelica allows developers to build up complex system models featuring both thermo-fluid system models (typically acausal) and realistic controls (typically causal) using a graphical, hierarchical approach, from individual equipment to subsystems and districts. These features allow models to be constructed with standard-form physical equations (i.e., not constrained to only input/output formulations) and may contain continuous, discrete, and hybrid differential equations, aiding flexibility for use cases and significantly reducing model development time [7]. Second, simulation code is generated automatically, and in contrast to many traditional tools, the model equations are separated from the simulation code. This code separation and the availability of several, selectable numerical solvers improves the developer's ability to successfully simulate complex system models, particularly "stiff" models with both fast and slow dynamics, which are common in DC applications. For applied energy projects, the development time reductions are particularly blaring when making small design modifications in an existing model to explore different energy efficiency measures. Lastly, the Modelica community has rich open-source libraries that span multiple domains, which enable users to construct their case study models from like systems and share resources among transdisciplinary, external groups.
The case study by Zabala et al. [15] leveraged the benefits of Modelica for DC applications, creating reduced-order models of a chiller plant based on high-fidelity Modelica models for real-time, model predictive control (MPC) applications. As part of this work, the District Cooling Open Source Library was publicly released [16]. While Zabala et al. [15] developed their chiller models from Modelica, they simulated the entire DC system externally using Python. This approach was highly effective for MPC, but numerous benefits of Modelica for DC analysis remain unexplored. Beyond this work, there is surprisingly very little Modelica-based research for DC analysis.
Leveraging the benefits of Modelica, this is one of the first DC system modeling and simulation studies for detailed energy analysis that considers hydraulics in the piping networks and realistic plant system design and control. To support future sci-entific endeavors, we developed and open-source released new DC models at the Modelica Buildings Library [17] (version 8.1.0), documented a systematic methodology for case studies of existing DC systems, and adopted the new models and proposed methodology for a real-world case study at the University of Colorado Boulder. The selected case study contains a chiller plant with a waterside economizer (WSE), a popular and effective way to improve chiller service life and reduce plant energy consumption [18,19,9] that is increasingly being required by code. Despite their importance, to our best knowledge, others have yet to model DC systems featuring a WSE. Further, we also evaluate and present the numerical performance (i.e., the scalability) of a detailed plant model connected to complete districts of several sizes, which is important to understand the capabilities and limitations of large Modelica models for DC analysis. Lastly, this paper's case study extends our previous work on the campus's central plant [20] by presenting a systematic methodology in modeling and simulation of an existing DC system, adding models for the distribution network with individual building loads, improving the central plant models, and extending the impact analysis to include financial and carbon savings.
It should be noted that there are many similarities between heating and cooling at district-scale [21]. In addition, Modelica-based modeling of other DHC types such as combined heat and power or ambient water networks [10,22,23] and chiller plants in non-DC applications [24,25,26,27,28] are also extensive. However, there is a general lack of publications on DC modeling and simulation [29], which we address herein through the contribution of new open-source models, a well-structured methodology, and an exemplifying case study. For demonstration, we select several energy efficiency measures that do not require capital investment, including control setpoint optimization, assessing the effectiveness of the nonintegrated WSE, and adjusting pump setpoints, but the opensource models and proposed methodology can be used for other analysis as well. This includes the evaluation of energy and financial impacts when replacing chillers or adding thermal storage, or system flexibility analysis made possible by the inherent thermal inertia of the distribution piping. Although out of the scope for this case study, other approaches such as these may require capital investment but can be highly lucrative.
The rest of this paper is organized as follows. To support future Modelica-based DC system studies, we summarize our systematic methodology in section 2. The case study's DC system is presented in section 3, followed by the corresponding Modelica models (section 4). After verification and validation of the DC baseline models in section 5, section 6 presents the baseline performance, energy efficiency improvements including quantification of the financial and environmental impacts, as well as the numerical performance and scalability analysis. Finally, section 7 presents our conclusions.

Methodology
As shown in Figure 1, we followed a systematic methodology for improving the energy performance of the DC system using Modelica, from problem definition to simulation-based impact evaluation. While the fundamental steps can be considered typical best practices when applying modeling and simulation for energy and control analyses, the details within each step are uniquely proposed for existing DC system case studies. In principle, this methodology is suitable for any existing DC systems with some degree of measured data available. The following sections describe each of the eight steps in this methodology.

Step 1: Understand the Problem
Before collecting data or beginning modeling, it is critical to first understand the problem. This includes defining the scientific experiment and understanding what question the models should be able to answer. The goal is not to have the models represent all aspects of the real system, but only the features that are necessary per the problem objective. For example, this work aimed to study the thermo-fluid and control performance of the DC system during normal operation. Emergency protocols and atypical operating conditions were not included. Similarly, any equipment pertaining to chemical treatment was excluded. After formulating the problem, we then draw system schematics and define control inputs and loops. These schematics are invaluable blueprints for successful model development.

Step 2: Collect Physical System Information and Operational Data
Physical system information and operational data can be used to design the system models, provide input data for the models directly, and verify and validate model accuracy. With existing DC systems, it is best to gather information and data from a variety of sources. District-scale systems often change over time as new buildings are added or plant equipment is exchanged. Thus, collecting redundant information helps verify sizing and performance specifications are up-to-date. The variety of sources includes mechanical and control specification documents, site visits (e.g., operator interviews, nameplates, control panels), and other manufacturer performance files.

Electrical System
For electrical measurement data, logged energy and power consumption of buildings and major equipment should be collected if available. Collecting real power data at 15-minute intervals is preferred for validation since it is also commonly used by electric utility companies when determining peak demand charges [30]. While equipment submetered data is the best way to validate equipment electrical performance, this is often unavailable. Thus, the electrical performance validation with measured data might only be possible at the system level, while equipment need to be validated by other means. If this is the case, equipment electrical performance can be verified with the available electrical specification documents to ensure the power is always operating within the expected range.

Thermo-Fluid System
For thermo-fluid measurement data, the plant and buildings contain sensors that log data at given intervals. We recommend the interval should preferably be no bigger than 15-minutes to capture transient dynamics of the cooling equipment as well as pumps and valves. However, to identify potential control instabilities, spot measurements may be done at shorter frequencies such as 1 minute. For chiller-only plants, it is important that the data collected span full and part load conditions at a minimum. Further, data for chiller plants with WSEs also need to include Free Cooling (FC), Mechanical Cooling (MC), and Partial MC modes, if present. In general, we recommend collecting at least two years of historical data if available to account for inevitable gaps and errors in the data sets as well as the variation of weather condition and cooling load profiles. Data types can include supply and return temperatures, mass flow rates, pump head pressure, and heat flow rates for major systems (e.g., plant, buildings) and equipment (e.g., chillers, WSE) on both the chilled water (CHW) and CW sides.

Step 3: Pre-Process Measured Data
Before using the data, it is necessary to pre-process the datasets so that they are suitable for use as model inputs and validation. While there are not data size restrictions on measured data for calibration purposes, data used as model inputs needs further refinement to reduce file sizes and verify expected performance. To pre-process the data, we first inspect the entire dataset for each building, plant, and equipment to see if they follow expected trends and fall within reasonable limits. Data points containing errors can then be eliminated. For example, the water in the building pipes stabilizes at the ambient room temperature (20 • C) when no water is flowing, which at times produces negative temperature drops due to variability across the sensor locations and sensitivity levels. These data points should be cleaned to correctly indicate zero cooling. If necessary, the datasets can be smoothed by a moving average, making it suitable for use as model inputs. For instance, CHW supply and return temperature data at the buildings can be smoothed to eliminate sensor noise from the measurements. Finally, we can round the building measured data to an appropriate decimal point to further reduce the file size for model inputs. For example, temperature measurements can be rounded to the nearest 0.1 • C and heat flow rates to the nearest 1 W.

Step 4: Develop the Models
There are numerous Modelica modeling resources available for DC system energy analysis that include hydraulic system dynamics and typical chiller plants. Open-source libraries contain packaged models from equipment through district-scale systems that can be adapted directly or modified for case study projects. We developed new models for DC applications that we open-source released in the Modelica Buildings Library version 8.1.0. These models include a DC plant with parallel electric chillers, parallel cooling towers with a bypass, and representative controls; a consumer connection model that can interchange the pipe models; and a vector-based distribution network model to allow various configurations. In addition to these new models we developed, the Data Center package of the Buildings Library [9] provides many equipment, subsystem, and control models for chiller plants with and without WSE. Further, AixLib [31], IDEAS [32], and Building Systems [33] are all open-source libraries that contain resources for DC system model development, in addition to the District Cooling Open Source Library [16] previously mentioned. While models from these open-source libraries can often be used directly, it is common that custom system-level models need to be developed for the case study application. As such, we present in this work the customized models we created for a real-world case study based on the new models we open-source released in the Buildings library.
For successful model development of large-scale thermalfluid systems, compartmentalizing the system into subcomponents is key. This allows the sub-components to be tested in isolation and errors to be detected. Further, adding unit tests helps to verify that the sub-components continue to perform as expected as the model evolves. Once all sub-components have been developed and tested, then the full system model can be assembled hierarchically.

Step 5: Run Numerical Simulations
Modelica simulations can run in multiple environments. Dymola [34] is a popular commercial tool that provides a userfriendly graphic interface and several numerical solvers suitable for DC system analysis. The OPTIMICA Compiler Toolkit [35], with its a graphical user interface IMPACT, is also available commercially. OpenModelica [36] and JModelica [37] provide open-source environments for Modelica-based modeling and simulation. While there are many suitable numerical solvers for solving DC system models, we frequently adopt CVODE [38] for its suitability for solving stiff problems. In our experience of adopting CVODE in Dymola, it typically simulates thermo-fluid systems quickly and robustly. For this case study, all models were simulated in Dymola using CVODE as the numerical solver with a tolerance of 10 −6 .
When simulating complex thermo-system models, debugging models is often required to improve numerical performance. From our experiences, applying second-order filters to step functions (e.g., for valve opening or pump signals) often helps to improve numerical stability. Such a filter is built-in to HVAC models of many Modelica libraries. Furthermore, much attention should be paid towards PI controller gains to avoid control instability. While out of the scope of this paper, there are many resources available to help first design numericallyefficient models [39], and then simulate them effectively [40].

Step 6: Validate the Baseline Model
District-scale heating and cooling systems can have larger uncertainties in measured data than for individual building sys-tems, particularly due to the lack of sufficient quantitative data at such scale [41]. While the uncertainty can have an averaging effect at district-scales for some features, it is potentially compounding for others [42]. Thus, it is important for modelers to recognize the impacts of uncertainly in measured data when validating complete DC systems. To address this challenge, we suggest verifying modeling assumptions and validating simulation results for DHC energy-based studies across two primary dimensions: (1) electrical and thermo-fluid systems; and (2) equipment and system levels.

Electrical System
For electrical validation, it is common that electrical data is available at the plant and buildings' main meters, but submetering is not available for individual equipment. Further, buildings' main meters typically have equipment beyond the cooling system. For example, central plants frequently contain both heating and cooling equipment that serve the district, as well as heating, ventilation, and air conditioning (HVAC) equipment serving the plant building itself and other electrical loads (e.g., interior and site lighting, plug loads). To address this uncertainty, we suggest to validate system level DC power consumption from simulations with respect to the range of expected DC real power from measured data at the building's main electrical meter. An upper limit to DC plant power range P Pla is given as where P Mea is the measured real power data from the electrical meter. To represent the lower limit P Pla , we adjust P Mea by scaling factor σ using In this formulation, σ represents the ratio of the peak demands from all individual cooling plant equipment D Pla,i to the total peak demand of the main meter, including other loads D Oth, j , which are not apart of the DC system. An additional safety factor f of 1.25 is included to account for the likelihood of equipment operating at different part load ratios. In this case study, σ was 0.42. An example of such analysis is shown in Figure 9.
For equipment-level validation, submetered data should be used if available. If not available, modelers can verify that the power consumption of each equipment is always at or under the stated nominal value during the annual simulation. An example of this alternate approach is shown in Figure 10, where box plots depict the annual distribution of power required for major equipment, normalized by their rated nominal powers.

Thermo-Fluid System
For thermo-fluid validation, on-site sensors provide temperature and flow rate data from various points throughout the system. Following ASHRAE Guideline 14 [43], we evaluate both the Coefficient of Variation of the Root Mean Square Error and the Normalized Mean Bias Error where y i is the individual measured data,ŷ i is the corresponding simulation data,ȳ is the mean of the measured dataset, and N is the total number of data points. With hourly data, the CVRMSE and NMBE need to be within 30% and 10% respectively for the model to be considered validated [43].
It is important to note that ASHRAE Guideline 14 is intended for measurement of energy and demand savings. Thus, they do not provide specific guidance on temperature and mass flow measurements. However, the same standard was accepted herein across all thermo-fluid data types. For temperature, units of Kelvin were adopted, while mass flow was represented in kg/s. Furthermore, it is possible that after pre-processing and removing sensor errors, a full year of clean measured data may not be available. In this case, both system and equipment level thermo-fluid conditions can be validated using representative time periods. For instance, this study selected two representative weeks in both summer and winter to validate the system models.

Step 7: Identify System Improvements
With the validated baseline model, several energy efficiency improvements for the DC system can be simulated to quantify their potential benefits. While providing guidance on which energy efficiency measures to select is out of the scope of this work, other literature provide comprehensive methodologies for identifying optimal energy efficiency measures for buildings [44] and control strategies for central plants [45]. For demonstration, we selected three measures that require little to no financial investment. This includes optimizing the CW supply temperature setpoint (T CW,set ), assessing the effectiveness of the WSE, and adjusting CW pump flow rate settings. These were selected after analyzing the baseline energy results and identifying high potential areas for improvement that can readily be adopted by the system operators.
As an example for the T CW,set optimization, we can formulate a single objective optimization problem to minimize the plant's annual energy consumption as where x represents the free parameters to be tuned during the optimization process, E Pla is the total plant energy during the optimization period t ∈ [t 1 , t 2 ), P CH is the power of the chillers, P CW is the power of the CW pumps, P CHW is the power of the CHW pumps, P CT is the power of the cooling towers, x contains the tuner parameter lower limits, and x contains the tuner parameter upper limits.
Three separate optimization problems are solved with T CW,set prescribed as where in Equation 9, T CW,set is a fixed temperature (Fixed T CW,set ) equal to constant x 1 ; in Equation 10, T CW,set is a fixed approach temperature (Fixed T app ) equal to x 1 offset from the measured ambient wetbulb temperature T wb ; and in Equation 11, T CW,set is an adjusted approach temperature (Adjusted T app ) offset x 1 from T wb and scaled by x 2 based on the central plant's part load ratio (PLR). In this case study, optimization problems are solved for t ∈ [0, 365) days with the Optimization library version 2.2.4 [46], a Modelica library that is released alongside Dymola that can solve nonlinear multicriteria parameter and trajectory optimization problems involving Modelica model simulations. There are numerous numerical optimization algorithms available in this library, and more information is available in [47]. For this case study, we adopt the local simplex method, and solve the optimization problems with optimization and simulation tolerances of 10 −5 and 10 −6 , respectively.

Step 8: Evaluate Impacts
Beyond energy consumption alone, evaluating the impacts that cooling infrastructure has on humans and the environment require higher-level evaluation metrics. Our approaches are as follows.

Energy Impacts
There are several metrics available to evaluate energy impacts for DC systems. Annual site and source energy consumption are two common metrics. In addition, peak power is a helpful indicator to understand the instantaneous rate of energy required at a given time. For cooling plants, the plant efficiency is often represented by kW/ton. This is defined as the ratio of total cooling plant power (including chillers, pumps, and cooling tower fans) to the total cooling load served by the plant. For chillers, the kW/ton is also a common metric to represent the equipment energy efficiency ratio. Similarly, the coefficient of performance is used as an equipment or system efficiency indicator. This is a unitless value representing the ratio of the cooling load to the electrical power. For buildings, common energy metrics also include the energy use intensity (ratio of cooling energy to total floor area served). In this study, we adopt the annual site energy consumption as the primary energy impact metric but also assess the kW/ton and peak power.

Economic Impacts
We evaluate the annual cost of electricity for operating the DC system based on the electric rate schedule from the local utility. While many electric utilities employ static pricing programs, there are increasing trends towards dynamic pricing programs (e.g., critical peak pricing, time of use pricing). Because the return on investment of different energy efficiency measures is often sensitive to the pricing program [48], it is important to represent the pricing scheme accurately. As an example for this case study, we adopt the static pricing program assigned to the central plant. In this case, electricity is charged on a monthly basis $0.00458 per kWh of energy used and $3.86 per kW of billing demand.

Environmental Impacts
Representative environmental indicators can be selected based on previous knowledge of the life cycle impacts of building systems. Operational phase energy consumption dominates typical building life cycle impacts and its greenhouse gas emissions represent the largest impact category [49]. This is particularly true for cooling plants where their primary function is to produce CHW for the district. Although calculating greenhouse gas emissions would be more comprehensive for assessing climate impacts, carbon dioxide (CO 2 ) is the dominating byproduct of fossil fuel-based electricity generation. Thus, it is reasonable to use CO 2 as a proxy indicator. To calculate CO 2 , Cambium's datasets provide hourly cost and operational data for the U.S. electricity sector [50]. For this case study, we adopt Cambium's 2018 Midcase Scenario data for the average CO 2 emission rate (kg/MWh) of generation induced by Colorado's regional load. This dataset provides hourly carbon emissions based on the dynamic electricity mix for the state, including intermittent renewable generation, while also including the effects of imported and exported power as well as distribution losses.

System Description
A case study of DC system modeling and optimization is conducted using the methodology mentioned in section 2. The case study site is a satellite campus of the University of Colorado Boulder with DC services provided to six buildings. This section describes the analyzed system, including the distribution network, the central plant, and the connected buildings.

Distribution Network
The distribution network is a radial layout ( Figure 2) connecting six buildings with a central plant. The total length of supply and return pipes is around 1.5 km. The pipes can be divided into 11 segments. Some segments are insulated pipes located in underground tunnels while others are direct buried. Additional information for the district pipes can be found in the Appendix (Table A.

Central Plant
As shown in Figure 3, the primary-only central plant includes two water-cooled chillers with a parallel WSE on both the CHW and CW legs. This WSE configuration is commonly referred as nonintegrated, and the WSE can only run when it is capable of meeting the entire cooling load. In contrast, an integrated WSE is in series with the chillers on the CW or CHW legs, such that the WSE can share the cooling load with the chillers when it cannot meet the entire demand. At this plant, CW is cooled by two open-circuit cooling towers before returning to the chillers/WSE. The nominal power for each chiller compressor (2 chillers total, 1 compressor per chiller), CW pump (2 total), CHW pump (3 total), and cooling tower fan (4 total) is 366 kW, 56 kW, 30 kW, and 22 kW, respectively. The plant's hierarchical control is implemented through a Supervisory Control and Data Acquisition system; at top Master Control level, the plant operates in either MC mode (running the chillers only) or FC mode (running the WSE only). Further details about the mechanical and control system configurations, including system schematics and nominal information for the plant equipment, can be found in [20].

Connected Buildings
The serviced buildings are primarily student dormitories with some academic and dining functions (Table 1). Each building is connected to the district via an energy transfer station (ETS) with a direct piping connection (no heat exchanger). A bypass loop at the ETS maintains minimum CHW return temperatures to the district. The six buildings account for a total of 93,990 m 2 gross floor area, while the cumulative peak load across all buildings is 1,772 kW. Like many district energy projects, the campus originally anticipated more buildings being added to the site than exists currently. As such, the sizing of the central plant accounted for an additional twelve buildings which have not been constructed yet. The measured load duration curves from this site (Figure 4) show how the cooling needs vary from building to building in 2018. The CHW heat flow rate is calculated using the supply and return temperature and mass flow rate sensors at each building's ETS. For example, the peak cooling demands range from 248 kW (building 1) to 653 kW (building 2). During 2018, buildings 2 and 4 had cooling demand year-round, while

Modelica Implementation
As part of this work, we developed and open-source released subsystem-and system-level models needed for DC applications in the Buildings library. The version 8.1.0 release includes a detailed central plant featuring both mechanical and control systems as well as several distribution piping networks. In addition, we are continuing to improve and expand the DHC package of the Buildings library with more system types and configurations, control implementations, and interconnected building features, to be included in future releases. Leveraging the new models to date, we developed the case study at the University of Colorado Boulder. In this section, we present the customized models pertinent to the case study in a top-down approach.
At the top level, the Modelica DC system model ( Figure 5) includes a central plant, a distribution network, and interconnected buildings. With a vector of fluid connectors in the Distribution Network model, a district with any number of buildings can quickly be represented, given that the computer has sufficient resources to solve the problem. Through packaged hierarchical modeling, each major DC component can seamlessly be replaced; for example, high-order or reduced-order building models can be integrated into the flexible modeling framework.

Boundary Conditions
Modeled weather and ground conditions are consistent with the local environment. Historical 2018 weather data from onsite sensors were collected and averaged on an hourly basis. We then overwrote the ambient drybulb temperature, dew point, relative humidity, and barometric pressure data from a typical meteorological year 3 (TMY3) data set corresponding to a regional weather station with the on-site data. The TMY3 file was then imported via the Weather Data block from the Buildings library, in which T wb is calculated. We assumed the tunnel temperature is the same as the ground three meters below the surface, represented by a sinusoid ranging from 4.6 • C (March) to 15.9 • C (September). This is consistent with average regional ground temperatures. Since we used measured data for the building loads, the model is not using solar irradiation data. Hence, overriding only the temperature data still provides consistent boundary conditions.

Distribution Network
A vector-style distribution network model was developed to allow interconnection with any number of buildings. Shown in Figure 6, the Distribution Network model consists of a series of Connection models. These distribution and connection models can adopt several different pipe models and are included in the new open-source release of the Buildings library. Supply and return pipes between each connection are modeled as plug flow pipes (using the model Buildings.Fluid.FixedResistances.PlugFlowPipe). Developed with DHC in mind, this pipe model efficiently and accurately evaluates pipe pressure drop, heat losses through the distribution network, and fluid transport delay. It also produced the smallest errors in comprehensive DHC experiments compared to other available models. While more information about this pipe model can be found in [51], the fundamental equations are as follows. The transfer of energy along the pipe and heat lost to the external environment can be expressed with a combined energy and continuity equation as where positiveq e is the heat loss from the pipe to the surroundings, x is the axial position with time t, ρ is the density, c v is the specific heat at constant volume, T is temperature, p is pressure, ν is the flow velocity, A is the pipe's cross sectional area, S is the pipe circumference, and f D is the Darcy friction coefficient, and k is the thermal conductivity. Eliminating negligible terms (e.g., diffusive heat transfer, heat loss due to pressure drops and wall friction), taking a Lagrangian approach, and integrating the simplified energy balance equation with respect to dT and  dt yields where subscripts in, out, and b signify the pipe inlet, pipe outlet, and surroundings, respectively; R is the fluid thermal resistance; and C is the thermal capacitance. The thermal capacity of the pipe wall is also added to the model to account for its thermal inertia, represented by a single capacitance at the outlet of each pipe segment. This thermal model is comparable to other works [52,53,54], reinforcing its suitability for DHC applications.
In addition to the heat transport and losses, the fluid delay time is described in the PlugFlowPipe model using the onedimensional wave equation where z(x, t) is the transported fluid quantity. This equation is used to describe the delay time of fluid parcels through the pipe and their associated properties (i.e., enthalpy), and produces similar results to that by Velut and Tummescheit [54].

Central Plant
The detailed central plant model represents the physics of the real mechanical and control systems, following the system schematics and control specifications of the plant. Even though this exact plant model cannot be publicly released for confidentially reasons, the DC plant model available in the Buildings library represents a typical configuration in both equipment and control. Following the proposed systematic methodology in this research, others can tailor this generic model for their own case studies, as done herein. While Modelica diagrams and model details about the case study's central plant can be found in [20], the mathematical models for key equipment (chiller, cooling tower, and pump) are as follows.
Chillers are modeled with the ElectricEIR model in the Buildings library, which is based on the DOE-2 electric chiller [55]. This model consists of three curves for the available capacity (CAPFT), full-load efficiency (EIRFT), and efficiency as a function of chiller PLR (EIRFPLR), defined as where T CHW is the chilled water supply temperature, T CW is the condenser water supply temperature, a i , b i , c i , e i , and f i are regression coefficients, and the PLR of the chiller is defined as whereQ is the chiller capacity andQ re f is the capacity at the reference evaporator and condenser temperatures where the curves come to unity. Cooling towers are modeled based on the variable speed Merkel model in EnergyPlus version 8.9.0 [56], which determines the total heat transfer between air and water entering the tower based on Merkel's theory. The fundamental basis for Merkel's theory is that the steady-state total heat transfer dQ tot is proportional to the difference between the enthalpy of air in the free stream h a and the enthalpy of saturated air at the wetted-surface temperature h s , represented by where c p is the specific heat of moist air at constant pressure, U is the cooling tower overall heat transfer coefficient, and A is the heat transfer surface area. For off-design conditions, Sheier's adjustment factors are included to modify the current UA value with respect to the current wetbulb temperature, air flow rate, and water flow rate. With these bases, the cooling tower performance is modeled using the effectiveness-NTU relationships for a counter-flow regime. Contrary to many past DC modeling projects, this work included pump models at the central plant that reflect the performance curves of the installed equipment. For example, the case study by Zabala et al. [15] assumed the pumping power was negligible, since their energy consumption was much smaller than the chillers. In contrast, the pumping power in this central plant is far from negligible, as is shown in Section 6.1 below. Indeed, pumping energy for large DC systems is often a major contributor to the overall energy consumption, and others have focused solely on different pumping configurations in order to save energy [57].
Because pump models that use affinity laws -change in pressure ∆p ∝ r(t) 2 , with r(t) being the rotational speed, and volumetric flow rateV(t) ∝ r(t) -can lead to singular sets of equations, we implement a composite pump model from the Modelica Buildings Library that is consistent with affinity laws. To ensure that solutions to the differential algebraic system of equations posed by the thermo-fluid model can be computed robustly and efficiently by Newton-based solvers, the pump is formulated in such a way that the resulting equations of the fluid flow network has a unique solution in each operable region, and is differentiable in all inputs. While complete details for the pump formulation are available in [58], the fundamental formulation when the pump operates far from the origin is as follows. Let δ = 0.05 be a small number that is below the typical normalized pump speed and S n = {(V i , ∆p i )} n i=1 be the user-supplied performance data at full speed r = 1, withV i ≥ 0 and ∆p i ≥ 0 for all i ∈ {1, . . . , n}. Here, n represents the total number of operating points. For conditions r > δ (i.e., far from the origin), the affinity laws are satisfied while the maximum volumetric flow rateV max and maximum pressure change ∆p max are linearly extrapolated aṡ ∆p n and (20) The fan performance curve for r(t) > δ is defined as with the curve end points represented by ∆p + (1,V max ) = 0 and ∆p + (1, 0) = ∆p max , while h(·, S n ) is a cubic hermite spline that mapsV to ∆p, and ∆p(V(t)) represents the model of flow resistance approximated by a linear function as In addition to conditions r > δ, Wetter [58] also defines formulations for near origin (r < δ/2) and composite (r ∈ [δ/2, δ]) conditions to complete the pump model. Lastly, a set of data points is defined for electrical power consumption P n = {(V i , P i )} n i=1 per the pump's performance documentation. Total efficiency η is then computed as from the flow work W f lo and electrical power input P, with W f lo defined per the first law as In Figure 7, the CW pump curves used to define S n and P n are shown as an example. When looking to revise flow setpoints during MC and/or FC modes, it is important to ensure that the minimum flow does not surpass the minimum continuous stable flow (MCSF) point. As an upper limit, the pump is restricted by the maximum power of the motor. Further, there are system-level limitations on the flow rate as well, including the minimum flow required by the cooling towers to keep the surfaces wetted, the minimum condenser flow required by the chillers, and the minimum pump speed required to meet the cooling towers' static lift [59].

Connected Buildings
Because the CHW supply and return temperatures and the mass flow rate are metered at the ETS, building models were created to accept these tabulated data files directly. The building model (Figure 8) includes the district-side CHW piping and equipment. The rejected heat from the building's cooling equipment is transferred to the district's CHW at the ETS (using the model Buildings.Fluid.HeatExchangers.Heater T). This model measures the CHW supply temperature, and returns the CHW based on the tabulated temperature difference. The control valve modulates the flow to match the tabulated mass flow rate for each building.
As seen in Figure 8, the building-side HVAC systems (including the primary CHW pumps), piping, and thermal zones are not modeled. Therefore, the heat transferred at each building represents the metered heat flow rather than the cooling loads at the thermal zones. Some limitations of this approach are that (1) some design optimization between building and district operation cannot be tested and (2) the thermal comfort at the rooms cannot be verified. However, the simplification does provide an accurate and numerically efficient implementation for time series loads metered at the ETS, which is commonly available among existing DC systems. Models to be added to the library will include lumped thermal zone load models to increase the flexibility of the modeling platform.

Electrical System
The central plant's electric power was validated with historical metered data from the plant's main electrical meter (see The upper and lower limits of this expected power range were determined using the methodology described in section 2.6.1. While this approach was necessary in this case study based on the available measured data, the validation could be improved by installing submeters on major cooling equipment.  For equipment-level validation, the electrical power of each equipment operated at or below their nominal rated power throughout the entire annual simulation (Figure 10). The box plot effectively depicts the distribution characteristics of each equipment's power consumption, normalized by the nominal rated power of each equipment. In Figure 10, the thick middle line represents the median, the upper and lower ends of the box represent the upper and lower quartiles, the whiskers represent power rates outside of the middle 50%, and the circles depict the outliers. In this simulation, while the chillers, CW pumps, and CHW pumps tended to run at low part load ratios (with median normalized powers less than 0.4), the cooling tower fans frequently ran at full power, as exemplified by the median normalized power near 1. When on, the cooling towers ran at full power for 50% of the time, while they ran at part loads between 0.1 and 1 for 25% of the time.

Thermo-Fluid System
Two typical periods from summer (Aug. 1-14, 2018) and winter (Jan. 28 -Feb. 11, 2018) seasons were selected for thermo-fluid validation. The results are summarized in Table 2 for the summer period and Table 3 for the winter period. On an hourly basis, the simulations fell within the acceptable ranges of modeling uncertainty for all locations. Please note that the acceptable ranges per ASHRAE Guideline 14 [43] depend on the measurement frequency (e.g., hourly, monthly) and other factors. The acceptable ranges quoted are specifically for hourly calibration with whole building simulation. During the summer period, the plant operated in MC mode with the chiller meeting the entire cooling demand. The buildings were well controlled across all data types, and temperature errors were small across all locations.
During the winter period, the plant operated in FC mode, with the WSE meeting the entire cooling demand. Two buildings (1 and 3) did not have a cooling load during this period, and thus were excluded from winter validation. Similar to the summer period, all simulation test locations fell within the acceptable ranges of modeling uncertainty, as seen in Table 3. By validating the DC model across multiple layers and system types, we can gain confidence in the model's representation of the real system for further energy analyses.

Baseline Energy Performance
In the annual simulation, the central plant consumed 552 MWh. As seen in Figure 11, the CW pumps represented the major energy consumer in the cold months, while In colder seasons, the CW pumps dominated the plant's energy consumption for three main reasons: (1) two pumps run during FC mode while only one pump runs during MC mode when one chiller is on; (2) the CW pumps operate at a lower efficiency during FC than MC (as shown in Figure 7); and (3) the CW side of the WSE has higher-pressure resistance compared to the chiller's condenser. Figure 12 shows the detailed analysis of the central plant operation for two sample windows during FC (Feb. [13][14][15][16][17][18]2018) and MC (Aug. [19][20][21][22][23][24]2018) modes. During FC mode, the plant's overall kW/ton was higher than that during MC mode. This surprising result is primarily due to (1) the lower cooling load during FC (260 kW average during the mid-February week compared to 2,700 kW during the late August week) and (2) the higher power required for the CW pumps during FC mode. Even though the cooling load is small during FC mode, the plant still needs to run equipment at their minimum limits. Conversely, the high-efficiency, variable speed, variable flow chillers during MC mode allow the chillers to operate at higher efficiency while serving larger cooling loads. During FC mode, the cooling tower fans run periodically, but the plant often bypasses the cooling towers to maintain T CW,set . The cooling tower bypass is used more consistently during MC mode, largely because T CW,set is higher in MC than in FC mode. Further, the CW mass flow rate through the WSE is consistently and significantly greater than the CHW mass flow for both FC and MC modes. This results in a very small ∆T on the CW side compared to the CHW side for the WSE (Figure 12c) and chillers (Figure 12f).

Energy Efficiency Improvements
Several energy efficiency measures were simulated based on the validated baseline results. As demonstration cases, three measures that do not require additional financial investments were pursued: (1) optimizing T CW,set during MC mode, (2) assessing the effectiveness of the WSE, and (3) revising the CW pump flow rate setpoints during FC mode. The results and discussion from these cases are as follows.

Condenser Water Supply Temperature Optimization
As seen in Table 4, all three T CW,set optimizations reduced the plant's annual energy consumption by 2.5% to 4.4%. Both the Baseline and the Fixed T CW,set optimization case (Equation 9) have fixed temperature setpoints. Simply changing the fixed T CW,set from 15.6 • C to 18.7 • C can save 2.5% energy. While some engineers reduce T CW,set as low as possible (around 15 • C) to allow the chillers to operate at a higher COP [60], there is often a trade-off between chiller energy and cooling tower fan energy [61]. This trade-off was present in this DC simulation. Raising T CW,set to a higher level caused chiller energy to increase 5.5%, but saved 42.8% in cooling tower fan energy. These results reinforce the need to consider the entire plant's operation when optimizing control setpoints.
Interestingly, the Fixed T app (Equation 10) and Adjusted T app (Equation 11) optimizations produced equivalent energy savings, reducing the annual energy consumption by 4.4%. Other researchers similarly studied Fixed and Adjusted T app based on wetbulb and part load ratios; for example, Liu and Chuah [62] found that optimizing the Adjusted T app (based on wetbulb, chiller load ratio, and chiller and cooling tower fan performance characteristics) saved more than 4% annual energy for chiller plants without WSE, with the greatest savings occurring in climates with high seasonal variation of T wb . However, additional savings from Fixed to Adjusted T app were not seen in this case study, primarily due to the presence of the WSE and the central plant's control logics. In this DC system, the optimized T CW,set is only applicable during MC mode; thus, the seasonal variation in T wb and part load ratios were reduced to the warm MC seasons only. Based on these results, it is recommended that the plant control T CW,set based on a fixed T app of 1.9 K above T wb .  A few limitations of the T CW,set optimization should be noted. First, we used a local optimization method, and it is possible that only a local rather than the global optimum was found. Future cases can use a hybrid approach to reduce the search space with a global optimization method, before selecting a local optimum. Second, the energy savings listed assume that the measurements are ideal, but physical T wb sensors can be inaccurate and prone to drift [63]. Thus, the energy savings from the Fixed and Adjusted T app cases may be overstated. Despite this, T CW,set based on wetbulb approach temperatures are adapted in industry [18,64,25].

The Effectiveness of the Waterside Economizer
To evaluate the savings that the WSE can bring, we simulated the same plant model but without the WSE. The results are shown in Figure 13. Through this experiment, we determined that the nonintegrated WSE saves 6.9% of the plant's annual energy consumption. While these savings are notable, the largest benefit of the WSE is the reduction in chiller run time. With the nonintegrated WSE, the total run time of the two chillers decreases from 373 to 172 days per year, which is a 54% reduction. These large changes in run time extend the service life of the chillers and reduce maintenance costs. Conversely, the total run times of the CW pumps and cooling towers are higher by approximately 202 days each when the WSE is added. This is because of the control specification that two CW pumps run during FC mode, and a lower T CW,set is required to meet the T CHWS setpoint with the WSE than with the chillers, causing the cooling tower fans to run more often. While out of the scope for this study, further energy efficiency improvements are likely possible by reconfiguring the WSE to be connected in series with the chillers on the CHW side (e.g., integrated). An integrated WSE is often the most energy efficient solution because it allows the WSE and chillers to share the cooling load even when the WSE cannot handle the entire load [19,12,65].  Figure 13: Energy consumption and total running times of equipment for DC system models with (baseline) and without the WSE.

Condenser Water Flow Rate Adjustments
With high CW pumping energy from the baseline model, we tested various pump flow rates and staging configuration to reduce the annual pumping energy during FC mode. The tested cases include operating one CW pump near its maximum flow rate (1-CWP-126), one CW pump at 100 kg/s (1-CWP-100), and two CW pumps at 50 kg/s each (2-CWP-50). As seen in Table 5, the 1-CWP-126 case saved 5.2% energy, while the other two cases saved 10.2%. In addition to energy savings, it is beneficial to operate one pump instead of two when possible to improve the pump service life and reduce maintenance costs. Figure 14 provides further insight into the CW flow rate adjustments. Interestingly, the 2-CWP-50 case operates at the lowest efficiency during FC mode, around 34% on average. However, the pumps are also operating at a lower power and with higher temperature differences across the CW-side of the WSE, which results in higher energy savings. These results fall in line with modern industry recommendations to decrease CW flow rate while increasing the temperature difference to save operating costs [59]. Pump flow rate adjustments such as these present practical cost-free retrofit measures that the central plant can readily implement.

Impact Evaluation
By improving the energy efficiency of the DC system, additional cost savings and greenhouse gas reduction are achieved. Figure 15 summarizes the results across all simulation cases. All cases except for the case without the WSE reduced energy, operating costs, and carbon emissions. From individual measures, the greatest cost savings came from the Fixed and Adjusted T app optimizations, primarily due to the larger decrease in monthly peak powers. Meanwhile, the largest carbon savings from individual measures came from changing the CW pump flow settings during FC mode to 1-CWP-100 and 2-CWP-50. It is important to note that the magnitude of the energy savings from individual measures (2-7%) can be considered insufficient with respect to the model accuracy. Thus, these results do not suggest that the DC system operators should employ individual measures alone. However, the magnitude of energy savings from combined measures is greater than the 10% minimum threshold recommended in international standards [66], and thus were recommended to the system operators, who agree with the findings. With combined measures of the Fixed T app optimization and operating two CW pumps at 50 kg/s, the best savings were achieved, reducing annual energy consumption by 84.6 MWh (15.3%), costs by $930 (8.9%), and carbon emissions by 58.0 metric tons (15.0%). While the cost dollar savings appear low, this is because of the relatively cheap cost of energy per the pricing structure stated in section 2.8.2. While this work evaluated economic and climate impacts as byproducts of energy efficiency strategies, future work can optimize these impact factors directly.

Numerical Performance
For this case study, all simulations ran using Dymola 2021 on a Linux operating system. The computer contained 32 GB of RAM and a 3.60 GHz Intel ® Xeon ® CPU. With CVODE As a scaling test, the simulation was repeated with the same plant model connected to multiples of the six-building district network (Table 6 and Figure 16). Across all simulation cases, the total cooling demand at the plant was held constant to satisfy system sizing requirements. The numerical problem contained 13 nonlinear SOE, regardless of the number of buildings (n). While the total number of nonlinear SOE was constant with respect to n, the dimension of the largest nonlinear system increased in size proportional to (n+5), i.e., linear. The number of continuous time states also increased near linearly with respect to n (more specifically, O(n 0.85 )). These two factors typically have the greatest impact on computing time. Shown in Figure 16, the computing time scaled by O(n 2.1 ) and O(s 2.5 ), where s represents the number of continuous time states. Over the range of district sizes, the computing time scaled closer to linearly when n was less than 20, while the quadratic to cubic scaling was present for districts larger than 20. With the given computer configuration (primarily limited by RAM capacity of 32 GB), the maximum number of buildings possible was about 300.
While scalability studies are limited for DHC, a few previous studies involving small to medium sized districts present similar results. Schweiger et al. [67] found that simulation CPU time scaled linearly with the number of buildings; however, that study only evaluated small district sizes with two to nine end users. Our results indicated similar trends, where districts with one to 20 customers had linear scaling, and it was not until medium sized districts were included that the results scaled at higher orders. Similarly, Jorissen et al. [68] simulated districts from 1 to 20 buildings and found computing time scales quadratically with DASSL solver, while Euler integration scaled linearly.

Conclusion
There is a growing interest and need to include district-scale cooling in advanced modeling platforms. While Modelica has been successful in chiller plants and district heating applications, DC applications remain scarce. This case study is one of the first modeling and simulation studies to include detailed plant thermo-fluid and control system models as well as representational network piping hydraulics, both of which are critical for detailed energy analysis. This contribution was enabled by the Modelica language, which is often more adept than traditional tools to model and simulate complex, stiff thermo-fluid system models with realistic controls, particularly when considering small system modifications through energy retrofit actions. We developed and open-source released new models at the Buildings library to support DC analysis of this variety, which provided a foundation for the real-world case study exemplified in this paper. Further, we aimed to include plant configurations (e.g., the free cooling WSE) that are underrepresented in modeling studies but common in practice, while also addressing concerns for the suitability of large physics-based Modelica models for large-scale district analysis (e.g., the scalability tests).
The case study included detailed models of the central plant's mechanical and control systems, a complete network topology, and simplified building models to reflect the thermo-fluid dynamics of the real-world system. This enabled us to quantify practical energy efficiency improvements that the system operators can readily implement. However in addition, it is important to identify the impact of our engineering systems on humans and the environment beyond energy consumption alone. To this end, we evaluated the electricity costs, carbon emissions, and equipment run time reductions as byproducts of the existing system configuration and energy improvement strategies. The best savings were achieved through combined measures of the optimized Fixed T app and operating two CW pumps at 50 kg/s, reducing annual energy consumption by 84.6 MWh (15.3%), costs by $930 (8.9%), and carbon emissions by 58.0 metric tons (15.0%), while the waterside economizer reduced the total run time of chillers by 201 days/year. Future studies will expand on this preliminary work to optimize both costs and environmental impacts directly.

Acknowledgements
This material is based upon work supported by the U.S. Department of Energy's Office of Energy Efficiency and Renewable Energy (EERE) under the Advanced Manufacturing Office, Award Number DE-EE0009139. This work was also supported by the Building Technologies Office of the U.S. Department of Energy, under contract numbers DE-AC36-08GO28308 and DE-AC02-05CH11231. Further, this work emerged from the IBPSA Project 1, an international project conducted under the umbrella of the International Building Performance Simulation Association (IBPSA). Project 1 will develop and demonstrate a BIM/GIS and Modelica Framework for building and community energy system design and operation. The authors would also like to thank the campus facilities team for their assistance with data collection, expert advice, and overall support of this project.

Appendix A. Distribution Network Parameters
See